Soma genérica de variáveis aleatórias gama

35

Eu li que a soma das variáveis aleatórias Gamma com o mesmo parâmetro de escala é outra variável aleatória Gamma. Também vi o artigo de Moschopoulos descrevendo um método para a soma de um conjunto geral de variáveis aleatórias gama. Eu tentei implementar o método de Moschopoulos, mas ainda não obtive sucesso.

Como é a soma de um conjunto geral de variáveis aleatórias Gamma? Para tornar essa pergunta concreta, como ela é:

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

Se os parâmetros acima não forem particularmente reveladores, sugira outros.

probability distributions gamma-distribution summations saddlepoint-approximation OSE
fonte

4

Uma solução explícita para a soma de quaisquer duas distribuições Gamma foi publicada em stats.stackexchange.com/a/252192 .

whuber

Um exemplo especial disso, em que todas as distribuições Gamma têm o parâmetro de forma 1 (isto é, são exponenciais) é chamada de distribuição hipoxponencial (família) . No caso de apenas duas distribuições exponenciais, também existe uma fórmula explícita dada em stats.stackexchange.com/questions/412849 .

whuber

37

Primeiro, combine todas as somas com o mesmo fator de escala : a mais uma forma . $\Gamma(n, \beta)$ $\Gamma(m,\beta)$ $\Gamma(n+m,\beta)$

Em seguida, observe que a função característica (cf) de é , de onde o cf de uma soma dessas distribuições é o produto $\Gamma(n, \beta)$ $(1-i \beta t)^{-n}$

\prod_{j} \frac{1}{(1 - Eu β_{j} t)^{n_{j}}} .

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

Quando os são todos integrais, esse produto se expande como uma fração parcial para uma combinação linear de que os são números inteiros entre e . No exemplo com (da soma de e ) e , encontramos $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 - Eu t)^{8}} \frac{1}{(1 - 2 Eu t)^{4}} = \frac{1}{(x + Eu)^{8}} - \frac{8 Eu}{(x + Eu)^{7}} - \frac{40.}{(x + Eu)^{6}} + \frac{160 Eu}{(x + Eu)^{5}} + \frac{560}{(x + Eu)^{4}} - \frac{1792 Eu}{(x + Eu)^{3}} - \frac{5376}{(x + Eu)^{2}} + \frac{15360 Eu}{x + Eu} + \frac{256}{(2 x + Eu)^{4}} + \frac{2048 Eu}{(2 x + Eu)^{3}} - \frac{9216}{(2 x + Eu)^{2}} - \frac{30720 Eu}{2 x + Eu} .

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

O inverso de tomar o cf é a Transformada de Fourier inversa, que é linear : isso significa que podemos aplicá-lo termo a termo. Cada termo é reconhecível como um múltiplo do cf de uma distribuição gama e, portanto, é prontamente invertido para produzir o PDF . No exemplo, obtemos

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

para o PDF da soma.

Essa é uma mistura finita de distribuições Gama com fatores de escala iguais aos da soma e fatores de forma menores ou iguais aos da soma. Exceto em casos especiais (onde pode ocorrer algum cancelamento), o número de termos é dado pelo parâmetro de forma total (assumindo que todos os são diferentes). $n_1 + n_2 + \cdots$ $n_j$

Como teste, eis um histograma de resultados obtidos adicionando desenhos independentes das distribuições e . Sobrepõe-se o gráfico de vezes a função anterior. O ajuste é muito bom. $10^4$ $\Gamma(8,1)$ $\Gamma(4,2)$ $10^4$

Figura

Moschopoulos leva essa idéia um passo adiante ao expandir o cf da soma em uma série infinita de funções características Gamma sempre que um ou mais dos não são integrais e, em seguida, termina a série infinita em um ponto em que esteja razoavelmente bem aproximada. $n_i$

whuber
fonte

2

Comentário secundário : Normalmente, uma mistura finita significa um pdf da forma que e , ou seja, são probabilidades e o pdf podem ser interpretados como a soma ponderada (lei da probabilidade total) de pdfs condicionais, dadas as várias condições que ocorrem com as probabilidades . No entanto, na soma acima, alguns dos coeficientes são negativos e, portanto, a interpretação padrão da mistura não se aplica.

f (x) = \sum_{Eu = 1}^{n} {uma}_{Eu} f_{Eu} (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$

precisa saber é o seguinte

@Dilip Esse é um bom ponto. O que torna este caso interessante é que, embora alguns dos coeficientes possam ser negativos, essa combinação ainda é uma distribuição válida (por sua própria construção).

whuber

Essa abordagem pode ser estendida para incluir a adição de variáveis dependentes? Em particular, quero adicionar 6 distribuições, cada uma tendo alguma correlação com as outras.

masher

11

Mostrarei outra solução possível, que é amplamente aplicável e, com o software R de hoje, bastante fácil de implementar. Essa é a aproximação da densidade do ponto de sela, que deve ser mais conhecida!

Para terminologia sobre a distribuição gama, seguirei https://en.wikipedia.org/wiki/Gamma_distribution com a parametrização shape / scale, é o parâmetro shape e é scale. Para a aproximação do ponto de sela, seguirei Ronald W. Butler: "Aproximações do ponto de sela com aplicações" (Cambridge UP). A aproximação do ponto de sela é explicada aqui: Como funciona a aproximação do ponto de sela? aqui vou mostrar como é usado nesta aplicação. $k$ $\theta$

Seja uma variável aleatória com a função geradora de momento existente que deve existir por em algum intervalo aberto que contenha zero. Em seguida, defina a função geradora cumulante por Sabe-se que . A equação do ponto de sela é que define implicitamente como uma função de (que deve estar no intervalo de ). Nós escrevemos essa função implicitamente definida como . Observe que a equação do ponto de sela sempre tem exatamente uma solução, porque a função cumulante é convexa. $X$

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

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

s

$s$

K (s) = \log 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)$

Então a aproximação do ponto de sela à densidade de é dada por Não é garantido que esta função de densidade aproximada seja integrada a 1, assim como a aproximação não normalizada do ponto de sela. Poderíamos integrá-lo numericamente e renormalizar para obter uma melhor aproximação. Mas essa aproximação é garantida para não ser negativa. $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)$

Agora, sejam variáveis variáveis aleatórias gama independentes, onde possui a distribuição com parâmetros . Então a função geradora cumulante é definida para . A primeira derivada é e a segunda derivada é A seguir, darei algum código para calcular isso e usarei os valores de parâmetro , , $X_1, X_2, \dots, X_n$ $X_i$ $(k_i, \theta_i)$

K (s) = - \sum_{Eu = 1}^{n} k_{Eu} em (1 - θ_{Eu} 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) = \sum_{Eu = 1}^{n} \frac{k_{Eu} θ_{Eu}}{1 - θ_{Eu} s}

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

K^{″} (s) = \sum_{Eu = 1}^{n} \frac{k_{Eu} θ_{Eu}^{2}}{(1 - θ_{Eu} 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)$ . Observe que o Rcódigo a seguir usa um novo argumento na função uniroot introduzida no R 3.1, portanto, não será executado nos R's mais antigos.

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

resultando na seguinte plotagem: insira a descrição da imagem aqui

Vou deixar a aproximação do ponto de sela normalizada como um exercício.

kjetil b halvorsen
fonte

1

Isso é interessante, mas não posso fazer seu Rcódigo funcionar para comparar a aproximação com a resposta exata. Qualquer tentativa de invocação fhatgera erros, aparentemente no uso de uniroot.

whuber

3

Qual é a sua versão R? Os códigos usam um novo argumento para uniroot, extendInt, que foi introduzido no R versão 3.1. Se o seu R for mais antigo, tente removê-lo (e estenda o intervalo dado ao uniroot). Mas isso tornará o código menos robusto!

Kjetil b halvorsen

10

A equação de Welch-Satterthwaite pode ser usada para fornecer uma resposta aproximada na forma de uma distribuição gama. Isso tem a boa propriedade de nos deixar tratar as distribuições gama como sendo (aproximadamente) fechadas em adição. Essa é a aproximação no teste t de Welch comumente usado.

(A distribuição gama pode ser vista como uma distribuição qui-quadrado em escala e permitindo um parâmetro de forma não inteiro).

Eu adaptei a aproximação à parametrização da distribuição gama: $k, \theta$

k_{s você m} = \frac{(\sum_{Eu} θ_{Eu} k_{Eu})^{2}}{\sum_{Eu} θ_{Eu}^{2} k_{Eu}}

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

θ_{s você m} = \frac{\sum θ_{Eu} k_{Eu}}{k_{s você m}}

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

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

Então temos aproximadamente Gamma (10.666 ..., 1.5)

Vemos que o parâmetro de forma foi mais ou menos totalizado, mas um pouco menos porque os parâmetros da escala de entrada diferem. é tal que a soma tem o valor médio correto. $k$ $\theta_i$ $\theta$

Paul Harrison
fonte

6

Uma solução exata para a convolução (isto é, soma) de distribuições gama é dada como Eq. (1) no pdf vinculado por DiSalvo . Como isso é um pouco longo, levará algum tempo para copiá-lo aqui. Para apenas duas distribuições gama, sua soma exata na forma fechada é especificada pela Eq. (2) de DiSalvo e sem pesos pela Eq. (5) de Wesolowski et al. , que também aparece no site do CV como resposta a essa pergunta. Isso é, $n$

G D C (uma, b, α, β; τ) = {\begin{array}{cc} \frac{b^{uma} β^{α}}{Γ (uma + α)} e^{- b τ} {τ^{uma + α}}^{- 1}_{1} F_{1} [α, uma + α, (b - β) τ], & τ > 0 0 \\ 0 0, τ \leq 0 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.,$ onde a notação nas perguntas acima; , aqui. Ou seja, e são constantes de velocidade aqui e não escalares tempo.

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

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

b

$b$

β

$\beta$

Carl
fonte

Soma genérica de variáveis ​​aleatórias gama

Respostas:

Soma genérica de variáveis aleatórias gama