Generische Summe von Gamma-Zufallsvariablen

Ich habe gelesen, dass die Summe der Gamma-Zufallsvariablen mit demselben Skalenparameter eine andere Gamma-Zufallsvariable ist. Ich habe auch gesehen, dass der Artikel von Moschopoulos eine Methode zur Summierung einer allgemeinen Menge von Gamma-Zufallsvariablen beschreibt. Ich habe versucht, die Methode von Moschopoulos zu implementieren , habe aber noch keinen Erfolg.

Wie sieht die Summe einer allgemeinen Menge von Gamma-Zufallsvariablen aus? Um diese Frage konkret zu machen, wie sieht sie aus:

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

Wenn die oben genannten Parameter nicht besonders aussagekräftig sind, schlagen Sie bitte andere vor.

Zunächst verbinden alle Summen denselben Skalierungsfaktor aufweist : a sowie eine variate Form a Veränderlichen.Γ ( m , β ) Γ ( n + m , β $\Gamma(n, \beta)$$\Gamma(m,\beta)$$\Gamma(n+m,\beta)$

Als nächstes werden beobachten , dass die Kennlinienfunktion (CF) des ist , von wo aus der cf einer Summe dieser Verteilungen ist das Produkt( 1 - i β t ) - $\Gamma(n, \beta)$$(1-i \beta t)^{-n}$

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

Wenn alle ganzzahlig sind, expandiert dieses Produkt als Teilbruch zu einer linearen Kombination von wobei die ganze Zahlen zwischen und . In dem Beispiel mit (aus der Summe von und ) und finden wir , ( 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}.$$

Das Gegenteil von cf ist die inverse Fourier-Transformation, die linear ist . Das heißt, wir können sie termweise anwenden. Jeder Term ist als ein Vielfaches des cf einer Gamma-Verteilung erkennbar und kann daher leicht invertiert werden, um das PDF zu erhalten . Im Beispiel erhalten wir

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

für das PDF der Summe.

Dies ist eine endliche Mischung von Gamma-Verteilungen mit Skalierungsfaktoren, die denen in der Summe entsprechen, und Formfaktoren, die denen in der Summe entsprechen oder darunter liegen. Mit Ausnahme von Sonderfällen (in denen ein gewisser Widerruf auftreten kann) wird die Anzahl der Terme durch den Gesamtformparameter (vorausgesetzt, alle sind unterschiedlich).n $n_1 + n_2 + \cdots$$n_j$

---

Als Test sehen Sie hier ein Histogramm mit Ergebnissen, die durch Addition von unabhängigen Zügen aus den Verteilungen und . Darauf wird der Graph der fachen vorhergehenden Funktion überlagert . Die Passform ist sehr gut. ≤ ( 8 , 1 ) ≤ ( 4 , 2 ) 10 $10^4$$\Gamma(8,1)$$\Gamma(4,2)$$10^4$

---

Moschopoulos führt diese Idee noch einen Schritt weiter, indem er den cf der Summe zu einer unendlichen Reihe von Gamma-Kennlinienfunktionen erweitert, wenn eines oder mehrere der nicht ganzzahlig sind, und dann die unendliche Reihe an einem Punkt beendet, an dem sie einigermaßen gut angenähert ist.$n_i$

Ich werde eine andere mögliche Lösung aufzeigen, die ziemlich weit verbreitet und mit der heutigen R-Software ziemlich einfach zu implementieren ist. Das ist die Sattelpunktdichte-Näherung, die weiter bekannt sein sollte!

Für die Terminologie zur Gammaverteilung folge ich https://en.wikipedia.org/wiki/Gamma_distribution mit der Form- / Skalierungsparametrierung, ist der Formparameter und ist die Skalierung. Für die Sattelpunktnäherung werde ich Ronald W Butler folgen: "Sattelpunktnäherungen mit Anwendungen" (Cambridge UP). Die Sattelpunktnäherung wird hier erklärt: Wie funktioniert die Sattelpunktnäherung? hier werde ich zeigen, wie es in dieser Anwendung verwendet wird.$k$$\theta$

Sei eine Zufallsvariable mit der existierenden momenterzeugenden Funktion die für in einem offenen Intervall existieren muss , das Null enthält. Definieren Sie dann die kumulative Erzeugungsfunktion durch Es ist bekannt, dass . Die Sattelpunktgleichung ist was implizit als eine Funktion von (die im Bereich von ). Wir schreiben diese implizit definierte Funktion als . Beachten Sie, dass die Sattelpunktgleichung immer genau eine Lösung hat, da die kumulative Funktion konvex ist. M ( s ) = E e s X s K ( s ) = log M ( s ) E X = K ' ( 0 ) , Var ( x ) = K " ( 0 ) K ' ( s ) = x s x x s ( x $X$

$$M(s) = E e^{sX}$$ $s$ $$K(s) = \log M(s)$$ $E X = K'(0), \text{Var} (X) = K''(0)$ $$K'(\hat{s}) = x$$ $s$$x$$X$$\hat{s}(x)$

Dann ist die Sattelpunktnäherung an die Dichte von gegeben durch Es ist nicht garantiert, dass diese ungefähre Dichtefunktion zu 1 integriert wird, ebenso wie die nicht normalisierte Sattelpunktnäherung. Wir könnten es numerisch integrieren und renormieren, um eine bessere Annäherung zu erhalten. Diese Annäherung ist jedoch garantiert nicht negativ.X f ( x ) = $f$$X$

$$\hat{f}(x) = \frac1{\sqrt{2\pi K''(\hat{s})}} \exp(K(\hat{s}) - \hat{s} x)$$

Nun seien unabhängige Gamma-Zufallsvariablen, wobei die Verteilung mit Parametern . Dann ist die kumulativ erzeugende Funktion definiert für . Die erste Ableitung ist und die zweite Ableitung ist Im Folgenden werde ich einen Code angeben, der dies berechnet, und die Parameterwerte , ,$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)$$ $s<1/\max(\theta_1, \theta_2, \dots, \theta_n)$ $$K'(s) = \sum_{i=1}^n \frac{k_i \theta_i}{1-\theta_i s}$$ n=3k=(1,2,3)θ $$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)$. Beachten Sie, dass der folgende RCode ein neues Argument in der in R 3.1 eingeführten Uniroot-Funktion verwendet und daher nicht in älteren Rs ausgeführt wird.

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")

Daraus ergibt sich die folgende Handlung:

Ich werde die normalisierte Sattelpunktnäherung als Übung verlassen.

Die Welch-Satterthwaite-Gleichung könnte verwendet werden, um eine ungefähre Antwort in Form einer Gammaverteilung zu geben . Dies hat die nette Eigenschaft, dass wir Gammaverteilungen als (ungefähr) geschlossen behandeln, wenn sie hinzugefügt werden. Dies ist die Näherung im allgemein verwendeten Welch-T-Test.

(Die Gamma-Verteilung kann als skalierte Chi-Quadrat-Verteilung betrachtet werden und erlaubt nicht ganzzahlige Formparameter.)

Ich habe die Approximation an die Parametrisierung der Gammaverteilung angepasst:$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} }$$

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

Wir erhalten also ungefähr Gamma (10.666 ..., 1.5)

Wir sehen, dass der Formparameter mehr oder weniger summiert wurde, aber etwas weniger, weil sich die Eingangsskalenparameter unterscheiden. ist so, dass die Summe den richtigen Mittelwert hat.$k$$\theta_i$$\theta$

Eine genaue Lösung der Faltung (dh der Summe) von Gammaverteilungen ist gegeben als Gl. (1) im verlinkten pdf von DiSalvo . Da dies etwas lang ist, wird es einige Zeit dauern, es hierher zu kopieren. Für nur zwei Gammaverteilungen ist ihre exakte Summe in geschlossener Form durch Gl. (2) von DiSalvo und ohne Gewichte nach Gl. (5) von Wesolowski et al. , der auch auf der CV-Seite als Antwort auf diese Frage erscheint. Das ist,$n$Gamma(a,b)→Γ(a,1/b)b

$$\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.,$$ wobei die Notation in den obigen Fragen; hier. Das heißt, und sind hier Geschwindigkeitskonstanten und keine Zeitskalare.$Gamma(a,b) \rightarrow \Gamma(a,1/b)$$b$$\beta$