Программирование двумерного нормального CDF на R

Question

У меня есть вопрос, касающийся кодирования функции, которая содержит двумерный нормальный CDF в R. Функция, которую я пытаюсь кодировать, требует одного двумерного нормального CDF, который должен рассчитываться по-разному в зависимости от наблюдения.В частности, в зависимости от значения определенной переменной корреляция должна «переключаться» между положительным и отрицательным, но в вызове не должно быть различий.

Этот стиль функции был закодирован в LIMDEP, и я пытаюсь воспроизвести его, но не смог заставить его работать в R. Команда в LIMDEP для вычисления двумерного нормального CDF - "BVN (x1, x2, r) ", что явно требует двух переменных, используемых для расчета (x1, x2) и корреляции (r).LIMDEP использует 15-точечную квадратуру Гаусса-Лагерра для вычисления двумерного нормального CDF.

В R получается, что два пакета вычисляют многомерный нормальный CDF.Я пробовал пакет mnormt (хотя есть и пакет mvtnorm - я не вижу существенной разницы), который использует метод Genz, который кажется похожим, но более общим, чем используемый квадратурный метод Gauss-Laguerre 15в LIMDEP (ссылаясь на документы в? pmnorm).

Каждый раз, когда я пытался использовать пакет mnormt, команда pmnorm () требует форму: pmnorm (data, mean, varcov), которую я не смог закодировать для корреляционного переключения.

Любые идеи, как заставить это работать ??

Вот пример некоторого тривиального кода, объясняющего, что я говорю о том, что я хотел бы сделать (кроме как внутри функции безцикл):

 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?Это не будет сильно отличаться.