Jumlah generik dari variabel acak Gamma

35

Saya telah membaca bahwa jumlah variabel acak Gamma dengan parameter skala yang sama adalah variabel acak Gamma lainnya. Saya juga telah melihat makalah oleh Moschopoulos yang menggambarkan metode untuk penjumlahan set umum variabel acak Gamma. Saya telah mencoba menerapkan metode Moschopoulos tetapi belum berhasil.

Seperti apa penjumlahan set umum variabel acak Gamma? Untuk membuat pertanyaan ini konkret, seperti apa bentuknya:

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

Jika parameter di atas tidak terlalu terbuka, harap sarankan orang lain.

OSE
sumber
4
Solusi eksplisit untuk jumlah dari setiap dua distribusi Gamma telah diposting di stats.stackexchange.com/a/252192 .
Whuber
Contoh khusus ini, di mana semua distribusi Gamma memiliki parameter bentuk 1 (yaitu, mereka eksponensial) disebut distribusi hypoexponential (keluarga) . Untuk kasus hanya dua distribusi eksponensial ada juga rumus eksplisit yang diberikan di stats.stackexchange.com/questions/412849 .
whuber

Jawaban:

37

Pertama, gabungkan setiap jumlah yang memiliki faktor skala yang sama : a ditambah variasi varian .Γ(n,β)Γ(m,β)Γ(n+m,β)

Selanjutnya, amati bahwa fungsi karakteristik (cf) dari adalah , di mana cf dari sejumlah distribusi ini adalah produkΓ(n,β)(1iβt)n

j1(1iβjt)nj.

Ketika semuanya integral, produk ini diperluas sebagai sebagian parsial menjadi kombinasi linear dari mana adalah bilangan bulat antara dan . Dalam contoh dengan (dari jumlah dan ) dan kita temukan ( 1 - i β j t ) - ν ν 1 n j β 1 = 1 , n 1 = 8 Γ ( 3 , 1 ) Γ ( 5 , 1 ) β 2 = 2 , n 2 = 4nj (1iβjt)νν1njβ1=1,n1=8Γ(3,1)Γ(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.

Kebalikan dari mengambil cf adalah invers Fourier Transform, yang linear : itu berarti kita dapat menerapkannya istilah demi istilah. Setiap istilah dikenali sebagai kelipatan cf dari distribusi Gamma dan karenanya mudah dibalik untuk menghasilkan PDF . Dalam contoh yang kita dapatkan

ett75040+190ett6+13ett5+203ett4+83et2t3+2803ett3128et2t2+896ett2+2304et2t+5376ett15360et2+15360et

untuk PDF dari jumlah tersebut.

Ini adalah campuran terbatas distribusi Gamma yang memiliki faktor skala sama dengan mereka yang ada dalam jumlah dan faktor bentuk kurang dari atau sama dengan mereka yang ada dalam jumlah. Kecuali dalam kasus khusus (di mana beberapa pembatalan mungkin terjadi), jumlah istilah diberikan oleh parameter bentuk total (dengan asumsi semua berbeda).n1+n2+nj


Sebagai tes, berikut adalah histogram hasil diperoleh dengan menambahkan undian independen dari distribusi dan . Di atasnya ditumpangkan grafik kali fungsi sebelumnya. Cocok sangat bagus.104Γ(8,1)Γ(4,2)104

Angka


Moschopoulos membawa ide ini selangkah lebih maju dengan memperluas cf dari jumlah menjadi deret fungsi fungsi Gamma yang tak terbatas setiap kali satu atau lebih dari adalah non-integral, dan kemudian mengakhiri deret infinite pada titik di mana ia cukup didekati dengan baik.ni

whuber
sumber
2
Minor komentar: Biasanya, campuran yang terbatas berarti pdf dari bentuk di mana sebuah i > 0 dan Σ i a i = 1 , yaitu, sebuah i adalah probabilitas dan pdf dapat diartikan sebagai (hukum probabilitas total) jumlah tertimbang bersyarat PDF diberikan berbagai kondisi yang terjadi dengan probabilitas sebuah i
f(x)=i=1naifi(x)
ai>0iai=1aiai. Namun, dalam jumlah di atas, beberapa koefisien negatif dan dengan demikian interpretasi standar dari campuran tidak berlaku.
Dilip Sarwate
@ Philip Itu poin yang bagus. Apa yang membuat kasus ini menarik adalah bahwa walaupun beberapa koefisien mungkin negatif, namun kombinasi ini masih merupakan distribusi yang valid (berdasarkan konstruksinya).
Whuber
Bisakah pendekatan ini diperluas untuk memperhitungkan penambahan variabel dependen? Secara khusus, saya ingin menambahkan hingga 6 distribusi dengan masing-masing memiliki beberapa korelasi dengan yang lain.
hidung belang
11

Saya akan menunjukkan solusi lain yang mungkin, yang cukup banyak diterapkan, dan dengan perangkat lunak R saat ini, cukup mudah untuk diterapkan. Itu adalah perkiraan kepadatan saddlepoint, yang seharusnya diketahui lebih luas!

Untuk terminologi tentang distribusi gamma, saya akan mengikuti https://en.wikipedia.org/wiki/Gamma_distribution dengan parametrization bentuk / skala, adalah parameter bentuk dan θ adalah skala. Untuk perkiraan saddlepoint saya akan mengikuti Ronald W Butler: "perkiraan Saddlepoint dengan aplikasi" (Cambridge UP). Perkiraan saddlepoint dijelaskan di sini: Bagaimana cara pendekatan saddlepoint bekerja? di sini saya akan menunjukkan bagaimana ini digunakan dalam aplikasi ini.kθ

Misalkan adalah variabel acak dengan fungsi penghasil momen yang ada M ( s ) = E e s X yang harus ada untuk s dalam beberapa interval terbuka yang berisi nol. Kemudian tentukan fungsi penghasil kumulant oleh K ( s ) = log M ( s ). Diketahui bahwa E X = K ( 0 ) , Var ( X ) = K ( 0 )X

M(s)=EesX
s
K(s)=logM(s)
EX=K(0),Var(X)=K(0). Persamaan saddlepoint adalah yang mendefinisikan secara implisit s sebagai fungsi dari x (yang harus dalam kisaran X ). Kami menulis fungsi ini didefinisikan secara implisit sebagai s ( x ) . Perhatikan bahwa persamaan saddlepoint selalu memiliki tepat satu solusi, karena fungsi kumulans adalah cembung.
K(s^)=x
sxXs^(x)

Kemudian saddlepoint pendekatan kepadatan dari X diberikan oleh f ( x ) = 1fX fungsi kepadatan perkiraan ini tidak dijamin untuk mengintegrasikan ke 1, sehingga adalah saddlepoint pendekatan unnormalized. Kita bisa mengintegrasikannya secara numerik dan renormalisasi untuk mendapatkan perkiraan yang lebih baik. Tetapi perkiraan ini dijamin tidak negatif.

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

Sekarang mari menjadi variabel acak gamma independen, di mana X i memiliki distribusi dengan parameter ( k i , θ i ) . Maka fungsi penghasil kumulans adalah K ( s ) = - n i = 1 k i ln ( 1 - θ i s ) yang ditentukan untuk s < 1 / maks (X1,X2,,XnXi(ki,θi)

K(s)=i=1nkiln(1θis)
. Derivatif pertama adalah K ( s ) = n i = 1 k i θ is<1/max(θ1,θ2,,θn) dan turunan kedua adalah K(s)= n i=1kiθ 2 i
K(s)=saya=1nksayaθsaya1-θsayas
Berikut ini saya akan memberikan beberapakode menghitung ini, dan akan menggunakan nilai parametern=3,k=(1,2,3),θ=(1,2,3). Perhatikan bahwakodeberikutmenggunakan argumen baru dalam fungsi uniroot yang diperkenalkan di R 3.1, jadi tidak akan berjalan di R yang lebih lama.
K(s)=saya=1nksayaθsaya2(1-θsayas)2.
Rn=3k=(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")

menghasilkan plot berikut: masukkan deskripsi gambar di sini

Saya akan meninggalkan pendekatan saddlepoint yang dinormalisasi sebagai latihan.

kjetil b halvorsen
sumber
1
Ini menarik, tetapi saya tidak dapat membuat Rkode Anda berfungsi untuk membandingkan perkiraan dengan jawaban yang tepat. Setiap upaya untuk memunculkan fhatmenghasilkan kesalahan, tampaknya dalam penggunaan uniroot.
Whuber
3
Apa versi R Anda? Kode menggunakan argumen baru untuk uniroot, extendedInt, yang diperkenalkan di R versi 3.1. Jika R Anda lebih tua, Anda dapat mencoba menghapusnya, (dan memperpanjang interval yang diberikan untuk uniroot). Tapi itu akan membuat kodenya kurang kuat!
kjetil b halvorsen
10

The Persamaan Welch-Satterthwaite dapat digunakan untuk memberikan perkiraan jawaban dalam bentuk distribusi gamma. Ini memiliki properti yang bagus untuk membiarkan kami memperlakukan distribusi gamma sebagai (kurang-lebih) ditutup dengan penambahan. Ini adalah perkiraan dalam uji-t Welch yang biasa digunakan.

(Distribusi gamma dapat dilihat sebagai distribusi chi-square diskalakan, dan memungkinkan parameter bentuk non-integer.)

k,θ

kskamum=(sayaθsayaksaya)2sayaθsaya2ksaya

θskamum=θsayaksayakskamum

k=(3,4,5)θ=(1,2,1)

Jadi kita mendapatkan sekitar Gamma (10.666 ..., 1.5)

kθsayaθ

Paul Harrison
sumber
6

n

GDC(a,b,α,β;τ)={baβαΓ(a+α)ebττa+α11F1[α,a+α,(bβ)τ],τ>00,τ0,
Gamma(a,b)Γ(a,1/b)bβ
Carl
sumber