我试图优化我的代码,它以四个向量X1,X2,X3,X4和一个相关矩阵R为参数。首先,对于每个观察I计算F(X3,X4)+F(X1,X2)-F(X2,X3)-F(X1,X4),其中F对应于以下函数:
F <- function(u,v,R){
pmvnorm(upper=c(qnorm(u,0,1), qnorm(v,0,1)),mean=c(0,0),lower=c(-Inf, -Inf),corr = R, sigma=NULL, algorithm = mvtnorm::GenzBretz(), keepAttr=FALSE))
}
然后,我将所有观察结果的对数(以某个delta为阈值)相加,并返回总和的相反值。
我希望我的代码更快,使用mapply
与使用for
循环相比没有显著差异。
目前,我正在使用下面的代码:
L_n <- function(R, X1, X2, X3, X4){
delta = 10**-9
mysummands <- mapply(F, X3, X4, MoreArgs=list(R=R)) + mapply(F, X1, X2, MoreArgs=list(R=R)) - mapply(F, X2, X3, MoreArgs=list(R=R)) - mapply(F, X1, X4, MoreArgs=list(R=R))
L <- sum(log(mapply(max,mysummands,MoreArgs=list(delta))))
return(-L)
}
如果我尝试通过指定u和v的向量来使用F作为向量化函数,我会得到一个错误,即相关矩阵的上界和对角线的大小不相同。(相关矩阵R对于所有观测保持相同)。
是否存在pmvnorm的矢量化等价物?
1条答案
按热度按时间3phpmpom1#
mnormt包中的
pmnorm
向量化: