Ogólna suma zmiennych losowych Gamma

35

Mam przeczytać , że suma gamma zmiennych losowych o tym samym parametrem skali jest inna zmienna losowa Gamma. Widziałem też artykuł Moschopoulosa opisujący metodę sumowania ogólnego zestawu zmiennych losowych Gamma. Próbowałem wdrożyć metodę Moschopoulosa, ale jeszcze nie odniosłem sukcesu.

Jak wygląda suma ogólnego zestawu zmiennych losowych Gamma? Aby sprecyzować to pytanie, jak to wygląda:

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

Jeśli powyższe parametry nie są szczególnie odkrywcze, sugeruj inne.

OSE
źródło
4
Wyraźne rozwiązanie dla sumy dowolnych dwóch dystrybucji gamma opublikowano na stronie stats.stackexchange.com/a/252192 .
whuber
Szczególny przykład tego, w którym wszystkie rozkłady gamma mają parametr kształtu 1 (czyli są wykładnicze), nazywa się rozkładem hipoeksponencjalnym (rodzina) . W przypadku tylko dwóch rozkładów wykładniczych istnieje również wyraźna formuła podana na stronie stats.stackexchange.com/questions/412849 .
whuber

Odpowiedzi:

37

Najpierw połącz dowolne sumy o tym samym współczynniku skali : a plus Γ ( m , β ) zmieniają się od a Γ ( n + m , β )Γ(n,β)Γ(m,β)Γ(n+m,β) .

Następnie zauważ, że funkcją charakterystyczną (cf) dla jest ( 1 - i β t ) - n , skąd cf sumy tych rozkładów jest iloczynemΓ(n,β)(1iβt)n

j1(1iβjt)nj.

Gdy integralną, produkt ten rozszerza się częściowo frakcji do liniowej kombinacji z ( 1 - i β J T ) - ν gdzie ν liczby całkowite od 1 i n j . W przykładzie z β 1 = 1 , n 1 = 8 (z sumy Γ ( 3 , 1 ) i Γ ( 5 )nj (1iβjt)νν1njβ1=1,n1=8Γ(3,1)Γ(5,1) ) i znajdujemyβ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.

The inverse of taking the cf is the inverse Fourier Transform, which is linear: that means we may apply it term by term. Each term is recognizable as a multiple of the cf of a Gamma distribution and so is readily inverted to yield the PDF. In the example we obtain

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

for the PDF of the sum.

This is a finite mixture of Gamma distributions having scale factors equal to those within the sum and shape factors less than or equal to those within the sum. Except in special cases (where some cancellation might occur), the number of terms is given by the total shape parameter n1+n2+ (assuming all the nj are different).


As a test, here is a histogram of 104 results obtained by adding independent draws from the Γ(8,1) and Γ(4,2) distributions. On it is superimposed the graph of 104 times the preceding function. The fit is very good.

Postać


Moschopoulos carries this idea one step further by expanding the cf of the sum into an infinite series of Gamma characteristic functions whenever one or more of the ni is non-integral, and then terminates the infinite series at a point where it is reasonably well approximated.

whuber
źródło
2
Minor comment: Typically, a finite mixture means a pdf of the form
f(x)=i=1naifi(x)
where ai>0 and iai=1, that is, the ai are probabilities and the pdf can be interpreted as the (law of total probability) weighted sum of conditional pdfs given various conditions that occur with probabilities ai. However, in the sum above, some of the coefficients are negative and thus the standard interpretation of the mixture does not apply.
Dilip Sarwate
@Dilip That's a good point. What makes this case interesting is that although some of the coefficients may be negative, nevertheless this combination is still a valid distribution (by its very construction).
whuber
Can this approach be extended to account for addition of dependent variables? In particular, I want to add up 6 distributions with each having some correlation with the others.
masher
11

I will show another possible solution, that is quite widely applicable, and with todays R software, quite easy to implement. That is the saddlepoint density approximation, which ought to be wider known!

For terminology about the gamma distribution, I will follow https://en.wikipedia.org/wiki/Gamma_distribution with the shape/scale parametrization, k is shape parameter and θ is scale. For the saddlepoint approximation I will follow Ronald W Butler: "Saddlepoint approximations with applications" (Cambridge UP). The saddlepoint approximation is explained here: How does saddlepoint approximation work? here I will show how it is used in this application.

X

M(s)=EesX
s
K(s)=logM(s)
EX=K(0),Var(X)=K(0). The saddlepoint equation is
K(s^)=x
which implicitely defines s as a function of x (which must be in the range of X). We write this implicitely defined function as s^(x). Note that the saddlepoint equation always has exactly one solution, because the cumulant function is convex.

Then the saddlepoint approximation to the density f of X is given by

f^(x)=12πK(s^)exp(K(s^)s^x)
This approximate density function is not guaranteed to integrate to 1, so is the unnormalized saddlepoint approximation. We could integrate it numerically and the renormalize to get a better approximation. But this approximation is guaranteed to be non-negative.

Now let X1,X2,,Xn be independent gamma random variables, where Xi has the distribution with parameters (ki,θi). Then the cumulant generating function is

K(s)=i=1nkiln(1θis)
defined for s<1/max(θ1,θ2,,θn). The first derivative is
K(s)=i=1nkiθi1θis
and the second derivative is
K(s)=i=1nkiθi2(1θis)2.
In the following I will give some R code calculating this, and will use the parameter values n=3, k=(1,2,3), θ=(1,2,3). Note that the following R code uses a new argument in the uniroot function introduced in R 3.1, so will not run in older R's.
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")

resulting in the following plot: enter image description here

I will leave the normalized saddlepoint approximation as an exercise.

kjetil b halvorsen
źródło
1
This is interesting, but I cannot make your R code work to compare the approximation to the exact answer. Any attempt to invoke fhat generates errors, apparently in the use of uniroot.
whuber
3
What is your R version? The codes uses a new argument to uniroot, extendInt, which was introduces in R version 3.1 If your R is older, you might try to remove that, (and extend the interval given to uniroot). But that will make the code less robust!
kjetil b halvorsen
10

The Welch–Satterthwaite equation could be used to give an approximate answer in the form of a gamma distribution. This has the nice property of letting us treat gamma distributions as being (approximately) closed under addition. This is the approximation in the commonly used Welch's t-test.

(The gamma distribution is can be viewed as a scaled chi-square distribution, and allowing non-integer shape parameter.)

I've adapted the approximation to the k,θ parametrization of the gamma distriubtion:

ksum=(iθiki)2iθi2ki

θsum=θikiksum

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

So we get approximately Gamma(10.666... ,1.5)

We see the shape parameter k has been more or less totalled, but slightly less because the input scale parameters θi differ. θ is such that the sum has the correct mean value.

Paul Harrison
źródło
6

An exact solution to the convolution (i.e., sum) of n gamma distributions is given as Eq. (1) in the linked pdf by DiSalvo. As this is a bit long, it will take some time to copy it over here. For only two gamma distributions, their exact sum in closed form is specified by Eq. (2) of DiSalvo and without weights by Eq. (5) of Wesolowski et al., which also appears on the CV site as an answer to that question. That is,

GDC(a,b,α,β;τ)={baβαΓ(a+α)ebττa+α11F1[α,a+α,(bβ)τ],τ>00,τ0,
where the notation in the questions above; Gamma(a,b)Γ(a,1/b), here. That is, b and β are rate constants here and not time scalars.
Carl
źródło