!sw_vers
!python -V

import sympy
from sympy import (
    symbols, Symbol, S, oo, Sum, limit, exp, log, sqrt, pi,
    simplify, binomial, factorial, diff, solve, Eq, Matrix,
    hessian, Function, init_printing, latex, gamma
)
from sympy.stats import Normal, Binomial, density, E, variance
import matplotlib.pyplot as plt
import numpy as np
import math
from IPython.display import display # 数式表示用

# 数式をLaTeX形式で表示するための設定
init_printing()

# グラフのスタイル設定
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['font.family'] = 'DejaVu Sans'

print("sympy version :", sympy.__version__)

# シンボルの定義
n = symbols('n', integer=True, positive=True)
k = symbols('k', integer=True, nonnegative=True)
p = symbols('p', real=True, positive=True)
lam = symbols('lambda', real=True, positive=True)

# 二項分布のPMF
pmf_binom = binomial(n, k) * p**k * (1 - p)**(n - k)

print("二項分布の確率質量関数 (PMF):")
display(pmf_binom)

# p -> lambda/n の置換
pmf_sub = pmf_binom.subs(p, lam / n)

# 項の分解
# Term A: nCk * k! / n^k
# binomial(n, k) = gamma(n+1) / (gamma(k+1) * gamma(n-k+1))
term_A = (gamma(n + 1) / gamma(n - k + 1)) / n**k

# Term B: lambda^k / k!
term_B = lam**k / factorial(k)

# Term C: (1 - lambda/n)^n
term_C = (1 - lam / n)**n

# Term D: (1 - lambda/n)^(-k)
term_D = (1 - lam / n)**(-k)

print("分解した項:")
print("Term A:")
display(term_A)
print("Term C:")
display(term_C)

# Term A の極限確認 (kを具体的な値にして確認)
print("Term A の極限値確認:")
for k_val in [0, 1, 2, 5]:
    term_A_concrete = term_A.subs(k, k_val)
    lim_val = limit(term_A_concrete, n, oo)
    print(f"k={k_val}: limit = {lim_val}")

# 従って、一般の有限整数kにおいて Term A -> 1 とみなす
lim_A = 1

# Term C, D の極限
lim_C = limit(term_C, n, oo)
lim_D = limit(term_D, n, oo)

print("Term C の極限:", lim_C)
print("Term D の極限:", lim_D)

# 全体の計算結果
poisson_pmf_derived = lim_A * term_B * lim_C * lim_D

print("導出されたポアソン分布:")
display(poisson_pmf_derived)

# 理論式との比較
expected_poisson_pmf = (lam**k * exp(-lam)) / factorial(k)
print("数式の一致確認:", simplify(poisson_pmf_derived - expected_poisson_pmf) == 0)

def plot_convergence(lam_val=5, n_list=[10, 20, 50, 100]):
    k_vals = np.arange(0, 20)

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

    # ポアソン分布（理論値）
    # np.math は非推奨のため math を使用
    poisson_probs = [np.exp(-lam_val) * (lam_val**k_i) / math.factorial(k_i) for k_i in k_vals]
    plt.plot(k_vals, poisson_probs, 'k--', linewidth=2, label=rf'Poisson($\lambda={lam_val}$)')

    # 二項分布
    markers = ['o', 's', '^', 'D']
    for i, n_val in enumerate(n_list):
        p_val = lam_val / n_val
        binom_probs = []
        for k_i in k_vals:
            if k_i > n_val:
                prob = 0
            else:
                coef = math.comb(n_val, k_i)
                prob = coef * (p_val**k_i) * ((1 - p_val)**(n_val - k_i))
            binom_probs.append(prob)

        plt.plot(k_vals, binom_probs, marker=markers[i % len(markers)], linestyle='-', alpha=0.7, label=f'Binomial($n={n_val}, p={p_val:.3g}$)')

    plt.title(rf'Convergence of Binomial to Poisson ($\lambda={lam_val}$)', fontsize=16)
    plt.xlabel('k', fontsize=12)
    plt.ylabel('P(X=k)', fontsize=12)
    plt.legend(fontsize=12)
    plt.grid(True, alpha=0.3)
    plt.show()

plot_convergence(lam_val=5)

# シンボルの定義
mu = symbols('mu', real=True)
sigma = symbols('sigma', real=True, positive=True)
x = symbols('x', real=True)
n_samples = symbols('n', integer=True, positive=True)

# 正規分布のPDF
pdf_normal = 1 / sqrt(2 * pi * sigma**2) * exp(-(x - mu)**2 / (2 * sigma**2))

print("正規分布のPDF:")
display(pdf_normal)

# 定数項と対数項
term1 = -n_samples / 2 * log(2 * pi)
term2 = -n_samples * log(sigma)

# 二乗和の展開項
s_1 = symbols('s_1', real=True) # sum(x_i)
s_2 = symbols('s_2', real=True) # sum(x_i^2)
sum_sq_expanded = s_2 - 2*mu*s_1 + n_samples * mu**2

# 対数尤度関数
log_L_calc = term1 + term2 - 1/(2*sigma**2) * sum_sq_expanded

print("対数尤度関数:")
display(log_L_calc)

# mu で偏微分
d_logL_d_mu = diff(log_L_calc, mu)

print("スコア関数 (mu):")
display(simplify(d_logL_d_mu))

# 方程式の解
mu_solution = solve(Eq(d_logL_d_mu, 0), mu)

print("mu の最尤推定量:")
display(mu_solution[0])

# sigma で偏微分
d_logL_d_sigma = diff(log_L_calc, sigma)

print("スコア関数 (sigma):")
display(simplify(d_logL_d_sigma))

# mu に解を代入して解く
d_logL_d_sigma_subs = d_logL_d_sigma.subs(mu, mu_solution[0])
sigma_solutions = solve(Eq(d_logL_d_sigma_subs, 0), sigma)

# 正の解を選択
sigma_mle = [s for s in sigma_solutions if not s.is_negative][0]

print("sigma の最尤推定量:")
display(sigma_mle)

# ヘッセ行列の計算
vars_list = [mu, sigma]
hessian_mat = hessian(log_L_calc, vars_list)

print("ヘッセ行列:")
display(hessian_mat)

# 推定値における評価 (1,1)成分
H_at_mle = hessian_mat.subs({
    mu: mu_solution[0],
    sigma: sigma_mle
})

h11_val = simplify(H_at_mle[0, 0])
print("(1,1)成分の値:")
display(h11_val)

# 期待値計算
# E[x_i] = mu, E[(x_i-mu)^2] = sigma^2
i11 = -hessian_mat[0, 0]
i12 = 0
term_sum_sq = n_samples * sigma**2
val_h22 = hessian_mat[1, 1].subs(sum_sq_expanded, term_sum_sq)
i22 = -(-2 * n_samples / sigma**2)

Fisher_Info = Matrix([
    [i11, i12],
    [i12, i22]
])

print("フィッシャー情報行列:")
display(Fisher_Info)

def plot_log_likelihood(n=100, true_mu=5.0, true_sigma=2.0):
    np.random.seed(42)
    data = np.random.normal(true_mu, true_sigma, n)

    sample_mean = np.mean(data)
    sample_std = np.std(data)

    mu_range = np.linspace(true_mu - 1.0, true_mu + 1.0, 100)
    sigma_range = np.linspace(true_sigma - 1.0, true_sigma + 1.0, 100)
    MU, SIGMA = np.meshgrid(mu_range, sigma_range)

    term1 = -n/2 * np.log(2 * np.pi)
    term2 = -n * np.log(SIGMA)
    SSD = np.sum((data[:, np.newaxis, np.newaxis] - MU)**2, axis=0)
    LOG_L = term1 + term2 - SSD / (2 * SIGMA**2)

    plt.figure(figsize=(10, 8))
    contours = plt.contour(MU, SIGMA, LOG_L, levels=20, cmap='viridis')
    plt.clabel(contours, inline=True, fontsize=8)

    plt.plot(sample_mean, sample_std, 'rx', markersize=12, markeredgewidth=2, label='MLE')
    plt.plot(true_mu, true_sigma, 'wo', markerfacecolor='none', markeredgewidth=2, label='True Param')

    plt.title(rf'Log-Likelihood Contour (n={n})', fontsize=16)
    plt.xlabel(rf'$\mu$', fontsize=14)
    plt.ylabel(rf'$\sigma$', fontsize=14)
    plt.legend()
    plt.colorbar(label='Log-Likelihood')
    plt.show()

plot_log_likelihood()