R Error : non-conformable arguments when doing double integration-CodePudding

I have the following data. The data has 2 groups (G) each with 2 individuals (n_i).

library(cubature)
set.seed(1)
G=2 ; # Suppose 2 groups
n_i<-2 # There are two individuals per group
nTot<-4 # In total we have 4 individuals 
z_j<-rbinom(nTot, 1, 0.5)
T_j<- runif(nTot, 5, 10)
Data<-round(data.frame(id=rep(1:nTot), group=rep(1:G, rep(2,G) ),z_j, T_j),0)
Data
  id group z_j T_j
1  1     1    0    6
2  2     1    0    9
3  3     2    1   10
4  4     2    1    8

For every individual(row), I have the following function f(x,y)= (x*exp(2y))^Zj[id] * exp(-x*exp(2y)*Tj[id]).

I am trying to evaluate the following double integral in R as follows.

Group 1

$\int_{-20}^{20} \int_{0}^{20} \left( \prod_{j=1}^{2} [x\exp(2y)]^{z_{j}} \exp[ - x\exp(2y) T_{j} ] \right)dx dy = \int_{-20}^{20} \int_{0}^{20} [x\exp(2y)]^{0} \exp[ -6 x\exp(2y) ] *[x\exp(2y)]^{0} \exp[ - 9 x\exp(2y) ] dx dy$

Group 2

$\int_{-20}^{20} \int_{0}^{20} \left( \prod_{j=1}^{2} [x\exp(2y)]^{z_{j}} \exp[ - x\exp(2y) T_{j} ] \right)dx dy = \int_{-20}^{20} \int_{0}^{20} [x\exp(2y)]^{1} \exp[ -10 x\exp(2y) ] *[x\exp(2y)]^{1} \exp[ - 8 x\exp(2y) ] dx dy$

I would like to use the quad2d function in the pracma package for double integration. Attempt code for the double integral yields Error in wx %*% Z : non-conformable arguments : as shown below.

fInt<- function(df) {
    function(x,y) {
    prod(
      sapply(
        1:nrow(df),
        function(i)  (x*exp(2*y) )^df$z_j[i]*exp( - x *exp(2*y)* df$T_j[i])
      )
    )
  }
}

GroupInt <- sapply(1:G, function(grp) quad2d(  fInt(subset(Data, group == grp, select = c("z_j", "T_j"))), 0, 20, -20, 20))
Error in wx %*% Z : non-conformable arguments

When you simplify manually the integral a solution exists. For example group 1 we have $\int_{-20}^{20} \int_{0}^{20} \exp(-15x \exp(2y)) dx dy$

quad2d(  function(x,y) exp(-15*x *exp(2*y)),  0, 20, -20, 20)
[1] 348.5557

CodePudding user response：

I believe that's because the integrand is not vectorized. Could you try:

fInt <- function(df) {
    Vectorize(function(x,y) {
    prod(
      sapply(
        1:nrow(df),
        function(i)  (x*exp(2*y) )^df$z_j[i]*exp( - x *exp(2*y)* df$T_j[i])
      )
    )
  })
}

EDIT

I took a look at the code of quad2d and I observed that the integrand is applied to two matrices having the same shape. So this should work:

fInt<- function(df) {
  function(x,y) {
    arr <- abind::abind(x, y, along = 3)
    apply(arr, c(1, 2), function(xy){
      x <- xy[1]; y <- xy[2]
      prod(
          sapply(
            1:nrow(df),
            function(i)  (x*exp(2*y) )^df$z_j[i]*exp( - x *exp(2*y)* df$T_j[i])
          )
      )
    })
  }
}