ホーム » 2024 » 1月 (ページ 3)

月別アーカイブ: 1月 2024

2021 Q3(4)

ポアソン分布のパラメータλの区間推定の中点がその最尤推定値より大きくなることを学びました。面白い。

コード

ポアソン分布の確率質量関数にλに3,5,10を与えてプロットします

#2021 Q3(4)  2024.8.30

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

# パラメータ設定
lambda_values = [3, 5, 10]  # 異なる λ の値
max_x = 25  # x軸の最大値(PMFの表示範囲)

# プロットの設定
plt.figure(figsize=(12, 6))

# 各 λ について PMF を計算し、プロット
for lambda_value in lambda_values:
    x = np.arange(0, max_x + 1)  # x軸の範囲
    pmf = poisson.pmf(x, lambda_value)  # PMF を計算
    plt.plot(x, pmf, marker='o', label=f'λ = {lambda_value}')

# グラフの装飾
plt.title('ポアソン分布の確率質量関数(PMF)の比較')
plt.xlabel('x')
plt.ylabel('P(X=x)')
plt.legend()
plt.grid(True)
plt.show()

λが大きくなるほど分散も大きくなります

λの信頼区間と最尤推定値の関係を見てみます

#2021 Q3(4)  2024.8.30

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
lambda_true = 5  # 真の λ の値
n = 10  # サンプル数
num_simulations = 1000  # シミュレーション回数

# 信頼区間と最尤推定値の保存リスト
lambda_L_values = []
lambda_U_values = []
lambda_hat_values = []

# シミュレーション開始
for _ in range(num_simulations):
    # ポアソン分布に従う n 個のサンプルを生成
    samples = np.random.poisson(lambda_true, n)
    T = np.sum(samples)  # T を計算
    
    # 最尤推定値を計算
    lambda_hat = T / n
    
    # 信頼区間の計算
    lambda_L = (T - 2*np.sqrt(T + 1) + 2) / n
    lambda_U = (T + 2*np.sqrt(T + 1) + 2) / n
    
    # 結果を保存
    lambda_L_values.append(lambda_L)
    lambda_U_values.append(lambda_U)
    lambda_hat_values.append(lambda_hat)

# 信頼区間の中点の計算
lambda_M_values = [(lambda_L_values[i] + lambda_U_values[i]) / 2 for i in range(num_simulations)]
lambda_M_mean = np.mean(lambda_M_values)

# 最尤推定値の平均を計算
lambda_hat_mean = np.mean(lambda_hat_values)

# 信頼区間のプロット
plt.figure(figsize=(12, 6))
for i in range(num_simulations):
    plt.plot([lambda_L_values[i], lambda_U_values[i]], [i, i], color='blue')
plt.axvline(lambda_hat_mean, color='purple', linestyle='--', label='最尤推定値の平均')
plt.axvline(lambda_M_mean, color='orange', linestyle='-', label='信頼区間の中点の平均')
plt.title('信頼区間と最尤推定値の関係')
plt.xlabel('λ の値')
plt.ylabel('シミュレーション回数')
plt.legend()
plt.grid(True)
plt.show()

# 中点の平均値を表示
print(f'信頼区間の中点の平均: {lambda_M_mean:.4f}')
print(f'最尤推定値の平均: {lambda_hat_mean:.4f}')
print(f'真の λ: {lambda_true}')
信頼区間の中点の平均: 5.2235
最尤推定値の平均: 5.0235
真の λ: 5

λの信頼区間の中点は最尤推定値よりも右にずれます。λが大きくなると分散が大きくなるためです。

X=5に固定したポアソン分布の確率質量関数P(X=5|λ)のグラフを見てみます

#2021 Q3(4)  2024.8.30

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

# パラメータ設定
x_value = 5  # 固定された x の値
n = 10  # サンプル数(ここでは仮に10個の観測データがあると仮定)
T = x_value * n  # 合計値 T は x * n として計算

# λの範囲を指定
lambda_range = np.linspace(1, 20, 200)

# λに対するPMFを計算
pmf_values = poisson.pmf(x_value, lambda_range)

# 最尤推定値(λ_hat)を計算
lambda_hat = T / n

# 信頼区間の計算
lambda_L = (T - 2*np.sqrt(T + 1) + 2) / n
lambda_U = (T + 2*np.sqrt(T + 1) + 2) / n
lambda_M = (lambda_L + lambda_U) / 2  # 信頼区間の中点

# グラフの描画
plt.figure(figsize=(12, 6))
plt.plot(lambda_range, pmf_values, label=f'P(X={x_value})')

# 最尤推定値、信頼区間の上下限、中点を表示
plt.axvline(lambda_hat, color='red', linestyle='--', label=f'最尤推定値 λ_hat = {lambda_hat:.2f}')
plt.axvline(lambda_L, color='blue', linestyle='--', label=f'信頼区間下限 λ_L = {lambda_L:.2f}')
plt.axvline(lambda_U, color='green', linestyle='--', label=f'信頼区間上限 λ_U = {lambda_U:.2f}')
plt.axvline(lambda_M, color='orange', linestyle='-', label=f'信頼区間中点 λ_M = {lambda_M:.2f}')

# グラフの装飾
plt.title(f'ポアソン分布のPMFと信頼区間の中点(x={x_value}固定)')
plt.xlabel('λ')
plt.ylabel(f'P(X={x_value})')
plt.grid(True)
plt.legend()
plt.show()

右に裾が長く、λの信頼区間の中点は最尤推定値よりも右にズレていることが分かります。

なお、このグラフは、ボソン分布の確率質量関数

より

となります。これは、k=6 , θ=1 としたときガンマ分布に一致します

2021 Q3(3)

ポアソン分布に従う複数の確率変数の和が正規分布に近似するとした場合の信頼区間を求める問題を学びました。

コード

二次方程式を解いて信頼区間を求めます

# 2021 Q3(3)  2024.8.29

import sympy as sp

# 変数の定義
n, lambda_, T = sp.symbols('n lambda_ T', real=True, positive=True)
c = 2

# 不等式を通常の方程式として変換
# まず、二次方程式の形を導出
lhs = sp.Eq((T - n*lambda_)**2, (c**2) * n * lambda_)

# λに関する二次方程式を解く
solutions = sp.solve(lhs, lambda_)

# 解の表示
lambda_L = solutions[0].simplify()
lambda_U = solutions[1].simplify()

# 結果の表示
display(lambda_L, lambda_U)

シミュレーションで、真のλが信頼区間に含まれる割合を求めます

# 2021 Q3(3)  2024.8.29

import numpy as np

# パラメータの設定
lambda_true = 5  # 真の λ の値
n = 10  # サンプル数
num_simulations = 1000  # シミュレーションの回数

# 結果を保存するカウンター
within_confidence = 0

# シミュレーション開始
for _ in range(num_simulations):
    # ポアソン分布に従う n 個のサンプルを生成
    samples = np.random.poisson(lambda_true, n)
    T = np.sum(samples)  # T を計算
    
    # 信頼区間の計算
    lambda_L = (T - 2*np.sqrt(T + 1) + 2) / n
    lambda_U = (T + 2*np.sqrt(T + 1) + 2) / n
    
    # 真のλが信頼区間内にあるか確認
    if lambda_L <= lambda_true <= lambda_U:
        within_confidence += 1

# 信頼区間内に真のλが含まれる割合を計算
coverage = within_confidence / num_simulations

# 結果の表示
print(f"{num_simulations}回のシミュレーションのうち、真のλが信頼区間内に含まれていた割合: {coverage:.4f}")
1000回のシミュレーションのうち、真のλが信頼区間内に含まれていた割合: 0.9550

c=2のとき信頼区間の信頼率は95.45%となり、シミュレーションと概ね一致します

上のシミュレーションを視覚化します

# 2021 Q3(3)  2024.8.29

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
lambda_true = 5  # 真の λ の値
n = 10  # サンプル数
num_simulations = 1000  # シミュレーションの回数

# 結果を保存するリスト
lower_limits = []
upper_limits = []
within_confidence = 0

# シミュレーション開始
for _ in range(num_simulations):
    # ポアソン分布に従う n 個のサンプルを生成
    samples = np.random.poisson(lambda_true, n)
    T = np.sum(samples)  # T を計算
    
    # 信頼区間の計算
    lambda_L = (T - 2*np.sqrt(T + 1) + 2) / n
    lambda_U = (T + 2*np.sqrt(T + 1) + 2) / n
    
    # 信頼区間の記録
    lower_limits.append(lambda_L)
    upper_limits.append(lambda_U)
    
    # 真のλが信頼区間内にあるか確認
    if lambda_L <= lambda_true <= lambda_U:
        within_confidence += 1

# 信頼区間内に真のλが含まれる割合を計算
coverage = within_confidence / num_simulations
print(f"{num_simulations}回のシミュレーションのうち、真のλが信頼区間内に含まれていた割合: {coverage:.4f}")

# ヒストグラムのプロット
plt.figure(figsize=(12, 6))
plt.hist(lower_limits, bins=30, alpha=0.5, label='信頼区間の下限', color='blue')
plt.hist(upper_limits, bins=30, alpha=0.5, label='信頼区間の上限', color='red')
plt.axvline(lambda_true, color='green', linestyle='--', label='真の λ')
plt.title('信頼区間の分布')
plt.xlabel('λ の値')
plt.ylabel('頻度')
plt.legend()
plt.grid(True)
plt.show()

# 信頼区間のプロット
plt.figure(figsize=(12, 6))
for i in range(num_simulations):
    plt.plot([lower_limits[i], upper_limits[i]], [i, i], color='blue')
plt.axvline(lambda_true, color='red', linestyle='--', label='真の λ')
plt.title('シミュレーションごとの信頼区間')
plt.xlabel('λ の値')
plt.ylabel('シミュレーション回数')
plt.legend()
plt.grid(True)
plt.show()

真のλが95.45%の確率で信頼区間内に含まれている様子が確認できました。

2021 Q3(2)[2-3]

ポアソン分布に従う確率変数の和Tからλの最尤推定量を求める問題をやりました。

コード

scipyのminimizeで、最尤推定量を探索する。

# 2021 Q3(2)[2-3]  2024.8.28

import numpy as np
from scipy.optimize import minimize
from scipy.stats import poisson

# パラメータの設定
lambda_value = 5  # 真のパラメータ(実際のλ)
n = 10  # サンプル数
num_simulations = 1000  # シミュレーション回数

# ポアソン分布に従う n 個のサンプルを生成
samples = np.random.poisson(lambda_value, n)

# 尤度関数(負の対数尤度関数)
def negative_log_likelihood(lambda_hat, data):
    # ポアソン分布に基づく対数尤度を計算。minimizeを用いるためマイナス符号をつける
    return -np.sum(poisson.logpmf(data, lambda_hat))

# 最尤推定量の推定
initial_guess = np.mean(samples)  # 初期値としてサンプルの平均を使用
result = minimize(negative_log_likelihood, x0=initial_guess, args=(samples,), bounds=[(0.001, None)])
lambda_mle = result.x[0]

# 結果を表示
print(f"実際のλ: {lambda_value}")
print(f"推定されたλの最尤推定値: {lambda_mle}")
実際のλ: 5
推定されたλの最尤推定値: 4.9

尤度関数のグラフを確認

# 2021 Q3(2)[2-3]  2024.8.28

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

# パラメータの設定
lambda_value = 5  # 真のパラメータ(実際のλ)
n = 10  # サンプル数

# ポアソン分布に従う n 個のサンプルを生成
samples = np.random.poisson(lambda_value, n)

# 尤度関数を定義
def likelihood(lambda_hat, data):
    return np.prod(poisson.pmf(data, lambda_hat))

# 尤度関数と対数尤度関数の計算
lambda_range = np.linspace(0.1, 10, 500)  # すべてのグラフで共通の範囲
likelihood_values = [likelihood(l, samples) for l in lambda_range]
log_likelihood_values = np.log(likelihood_values)

# 最大の尤度を与えるλを求める
max_likelihood_lambda = lambda_range[np.argmax(likelihood_values)]

# グラフのプロット
fig, ax = plt.subplots(2, 1, figsize=(10, 6))

# 尤度関数のプロット
ax[0].plot(lambda_range, likelihood_values, label='尤度', color='blue')
ax[0].axvline(x=max_likelihood_lambda, color='red', linestyle='--', label=f'最大尤度 λ = {max_likelihood_lambda:.2f}')
ax[0].set_title("尤度関数")
ax[0].set_xlabel("λ の値")
ax[0].set_ylabel("尤度")
ax[0].set_xlim(lambda_range[0], lambda_range[-1])  # 横軸を合わせる
ax[0].legend()
ax[0].grid(True)

# 対数尤度関数のプロット
ax[1].plot(lambda_range, log_likelihood_values, label='対数尤度', color='green')
ax[1].axvline(x=max_likelihood_lambda, color='red', linestyle='--', label=f'最大尤度 λ = {max_likelihood_lambda:.2f}')
ax[1].set_title("対数尤度関数")
ax[1].set_xlabel("λ の値")
ax[1].set_ylabel("対数尤度")
ax[1].set_xlim(lambda_range[0], lambda_range[-1])  # 横軸を合わせる
ax[1].legend()
ax[1].grid(True)

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

# 最終結果の表示
print(f"実際のλ: {lambda_value}")
print(f"推定されたλの最尤推定値: {max_likelihood_lambda}")

尤度関数でも対数尤度関数でも同じλでピークになる。

nが小さいとλの推定量の分散が大きくなるので誤差が大きくなる。上のグラフはたまたま5.00になった。n=100ぐらいにすると安定する。

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

# 2021 Q3(2)[2-3]  2024.8.28

import numpy as np

# パラメータの設定
lambda_value = 5  # 真のパラメータ(実際のλ)
n = 10  # サンプル数
num_simulations = 1000  # シミュレーション回数

# ポアソン分布に従う n 個のサンプルを num_simulations 回生成
samples = np.random.poisson(lambda_value, (num_simulations, n))

# 各サンプルの和 T を計算
T = np.sum(samples, axis=1)

# λ の最尤推定量を計算
lambda_hat = T / n

# 結果を表示
print(f"実際のλ: {lambda_value}")
print(f"推定されたλの平均: {np.mean(lambda_hat)}")
print(f"推定されたλの分散: {np.var(lambda_hat)}")
実際のλ: 5
推定されたλの平均: 5.0058
推定されたλの分散: 0.49292636

2021 Q3(2)[2-2]

ポアソン分布に従う独立な変数の和はパラメータλの十分統計量であることを学びました。

コード

ポアソン分布に従う10個の乱数の和 Tの分布を見てみます。

# 2024 Q3(2)[2-2]  2024.8.27

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
lambda_value = 5  # 真のパラメータ
n = 10  # サンプル数
num_simulations = 1000000  # シミュレーション回数

# ポアソン分布に従う n 個のサンプルを生成
samples = np.random.poisson(lambda_value, (num_simulations, n))

# 和 T を計算
T = np.sum(samples, axis=1)

# 和 T のヒストグラムをプロット
plt.hist(T, bins=30, density=True, alpha=0.75, color='blue')
plt.title(f'Tの分布 (n 個のポアソン変数の和)')
plt.xlabel('T')
plt.ylabel('確率密度')
plt.grid(True)
plt.show()

T~Po(nλ)になっています。

つぎにT=50の条件下での確率変数X1の分布を見てます。

# 2024 Q3(2)[2-2]  2024.8.27

import numpy as np
import matplotlib.pyplot as plt

# 条件付き分布の確認
# 例えば、T = 50 の場合におけるサンプルの分布を確認

T_value = 50  # T の特定の値
samples_given_T = samples[np.sum(samples, axis=1) == T_value]

# 確率変数の1つ (X1) の分布をプロット
plt.hist(samples_given_T[:, 0], bins=15, density=True, alpha=0.75, color='green')
plt.title(f'T={T_value} の条件付き X1 の分布')
plt.xlabel('X1')
plt.ylabel('確率密度')
plt.grid(True)
plt.show()

二項分布に従っているようです。

異なるλに対して、和T=50の条件下でのX1の分布を重ねて見てます。

# 2024 Q3(2)[2-2]  2024.8.27

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

# パラメータ設定
lambda_values = [3, 5, 7]  # 異なる λ の値
n = 10  # サンプル数
num_simulations = 1000000  # シミュレーション回数
T_value = 50  # T の特定の値

plt.figure(figsize=(12, 6))

for lambda_value in lambda_values:
    # ポアソン分布に従う n 個のサンプルを生成
    samples = np.random.poisson(lambda_value, (num_simulations, n))
    
    # T = T_value のサンプルを抽出
    samples_given_T = samples[np.sum(samples, axis=1) == T_value]
    
    # X1 の分布をプロット
    counts, bin_edges, _ = plt.hist(samples_given_T[:, 0], bins=15, density=True, alpha=0.5, label=f'λ = {lambda_value}')

# 理論二項分布を描画
x1_min, x1_max = np.min(samples_given_T[:, 0]), np.max(samples_given_T[:, 0])
x = np.arange(x1_min, x1_max + 1)
estimated_p = 1 / n  # 二項分布の p は 1/n で推定
binom_pmf = binom.pmf(x, T_value, estimated_p)
plt.plot(x, binom_pmf, 'k--', label=f'理論二項分布 Bin(T={T_value}, p=1/{n})')

plt.title(f'X1の分布 (T={T_value} 条件付き) と理論二項分布')
plt.xlabel('X1')
plt.ylabel('確率密度')
plt.legend()
plt.grid(True)
plt.show()

λに依存せず同じ二項分布に従います。よって統計量Tはλに対する十分統計量であることが示されました。

なお、X1~Xnの組は多項分布に従うが、簡単のためX1のみの分布を確認しました。X1は二項分布に従います。

2021 Q3(2)[2-1]

ポアソン分布の再生性をモーメント母関数の積から示す問題をやりました。

コード

数式を使った計算

# 2021 Q3(2)[2-1]  2024.8.26

import sympy as sp

# 変数の定義
s, lambda_, n = sp.symbols('s lambda_ n', real=True, positive=True)

# 各 X_i のモーメント母関数
M_X = sp.exp(lambda_ * (sp.exp(s) - 1))

# 和 T のモーメント母関数 M_T(s)
M_T = M_X**n

# M_T(s) を簡略化
M_T_simplified = sp.simplify(M_T)

# 結果を表示
display(M_T_simplified)

Po(nλ)のモーメント母関数に等しいです。

モーメント母関数を使って期待値と分散を求めます。

# 2021 Q3(2)[2-1]  2024.8.26

import sympy as sp

# 変数の定義
s, lambda_, n = sp.symbols('s lambda_ n', real=True, positive=True)

# モーメント母関数 M_T(s)
M_T = sp.exp(-lambda_ * n * (1 - sp.exp(s)))

# 1次モーメント(期待値)の計算
M_T_prime = sp.diff(M_T, s)
expectation = M_T_prime.subs(s, 0)

# 2次モーメントの計算
M_T_double_prime = sp.diff(M_T_prime, s)
second_moment = M_T_double_prime.subs(s, 0)

# 分散の計算
variance = second_moment - expectation**2

# 結果を表示
display(expectation, variance)

Po(nλ)の期待値と分散に一致します。