[PYTHON] Zusammenfassung der Spezifikationsmethode für Gammaverteilungsparameter

Gammaverteilung

Die Gammaverteilung wird verwendet, um Werte zu modellieren, die positive kontinuierliche Werte annehmen. Es gibt zwei Parameter, aber die Definitionen unterscheiden sich je nach Bibliothek geringfügig, und es ist schwierig, den Durchschnittswert und die Varianz zu verstehen. Ich habe das Gefühl, dass ich sie immer recherchiere, also werde ich sie zusammenfassen.

Wahrscheinlichkeitsdichtefunktion

\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 Beides sind Parameter mit positiven Werten. Abhängig von der Bibliothek kann dies jedoch mit $ \ lambda = \ frac {1} {\ theta} $ ausgedrückt werden.

Es sieht so aus, wenn es mit einigen Parametern gezeichnet wird.

Statistiken

Durchschnittlich $ \ mu = k \ theta = \ frac {k} {\ lambda} $ Verteilung $ v = k \ theta ^ 2 = \ frac {k} {\ lambda ^ 2} $

Im Gegenteil, wenn Sie die Parameter aus dem Mittelwert und der Varianz bestimmen möchten, verwenden Sie diese. \theta = \frac{v}{\mu} \ (\lambda = \frac{\mu}{v}) k = \frac{\mu^2}{v}

So legen Sie Parameter in jeder Bibliothek fest

Zusammenfassend sieht es so aus. Ich wollte diese Tabelle schreiben!

numpy und scipy verwenden den Skalierungsparameter $ \ theta $, aber in sogenannten probabilistischen Programmiersprachen wie PyMC3 und Stan TFP werden sie als Inverse von $ \ theta $ angegeben.

Die Umkehrung von $ \ theta $, $ \ lambda $, wird als Ratenparameter bezeichnet, und Bibliotheken, die die Parameter $ \ alpha, \ beta $ aufrufen, scheinen die Definition von $ \ lambda $ zu übernehmen.

In PyMC3 ist es auch möglich, die Gammaverteilung durch Mittelwert (mu) und Standardabweichung (Sigma) anzugeben. Außerdem scheint in R entweder Form oder Rate angegeben werden zu können.

Bestätigung der Umsetzung

Um zu bestätigen, ob die obige Parameterliste korrekt ist, habe ich 10.000 Zufallszahlen aus $ Gamma (2, 2) $ in jeder Bibliothek erhalten, die Wahrscheinlichkeitsdichtefunktion geschätzt und verglichen.

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)

Die Ergebnisse sind in der folgenden Abbildung dargestellt, und es wurde bestätigt, dass alle Bibliotheken korrekt implementiert wurden.

Nur Stan wusste nicht, wie man Zufallszahlen direkt aus der Wahrscheinlichkeitsverteilung erhält, also versuchte ich stattdessen, die Parameter der Gammaverteilung aus den oben erhaltenen Zufallszahlen zu schätzen.

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)

Die geschätzten Werte der Form- und Geschwindigkeitsparameter betragen 1,98 bzw. 0,49, was ebenfalls die erwarteten Ergebnisse sind.

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

Bezeichnung der Verteilung durch Mittelwert und Standardabweichung

Es ist zweckmäßig, eine Funktion vorzubereiten, um die Form und den Maßstab der entsprechenden Gammaverteilung aus dem Mittelwert und der Standardabweichung zu berechnen.

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

Gammaverteilung mit einem Durchschnitt von 4 und einer Standardabweichung von 2. Beachten Sie, dass der häufigste Wert (der wahrscheinlichste Wert) kleiner als der Durchschnittswert ist, da die Verteilung rechts einen langen Schwanz hat.

Bonus: Beziehung zu anderen Wahrscheinlichkeitsverteilungen

Wenn $ k = 1 $ ist, stimmt die Gammaverteilung mit der Exponentialverteilung des Parameters $ \ theta $ überein. Wenn $ k = \ frac {n} {2} (n = 1,2, \ dots), \ \ theta = 2 $, stimmt die Gammaverteilung mit der Chi-Quadrat-Verteilung mit $ n $ der Freiheit überein.

Klicken Sie hier für den 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')