library (reprex)
#> Warning: package 'reprex' was built under R version 3.4.3
I want to calculate (1) absolute of determinant of Jacobean of a two functions with respect to a and b. Then integrate as in step (2)
(2) I calculated the integral over a and b such that the integral from zero to one. The problem that Jacobean should be calculated numerically at specific points for a and b. To avoid that I calculated partial derivative for each function at aand b and combine them in a matrix "d". I defined a function of a and b; "func". When I calculated integral over a and b for the function "func", I got the following error message:
Error in determinant.matrix(x, logarithm = TRUE, ...) : (list) object cannot be coerced to type 'double'
How I can fix this problem. Also, how can the boundary points; zero and one not included in the integral. i.e. the integrals started above zero and reach one (since the functions will not be defined at zero and one).
The following is the codes and results. Thank you.
library (reprex)
Warning message:
package ‘reprex’ was built under R version 3.4.3
> reprex()
Rendered reprex is on the clipboard.
> library("Deriv")
Warning message:
package ‘Deriv’ was built under R version 3.4.3
> c1<-list()
> c1[[1]]<-Deriv(~(choose(10,0)*beta((a*(1-b))/b+0,(1-a)*(1-b)/b+10)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,1)*beta((a*(1-b))/b+1,(1-a)*(1-b)/b+9)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,2)*beta((a*(1-b))/b+2,(1-a)*(1-b)/b+8)/beta((a*(1-b))/b,(1-a)*(1-b)/b)),"a")
> c1[[2]]<-Deriv(~(choose(10,0)*beta((a*(1-b))/b+0,(1-a)*(1-b)/b+10)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,1)*beta((a*(1-b))/b+1,(1-a)*(1-b)/b+9)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,2)*beta((a*(1-b))/b+2,(1-a)*(1-b)/b+8)/beta((a*(1-b))/b,(1-a)*(1-b)/b)),"b")
> c1
[[1]]
{
.e1 <- 1 - b
.e3 <- (1 - a) * .e1/b
.e5 <- a * .e1/b
.e6 <- digamma(.e3)
.e7 <- .e3 + 10
.e8 <- .e3 + 8
.e9 <- .e3 + 9
.e10 <- 1 + .e5
.e11 <- 2 + .e5
.e13 <- digamma(.e5) - .e6
((10 * (digamma(.e10) - digamma(.e9)) - 10 * .e13) * beta(.e10,
.e9) + (45 * (digamma(.e11) - digamma(.e8)) - 45 * .e13) *
beta(.e11, .e8) + beta(.e5, .e7) * (.e6 - digamma(.e7))) *
.e1/(b * beta(.e5, .e3))
}
[[2]]
{
.e1 <- 1 - b
.e2 <- 1 - a
.e3 <- .e1/b
.e5 <- .e2 * .e1/b
.e7 <- a * .e1/b
.e9 <- digamma(.e3 + 10)
.e10 <- digamma(.e3)
.e11 <- digamma(.e5)
.e12 <- .e5 + 10
.e13 <- .e5 + 8
.e14 <- .e5 + 9
.e16 <- .e2 * (.e11 - .e10) + a * (digamma(.e7) - .e10)
.e17 <- 1 + .e7
.e18 <- 2 + .e7
((.e2 * (.e11 + .e9 - (digamma(.e12) + .e10)) + a * (.e9 -
.e10)) * beta(.e7, .e12) + (10 * .e16 - 10 * (.e2 * (digamma(.e14) -
.e9) + a * (digamma(.e17) - .e9))) * beta(.e17, .e14) +
(45 * .e16 - 45 * (.e2 * (digamma(.e13) - .e9) + a *
(digamma(.e18) - .e9))) * beta(.e18, .e13)) * (.e3 +
1)/(b * beta(.e7, .e5))
}
> d1<-matrix(unlist(c1),ncol=2)
> d1
[,1] [,2]
[1,] Expression Expression
> library("Deriv")
> c2<-list()
> c2[[1]]<-Deriv(~(choose(10,0)*beta((a*(1-b))/b+0,(1-a)*(1-b)/b+10)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,1)*beta((a*(1-b))/b+1,(1-a)*(1-b)/b+9)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,2)*beta((a*(1-b))/b+2,(1-a)*(1-b)/b+8)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,3)*beta((a*(1-b))/b+3,(1-a)*(1-b)/b+7)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,4)*beta((a*(1-b))/b+4,(1-a)*(1-b)/b+6)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,5)*beta((a*(1-b))/b+5,(1-a)*(1-b)/b+5)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,6)*beta((a*(1-b))/b+6,(1-a)*(1-b)/b+4)/beta((a*(1-b))/b,(1-a)*(1-b)/b)),"a")
> c2[[2]]<-Deriv(~(choose(10,0)*beta((a*(1-b))/b+0,(1-a)*(1-b)/b+10)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,1)*beta((a*(1-b))/b+1,(1-a)*(1-b)/b+9)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,2)*beta((a*(1-b))/b+2,(1-a)*(1-b)/b+8)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,3)*beta((a*(1-b))/b+3,(1-a)*(1-b)/b+7)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,4)*beta((a*(1-b))/b+4,(1-a)*(1-b)/b+6)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,5)*beta((a*(1-b))/b+5,(1-a)*(1-b)/b+5)/beta((a*(1-b))/b,(1-a)*(1-b)/b)+
+ choose(10,6)*beta((a*(1-b))/b+6,(1-a)*(1-b)/b+4)/beta((a*(1-b))/b,(1-a)*(1-b)/b)),"b")
> c2
[[1]]
{
.e1 <- 1 - b
.e3 <- (1 - a) * .e1/b
.e5 <- a * .e1/b
.e6 <- digamma(.e3)
.e8 <- digamma(.e5) - .e6
.e9 <- .e3 + 10
.e10 <- .e3 + 4
.e11 <- .e3 + 5
.e12 <- .e3 + 6
.e13 <- .e3 + 7
.e14 <- .e3 + 8
.e15 <- .e3 + 9
.e16 <- 1 + .e5
.e17 <- 2 + .e5
.e18 <- 210 * .e8
.e19 <- 3 + .e5
.e20 <- 4 + .e5
.e21 <- 5 + .e5
.e22 <- 6 + .e5
((10 * (digamma(.e16) - digamma(.e15)) - 10 * .e8) * beta(.e16,
.e15) + (120 * (digamma(.e19) - digamma(.e13)) - 120 *
.e8) * beta(.e19, .e13) + (210 * (digamma(.e20) - digamma(.e12)) -
.e18) * beta(.e20, .e12) + (210 * (digamma(.e22) - digamma(.e10)) -
.e18) * beta(.e22, .e10) + (252 * (digamma(.e21) - digamma(.e11)) -
252 * .e8) * beta(.e21, .e11) + (45 * (digamma(.e17) -
digamma(.e14)) - 45 * .e8) * beta(.e17, .e14) + beta(.e5,
.e9) * (.e6 - digamma(.e9))) * .e1/(b * beta(.e5, .e3))
}
[[2]]
{
.e1 <- 1 - b
.e2 <- 1 - a
.e3 <- .e1/b
.e5 <- .e2 * .e1/b
.e7 <- a * .e1/b
.e9 <- digamma(.e3 + 10)
.e10 <- digamma(.e3)
.e11 <- digamma(.e5)
.e13 <- .e2 * (.e11 - .e10) + a * (digamma(.e7) - .e10)
.e14 <- .e5 + 10
.e15 <- .e5 + 4
.e16 <- .e5 + 5
.e17 <- .e5 + 6
.e18 <- .e5 + 7
.e19 <- .e5 + 8
.e20 <- .e5 + 9
.e21 <- 1 + .e7
.e22 <- 2 + .e7
.e23 <- 210 * .e13
.e24 <- 3 + .e7
.e25 <- 4 + .e7
.e26 <- 5 + .e7
.e27 <- 6 + .e7
((.e2 * (.e11 + .e9 - (digamma(.e14) + .e10)) + a * (.e9 -
.e10)) * beta(.e7, .e14) + (10 * .e13 - 10 * (.e2 * (digamma(.e20) -
.e9) + a * (digamma(.e21) - .e9))) * beta(.e21, .e20) +
(120 * .e13 - 120 * (.e2 * (digamma(.e18) - .e9) + a *
(digamma(.e24) - .e9))) * beta(.e24, .e18) + (.e23 -
210 * (.e2 * (digamma(.e15) - .e9) + a * (digamma(.e27) -
.e9))) * beta(.e27, .e15) + (.e23 - 210 * (.e2 *
(digamma(.e17) - .e9) + a * (digamma(.e25) - .e9))) *
beta(.e25, .e17) + (252 * .e13 - 252 * (.e2 * (digamma(.e16) -
.e9) + a * (digamma(.e26) - .e9))) * beta(.e26, .e16) +
(45 * .e13 - 45 * (.e2 * (digamma(.e19) - .e9) + a *
(digamma(.e22) - .e9))) * beta(.e22, .e19)) * (.e3 +
1)/(b * beta(.e7, .e5))
}
> d2<-matrix(unlist(c2),ncol=2)
> d=rbind(d1,d2)
> d
[,1] [,2]
[1,] Expression Expression
[2,] Expression Expression
> func=function(a,b){abs(det(d))}
> #hcubature(func,lowerLimit = c(0,0),upperLimit=c(1,1))
> #install.Packages("pracma",lib="/my/own/R-packages/")
>
> #Another way to integrate
> library(rmutil)
Attaching package: ‘rmutil’
The following object is masked from ‘package:stats’:
nobs
> int2(func,a=c(0,0),b=c(1,1))
Show Traceback
Rerun with Debug
Error in determinant.matrix(x, logarithm = TRUE, ...) :
(list) object cannot be coerced to type 'double'
> library (reprex)
> library (reprex)
> reprex:::reprex_addin()
Loading required package: shiny
Warning: package ‘shiny’ was built under R version 3.4.3
Listening on http://127.0.0.1:4467
Rendered reprex is on the clipboard.
>