Блог

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

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

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

Каждый раз, когда я пытался использовать пакет gnormt, команде pmnorm() требуется форма: pmnorm(данные, среднее значение, 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
  

Спасибо за вашу помощь!

Похоже, все, что вам нужно, это 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? Это не сильно отличалось бы.