Условная вероятность на основе копулы - PullRequest
Новая
       14

Условная вероятность на основе копулы

0 голосов
/ 11 мая 2018

Чтение статьи «Многомерная модель морских штормов с использованием копул» (De Michele et al., 2007) https://www.sciencedirect.com/science/article/pii/S0378383907000592 Я застрял на вычислении частных производных в R.

Математический фон:

Если у вас есть три переменные H, D, I (с U1 = F (H) и т. Д.), И вам необходимо оценить условную вероятность

P (U3 | U1, U2)

сначала вы должны оценить соотношение частных производных, таких как

(∂C (u1, u2, u3) / ∂u1 ∂u2) / (∂C (u1, u2) / 1u1 ∂u2)

(подробнее см. На изображении), According to De Michele et al., 2007

код:

После пакета R "VineCopula" (https://cran.r -project.org / web / packages / VineCopula / VineCopula.pdf ) функция dduCopula вычисляет частную производную

∂C (u1, u2) / ∂u1

, поэтому эта процедура проста для простых частных производных (на рисунке они обозначены k, m и n).

Но как я могу оценить частную производную, которая относится к другой переменной "u1" по сравнению с связкой C (k, m).

∂C (к, м) / 1u1?

Для частных производных k и n: 

library(VineCopula)
u1<-pobs(H)
u2<-pobs(D) 
library(copula)
C_hd <- BB1Copula()
k <- ddvCopula(cbind(u1,u2), C_hd)
n <- dduCopula(cbind(u1,u2), C_hd)

У вас есть идеи?

1 Ответ

0 голосов
/ 14 мая 2018

Принцип построения, предложенный De Michele et al. (2017) в случае тривариата широко известен как лоза связка или конструкция парная связка . Если вы хотите работать с такими моделями, я предлагаю сначала прочитать о них немного больше. Хорошей отправной точкой является статья Aas et al. (2009) и вы найдете много других на http://www.vine -copula.org или google scholar .

Теперь к вашему вопросу: по правилу цепочки получаем это и, следовательно, , что .

Пример в R: 

# required package
library(VineCopula)

## simulate data
H <- rnorm(100)
D <- rnorm(100)
I <- rnorm(100)

## probability integral transforms
u_H <- pobs(H)
u_D <- pobs(D) 
u_I <- pobs(I)

## define dummy pair-copulas
C_HD <- BB1Copula()
C_DI <- BB1Copula()
C_HIgivenD <- BB1Copula()  # this is conditional dependence!

## calculate the conditional probability
k <- ddvCopula(cbind(u_H, u_D), C_HD)
m <- dduCopula(cbind(u_D, u_I), C_DI)
cond_prob <- dduCopula(cbind(k, m), C_HIgivenD)
...