Somme générique de variables aléatoires Gamma

35

J'ai lu que la somme des variables aléatoires Gamma ayant le même paramètre d'échelle est une autre variable aléatoire Gamma. J'ai également vu le document de Moschopoulos décrivant une méthode de sommation d'un ensemble général de variables aléatoires Gamma. J'ai essayé d'implémenter la méthode de Moschopoulos mais je n'ai pas encore réussi.

À quoi ressemble la somme d'un ensemble général de variables aléatoires Gamma? Pour concrétiser cette question, à quoi ça ressemble:

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

Si les paramètres ci-dessus ne sont pas particulièrement révélateurs, veuillez en suggérer d'autres.

probability distributions gamma-distribution summations saddlepoint-approximation OSE
la source

4

Une solution explicite pour la somme de deux distributions Gamma quelconques a été publiée à l' adresse stats.stackexchange.com/a/252192 .

whuber

Un exemple particulier, où toutes les distributions Gamma ont le paramètre de forme 1 (c'est-à-dire qu'elles sont exponentielles), est appelé distribution hypoexponentielle (famille) . Pour le cas de deux distributions exponentielles uniquement, une formule explicite est également donnée à l' adresse stats.stackexchange.com/questions/412849 .

whuber

37

Commencez par combiner les sommes ayant le même facteur d'échelle : a plus une variable sous une forme . $\Gamma(n, \beta)$ $\Gamma(m,\beta)$ $\Gamma(n+m,\beta)$

Ensuite, observons que la fonction caractéristique (cf) de est , d’où le cf d’une somme de ces distributions est le produit $\Gamma(n, \beta)$ $(1-i \beta t)^{-n}$

\underset{j}{Π} \frac{1}{(1 - je β_{j} t)^{n_{j}}} .

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

Lorsque les sont tous , ce produit se transforme en fraction partielle en une combinaison linéaire de où les sont des entiers compris entre et . Dans l'exemple avec (de la somme de et de ) et de on trouve $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 - 2 je t)^{4}} = \frac{1}{(X + je)^{8}} - \frac{8 je}{(X + je)^{7}} - \frac{40}{(X + je)^{6}} + \frac{160 je}{(X + je)^{5}} + \frac{560}{(X + je)^{4}} - \frac{1792 je}{(X + je)^{3}} - \frac{5376}{(X + je)^{2}} + \frac{15360 je}{X + je} + \frac{256}{(2 X + je)^{4}} + \frac{2048 je}{(2 X + je)^{3}} - \frac{9216}{(2 X + je)^{2}} - \frac{30720 je}{2 X + je} .

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

L'inverse de prendre le cf est l'inverse de la transformée de Fourier, qui est linéaire : cela signifie que nous pouvons l'appliquer terme à terme. Chaque terme est reconnaissable comme un multiple du cf d'une distribution gamma et est donc facilement inversé pour produire le PDF . Dans l'exemple on obtient

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

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

pour le PDF de la somme.

Il s'agit d'un mélange fini de distributions gamma ayant des facteurs d'échelle égaux à ceux de la somme et des facteurs de forme inférieurs ou égaux à ceux de la somme. Sauf dans des cas particuliers (où une certaine annulation peut survenir), le nombre de termes est donné par le paramètre de forme totale (en supposant que tous les sont différents). $n_1 + n_2 + \cdots$ $n_j$

À titre de test, voici un histogramme de résultats obtenus en ajoutant des tirages indépendants à partir des distributions et . Sur elle se superpose le graphe de fois la fonction précédente. La coupe est très bonne. $10^4$ $\Gamma(8,1)$ $\Gamma(4,2)$ $10^4$

Moschopoulos va encore plus loin dans cette idée en élargissant le cf de la somme en une série infinie de fonctions caractéristiques gamma chaque fois qu'un ou plusieurs des sont pas entiers, puis termine la série infinie en un point où elle est raisonnablement approximée. $n_i$

whuber
la source

2

Commentaire Minor: En règle générale, un mélange fini signifie un pdf de la forme où et , qui est, les sont les probabilités et le pdf peuvent être interprétés comme la somme pondérée (loi de la probabilité totale) des fichiers PDF conditionnels, compte tenu de diverses conditions qui se produisent avec des probabilités . Cependant, dans la somme ci-dessus, certains des coefficients sont négatifs et l'interprétation standard du mélange ne s'applique donc pas.

f (x) = \sum_{i = 1}^{n} a_{i} f_{i} (x)

$f(x) = \sum_{i=1}^n a_i f_i(x)$

a_{i} > 0

$a_i > 0$

\sum_{i} a_{i} = 1

$\sum_i a_i = 1$

a_{i}

$a_i$

a_{i}

$a_i$

Dilip Sarwate

@Dilip C'est un bon point. Ce qui rend ce cas intéressant est que, bien que certains des coefficients puissent être négatifs, cette combinaison reste néanmoins une distribution valide (de par sa construction même).

whuber

Cette approche peut-elle être étendue pour tenir compte de l'ajout de variables dépendantes? En particulier, je veux additionner 6 distributions, chacune ayant une corrélation avec les autres.

Masher

11

Je vais montrer une autre solution possible, qui est assez largement applicable, et avec le logiciel actuel R, assez facile à mettre en œuvre. C’est l’approximation de la densité de la pointe de la selle, qui devrait être plus largement connue!

Pour la terminologie relative à la distribution gamma, je suivrai https://en.wikipedia.org/wiki/Gamma_distribution avec la paramétrisation shape / scale, est le paramètre de forme et est scale. Pour l’approximation du point à cheval, je vais suivre Ronald W Butler: "Approximations à l’aide des applications" (Cambridge UP). L’approximation du point de selle est expliquée ici: Comment fonctionne l’approximation du point de selle? Ici, je vais montrer comment il est utilisé dans cette application. $k$ $\theta$

Soit une variable aléatoire avec une fonction génératrice de moment existante qui doit exister pour dans un intervalle ouvert contenant zéro. Puis définissez la fonction génératrice du cumulant par On sait que . L'équation du point de selle est qui définit implicitement en fonction de (qui doit être compris dans la plage de ). Nous écrivons cette fonction définie implicitement comme . Notez que l'équation du point d'équilibre a toujours exactement une solution, car la fonction de cumulant est convexe. $X$

M (s) = E e^{s X}

$M(s) = E e^{sX}$

s

$s$

K (s) = bûche M (s)

$K(s) = \log M(s)$

E X = K^{'} (0), Var (X) = K^{″} (0)

$E X = K'(0), \text{Var} (X) = K''(0)$

K^{'} (\hat{s}) = X

$K'(\hat{s}) = x$

s

$s$

x

$x$

X

$X$

\hat{s} (x)

$\hat{s}(x)$

Ensuite, l’approximation du point-selle de la densité de est donnée par L'intégration de la fonction de densité approximative à 1 n'est pas garantie, pas plus que l'approximation du point de selle non normalisée. Nous pourrions l’intégrer numériquement et le renormaliser pour obtenir une meilleure approximation. Mais cette approximation est garantie non négative. $f$ $X$

\hat{F} (X) = \frac{1}{\sqrt{2 π K^{"} (\hat{s})}} \exp (K (\hat{s}) - \hat{s} X)

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

Soit maintenant des variables aléatoires gamma indépendantes, où a la distribution avec les paramètres . La fonction génératrice du cumulant est alors défini pour . La première dérivée est et la dérivée seconde est Dans ce qui suit, je donnerai du code le calculant et utiliserai les valeurs de paramètre , , $X_1, X_2, \dots, X_n$ $X_i$ $(k_i, \theta_i)$

K (s) = - Σ_{je = 1}^{n} k_{je} dans (1 - θ_{je} s)

$K(s) = -\sum_{i=1}^n k_i \ln(1-\theta_i s)$

s < 1 / max (θ_{1}, θ_{2}, \dots, θ_{n})

$s<1/\max(\theta_1, \theta_2, \dots, \theta_n)$

K^{'} (s) = Σ_{je = 1}^{n} \frac{k_{je} θ_{je}}{1 - θ_{je} s}

$K'(s) = \sum_{i=1}^n \frac{k_i \theta_i}{1-\theta_i s}$

K^{"} (s) = Σ_{je = 1}^{n} \frac{k_{je} θ_{je}^{2}}{(1 - θ_{je} s)^{2}} .

$K''(s) = \sum_{i=1}^n \frac{k_i \theta_i^2}{(1-\theta_i s)^2}.$ R

n = 3

$n=3$

k = (1, 2, 3)

$k=(1,2,3)$

θ = (1, 2, 3)

$\theta=(1,2,3)$ . Notez que le Rcode suivant utilise un nouvel argument dans la fonction uniroot introduite dans R 3.1 et ne s'exécutera donc pas dans les anciens 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")

résultant en l'intrigue suivante: entrez la description de l'image ici

Je laisserai l'approximation du point de selle normalisée comme un exercice.

kjetil b halvorsen
la source

1

C'est intéressant, mais je ne peux pas faire fonctionner votre Rcode pour comparer l'approximation à la réponse exacte. Toute tentative d' fhatappel génère des erreurs, apparemment dans l'utilisation de uniroot.

whuber

3

Quelle est votre version R? Les codes utilisent un nouvel argument pour uniroot, extendInt, qui a été introduit dans la version 3.1 de R Si votre R est plus ancien, vous pouvez essayer de le supprimer (et d'étendre l'intervalle donné à uniroot). Mais cela rendra le code moins robuste!

kjetil b halvorsen

10

L' équation de Welch – Satterthwaite pourrait être utilisée pour donner une réponse approximative sous la forme d'une distribution gamma. Ceci a la belle propriété de nous permettre de traiter les distributions gamma comme étant (approximativement) fermées sous addition. C'est l'approximation du test t de Welch couramment utilisé.

(La distribution gamma peut être vue comme une distribution chi-carré mise à l'échelle, et permettant un paramètre de forme non entier.)

J'ai adapté l'approximation à la paramétrisation de la distribution gamma: $k, \theta$

k_{s vous m} = \frac{(\underset{je}{Σ} θ_{je} k_{je})^{2}}{\underset{je}{Σ} θ_{je}^{2} k_{je}}

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

θ_{s vous m} = \frac{Σ θ_{je} k_{je}}{k_{s vous m}}

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

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

Nous obtenons donc environ Gamma (10.666 ..., 1.5)

Nous voyons que le paramètre de forme a été plus ou moins totalisé, mais un peu moins parce que les paramètres d’échelle en entrée diffèrent. est tel que la somme a la valeur moyenne correcte. $k$ $\theta_i$ $\theta$

Paul Harrison
la source

6

Une solution exacte à la convolution (c'est-à-dire la somme) de distributions gamma est donnée par Eq. (1) dans le pdf lié par DiSalvo . Comme c'est un peu long, il faudra du temps pour le copier ici. Pour seulement deux distributions gamma, leur somme exacte sous forme fermée est spécifiée par Eq. (2) de DiSalvo et sans poids par Eq. (5) de Wesolowski et al. , qui figure également sur le site du CV en réponse à cette question. C'est, $n$

g ré C (une, b, α, β; τ) = {\begin{array}{cc} \frac{b^{une} β^{α}}{Γ (une + α)} e^{- b τ} {τ^{une + α}}^{- 1}_{1} F_{1} [α, une + α, (b - β) τ], & τ > 0 \\ 0, τ \leq 0 \end{array},

$\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.,$ où la notation dans les questions ci-dessus; , ici. C'est-à-dire que et sont des constantes de vitesse ici et non des scalaires temporels.

G a m m a (a, b) \to Γ (a, 1 / b)

$Gamma(a,b) \rightarrow \Gamma(a,1/b)$

b

$b$

β

$\beta$

Carl
la source

Somme générique de variables aléatoires Gamma

Réponses: