[PYTHON] Résumé de la méthode de spécification des paramètres de distribution gamma

Distribution gamma

La distribution gamma est utilisée pour modéliser des valeurs qui prennent des valeurs continues positives. Il y a deux paramètres, mais les définitions sont légèrement différentes selon la bibliothèque, et il est difficile de comprendre la valeur moyenne et la variance, donc j'ai l'impression de toujours les rechercher, donc je vais les résumer.

Fonction de densité de probabilité

\begin{align}
f(x) &= \frac{1}{\Gamma(k) \theta^k} x^{k-1} e^{- \frac{x}{\theta}} \\
&= \frac{\lambda^k}{\Gamma(k)} x^{k-1} e^{- \lambda x}
\end{align}

k: shape parameter \theta: scale parameter Les deux sont des paramètres avec des valeurs positives. Cependant, selon la bibliothèque, il peut être exprimé en utilisant $ \ lambda = \ frac {1} {\ theta} $.

Cela ressemble à ceci lorsqu'il est tracé avec certains paramètres.

Statistiques

Moyenne $ \ mu = k \ theta = \ frac {k} {\ lambda} $ Distribution $ v = k \ theta ^ 2 = \ frac {k} {\ lambda ^ 2} $

Au contraire, si vous souhaitez déterminer les paramètres à partir de la moyenne et de la variance, utilisez ceci. \theta = \frac{v}{\mu} \ (\lambda = \frac{\mu}{v}) k = \frac{\mu^2}{v}

Comment spécifier des paramètres dans chaque bibliothèque

En résumé, ça ressemble à ça. Je voulais écrire ce tableau!

numpy et scipy utilisent le paramètre d'échelle $ \ theta $, mais dans les langages de programmation dits probabilistes tels que PyMC3 et Stan TFP, ils sont spécifiés comme l'inverse de $ \ theta $.

L'inverse de $ \ theta $, $ \ lambda $, est appelé paramètre de taux, et les bibliothèques qui appellent les paramètres $ \ alpha, \ beta $ semblent adopter la définition de $ \ lambda $.

Dans PyMC3, il est également possible de spécifier la distribution gamma par la moyenne (mu) et l'écart type (sigma). En outre, dans R, il semble que la forme ou le taux puissent être spécifiés.

Confirmation de la mise en œuvre

Afin de confirmer si la liste de paramètres ci-dessus est correcte, j'ai obtenu 10000 nombres aléatoires de $ Gamma (2, 2) $ dans chaque bibliothèque, estimé la fonction de densité de probabilité et les ai comparés.

import numpy as np
import scipy as sp
import pymc3 as pm
import tensorflow_probability as tfp
import pystan
import matplotlib.pyplot as plt
import seaborn as sns

shape = 2
scale = 2
rate = 1 / scale
n_sample = 10000

xx = np.linspace(0, 20)
# ground truth
gamma_pdf = sp.stats.gamma(a=shape, scale=scale).pdf

s_np = np.random.gamma(shape=shape, scale=scale, size=n_sample)
s_sp = sp.stats.gamma(a=shape, scale=scale).rvs(size=n_sample)
s_tfp = tfp.distributions.Gamma(concentration=shape, rate=rate).sample(n_sample).numpy()
s_pm = pm.Gamma.dist(alpha=shape, beta=rate).random(size=n_sample)

fig, ax = plt.subplots()
ax.plot(xx, gamma_pdf(xx), label='truth', lw=2, c='k')
sns.kdeplot(s_np, ax=ax, label='numpy', alpha=0.8)
sns.kdeplot(s_sp, ax=ax, label='scipy', alpha=0.8)
sns.kdeplot(s_tfp, ax=ax, label='TFP', alpha=0.8)
sns.kdeplot(s_pm, ax=ax, label='PyMC3', alpha=0.8)

Les résultats sont présentés dans la figure ci-dessous et il a été confirmé que toutes les bibliothèques étaient correctement implémentées.

Seul Stan ne savait pas comment obtenir des nombres aléatoires directement à partir de la distribution de probabilité, alors j'ai plutôt essayé d'estimer les paramètres de la distribution gamma à partir des nombres aléatoires obtenus ci-dessus.

stan_code = '''

data {
  int N;
  vector<lower=0>[N] Y;
}

parameters {
  real<lower=0> shape;
  real<lower=0> rate;
}

model {
  Y ~ gamma(shape, rate);
}

'''

data = dict(N=n_sample, Y=s_np)
stan_model = pystan.StanModel(model_code=stan_code)
fit = stan_model.sampling(data=data)
print(fit)

Les valeurs estimées des paramètres de forme et de débit sont respectivement de 1,98 et 0,49, qui sont également les résultats attendus.

Inference for Stan model: anon_model_6a5d60bed963727c801dc434b96a49a1.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
shape   1.98  8.3e-4   0.03   1.93   1.97   1.98    2.0   2.04   1016    1.0
rate    0.49  2.3e-4 7.4e-3   0.47   0.48   0.49   0.49    0.5   1020    1.0
lp__  -2.3e4    0.03   1.01 -2.3e4 -2.3e4 -2.3e4 -2.3e4 -2.3e4   1192    1.0

Désignation de la distribution par moyenne et écart type

Il est pratique de préparer une fonction pour calculer la forme et l'échelle de la distribution gamma correspondante à partir de la moyenne et de l'écart type.

def calc_gamma_param(mu, sigma):
    return (mu / sigma)**2, sigma**2 / mu

mu, sigma = 4, 2
shape, scale = calc_gamma_param(mu, sigma)

def plot_gamma(xx, shape, scale):
    plt.plot(xx, sp.stats.gamma(a=shape, scale=scale).pdf(xx), label=f'shape={shape}, scale={scale}')
    
xx = np.linspace(0, 10)
plot_gamma(xx, shape, scale)
plt.legend()

Distribution gamma avec une moyenne de 4 et un écart type de 2. Notez que la valeur la plus fréquente (la valeur la plus probable) est inférieure à la valeur moyenne car la distribution a une longue queue vers la droite.

Bonus: relation avec d'autres distributions de probabilité

Lorsque $ k = 1 $, la distribution gamma correspond à la distribution exponentielle du paramètre $ \ theta $. Quand $ k = \ frac {n} {2} (n = 1,2, \ points), \ \ theta = 2 $, la distribution gamma correspond à la distribution du chi carré avec $ n $ de liberté.

Cliquez-ici pour le code.

from scipy import stats

xx = np.linspace(0, 10)
fig, ax = plt.subplots(1, 2, figsize=(10, 3))

shape, scale = 1, 3
ax[0].plot(xx, stats.gamma(a=shape, scale=scale).pdf(xx), label=f'Gamma({shape}, {scale})')
ax[0].plot(xx, stats.expon(scale=scale).pdf(xx), label=f'Exp({scale})')
ax[0].legend()

shape, scale = 3/2, 2 
ax[1].plot(xx, stats.gamma(a=shape, scale=scale).pdf(xx), label=f'Gamma({shape}, {scale})')
ax[1].plot(xx, stats.chi2(df=2*shape).pdf(xx), label=f'Chi2({int(2*shape)})')
ax[1].legend()

plt.savefig('gamma_exp_chi2.png')