# 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

# 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

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)

## 2022 Q3(2)  2024.8.4

import sympy as sp

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

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

# 期待値 E[Λ] の計算
expected_value = sp.integrate(lambda_var * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

# E[Λ^2] の計算
expected_value_Lambda2 = sp.integrate(lambda_var**2 * gamma_pdf, (lambda_var, 0, sp.oo)).simplify()

# 分散 V[Λ] の計算
variance = (expected_value_Lambda2 - expected_value**2).simplify()

# 結果を辞書形式で表示
{
    "期待値 E[Λ]": expected_value,
    "分散 V[Λ]": variance
}

{'期待値 E[Λ]': alpha/beta, '分散 V[Λ]': alpha/beta**2}

## 2022 Q3(2)  2024.8.4

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

# パラメータ α と β の設定
alpha = 3.0
beta = 2.0
# サンプルサイズ
sample_size = 10000

# ガンマ分布に従う乱数を生成
samples = np.random.gamma(alpha, 1/beta, sample_size)

# 期待値と分散の計算
expected_value = np.mean(samples)
variance = np.var(samples)

# 結果の表示
print(f"期待値のシミュレーション結果: {expected_value}")
print(f"分散のシミュレーション結果: {variance}")

# ヒストグラムの描画
plt.hist(samples, bins=50, density=True, alpha=0.75, color='blue', edgecolor='black')

# 理論的なガンマ分布の確率密度関数をプロット
x = np.linspace(0, max(samples), 100)
gamma_pdf = gamma.pdf(x, alpha, scale=1/beta)
plt.plot(x, gamma_pdf, 'r', linestyle='-', label='理論的なガンマ分布')

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

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

期待値のシミュレーション結果: 1.4947527551017439
分散のシミュレーション結果: 0.7452822182213243

# パラメータ α と β の設定
alpha = 3.0
beta = 2.0

# 理論値の計算
theoretical_expected_value = alpha / beta
theoretical_variance = alpha / beta**2

# 結果の表示
print(f"理論値 - 期待値: {theoretical_expected_value}")
print(f"理論値 - 分散: {theoretical_variance}")

理論値 - 期待値: 1.5
理論値 - 分散: 0.75

# 2022 Q3(1)  2024.8.3

import sympy as sp

# 定義
k = sp.symbols('k')
lambda_param = sp.symbols('lambda')

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

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

# E[X^2] の計算
expected_value_X2 = sp.summation(k**2 * poisson_pmf, (k, 0, sp.oo)).simplify()

# 分散 V[X] の計算
variance = (expected_value_X2 - expected_value**2).simplify()

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

{'期待値 E[X]': lambda, '分散 V[X]': lambda}

import numpy as np

# パラメータ λ の設定
lambda_param = 5
# シミュレーションの回数
num_simulations = 10000
# 各シミュレーションで生成するサンプルの数
sample_size = 1000

# シミュレーション結果の保存用
expected_values = []
variances = []

for _ in range(num_simulations):
    # ポアソン分布に従う乱数を生成
    samples = np.random.poisson(lambda_param, sample_size)
    # 期待値を計算
    expected_values.append(np.mean(samples))
    # 分散を計算
    variances.append(np.var(samples))

# シミュレーション結果の平均を計算
average_expected_value = np.mean(expected_values)
average_variance = np.mean(variances)

print(f"期待値のシミュレーション結果: {average_expected_value}")
print(f"分散のシミュレーション結果: {average_variance}")

期待値のシミュレーション結果: 4.999602
分散のシミュレーション結果: 4.9992533574

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial

# パラメータ λ の設定
lambda_param = 5
# サンプルサイズ
sample_size = 10000

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

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

# 理論的なポアソン分布の確率質量関数をプロット
x = np.arange(0, max(samples) + 1)
poisson_pmf = (lambda_param**x * np.exp(-lambda_param)) / factorial(x)
plt.plot(x, poisson_pmf, 'r', marker='o', linestyle='-', label='理論的なポアソン分布')

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

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