Generische Summe von Gamma-Zufallsvariablen

35

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:

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

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

OSE
quelle
4
Eine explizite Lösung für die Summe von zwei beliebigen Gamma-Verteilungen wurde unter stats.stackexchange.com/a/252192 veröffentlicht .
whuber
Ein spezielles Beispiel hierfür, bei dem alle Gamma-Verteilungen den Formparameter 1 haben ( dh exponentiell sind), ist die hypoexponentielle Verteilung (Familie) . Für den Fall von nur zwei Exponentialverteilungen gibt es auch eine explizite Formel unter stats.stackexchange.com/questions/412849 .
Whuber

Antworten:

37

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

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 ) - nΓ(n,β)(1iβt)n

j1(1iβjt)nj.

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 = 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.

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

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

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 jn1+n2+nj


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 4104Γ(8,1)Γ(4,2)104

Zahl


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.ni

whuber
quelle
2
Kleiner Kommentar: Typischerweise bedeutet eine endliche Mischung ein PDF der Form wobei und , das heißt, die sind Wahrscheinlichkeiten und das pdf können als die (Gesetz der Gesamtwahrscheinlichkeit) gewichtete Summe von bedingten pdfs unter verschiedenen Bedingungen interpretiert werden, die mit Wahrscheinlichkeiten . In der obigen Summe sind jedoch einige der Koeffizienten negativ, und daher findet die Standardinterpretation der Mischung keine Anwendung. a i > 0 Σ i a i = 1 a i a i
f(x)=i=1naifi(x)
ai>0iai=1aiai
Dilip Sarwate
@ Dilip Das ist ein guter Punkt. Was diesen Fall interessant macht, ist, dass, obwohl einige der Koeffizienten negativ sein können, diese Kombination dennoch eine gültige Verteilung ist (aufgrund ihrer Konstruktion).
Whuber
Kann dieser Ansatz dahingehend erweitert werden, dass abhängige Variablen hinzugefügt werden? Insbesondere möchte ich 6 Verteilungen addieren, von denen jede eine gewisse Korrelation mit den anderen aufweist.
Stampfer
11

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θ

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)=EesX
s
K(s)=logM(s)
EX=K(0),Var(X)=K(0)
K(s^)=x
sxXs^(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 ) = 1fX

f^(x)=12πK(s^)exp(K(s^)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 , ,X1,X2,,XnXi(ki,θi)

K(s)=i=1nkiln(1θis)
s<1/max(θ1,θ2,,θn)
K(s)=i=1nkiθi1θis
n=3k=(1,2,3)θ=
K(s)=i=1nkiθi2(1θis)2.
Rn=3k=(1,2,3)θ=(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: Bildbeschreibung hier eingeben

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

kjetil b halvorsen
quelle
1
Das ist interessant, aber ich kann nicht dafür sorgen, dass Ihr RCode funktioniert, um die Annäherung an die exakte Antwort zu vergleichen. Jeder Aufrufversuch fhaterzeugt Fehler, anscheinend bei der Verwendung von uniroot.
whuber
3
Was ist deine R-Version? In den Codes wird ein neues Argument zum Unirooten verwendet, "extendInt", das in R Version 3.1 eingeführt wurde. Wenn Ihr R älter ist, können Sie versuchen, dieses zu entfernen (und das Intervall für "uniroot" zu verlängern). Aber das macht den Code weniger robust!
kjetil b halvorsen
10

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,θ

ksum=(iθiki)2iθi2ki

θsum=θikiksum

Sei ,k=(3,4,5)θ=(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θiθ

Paul Harrison
quelle
6

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,nGamma(a,b)Γ(a,1/b)bβ

GDC(a,b,α,β;τ)={baβαΓ(a+α)ebττa+α11F1[α,a+α,(bβ)τ],τ>00,τ0,
wobei die Notation in den obigen Fragen; hier. Das heißt, und sind hier Geschwindigkeitskonstanten und keine Zeitskalare.Gamma(a,b)Γ(a,1/b)bβ
Carl
quelle