R 中的积分误差:极限是 NA 或 NaN
Integral error in R : a limit is NA or NaN
fInt1
和fInt2
下面两个函数的区别在于增加了乘法项df$ui[i]
。集成 fInt1
有效并给出了解决方案。但是,积分 fInt2
得到 limit is NA or NAN error
。我哪里会出错?这是
的后续问题
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
错误是 sapply
中的索引。它应该使用 1:nrow(df)
而不是 seq_along(df)
,即 column-wise.
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
fInt1
和fInt2
下面两个函数的区别在于增加了乘法项df$ui[i]
。集成 fInt1
有效并给出了解决方案。但是,积分 fInt2
得到 limit is NA or NAN error
。我哪里会出错?这是
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
错误是 sapply
中的索引。它应该使用 1:nrow(df)
而不是 seq_along(df)
,即 column-wise.
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