квантили сумм, использующие распределения копулы, слишком медленные - PullRequest
Новая
       6

квантили сумм, использующие распределения копулы, слишком медленные

6 голосов
/ 12 ноября 2011

Попытка создать таблицу для квантилей суммы двух зависимых случайных величин, используя встроенные распределения связок (Clayton, Frank, Gumbel) с бета-маргиналами. Пробовал NProbability и FindRoot различными методами - недостаточно быстро. Пример комбинаций копула-маргинал, которые мне нужно исследовать, следующий: 

nProbClayton[t_?NumericQ, c_?NumericQ] := 
        NProbability[  x + y <= t, {x, y}  \[Distributed]    
               CopulaDistribution[{"Clayton", c}, {BetaDistribution[8, 2], 
                                                   BetaDistribution[8, 2]}]]

Для однократной оценки числовой вероятности используется 

nProbClayton[1.9, 1/10] // Timing // Quiet

Я получаю 

{4.914, 0.939718}

на компьютере с 64-битной версией Core2 Duo T9600 2,80 ГГц (MMA 8.0.4)

Чтобы получить квантиль суммы, используйте 

FindRoot[nProbClayton[q, 1/10] == 1/100, {q, 1, 0, 2}// Timing // Quiet

различными методами 

( `Method -> Automatic`, `Method -> "Brent"`, `Method -> "Secant"` )

требуется около минуты, чтобы найти один квантиль: время равно 

{48.781, {q -> 0.918646}}
{50.045, {q -> 0.918646}}
{65.396, {q -> 0.918646}}

Для других комбинаций копула-маргинальное время немного лучше.

Потребность: любые хитрости / методы для улучшения времени.

1 Ответ

7 голосов
/ 12 ноября 2011

CDF копулы Клейтона-Парето с параметром c можно рассчитать в соответствии с 

cdf[c_] := Module[{c1 = CDF[BetaDistribution[8, 2]]}, 
   (c1[#1]^(-1/c) + c1[#2]^(-1/c) - 1)^(-c) &]

Тогда cdf[c][t1,t2] - это вероятность того, что x<=t1 и y<=t2. Это означает, что вы можете рассчитать вероятность того, что x+y<=t в соответствии с 

prob[t_?NumericQ, c_?NumericQ] := 
   NIntegrate[Derivative[1, 0][cdf[c]][x, t - x], {x, 0, t}]

Время, которое я получаю на своей машине, 

prob[1.9, .1] // Timing

(* ==> {0.087518, 0.939825} *)

Обратите внимание, что я получаю другое значение для вероятности, чем в оригинальном сообщении. Однако выполнение nProbClayton[1.9,0.1] выдает предупреждение о медленной конвергенции, что может означать, что результат в исходном сообщении отключен. Кроме того, если я изменю x+y<=t на x+y>t в исходном определении nProbClayton и вычислю 1-nProbClayton[1.9,0.1], я получу 0.939825 (без предупреждений), что соответствует результату, описанному выше.

За квантиль суммы я получаю 

FindRoot[prob[q, .1] == .01, {q, 1, 0, 2}] // Timing

(* ==> {1.19123, {q -> 0.912486}} *)

Опять же, я получаю результат, отличный от того, что был в оригинальном сообщении, но аналогично предыдущему, меняя x+y<=t на x+y>t и вычисляя FindRoot[nProbClayton[q, 1/10] == 1-1/100, {q, 1, 0, 2}], возвращает то же значение для q, что и выше.

...