2021 Q3(1)
ポアソン分布のモーメント母関数を求めました。
![](https://statistics.blue/wp-content/uploads/2024/07/q_note_2021_2_ページ_8-728x1024.jpg)
コード
数式を使った計算
# 2021 Q3(1) 2024.8.25
import sympy as sp
# 変数の定義
s, x, lambda_ = sp.symbols('s x lambda_', real=True, positive=True)
# ポアソン分布の確率質量関数 (PMF)
poisson_pmf = (lambda_**x * sp.exp(-lambda_)) / sp.factorial(x)
# モーメント母関数 M_X(s) の定義
M_X = sp.summation(sp.exp(s*x) * poisson_pmf, (x, 0, sp.oo))
# 結果を簡略化
M_X_simplified = sp.simplify(M_X)
# 結果を表示
M_X_simplified
![](https://statistics.blue/wp-content/uploads/2024/02/image-11.png)
無限級数がe^xの形に変換されないようです。
exp(s) を t に置き換えて、最後に戻してみます。
# 2021 Q3(1) 2024.8.25
import sympy as sp
# 変数の定義
s, x, lambda_, t = sp.symbols('s x lambda_ t', real=True, positive=True)
# exp(s) を t に置き換え
summand = (lambda_ * t)**x / sp.factorial(x)
# モーメント母関数 M_X(s) の定義
M_X = sp.exp(-lambda_) * sp.summation(summand, (x, 0, sp.oo))
# 簡略化
M_X_simplified = sp.simplify(M_X)
# 最後に t を exp(s) に戻す
M_X_final = M_X_simplified.subs(t, sp.exp(s))
# 結果を表示
M_X_final
![](https://statistics.blue/wp-content/uploads/2024/02/image-13.png)
手計算と同じ形になりました。
モーメント母関数を使って期待値と分散を求めます。
import sympy as sp
# 変数の定義
s, lambda_ = sp.symbols('s lambda_', real=True, positive=True)
# モーメント母関数 M_X(s) の定義
M_X = sp.exp(lambda_ * (sp.exp(s) - 1))
# 1次モーメント(期待値)の計算
M_X_prime = sp.diff(M_X, s)
expectation = M_X_prime.subs(s, 0)
# 2次モーメントの計算
M_X_double_prime = sp.diff(M_X_prime, s)
second_moment = M_X_double_prime.subs(s, 0)
# 分散の計算
variance = second_moment - expectation**2
# 結果を表示
display(expectation, variance)
![](https://statistics.blue/wp-content/uploads/2024/02/image-14.png)
2021 Q2(4)
事前分布と尤度から事後確率が最大となるパラメータの推定値を求めました。
![](https://statistics.blue/wp-content/uploads/2024/07/q_note_2021_2_ページ_7-728x1024.jpg)
コード
数式を使った計算
# 2021 Q2(4) 2024.8.24
import numpy as np
from scipy.special import comb
# 事前分布 P(N_A)
def prior(NA):
return NA + 1
# 尤度関数 P(X = 4 | N_A)
def likelihood(NA):
return comb(NA, 4) * comb(100 - NA, 11) / comb(100, 15)
# 事後分布 P(N_A | X = 4) (正規化定数は省略)
def posterior(NA):
return prior(NA) * likelihood(NA)
# N_A の範囲
NA_values = np.arange(4, 90)
# 事後分布の計算
posterior_values = [posterior(NA) for NA in NA_values]
# 最大値を取る N_A (事後モード) を探索
NA_mode = NA_values[np.argmax(posterior_values)]
NA_mode
30
プロット
# 2021 Q2(4) 2024.8.24
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb
# 事前分布 P(N_A)
def prior(NA):
return NA + 1
# 正規化定数
C = 1 / 5151
# 尤度関数 P(X = 4 | N_A)
def likelihood(NA):
return comb(NA, 4) * comb(100 - NA, 11) / comb(100, 15)
# 事後分布 P(N_A | X = 4)
def posterior(NA):
return prior(NA) * likelihood(NA)
# N_A の範囲を0から100まで(事前分布用)
NA_values_full = np.arange(0, 101)
# 事前分布を計算(0から100まで)
prior_values_full = [C * prior(NA) for NA in NA_values_full]
# N_A の範囲を4から89まで(事後分布用)
NA_values = np.arange(4, 90)
# 事後分布の計算
posterior_values = [posterior(NA) for NA in NA_values]
# 事後分布の正規化
posterior_sum = sum(posterior_values)
posterior_values_normalized = [value / posterior_sum for value in posterior_values]
# 事前分布が存在しない場合の事後分布(尤度のみを正規化)
likelihood_values = [likelihood(NA) for NA in NA_values]
likelihood_sum = sum(likelihood_values)
likelihood_values_normalized = [value / likelihood_sum for value in likelihood_values]
# 3つの分布を重ねて表示
plt.figure(figsize=(10, 6))
# 事前分布のプロット(0から100まで)
plt.plot(NA_values_full, prior_values_full, 'g-', label=r'事前分布 $P(N_A)$', markersize=5)
# 事後分布のプロット(4から89まで)
plt.plot(NA_values, posterior_values_normalized, 'bo-', label=r'事後分布 $P(N_A | X = 4)$', markersize=5)
# 事前分布が存在しない場合の事後分布(尤度のみ)
plt.plot(NA_values, likelihood_values_normalized, 'r-', label=r'事前分布なし $P(N_A | X = 4)$', markersize=5)
# 事後モードのプロット
NA_mode = NA_values[np.argmax(posterior_values_normalized)]
plt.axvline(NA_mode, color='blue', linestyle='--', label=f'事後モード: {NA_mode}')
# 事前分布なしの場合の最尤推定値のプロット
NA_mode_likelihood = NA_values[np.argmax(likelihood_values_normalized)]
plt.axvline(NA_mode_likelihood, color='red', linestyle=':', label=f'最尤推定値 (事前分布なし): {NA_mode_likelihood}')
# グラフの設定
plt.xlabel(r'$N_A$')
plt.ylabel(r'確率')
plt.title(r'事前分布ありとなしの比較')
plt.grid(True)
plt.legend()
plt.show()
![](https://statistics.blue/wp-content/uploads/2024/02/Untitled-14.png)