# 2019 Q1(1)  2024.9.9

import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import derivative

# パラメータ設定
lambda_val = 3  # ポアソン分布のパラメータ λ
sample_size = 10000  # サンプル数

# 1. ポアソン分布に従う乱数データを生成
samples = np.random.poisson(lambda_val, sample_size)

# 2. P(X=k) の推定（サンプル度数から確率質量関数の推定）
max_k = np.max(samples)
pmf_estimate = np.zeros(max_k + 1)
for k in range(max_k + 1):
    pmf_estimate[k] = np.sum(samples == k) / sample_size

# 3. G_X(t) の推定
def G_X(t):
    return np.sum([pmf_estimate[k] * t**k for k in range(max_k + 1)])

# 4. 理論的な G_X(t) (ポアソン分布の確率母関数)
def G_X_theoretical(t, lam):
    return np.exp(lam * (t - 1))

# 5. G_X(t) のプロット
t_values = np.linspace(0, 1, 100)
G_X_values = [G_X(t) for t in t_values]
G_X_theoretical_values = [G_X_theoretical(t, lambda_val) for t in t_values]

plt.plot(t_values, G_X_values, label="推定された G_X(t)", color="blue")
plt.plot(t_values, G_X_theoretical_values, label="理論的な G_X(t)", linestyle='--', color="red")
plt.title("推定された確率母関数と理論的な G_X(t) の比較")
plt.xlabel("t")
plt.ylabel("G_X(t)")
plt.legend()
plt.show()

# 6. G_X(t) の微分と2階微分
def G_X_diff_1(t):
    return derivative(G_X, t, dx=1e-6)

def G_X_diff_2(t):
    return derivative(G_X_diff_1, t, dx=1e-6)

# 7. t = 1 における期待値と分散の推定
E_X_estimated = G_X_diff_1(1)  # 期待値
Var_X_estimated = G_X_diff_2(1) + E_X_estimated - E_X_estimated**2  # 分散

# 結果の出力
print(f"推定された期待値 E[X]: {E_X_estimated}")
print(f"推定された分散 Var(X): {Var_X_estimated}")

# 真の λ と比較
print(f"真の期待値 E[X]: {lambda_val}")
print(f"真の分散 Var(X): {lambda_val}")

推定された期待値 E[X]: 3.0239000000098493
推定された分散 Var(X): 3.0824841743413582
真の期待値 E[X]: 3
真の分散 Var(X): 3