R: Интеграционный частный производный тривариат гауссова копула

Question

У меня есть небольшая проблема, пытаясь вычислить код ниже, используя R.

Немного контекста: я пытаюсь вычислить частную производную второго порядка, которая получена из трехвариантной гауссовой связки с корреляционной матрицей Sigma. (Извините, я не могу найти, как включить код LaTeX, чтобы помочь вам понять, какую формулу я точно вычисляю, в LaTeX это

\int_{-\infty}^{\Phi^{-1}(u_1)} \frac{1}{(2\pi)^{\frac{3}{2}}|\Sigma|^{\frac{1}{2}}} \exp \Big(-\frac{1}{2}(w_1 \, \Phi^{-1}(u_2) \, \Phi^{-1}(u_3)) \ Сигма ^{-1} (w_1 \,\Phi^{-1}(u_2) \, \Phi^{-1}(u_3))^t \Big) \mathrm dw_1 \, \frac{1}{\phi(\Phi^{-1}(u_2))} \, \frac{1}{\phi(\Phi^{-1}(u_3))})

Поскольку я являюсь относительно новым пользователем R, мой код для вычисления частной производной, учитывая значения u1, u2, u3 и матрицу корреляции, может быть неоптимальным, но то, что я написал, ниже (мне пришлось самому выполнять операции с матрицами "msgstr "по какой-то причине я не понял, что это не будет вычислять интеграцию, когда внутренняя функция была написана в матричной форме).

Моя проблема заключается в следующем: при использовании примера я использовал u1 = 0,1695729, u2 = 0,7357533, u3 = 0,4385701 вместе с Sigma <- c (0,3,0,5,0,4). По той причине, по которой я не понимаю, вычисление интеграла кажется невозможным, хотя при более тщательной проверке я построил внутреннюю функцию для различных значений u1 и получил функцию, которая кажется довольно "простой" для интеграции.

Кто-нибудь есть идеи, почему интеграция в неосуществимой? Любая помощь будет очень кстати!

gaus.der.part.u2u3<-function(u1,u2,u3,sigma){
  gaus.triv.2<-function(s,u2,u3,sigma){
    sigma.mat<-matrix(c(1,sigma[1],sigma[2],sigma[1],1,sigma[3],sigma[2],sigma[3],1),ncol=3)
    vec.1 <- matrix(c(s,qnorm(u2),qnorm(u3)),nrow=1)
    mat.inv <- solve(sigma.mat)
    vec.2 <- c(s*mat.inv[1,1]+qnorm(u2)*mat.inv[2,1]+qnorm(u3)*mat.inv[3,1],
               s*mat.inv[1,2]+qnorm(u2)*mat.inv[2,2]+qnorm(u3)*mat.inv[3,2],
               s*mat.inv[1,3]+qnorm(u2)*mat.inv[2,3]+qnorm(u3)*mat.inv[3,3])
    vec.3 <- vec.2[1]*s + vec.2[2]*qnorm(u2) + vec.2[3]*qnorm(u3)
    res <- 1/((2*pi)^(1.5)*abs(det(sigma.mat))^0.5)*exp(-0.5*vec.3)*1/dnorm(qnorm(u2))*1/dnorm(qnorm(u3))
    return(as.numeric(res))
  }  
  res.u2u3<-integrate(gaus.triv.2,-Inf,qnorm(u1), u2=u2, u3=u3, sigma=sigma,stop.on.error = FALSE)$value
  return(res.u2u3)
}