计算函数的积分，该函数是 R 中另外两个积分函数的乘积

Question

我有以下数据。数据有 2 个组 (G)，每个组有 2 个个体 (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

对于每个人（行），我有以下函数 f(x,y)= xy\exp(x+y+z[id]) 其中 x 的积分上限是 Tij[id].

我想要一个可以计算以下积分的 R 代码

第 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$

第 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$

我的尝试出错了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
}

Answer 1

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(
        1:nrow(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

计算函数的积分，该函数是 R 中另外两个积分函数的乘积

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

r

function

dplyr

integral