モーメント法によりガンマポアソン分布のパラメータの推定量を求めました。
コード
数式を使った計算
# 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
アルゴリズム
シミュレーションによる計算
パラメータの設定:
- 真のパラメータ と を設定する。
負の二項分布のパラメータ計算:
- パラメータ
- パラメータ
負の二項分布から乱数を生成:
- サンプルサイズ を設定し、負の二項分布から乱数を生成する。
サンプル統計量の計算:
- サンプル平均 と分散 を計算する。
モーメント法によるパラメータ推定:
- 推定された を で計算する。
- 推定された を で計算する。
推定結果の表示:
- 推定された と を表示する。
- 真の と を表示する。