R中的积分误差：限制为Na或NAN

Question

fInt1 和 fInt2 下面的两个函数之间的区别是附加的乘法项 df$ui[i]。积分fInt1 可以工作并给出解决方案。但是，集成fInt2会产生限制为NA或NAN错误。我哪里可能出错了？这是 计算一个函数的积分，该函数是 R 中其他两个积分函数的乘积

set.seed(1234)
G <- 5# Suppose 5 groups
theta<-0.5
n_i <- 2 # There are two individuals per group
nTot <- n_i*G # In total we have 4 individuals 
z_ij <- rnorm(nTot, 0, 0.1)
ui<-rnorm(nTot, 0, 0.2)
T_ij <- runif(nTot, 0, 15)
Data <- round(data.frame(id = rep(1:nTot), group = rep(1:G, rep(2,G)), ui,z_ij, T_ij=round(T_ij,1)) , 3)
head(Data)
  id group     ui   z_ij T_ij
1  1     1 -0.095 -0.121  8.3
2  2     1 -0.200  0.028  9.7
3  3     2 -0.155  0.108  4.7
4  4     2  0.013 -0.235  9.3
5  5     3  0.192  0.043  4.9
6  6     3 -0.022  0.051  7.5

函数作为内部积分

fInt1 <- 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
      )
    )
  }})
}

fInt2 <- function(df) {
  Vectorize({function(y) {
    prod(
      sapply(
        seq_along(df),
        function(i) integrate(function(x) x*df$ui[i]*y*exp(x + y + df$z_ij[i]), 0, df$T_ij[i])$value
      )
    )
  }})
}

GroupInt1 <- sapply(1:G, function(grp) integrate(fInt1(subset(Data, group == grp, select = c("z_ij","T_ij"))), -5, 5)$value)
GroupInt1
[1] [1] 8.579064e+14 7.361849e+12 1.529633e+12 4.659699e+14 2.230921e+13

函数产生错误

GroupInt2 <- sapply(1:G, function(grp) integrate(fInt2(subset(Data, group == grp, select = c("z_ij","ui", "T_ij"))), -5, 5)$value)
Error in integrate(function(x) x * df$ui[i] * y * exp(x + y + df$z_ij[i]), : a limit is NA or NaN

Answer 1

错误是sapply中的索引。它应该使用1：nrow（df），而不是seq_along（df），是列的。

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

fInt1 <- function(df) {
  Vectorize({function(y) {
    prod(
      sapply(
        1:nrow(df),
        function(i) y*exp(y + df$z_ij[i])*integrate(function(x) x*exp(x), 0, df$T_ij[i])$value
      )
    )
  }})
}

fInt2 <- function(df) {
  Vectorize({function(y) {
    prod(
      sapply(
        1:nrow(df),
        function(i) df$ui[i]*y*exp(y + df$z_ij[i])*integrate(function(x) x*exp(x), 0, df$T_ij[i])$value
      )
    )
  }})
}

GroupInt1 <- sapply(1:G, function(grp) integrate(fInt1(subset(Data, group == grp, select = c("z_ij","T_ij"))), -5, 5)$value)
GroupInt1
#> [1] 8.579064e+14 7.361849e+12 1.529633e+12 4.659699e+14 2.230921e+13

GroupInt2 <- sapply(1:G, function(grp) integrate(fInt2(subset(Data, group == grp, select = c("z_ij", "ui", "T_ij"))), -5, 5)$value)
GroupInt2
#> [1]  1.630022e+13 -1.483413e+10 -6.461171e+09  8.650265e+12 -1.799484e+12