# 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

# 2022 Q3(4)  2024.8.7

import sympy as sp

# 定義
k, lambda_var = sp.symbols('k lambda')
alpha, beta = sp.symbols('alpha beta', positive=True)

# ポアソン分布の確率質量関数
poisson_pmf = (lambda_var**k * sp.exp(-lambda_var)) / sp.factorial(k)

# ガンマ分布の確率密度関数
gamma_pdf = (beta**alpha / sp.gamma(alpha)) * lambda_var**(alpha - 1) * sp.exp(-beta * lambda_var)

# 周辺分布の計算
marginal_distribution = sp.integrate(poisson_pmf * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

# 期待値 E[X] の計算
expected_value = sp.summation(k * marginal_distribution, (k, 0, sp.oo)).simplify()

# 分散 V[X] の計算
expected_value_of_square = sp.summation(k**2 * marginal_distribution, (k, 0, sp.oo)).simplify()
variance = (expected_value_of_square - expected_value**2).simplify()

# 結果を表示
results = {
    "期待値 E[X]": expected_value,
    "分散 V[X]": variance
}

results

{'期待値 E[X]': alpha/beta, '分散 V[X]': alpha*(beta + 1)/beta**2}

import numpy as np

# パラメータの設定
alpha = 3.0
beta = 2.0
sample_size = 10000

# ガンマ分布から λ を生成
lambda_samples = np.random.gamma(alpha, 1/beta, sample_size)

# 生成された λ を使ってポアソン分布から X を生成
poisson_samples = [np.random.poisson(lam) for lam in lambda_samples]

# サンプルの期待値と分散を計算
sample_mean = np.mean(poisson_samples)
sample_variance = np.var(poisson_samples)

# 理論値を計算
theoretical_mean = alpha / beta
theoretical_variance = alpha * (beta + 1) / beta**2

# 結果を表示
results_with_caption = {
    "シミュレーションによる期待値": sample_mean,
    "理論値による期待値": theoretical_mean,
    "シミュレーションによる分散": sample_variance,
    "理論値による分散": theoretical_variance
}

results_with_caption

{'シミュレーションによる期待値': 1.5027,
 '理論値による期待値': 1.5,
 'シミュレーションによる分散': 2.31159271,
 '理論値による分散': 2.25}

# 2022 Q3(3)  2024.8.5

import sympy as sp

# 定義
k = sp.symbols('k', integer=True)
alpha, beta = sp.symbols('alpha beta', positive=True)
lambda_var = sp.symbols('lambda', positive=True)

# ポアソン分布の確率質量関数
poisson_pmf = (lambda_var**k * sp.exp(-lambda_var)) / sp.factorial(k)

# ガンマ分布の確率密度関数
gamma_pdf = (beta**alpha / sp.gamma(alpha)) * lambda_var**(alpha - 1) * sp.exp(-beta * lambda_var)

# 周辺分布の計算
marginal_distribution = sp.integrate(poisson_pmf * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

marginal_distribution

# 2022 Q3(3)  2024.8.5

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma, nbinom

# パラメータの設定
alpha = 3.0
beta = 2.0
sample_size = 10000

# ガンマ分布から λ を生成
lambda_samples = np.random.gamma(alpha, 1/beta, sample_size)

# 生成された λ を使ってポアソン分布から X を生成
poisson_samples = [np.random.poisson(lam) for lam in lambda_samples]

# 負の二項分布のパラメータを設定
r = alpha
p = beta / (beta + 1)

# 負の二項分布から理論的な確率質量関数を計算
x = np.arange(0, max(poisson_samples) + 1)
nbinom_pmf = nbinom.pmf(x, r, p)

# ヒストグラムの描画
plt.hist(poisson_samples, bins=np.arange(0, max(poisson_samples) + 1) - 0.5, density=True, alpha=0.75, color='blue', edgecolor='black', rwidth=0.8)

# 理論的な負の二項分布の確率質量関数をプロット
plt.plot(x, nbinom_pmf, 'r', linestyle='-', label='理論的な負の二項分布')

# グラフのタイトルとラベル
plt.title('ポアソン-ガンマ混合分布のシミュレーション結果')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

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