Computing integral of a function which is a multiplication of two other integral functions in R-CodePudding

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

rm(list=ls()); set.seed(1234)
G=2 ; # Suppose 2 groups
n_i<-2 # There are two individuals per group
nTot<-4 # In total we have 4 individuals 
z_ij<-runif(nTot, 0, 5)
T_ij<- runif(nTot, 5, 10)
Data<-round(data.frame(id=rep(1:nTot), group=rep(1:G, rep(2,G) ),z_ij, T_ij),1)
Data
  id group z_ij T_ij
1  1     1  0.6  9.3
2  2     1  3.1  8.2
3  3     2  3.0  5.0
4  4     2  3.1  6.2

For every individual(row), I have the following function f(x,y)= xy\exp(x y z[id]) where upper limit in integrating for x is Tij[id].

I would like a code in R that can compute the following integral

Group 1:

$\int_{0}^{5} \left( \int_{0}^{T_{ij}[1]} xy \exp( x y z[1]) dx * \int_{0}^{T_{ij}[2]} xy \exp( x y z[2]) dx \right)dy = \int_{0}^{5} \left( \int_{0}^{9.3} xy \exp( x y 0.6) dx * \int_{0}^{8.2} xy \exp( x y 3.1) dx \right)dy$

Group 2

$\int_{0}^{5} \left( \int_{0}^{T_{ij}[3]} xy \exp( x y z[3]) dx * \int_{0}^{T_{ij}[4]} xy \exp( x y z[4]) dx \right)dy = \int_{0}^{5} \left( \int_{0}^{5} xy \exp( x y 3) dx * \int_{0}^{6.2} xy \exp( x y 3.1) dx \right)dy$

My attempt which yields an error Error in integrate(function(y) as.numeric(sapply(split(sapply(1:nTot, : evaluation of function gave a result of wrong lengt

GroupInt<-rep(NA,G)
for(i in 1:G) ## s<-3
{
GroupInt[i]<- integrate( function(y) 
  as.numeric( sapply( split(     sapply( 1:nTot, function(id)  integrate(function(x)
    x*y*exp( x y Data$z_ij[id] ), 0, Data$T_ij[id] )$value   ), Data$group),prod) )[i] , 
0 , 5 )$value
}

CodePudding user response：

set.seed(1234)
G <- 2 ; # Suppose 2 groups
n_i <- 2 # There are two individuals per group
nTot <- 4 # In total we have 4 individuals 
z_ij <- runif(nTot, 0, 5)
T_ij <- runif(nTot, 5, 10)
Data <- round(data.frame(id = rep(1:nTot), group = rep(1:G, rep(2,G)), z_ij, T_ij), 1)

# function for the inner integral
fInner <- function(df) {
  Vectorize({function(y) {
    prod(
      sapply(
        seq_along(df),
        function(i) integrate(function(x) x*y*exp(x   y   df$z_ij[i]), 0, df$T_ij[i])$value
      )
    )
  }})
}

GroupInt <- sapply(1:G, function(grp) integrate(fInner(subset(Data, group == grp, select = c("z_ij", "T_ij"))), 0, 5)$value)
GroupInt
#> [1] 2.173417e 16 1.534357e 14

# compare to the analytical solution
fAnal <- function(df) {
  exp(sum(df$z_ij))*(41*exp(10) - 1)*prod(exp(df$T_ij)*(df$T_ij - 1)   1)/4
}
GroupInt2 <- sapply(1:G, function(grp) fAnal(subset(Data, group == grp, select = c("z_ij", "T_ij"))))
GroupInt2
#> [1] 2.173417e 16 1.534357e 14

all.equal(GroupInt, GroupInt2)
#> [1] TRUE