不确定是什么导致了 R 中的这个错误(pmvt - 包:mvtnorm)?
Unsure what's causing this error in R (pmvt - Package: mvtnorm)?
我有一个简单的冒险函数,导致错误的行被标记了。
h <- function(t,u) {
x <- 1 - Sa(t)
y <- 1 - Sm(u)
invx <- as.numeric(qt(x,df=d1))
invy <- as.numeric(qt(x,df=d1))
[ERROR LINE] copula <- pmvt(lower=as.numeric(cbind(-9999,-9999)),upper=cbind(invx,invy),df=d1,corr=matrix(cbind(1,d2,d2,1),byrow=T,ncol=2) )
density <- dmvt(cbind(invx,invy),sigma=matrix(cbind(1,d2,d2,1),byrow=T,ncol=2),df=d1)
num <- (sa(t)*sm(u))*density/dt(invx,df=d1)/dt(invy,df=d1)
den <- 1 - x - y + copula
hazard <- num/den
return(hazard)
}
然后由似然函数调用此风险函数:
# log Likelihood function for each individual car i
lli <- function(data) {
result <- 0;
# for all claims, evaluate hazard function at that point
if (nrow(data)> 2) {
for (k in 1:nrow(data)) {
if (data[k,3] == 1) {
result <- result + log(h(data[k,2],data[k,1]));
}
}
}
# integrate hazard function over areas between claims
for (k in 1:(nrow(data)-1)) {
integral <- quad2d(h,data[k,2],data[k+1,2],data[k,1],data[k+1,1]);
result <- result - integral;
}
return(result)
}
现在这个似然函数被第三个函数调用,用于我的整个数据集; 但是导致错误的是上面的函数,而不是下面的函数
# log Likelihood function over all vehicles
ll <- function(x) {
# Unpack parameters
d1 <<- x[1];
d2 <<- x[2];
total <- 0;
# Get log Likelihood for each vehicle
for (i in 1:length(alldata)) {
total <- total + lli(alldata[[i]]);
#print(sprintf("Found candidate solution %d value: %f",i,total));
}
#print(sprintf("Found candidate solution value: %f",total));
if (is.nan(total)) { #If it is undefined, make it a large negative number
total <- -2147483647 ;
}
return(-1*total); # Minimise instead of maximise
}
报错信息如下:
> ll(cbind(50,0.923))
Error in checkmvArgs(lower = lower, upper = upper, mean = delta, corr = corr, :
‘diag(corr)’ and ‘lower’ are of different length
我在使用 pmvnorm 时一直遇到同样的错误,最后不得不使用 pbivnorm 包来解决这个问题。不过,我找不到双变量 t 分布的替代包。我不明白问题出在哪里。当我自己调用函数 h(t,u) 时,它执行时没有问题。但是当 lli(data) 调用 h(t,u) 时,它不起作用。更奇怪的是,它们的长度是一样的。
> length(as.numeric(cbind(-9999,-9999)))
[1] 2
> length(diag(matrix(cbind(1,d2,d2,1),byrow=T,ncol=2)))
[1] 2
对于混乱的代码,我深表歉意。我不怎么用R。无论如何,这让我完全难住了。
数据文件在这里:https://files.fm/u/yx9pw2b3
我忘记包括的附加代码,基本上是一些常量和边际 CDF 函数:
Marginals.R:
p1 <- 0.4994485;
p2 <- 0.2344439;
p3 <- 0.1151654;
p4 <- 0.1509421;
b1 <- 0.7044292
t1 <- 1713.3170267
mu1 <- 7.014415
sig1 <- 1.394735
mu2 <- 6.926146
sig2 <- 1.056647
mu3 <- 6.7995896
sig3 <- 0.7212853
b2 <- 0.6444582
t2 <- 762.9962093
b3 <- 1.494303
t3 <- 410.828780
b1 <- 0.903
t1 <- 864.896
b2 <- 0.9109
t2 <- 314.2946
# Marginal survival distribution and density
Sa <- function(t) {return(exp(-(t / t1) ** b1))}
Sm <- function(u) {return(exp(-(u / t2) ** b2))}
sa <- function(t) {return((t / t1) ** b1 * b1 * exp(-(t / t1) ** b1) / t ) }
sm <- function(u) {return((u / t2) ** b2 * b2 * exp(-(u / t2) ** b2) / u ) }
总结:
问题是调用pvmt
时lower
和upper
的长度差,upper
的长度是2048,而lower
的长度是2.
推理:
1. pmvt
通过调用 mvtnorm
包中的 checkmvArgs
检查传入的参数。
2. 在checkmvArgs
中,lower
、upper
和mean
已经被rec <- cbind(lower, upper, mean)
放在一起了。这里的新数据 rec
有 2048 行而不是 2.
3. lower
然后被 lower <- rec[, "lower"]
替换,lower
现在的长度为 2048 而不是 2.
4.鉴于corr
仍然是一个2*2的矩阵,检查length(corr) != length(lower)
时出错
解法:
invx <- as.numeric(qt(x,df=d1))
invy <- as.numeric(qt(x,df=d1))
upper
表示长度为 2 的向量,因此 invx
和 invy
需要是单个数字。
由于不确定您要定义的上限是多少,我无法进一步解决。可能的一种是:
invx <- as.numeric(qt(x,df=d1))
invy <- as.numeric(qt(x,df=d1))
copula <- pmvt(lower=as.numeric(cbind(-9999,-9999)),upper=range(c(invx,invy)),df=d1,corr=matrix(c(1,d2,d2,1),byrow=T,ncol=2) )
它使用 invx
和 invy
的范围作为输入。因此 dmvt
不会受到影响。
注:
由于未提供值 a
,下面的下一行(调用 dmvt
)错误行失败。
编辑:
为了使问题更具体:
1. quad2d
将生成一个 Gauss-Legendre 正交,默认情况下将在给定范围内创建长度为 32 的正交。而且,
2. 然后使用此 Gauss-Legendre 积分中的 x 和 y 调用您的函数 h
。因此,h
中定义的t
和u
不是一个数字,而是一个向量。
我有一个简单的冒险函数,导致错误的行被标记了。
h <- function(t,u) {
x <- 1 - Sa(t)
y <- 1 - Sm(u)
invx <- as.numeric(qt(x,df=d1))
invy <- as.numeric(qt(x,df=d1))
[ERROR LINE] copula <- pmvt(lower=as.numeric(cbind(-9999,-9999)),upper=cbind(invx,invy),df=d1,corr=matrix(cbind(1,d2,d2,1),byrow=T,ncol=2) )
density <- dmvt(cbind(invx,invy),sigma=matrix(cbind(1,d2,d2,1),byrow=T,ncol=2),df=d1)
num <- (sa(t)*sm(u))*density/dt(invx,df=d1)/dt(invy,df=d1)
den <- 1 - x - y + copula
hazard <- num/den
return(hazard)
}
然后由似然函数调用此风险函数:
# log Likelihood function for each individual car i
lli <- function(data) {
result <- 0;
# for all claims, evaluate hazard function at that point
if (nrow(data)> 2) {
for (k in 1:nrow(data)) {
if (data[k,3] == 1) {
result <- result + log(h(data[k,2],data[k,1]));
}
}
}
# integrate hazard function over areas between claims
for (k in 1:(nrow(data)-1)) {
integral <- quad2d(h,data[k,2],data[k+1,2],data[k,1],data[k+1,1]);
result <- result - integral;
}
return(result)
}
现在这个似然函数被第三个函数调用,用于我的整个数据集; 但是导致错误的是上面的函数,而不是下面的函数
# log Likelihood function over all vehicles
ll <- function(x) {
# Unpack parameters
d1 <<- x[1];
d2 <<- x[2];
total <- 0;
# Get log Likelihood for each vehicle
for (i in 1:length(alldata)) {
total <- total + lli(alldata[[i]]);
#print(sprintf("Found candidate solution %d value: %f",i,total));
}
#print(sprintf("Found candidate solution value: %f",total));
if (is.nan(total)) { #If it is undefined, make it a large negative number
total <- -2147483647 ;
}
return(-1*total); # Minimise instead of maximise
}
报错信息如下:
> ll(cbind(50,0.923))
Error in checkmvArgs(lower = lower, upper = upper, mean = delta, corr = corr, :
‘diag(corr)’ and ‘lower’ are of different length
我在使用 pmvnorm 时一直遇到同样的错误,最后不得不使用 pbivnorm 包来解决这个问题。不过,我找不到双变量 t 分布的替代包。我不明白问题出在哪里。当我自己调用函数 h(t,u) 时,它执行时没有问题。但是当 lli(data) 调用 h(t,u) 时,它不起作用。更奇怪的是,它们的长度是一样的。
> length(as.numeric(cbind(-9999,-9999)))
[1] 2
> length(diag(matrix(cbind(1,d2,d2,1),byrow=T,ncol=2)))
[1] 2
对于混乱的代码,我深表歉意。我不怎么用R。无论如何,这让我完全难住了。
数据文件在这里:https://files.fm/u/yx9pw2b3
我忘记包括的附加代码,基本上是一些常量和边际 CDF 函数:
Marginals.R:
p1 <- 0.4994485;
p2 <- 0.2344439;
p3 <- 0.1151654;
p4 <- 0.1509421;
b1 <- 0.7044292
t1 <- 1713.3170267
mu1 <- 7.014415
sig1 <- 1.394735
mu2 <- 6.926146
sig2 <- 1.056647
mu3 <- 6.7995896
sig3 <- 0.7212853
b2 <- 0.6444582
t2 <- 762.9962093
b3 <- 1.494303
t3 <- 410.828780
b1 <- 0.903
t1 <- 864.896
b2 <- 0.9109
t2 <- 314.2946
# Marginal survival distribution and density
Sa <- function(t) {return(exp(-(t / t1) ** b1))}
Sm <- function(u) {return(exp(-(u / t2) ** b2))}
sa <- function(t) {return((t / t1) ** b1 * b1 * exp(-(t / t1) ** b1) / t ) }
sm <- function(u) {return((u / t2) ** b2 * b2 * exp(-(u / t2) ** b2) / u ) }
总结:
问题是调用pvmt
时lower
和upper
的长度差,upper
的长度是2048,而lower
的长度是2.
推理:
1. pmvt
通过调用 mvtnorm
包中的 checkmvArgs
检查传入的参数。
2. 在checkmvArgs
中,lower
、upper
和mean
已经被rec <- cbind(lower, upper, mean)
放在一起了。这里的新数据 rec
有 2048 行而不是 2.
3. lower
然后被 lower <- rec[, "lower"]
替换,lower
现在的长度为 2048 而不是 2.
4.鉴于corr
仍然是一个2*2的矩阵,检查length(corr) != length(lower)
解法:
invx <- as.numeric(qt(x,df=d1))
invy <- as.numeric(qt(x,df=d1))
upper
表示长度为 2 的向量,因此 invx
和 invy
需要是单个数字。
由于不确定您要定义的上限是多少,我无法进一步解决。可能的一种是:
invx <- as.numeric(qt(x,df=d1))
invy <- as.numeric(qt(x,df=d1))
copula <- pmvt(lower=as.numeric(cbind(-9999,-9999)),upper=range(c(invx,invy)),df=d1,corr=matrix(c(1,d2,d2,1),byrow=T,ncol=2) )
它使用 invx
和 invy
的范围作为输入。因此 dmvt
不会受到影响。
注:
由于未提供值 a
,下面的下一行(调用 dmvt
)错误行失败。
编辑:
为了使问题更具体:
1. quad2d
将生成一个 Gauss-Legendre 正交,默认情况下将在给定范围内创建长度为 32 的正交。而且,
2. 然后使用此 Gauss-Legendre 积分中的 x 和 y 调用您的函数 h
。因此,h
中定义的t
和u
不是一个数字,而是一个向量。