ホーム » 分布 » ガンマ分布

ガンマ分布」カテゴリーアーカイブ

2019 Q2(3)

(2)で求めた確率密度関数において確率変数Uの逆数の期待値を求めました。

 

コード

数式を使って1/Uの期待値を求めます。

# 2019 Q2(3)  2024.9.15

import sympy as sp

# 変数の定義
u, lambda_value = sp.symbols('u lambda', positive=True, real=True)

# ガンマ分布の確率密度関数 g(u) = λ^2 * u * exp(-λ * u)
g_u = lambda_value**2 * u * sp.exp(-lambda_value * u)

# 1/U の期待値の積分
integrand = (1/u) * g_u  # 1/u * g(u)

# 積分の実行
expectation = sp.integrate(integrand, (u, 0, sp.oo))

# 結果の簡略化
expectation_simplified = sp.simplify(expectation)

# 結果の表示
display(expectation_simplified)

手計算と一致します。

次に、数値シミュレーションを行います。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
lambda_value = 2  # λ=2
num_samples = 100000  # サンプルの数

# 指数分布に従う乱数 X1, X2 を生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)

# U = X1 + X2 の計算
U = X1 + X2

# 1/U の計算
inv_U = 1 / U

# シミュレーションの結果からの平均値
simulation_result = np.mean(inv_U)

# 理論値 E[1/U] = λ
theoretical_value = lambda_value

# 結果の表示
print(f'シミュレーションによる E[1/U]: {simulation_result}')
print(f'理論値 E[1/U]: {theoretical_value}')

# ヒストグラムの描画
plt.hist(inv_U, bins=50, density=True, alpha=0.7, label='シミュレーション結果')
plt.axvline(theoretical_value, color='r', linestyle='--', label=f'理論値 E[1/U] = {theoretical_value}')
plt.xlabel('1/U')
plt.ylabel('確率密度')
plt.title('シミュレーションによる 1/U の分布と理論値の比較')
plt.legend()
plt.show()
シミュレーションによる E[1/U]: 1.9928303711222939
理論値 E[1/U]: 2

シミュレーションによりE[1/U]は理論値に近づきました。分散は発散するためヒストグラムにすると大きなバラつきが見られます。

2019 Q2(2)

独立な指数分布に従う2変数の和の確率密度関数を求めました。

 

コード

数式を使ってUの確率密度関数を求めます。

# 2019 Q2(2)  2024.9.14

import sympy as sp

# 変数の定義
x, u, lambda_value = sp.symbols('x u lambda', positive=True, real=True)

# 指数分布の確率密度関数 (PDF) f(x)
f_x = lambda_value * sp.exp(-lambda_value * x)

# 畳み込み積分を計算して、g(u) = f(x) * f(u - x) の形式にする
g_u = sp.integrate(f_x * lambda_value * sp.exp(-lambda_value * (u - x)), (x, 0, u))

# 簡略化
g_u_simplified = sp.simplify(g_u)

# 結果の表示
display(g_u_simplified)

X1 (X2)と、Uの確率密度関数を重ねて描画してみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

# パラメータ設定
lambda_value = 2  # λ=2
x_range = np.linspace(0, 5, 1000)  # X1, X2 の範囲

# X1, X2 の指数分布の確率密度関数 (PDF) - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)

# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
pdf_U = gamma.pdf(u_range, a=2, scale=1/lambda_value)

# グラフのプロット
plt.figure(figsize=(8, 6))

# X1 (X2) の指数分布の描画
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')

# U のガンマ分布の描画
plt.plot(u_range, pdf_U, label='U = X1 + X2 (ガンマ分布)', color='red')

# グラフの装飾
plt.title('X1 (X2) と U の確率密度関数の比較')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

# グラフの表示
plt.show()

グラフの形状から、U = X1 + X2の関係は想像しづらいですね。

ガンマ分布の形状パラメータを1~2に変化させて、X1 (X2)からUへの変化を確認します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

# パラメータ設定
lambda_value = 2  # λ=2
x_range = np.linspace(0, 5, 1000)  # X1, X2, U の範囲を0~5に設定

# ガンマ分布の形状パラメータを1から2まで0.1ずつ変化させる
shape_params = np.arange(1, 2.1, 0.1)  # 1から2までの形状パラメータを0.1ステップで変化

# グラフのプロット
plt.figure(figsize=(8, 6))

# 形状パラメータが1から2に変化するガンマ分布の描画
for k in shape_params:
    gamma_pdf = gamma.pdf(x_range, a=k, scale=1/lambda_value)
    plt.plot(x_range, gamma_pdf, label=f'形状パラメータ k={k:.1f}')

# X1 (X2) の指数分布の描画 - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')

# グラフの装飾
plt.title('指数分布からガンマ分布への変化 (形状パラメータの0.1ステップ変化)')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

# グラフの表示
plt.show()

X1 (X2)からUへの変化が可視化できました。

次に、X1 (X2)と、Uの累積分布関数を重ねて描画してみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import expon, gamma

# パラメータ設定
lambda_value = 2  # λ=2
x_range = np.linspace(0, 5, 1000)  # X1, X2 の範囲

# X1, X2 の指数分布の累積分布関数 (CDF) - 同じなので1つにまとめる
cdf_X1_X2 = expon.cdf(x_range, scale=1/lambda_value)

# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
cdf_U = gamma.cdf(u_range, a=2, scale=1/lambda_value)

# グラフのプロット
plt.figure(figsize=(8, 6))

# X1 (X2) の指数分布のCDFの描画
plt.plot(x_range, cdf_X1_X2, label='X1 (X2) (指数分布 CDF)', color='blue', linestyle='--')

# U のガンマ分布のCDFの描画
plt.plot(u_range, cdf_U, label='U = X1 + X2 (ガンマ分布 CDF)', color='red')

# グラフの装飾
plt.title('X1 (X2) と U の累積分布関数 (CDF) の比較')
plt.xlabel('値')
plt.ylabel('累積確率')
plt.legend()

# グラフの表示
plt.show()

グラフの形状から、U = X1 + X2の関係を想像しやすくなりました。

2019 Q2(1)

独立な指数分布に従う2変数の和の期待値を求めました。

 

コード

数式を使ってUの期待値を求めます

import sympy as sp

# 変数の定義
x = sp.Symbol('x', positive=True)
lambda_value = sp.Symbol('lambda', positive=True)

# 指数分布の確率密度関数 (PDF)
f_x = lambda_value * sp.exp(-lambda_value * x)

# 期待値の式を計算 (積分)
E_X1 = sp.integrate(x * f_x, (x, 0, sp.oo))

# 期待値
E_X1_simplified = sp.simplify(E_X1)

# X1 + X2 の場合 (独立性より)
E_U = 2 * E_X1_simplified

# 結果の表示
display(E_U)

手計算と一致します

次にUに従う乱数を生成し分布を確認します。また形状パラメータ2のガンマ分布と重ねて描画します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma

# パラメータ設定
lambda_value = 2  # 例としてλ=2
num_samples = 1000  # サンプルの数

# 指数分布に従う乱数の生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)

# U = X1 + X2 の計算
U = X1 + X2

# 理論的な期待値
expected_value = 2 / lambda_value

# シミュレーション結果のヒストグラムのプロット (確率密度にスケール)
plt.hist(U, bins=30, density=True, alpha=0.7, label='シミュレーション結果')

# ガンマ分布 (形状パラメータ k=2, スケールパラメータ θ=1/λ) の確率密度関数 (PDF) を計算
k = 2  # 形状パラメータ (k = 2)
theta = 1 / lambda_value  # スケールパラメータ (θ = 1/λ)

# ガンマ分布に基づくx軸の範囲
x = np.linspace(0, max(U), 1000)

# ガンマ分布の確率密度関数 (PDF) を計算
gamma_pdf = gamma.pdf(x, a=k, scale=theta)

# ガンマ分布の折れ線グラフをプロット
plt.plot(x, gamma_pdf, 'r-', label=f'ガンマ分布 (k={k}, θ={theta:.2f})')

# グラフの装飾
plt.axvline(expected_value, color='g', linestyle='--', label=f'理論値 = {expected_value:.2f}')
plt.xlabel('U = X1 + X2')
plt.ylabel('確率密度')
plt.title('Uの分布とガンマ分布の比較')
plt.legend()
plt.show()

# 実際のサンプル平均値と理論的期待値の比較
print(f'理論的な期待値: {expected_value}')
print(f'シミュレーションからの平均値: {np.mean(U)}')
理論的な期待値: 1.0
シミュレーションからの平均値: 0.9902440726761544

Uの分布と形状パラメータ2のガンマ分布は重なりました。Uはガンマ分布に従うことが分かります。

2022 Q3(5)

モーメント法によりガンマポアソン分布のパラメータの推定量を求めました。

コード

数式を使った計算

# 2022 Q3(5)  2024.8.6

from sympy import init_printing
from sympy import symbols, Eq, solve

# これでTexの表示ができる
init_printing()

# 変数の定義
r, p, X_bar, S_squared = symbols('r p X_bar S_squared')
alpha, beta = symbols('alpha beta')

# モーメント法の方程式を定義
# 期待値 E[X] にサンプル平均を一致させる
eq1_correct = Eq(X_bar, alpha / beta)
# 分散 V[X] にサンプル分散を一致させる
eq2_correct = Eq(S_squared, alpha / beta**2 + alpha / beta)

# 方程式を解いて α と β を求める
solutions_correct = solve((eq1_correct, eq2_correct), (alpha, beta))
solutions_correct

シミュレーションによる計算

# 2022 Q3(5)  2024.8.6

import numpy as np
from scipy.stats import nbinom

# 1. パラメータの設定
alpha_true = 3.0
beta_true = 2.0

# 2. 負の二項分布のパラメータ計算
r = alpha_true
p = beta_true / (beta_true + 1)

# 3. 負の二項分布から乱数を生成
sample_size = 10000
nbinom_samples = nbinom.rvs(r, p, size=sample_size)

# 4. サンプル統計量の計算
X_bar = np.mean(nbinom_samples)
S_squared = np.var(nbinom_samples, ddof=1)

# 5. モーメント法によるパラメータ推定
alpha_estimated = X_bar**2 / (S_squared - X_bar)
beta_estimated = X_bar / (S_squared - X_bar)

# 6. 推定結果の表示
print(f"推定された α: {alpha_estimated}")
print(f"推定された β: {beta_estimated}")
print(f"理論値 α: {alpha_true}")
print(f"理論値 β: {beta_true}")
推定された α: 2.93536458612626
推定された β: 1.95157541794180
理論値 α: 3.0
理論値 β: 2.0

アルゴリズム

シミュレーションによる計算

パラメータの設定:

  • 真のパラメータ \alpha\beta を設定する。

負の二項分布のパラメータ計算:

  • パラメータ r = \alpha
  • パラメータ p = \frac{\beta}{\beta + 1}

負の二項分布から乱数を生成:

  • サンプルサイズ n を設定し、負の二項分布から乱数を生成する。

サンプル統計量の計算:

  • サンプル平均 \bar{X} と分散 S^2 を計算する。

モーメント法によるパラメータ推定:

  • 推定された \alpha\alpha = \frac{\bar{X}^2}{S^2 - \bar{X}} で計算する。
  • 推定された \beta\beta = \frac{\bar{X}}{S^2 - \bar{X}} で計算する。

推定結果の表示:

  • 推定された \alpha\beta を表示する。
  • 真の \alpha\beta を表示する。

2022 Q3(4)

条件付き期待値の公式と条件付き分散の公式を使ってガンマポアソン分布の期待値と分散を求めました。

コード

数式を使った計算

# 2022 Q3(4)  2024.8.7

import sympy as sp

# 定義
k, lambda_var = sp.symbols('k lambda')
alpha, beta = sp.symbols('alpha beta', positive=True)

# ポアソン分布の確率質量関数
poisson_pmf = (lambda_var**k * sp.exp(-lambda_var)) / sp.factorial(k)

# ガンマ分布の確率密度関数
gamma_pdf = (beta**alpha / sp.gamma(alpha)) * lambda_var**(alpha - 1) * sp.exp(-beta * lambda_var)

# 周辺分布の計算
marginal_distribution = sp.integrate(poisson_pmf * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

# 期待値 E[X] の計算
expected_value = sp.summation(k * marginal_distribution, (k, 0, sp.oo)).simplify()

# 分散 V[X] の計算
expected_value_of_square = sp.summation(k**2 * marginal_distribution, (k, 0, sp.oo)).simplify()
variance = (expected_value_of_square - expected_value**2).simplify()

# 結果を表示
results = {
    "期待値 E[X]": expected_value,
    "分散 V[X]": variance
}

results
{'期待値 E[X]': alpha/beta, '分散 V[X]': alpha*(beta + 1)/beta**2}

シミュレーションによる計算

import numpy as np

# パラメータの設定
alpha = 3.0
beta = 2.0
sample_size = 10000

# ガンマ分布から λ を生成
lambda_samples = np.random.gamma(alpha, 1/beta, sample_size)

# 生成された λ を使ってポアソン分布から X を生成
poisson_samples = [np.random.poisson(lam) for lam in lambda_samples]

# サンプルの期待値と分散を計算
sample_mean = np.mean(poisson_samples)
sample_variance = np.var(poisson_samples)

# 理論値を計算
theoretical_mean = alpha / beta
theoretical_variance = alpha * (beta + 1) / beta**2

# 結果を表示
results_with_caption = {
    "シミュレーションによる期待値": sample_mean,
    "理論値による期待値": theoretical_mean,
    "シミュレーションによる分散": sample_variance,
    "理論値による分散": theoretical_variance
}

results_with_caption
{'シミュレーションによる期待値': 1.5027,
 '理論値による期待値': 1.5,
 'シミュレーションによる分散': 2.31159271,
 '理論値による分散': 2.25}

2022 Q3(3)

ポアソン分布のパラメータλがガンマ分布に従うと負の二項分布になる。

コード

数式を使った計算

# 2022 Q3(3)  2024.8.5

import sympy as sp

# 定義
k = sp.symbols('k', integer=True)
alpha, beta = sp.symbols('alpha beta', positive=True)
lambda_var = sp.symbols('lambda', positive=True)

# ポアソン分布の確率質量関数
poisson_pmf = (lambda_var**k * sp.exp(-lambda_var)) / sp.factorial(k)

# ガンマ分布の確率密度関数
gamma_pdf = (beta**alpha / sp.gamma(alpha)) * lambda_var**(alpha - 1) * sp.exp(-beta * lambda_var)

# 周辺分布の計算
marginal_distribution = sp.integrate(poisson_pmf * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

marginal_distribution

次式と同じなので、負の二項分布となる。

シミュレーションによる計算

# 2022 Q3(3)  2024.8.5

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma, nbinom

# パラメータの設定
alpha = 3.0
beta = 2.0
sample_size = 10000

# ガンマ分布から λ を生成
lambda_samples = np.random.gamma(alpha, 1/beta, sample_size)

# 生成された λ を使ってポアソン分布から X を生成
poisson_samples = [np.random.poisson(lam) for lam in lambda_samples]

# 負の二項分布のパラメータを設定
r = alpha
p = beta / (beta + 1)

# 負の二項分布から理論的な確率質量関数を計算
x = np.arange(0, max(poisson_samples) + 1)
nbinom_pmf = nbinom.pmf(x, r, p)

# ヒストグラムの描画
plt.hist(poisson_samples, bins=np.arange(0, max(poisson_samples) + 1) - 0.5, density=True, alpha=0.75, color='blue', edgecolor='black', rwidth=0.8)

# 理論的な負の二項分布の確率質量関数をプロット
plt.plot(x, nbinom_pmf, 'r', linestyle='-', label='理論的な負の二項分布')

# グラフのタイトルとラベル
plt.title('ポアソン-ガンマ混合分布のシミュレーション結果')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

# グラフの表示
plt.grid(True)
plt.show()