Jumlah generik dari variabel acak Gamma

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:

$\text{Gamma}(3,1) + \text{Gamma}(4,2) + \text{Gamma}(5,1)$

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

Pertama, gabungkan setiap jumlah yang memiliki faktor skala yang sama : a ditambah variasi varian .$\Gamma(n, \beta)$$\Gamma(m,\beta)$$\Gamma(n+m,\beta)$

Selanjutnya, amati bahwa fungsi karakteristik (cf) dari adalah , di mana cf dari sejumlah distribusi ini adalah produk$\Gamma(n, \beta)$$(1-i \beta t)^{-n}$

$$\prod_{j} \frac{1}{(1-i \beta_j t)^{n_j}}.$$

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 = $n_j$ $(1-i \beta_j t)^{-\nu}$$\nu$$1$$n_j$$\beta_1 = 1, n_1=8$$\Gamma(3,1)$$\Gamma(5,1)$$\beta_2 = 2, n_2=4$

$$\frac{1}{(1-i t)^{8}}\frac{1}{(1- 2i t)^{4}} = \\ \frac{1}{(x+i)^8}-\frac{8 i}{(x+i)^7}-\frac{40}{(x+i)^6}+\frac{160 i}{(x+i)^5}+\frac{560}{(x+i)^4}-\frac{1792 i}{(x+i)^3}\\-\frac{5376}{(x+i)^2}+\frac{15360 i}{x+i}+\frac{256}{(2 x+i)^4}+\frac{2048 i}{(2 x+i)^3}-\frac{9216}{(2 x+i)^2}-\frac{30720 i}{2 x+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

$$\frac{e^{-t} t^7}{5040}+\frac{1}{90} e^{-t} t^6+\frac{1}{3} e^{-t} t^5+\frac{20}{3} e^{-t} t^4+\frac{8}{3} e^{-\frac{t}{2}} t^3+\frac{280}{3} e^{-t} t^3\\ -128 e^{-\frac{t}{2}} t^2+896 e^{-t} t^2+2304 e^{-\frac{t}{2}} t+5376 e^{-t} t-15360 e^{-\frac{t}{2}}+15360 e^{-t}$$

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).$n_1 + n_2 + \cdots$$n_j$

---

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

---

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.$n_i$

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$$\theta$

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) = E e^{sX}$$ $s$ $$K(s) = \log M(s)$$ $E X = K'(0), \text{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'(\hat{s}) = x$$ $s$$x$$X$$\hat{s}(x)$

Kemudian saddlepoint pendekatan kepadatan dari X diberikan oleh f ( x ) = $f$$X$ 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.

$$\hat{f}(x) = \frac1{\sqrt{2\pi K''(\hat{s})}} \exp(K(\hat{s}) - \hat{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 $X_1, X_2, \dots, X_n$$X_i$$(k_i, \theta_i)$

$$K(s) = -\sum_{i=1}^n k_i \ln(1-\theta_i s)$$ . Derivatif pertama adalah K ′ ( s ) = n ∑ i = 1 k i θ $s<1/\max(\theta_1, \theta_2, \dots, \theta_n)$ dan turunan kedua adalah K″(s)= n ∑ i=1kiθ 2 $$K'(s) = \sum_{i=1}^n \frac{k_i \theta_i}{1-\theta_i s}$$ 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) = \sum_{i=1}^n \frac{k_i \theta_i^2}{(1-\theta_i s)^2}.$$ R$n=3$$k=(1,2,3)$$\theta=(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:

Saya akan meninggalkan pendekatan saddlepoint yang dinormalisasi sebagai latihan.

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, \theta$

$$k_{sum} = { (\sum_i \theta_i k_i)^2 \over \sum_i \theta_i^2 k_i }$$

$$\theta_{sum} = { { \sum \theta_i k_i } \over k_{sum} }$$

$k=(3,4,5)$$\theta=(1,2,1)$

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

$k$$\theta_i$$\theta$

$n$

$$\mathrm{G}\mathrm{D}\mathrm{C}\left(\mathrm{a}\kern0.1em ,\mathrm{b}\kern0.1em ,\alpha, \beta; \tau \right)=\left\{\begin{array}{cc}\hfill \frac{{\mathrm{b}}^{\mathrm{a}}{\beta}^{\alpha }}{\Gamma \left(\mathrm{a}+\alpha \right)}{e}^{-\mathrm{b}\tau }{\tau^{\mathrm{a}+\alpha}}^{-1}{}_1F_1\left[\alpha, \mathrm{a}+\alpha, \left(\mathrm{b}-\beta \right)\tau \right],\hfill & \hfill \tau >0\hfill \\ {}\hfill \kern2em 0\kern6.6em ,\hfill \kern5.4em \tau \kern0.30em \le \kern0.30em 0\hfill \end{array}\right.,$$ $Gamma(a,b) \rightarrow \Gamma(a,1/b)$$b$$\beta$