在 R 中编程二元正态 CDF

Question

我有一个关于在 R 中编写包含二元正态 CDF 的函数的问题。我尝试编写的函数需要一个二元正态 CDF，应根据观察结果进行不同的计算。具体来说，根据某个变量的值，相关性应该在正负之间“切换”，但调用中应该没有差异。

这种风格的函数已在 LIMDEP 中编码，我正在尝试复制它，但无法让它在 R 中工作。LIMDEP 中计算二元正态 CDF 的命令是“BVN(x1, x2 , r)"，它明确需要用于计算的两个变量 (x1, x2) 和相关性 (r)。 LIMDEP 使用 Gauss-Laguerre 15 点求积来计算二元正态 CDF。

在 R 中，似乎有两个包计算多元正态 CDF。我一直在尝试使用 Genz 方法的 mnormt 包（虽然也有 mvtnorm 包——我没有看到主要区别），这似乎是相似的，但比使用的 Gauss-Laguerre 15 正交方法更通用在 LIMDEP 中（参考 ?pmnorm 下的论文）。

每次我尝试使用 mnormt 包时，命令 pmnorm() 都需要以下形式：pmnorm(data, mean, varcov)，我无法为相关切换编写代码。

任何想法如何让它工作？？

下面是一些简单代码的示例，用于解释我在说什么我想做的事情（除了没有 for 循环的函数内部）：

 library(mnormt)
 A <- c(0,1, 1, 1, 0, 1, 0, 1, 0, 1)
 q <- 2*A-1 
 set.seed(1234)
 x <- rnorm(10)
 y <- rnorm(10, 2, 2)

 #Need to return a value of the CDF for each row of data:
 cdf.results <- 0
 for(i in 1:length(A)){
 vc.mat <- matrix(c(1, q[i]*.7, q[i]*.7, 1.3), 2, 2)
 cdf.results[i] <- pmnorm(cbind(x[i], y[i]), c(0, 0), vc.mat) 
 }
 cdf.results

感谢您的帮助！

Answer 1

听起来您所需要的只是 1) 将您的脚本变成一个函数，以便它可以应用于任意 x、y 和 q 以及 2) 摆脱 for 循环。如果是这种情况，?function 和 ?apply 应该可以满足您的需求。

BVN=function(x,y,q) {

  cdf.results=apply(cbind(x,y,q),1,FUN=function(X) 
       {
      x=X[1]
      y=X[2]
      q=X[3]
          vc.mat <- matrix(c(1, q*.7, q*.7, 1.3), 2, 2)
          pmnorm(cbind(x, y), c(0, 0), vc.mat)     
                    # I think you may want c(0,2) but not sure
        }
                   )
 cdf.results
}



BVN(x,y,q)

这里的 x,y 和 q 是你上面写的。也许您希望函数采用矩阵 r，就像在 limdep 中一样？不会有太大的不同。