Общая сумма гамма-случайных величин

35

Я читал, что сумма гамма-случайных величин с тем же параметром масштаба является еще одной гамма-случайной величиной. Я также видел статью Moschopoulos, описывающую метод суммирования общего набора гамма-случайных величин. Я пытался реализовать метод Мосхопулоса, но пока не добился успеха.

Как выглядит суммирование общего набора гамма-случайных величин? Чтобы конкретизировать этот вопрос, как он выглядит:

Gamma(3,1)+Gamma(4,2)+Gamma(5,1)

Если приведенные выше параметры не особенно показательны, пожалуйста, предложите другие.

OSE
источник
4
Явное решение для суммы любых двух гамма-распределений было размещено на stats.stackexchange.com/a/252192 .
whuber
Особый пример этого, где все гамма-распределения имеют параметр формы 1 (то есть они экспоненциальные), называется гипоэкспоненциальным распределением (семейством) . Для случая только двух экспоненциальных распределений также есть явная формула, приведенная на stats.stackexchange.com/questions/412849 .
uber

Ответы:

37

Сначала объедините любые суммы, имеющие один и тот же масштабный коэффициент : переменная Γ(N,β) плюс a Γ(м,β) образуют переменную Γ(N+м,β) .

Далее заметим, что характеристическая функция (cf) для Γ(N,β) равна (1-яβT)-N , откуда cf суммы этих распределений является произведением

ΠJ1(1-яβJT)NJ,

Когда весь интеграл, этот продукт расширяется в виде частичной доли в линейную комбинацию из ( 1 - я β J т ) - ν , где ν целых числа от 1 и п J . В примере с β 1 = 1 , n 1 = 8 (из суммы Γ ( 3 , 1 ) и Γ ( 5 , 1)NJ (1-яβJT)-νν1NJβ1знак равно1,N1знак равно8Γ(3,1) ) и β 2 = 2 , n 2 = 4 находимΓ(5,1)β2знак равно2,N2знак равно4

1(1it)81(12it)4=1(x+i)88i(x+i)740(x+i)6+160i(x+i)5+560(x+i)41792i(x+i)35376(x+i)2+15360ix+i+256(2x+i)4+2048i(2x+i)39216(2x+i)230720i2x+i.

Обратное взятие cf является обратным преобразованием Фурье, которое является линейным : это означает, что мы можем применять его термин за термином. Каждый член распознается как кратное cf гамма-распределения и поэтому легко инвертируется для получения PDF . В примере мы получаем

е-TT75040+190е-TT6+13е-TT5+203е-TT4+83е-T2T3+2803е-TT3-128е-T2T2++896е-TT2+2304е-T2T+5376е-TT-15360е-T2+15360е-T

для PDF суммы.

Это конечная смесь гамма-распределений с масштабными коэффициентами, равными коэффициентам в сумме, и коэффициентами формы, меньшими или равными коэффициентам в сумме. За исключением особых случаев (где может произойти некоторое аннулирование), число слагаемых задается параметром общей формы (при условии, что все n j различны).N1+N2+NJ


В качестве теста приведем гистограмму из результатов, полученных путем добавления независимых отрисовок из распределений Γ ( 8 , 1 ) и Γ ( 4 , 2 ) . На него накладывается график, в 10 4 раза превышающий предыдущую функцию. Подгонка очень хорошая.104Γ(8,1)Γ(4,2)104

фигура


Мошопулос продвигает эту идею на один шаг вперед, расширяя cf суммы в бесконечный ряд гамма-характеристических функций всякий раз, когда один или несколько из нецелочислен, а затем завершает бесконечный ряд в точке, где он достаточно хорошо аппроксимируется. ,Nя

Whuber
источник
2
Незначительное комментарий: Как правило, конечное смесь означает PDF вида где я > 0 и Е я в я = 1 , то есть в виде I являются вероятностями, и PDF может быть интерпретирован как (закон полной вероятности) взвешенной суммы условных PDF, учитывая различные условия, которые возникают с вероятностями a i
f(x)=i=1naifi(x)
ai>0iai=1aiaя, Однако в приведенной выше сумме некоторые коэффициенты являются отрицательными, и поэтому стандартная интерпретация смеси не применяется.
Дилип Сарватэ
@Dilip Это хороший момент. Что делает этот случай интересным, так это то, что хотя некоторые коэффициенты могут быть отрицательными, тем не менее эта комбинация все еще является действительным распределением (по самой своей конструкции).
whuber
Можно ли расширить этот подход для учета добавления зависимых переменных? В частности, я хочу добавить 6 дистрибутивов, каждый из которых имеет некоторую корреляцию с другими.
masher
11

Я покажу другое возможное решение, которое довольно широко применимо, и с сегодняшним программным обеспечением R, довольно простое в реализации. Это приближение плотности седловой точки, которое должно быть более широко известным!

Для терминологии о гамма-распределении я буду следовать https://en.wikipedia.org/wiki/Gamma_distribution с параметризацией формы / масштаба, - параметр формы, а θ - масштаб. Для приближения седловой точки я буду следовать Рональду В. Батлеру: «Приближения седловой точки с приложениями» (Кембридж UP). Приближение седловой точки объясняется здесь: Как работает приближение седловой точки? здесь я покажу, как это используется в этом приложении.Кθ

Пусть - случайная величина с существующей порождающей момент функцией M ( s ) = E e s X, которая должна существовать для s в некотором открытом интервале, который содержит ноль. Затем определим производящую функцию кумулянта как K ( s ) = log M ( s ). Известно, что E X = K ( 0 ) , Var ( X ) = K ( 0 )Икс

M(s)знак равноЕеsИкс
s
К(s)знак равножурналM(s)
ЕИксзнак равноК'(0),Var(Икс)знак равноК"(0), Уравнение является перевал , который определяет неявный ы как функции х (которое должно быть в диапазоне X ). Обозначим эту функцию неявно определенную как з ( х ) . Обратите внимание, что уравнение седловой точки всегда имеет ровно одно решение, потому что кумулянтная функция является выпуклой.
K(s^)=x
sxXs^(x)

Тогда перевала приближение к плотности из X задается ф ( х ) = 1fX

f^(x)=12πK(s^)exp(K(s^)s^x)

X1,X2,,XnXi(ki,θi)

K(s)=i=1nkiln(1θis)
s<1/max(θ1,θ2,,θn)
K(s)=i=1nkiθi1θis
K(s)=i=1nkiθi2(1θis)2.
Rn=3Кзнак равно(1,2,3)θзнак равно(1,2,3)R
shape <- 1:3 #ki
scale <- 1:3 # thetai
# For this case,  we get expectation=14,  variance=36
make_cumgenfun  <-  function(shape, scale) {
      # we return list(shape, scale, K, K', K'')
      n  <-  length(shape)
      m <-   length(scale)
      stopifnot( n == m, shape > 0, scale > 0 )
      return( list( shape=shape,  scale=scale, 
                    Vectorize(function(s) {-sum(shape * log(1-scale * s) ) }),
                    Vectorize(function(s) {sum((shape*scale)/(1-s*scale))}) ,
                    Vectorize(function(s) { sum(shape*scale*scale/(1-s*scale)) }))    )
}

solve_speq  <-  function(x, cumgenfun) {
          # Returns saddle point!
          shape <- cumgenfun[[1]]
          scale <- cumgenfun[[2]]
          Kd  <-   cumgenfun[[4]]
          uniroot(function(s) Kd(s)-x,lower=-100,
                  upper = 0.3333, 
                  extendInt = "upX")$root
}

make_fhat <-  function(shape,  scale) {
    cgf1  <-  make_cumgenfun(shape, scale)
    K  <-  cgf1[[3]]
    Kd <-  cgf1[[4]]
    Kdd <- cgf1[[5]]
    # Function finding fhat for one specific x:
    fhat0  <- function(x) {
        # Solve saddlepoint equation:
        s  <-  solve_speq(x, cgf1)
        # Calculating saddlepoint density value:
        (1/sqrt(2*pi*Kdd(s)))*exp(K(s)-s*x)
    }
    # Returning a vectorized version:
    return(Vectorize(fhat0))
} #end make_fhat

 fhat  <-  make_fhat(shape, scale)
plot(fhat, from=0.01,  to=40, col="red", main="unnormalized saddlepoint approximation\nto sum of three gamma variables")

в результате на следующем участке: введите описание изображения здесь

Я оставлю нормализованное приближение седловой точки в качестве упражнения.

Къетил б Халворсен
источник
1
Это интересно, но я не могу заставить ваш Rкод работать, чтобы сравнить приближение с точным ответом. Любая попытка вызвать fhatприводит к ошибкам, по-видимому, при использовании uniroot.
whuber
3
Какая у тебя версия R? В кодах для uniroot используется новый аргумент extensionInt, который был представлен в версии 3.1 R. Если ваш R старше, вы можете попытаться удалить его (и увеличить интервал, заданный uniroot). Но это сделает код менее надежным!
kjetil b halvorsen
10

Уравнение Уэлча – Саттервейта может быть использовано для получения приблизительного ответа в форме гамма-распределения. Это имеет приятное свойство, позволяющее нам рассматривать гамма-распределения как (приблизительно) закрытые при добавлении. Это приближение в обычно используемом t-тесте Уэлча.

(Гамма-распределение можно рассматривать как масштабированное распределение хи-квадрат и допускает нецелочисленный параметр формы.)

К,θ

КsUмзнак равно(ΣяθяКя)2Σяθя2Кя

θsUмзнак равноΣθяКяКsUм

Кзнак равно(3,4,5)θзнак равно(1,2,1)

Таким образом, мы получаем примерно гамма (10,666 ..., 1,5)

Кθяθ

Пол Харрисон
источник
6

N

гDС(a,б,α,β;τ)знак равно{бaβαΓ(a+α)е-бττa+α-11F1[α,a+α,(б-β)τ],τ>00,τ0,
гaммa(a,б)Γ(a,1/б)бβ
деревенщина
источник