# 2016 Q2(4)  2024.11.19

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
lambda_value = 2  # 固定されたλ
c = 1.0  # 定数c
n_values = range(1, 21)  # nを1から20まで変化させる
sample_size = 1000  # シミュレーションのサンプルサイズ

simulated_expectations = []  # シミュレーション期待値を格納するリスト
theoretical_expectations = []  # 理論値を格納するリスト

# 修正版コード
for n in n_values:
    # ガンマ分布のサンプルを生成 (Y ~ Gamma(n, λ))
    Y_samples = np.random.gamma(shape=n, scale=1/lambda_value, size=sample_size)
    
    # 条件付きで (1 - c/Y)^(n-1) を計算
    values = np.where(Y_samples >= c, (1 - c / Y_samples)**(n - 1), 0)
    
    # シミュレーション期待値を計算
    simulated_expectations.append(np.mean(values))
    
    # 理論値 e^(-λc) を計算
    theoretical_expectations.append(np.exp(-lambda_value * c))

# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(n_values, theoretical_expectations, 'r-o', label='理論値: e^(-λc)', linewidth=2)
plt.plot(n_values, simulated_expectations, 'b-s', label='シミュレーション期待値', linewidth=2)

# グラフの設定
plt.title('理論値とシミュレーション期待値の比較 (n を変化)', fontsize=16)
plt.xlabel('n (形状パラメータ)', fontsize=14)
plt.ylabel('期待値', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True)

# プロットの表示
plt.show()

# 2016 Q2(4)  2024.11.18

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

# パラメータ設定
lambda_value = 2  # 真のλ
n = 5  # 指数分布の和の数 (ガンマ分布の形状パラメータ)
sample_size = 10000  # シミュレーションのサンプルサイズ

# 1. 指数分布の和をシミュレーション
sum_of_exponentials = np.sum(np.random.exponential(scale=1/lambda_value, size=(sample_size, n)), axis=1)

# 2. ガンマ分布の理論値を計算
x = np.linspace(0, max(sum_of_exponentials), 1000)
gamma_pdf = gamma.pdf(x, a=n, scale=1/lambda_value)

# 3. ヒストグラムと理論的なガンマ分布をプロット
plt.figure(figsize=(10, 6))
plt.hist(sum_of_exponentials, bins=50, density=True, alpha=0.6, color='skyblue', label='シミュレーション (指数分布の和)')
plt.plot(x, gamma_pdf, 'r-', label='理論的なガンマ分布', linewidth=2)

# グラフの設定
plt.title('指数分布の和 (n=5) とガンマ分布の比較', fontsize=16)
plt.xlabel('値', fontsize=14)
plt.ylabel('密度', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True)

# プロット表示
plt.show()

# 2016 Q2(4)  2024.11.18

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

# パラメータ設定
lambda_value = 2  # 真のλ
n_values = [2, 5, 10, 20]  # nの値を変化させる
sample_size = 10000  # シミュレーションのサンプルサイズ

# 2x2のグリッドでプロット
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

for i, n in enumerate(n_values):
    # 指数分布の和をシミュレーション
    sum_of_exponentials = np.sum(np.random.exponential(scale=1/lambda_value, size=(sample_size, n)), axis=1)
    
    # ガンマ分布の理論値を計算
    x = np.linspace(0, 20, 1000)  # 横軸の範囲を0～20に設定
    gamma_pdf = gamma.pdf(x, a=n, scale=1/lambda_value)
    
    # ヒストグラムと理論的なガンマ分布をプロット
    ax = axes[i]
    ax.hist(sum_of_exponentials, bins=50, range=(0, 20), density=True, alpha=0.6, color='skyblue', label='シミュレーション (指数分布の和)')
    ax.plot(x, gamma_pdf, 'r-', label='理論的なガンマ分布', linewidth=2)
    
    # グラフの設定
    ax.set_title(f'n={n}', fontsize=14)
    ax.set_xlabel('値', fontsize=12)
    ax.set_ylabel('密度', fontsize=12)
    ax.legend(fontsize=10)
    ax.grid(True)

# 全体のレイアウト調整
plt.suptitle('nを変化させた指数分布の和とガンマ分布の比較', fontsize=16)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

# 2019 Q2(3)  2024.9.15

import sympy as sp

# 変数の定義
u, lambda_value = sp.symbols('u lambda', positive=True, real=True)

# ガンマ分布の確率密度関数 g(u) = λ^2 * u * exp(-λ * u)
g_u = lambda_value**2 * u * sp.exp(-lambda_value * u)

# 1/U の期待値の積分
integrand = (1/u) * g_u  # 1/u * g(u)

# 積分の実行
expectation = sp.integrate(integrand, (u, 0, sp.oo))

# 結果の簡略化
expectation_simplified = sp.simplify(expectation)

# 結果の表示
display(expectation_simplified)

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
lambda_value = 2  # λ=2
num_samples = 100000  # サンプルの数

# 指数分布に従う乱数 X1, X2 を生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)

# U = X1 + X2 の計算
U = X1 + X2

# 1/U の計算
inv_U = 1 / U

# シミュレーションの結果からの平均値
simulation_result = np.mean(inv_U)

# 理論値 E[1/U] = λ
theoretical_value = lambda_value

# 結果の表示
print(f'シミュレーションによる E[1/U]: {simulation_result}')
print(f'理論値 E[1/U]: {theoretical_value}')

# ヒストグラムの描画
plt.hist(inv_U, bins=50, density=True, alpha=0.7, label='シミュレーション結果')
plt.axvline(theoretical_value, color='r', linestyle='--', label=f'理論値 E[1/U] = {theoretical_value}')
plt.xlabel('1/U')
plt.ylabel('確率密度')
plt.title('シミュレーションによる 1/U の分布と理論値の比較')
plt.legend()
plt.show()

シミュレーションによる E[1/U]: 1.9928303711222939
理論値 E[1/U]: 2

# 2019 Q2(2)  2024.9.14

import sympy as sp

# 変数の定義
x, u, lambda_value = sp.symbols('x u lambda', positive=True, real=True)

# 指数分布の確率密度関数 (PDF) f(x)
f_x = lambda_value * sp.exp(-lambda_value * x)

# 畳み込み積分を計算して、g(u) = f(x) * f(u - x) の形式にする
g_u = sp.integrate(f_x * lambda_value * sp.exp(-lambda_value * (u - x)), (x, 0, u))

# 簡略化
g_u_simplified = sp.simplify(g_u)

# 結果の表示
display(g_u_simplified)

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

# パラメータ設定
lambda_value = 2  # λ=2
x_range = np.linspace(0, 5, 1000)  # X1, X2 の範囲

# X1, X2 の指数分布の確率密度関数 (PDF) - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)

# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
pdf_U = gamma.pdf(u_range, a=2, scale=1/lambda_value)

# グラフのプロット
plt.figure(figsize=(8, 6))

# X1 (X2) の指数分布の描画
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')

# U のガンマ分布の描画
plt.plot(u_range, pdf_U, label='U = X1 + X2 (ガンマ分布)', color='red')

# グラフの装飾
plt.title('X1 (X2) と U の確率密度関数の比較')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

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

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

# パラメータ設定
lambda_value = 2  # λ=2
x_range = np.linspace(0, 5, 1000)  # X1, X2, U の範囲を0～5に設定

# ガンマ分布の形状パラメータを1から2まで0.1ずつ変化させる
shape_params = np.arange(1, 2.1, 0.1)  # 1から2までの形状パラメータを0.1ステップで変化

# グラフのプロット
plt.figure(figsize=(8, 6))

# 形状パラメータが1から2に変化するガンマ分布の描画
for k in shape_params:
    gamma_pdf = gamma.pdf(x_range, a=k, scale=1/lambda_value)
    plt.plot(x_range, gamma_pdf, label=f'形状パラメータ k={k:.1f}')

# X1 (X2) の指数分布の描画 - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')

# グラフの装飾
plt.title('指数分布からガンマ分布への変化 (形状パラメータの0.1ステップ変化)')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()

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

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

# パラメータ設定
lambda_value = 2  # λ=2
x_range = np.linspace(0, 5, 1000)  # X1, X2 の範囲

# X1, X2 の指数分布の累積分布関数 (CDF) - 同じなので1つにまとめる
cdf_X1_X2 = expon.cdf(x_range, scale=1/lambda_value)

# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
cdf_U = gamma.cdf(u_range, a=2, scale=1/lambda_value)

# グラフのプロット
plt.figure(figsize=(8, 6))

# X1 (X2) の指数分布のCDFの描画
plt.plot(x_range, cdf_X1_X2, label='X1 (X2) (指数分布 CDF)', color='blue', linestyle='--')

# U のガンマ分布のCDFの描画
plt.plot(u_range, cdf_U, label='U = X1 + X2 (ガンマ分布 CDF)', color='red')

# グラフの装飾
plt.title('X1 (X2) と U の累積分布関数 (CDF) の比較')
plt.xlabel('値')
plt.ylabel('累積確率')
plt.legend()

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

import sympy as sp

# 変数の定義
x = sp.Symbol('x', positive=True)
lambda_value = sp.Symbol('lambda', positive=True)

# 指数分布の確率密度関数 (PDF)
f_x = lambda_value * sp.exp(-lambda_value * x)

# 期待値の式を計算 (積分)
E_X1 = sp.integrate(x * f_x, (x, 0, sp.oo))

# 期待値
E_X1_simplified = sp.simplify(E_X1)

# X1 + X2 の場合 (独立性より)
E_U = 2 * E_X1_simplified

# 結果の表示
display(E_U)

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

# パラメータ設定
lambda_value = 2  # 例としてλ=2
num_samples = 1000  # サンプルの数

# 指数分布に従う乱数の生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)

# U = X1 + X2 の計算
U = X1 + X2

# 理論的な期待値
expected_value = 2 / lambda_value

# シミュレーション結果のヒストグラムのプロット (確率密度にスケール)
plt.hist(U, bins=30, density=True, alpha=0.7, label='シミュレーション結果')

# ガンマ分布 (形状パラメータ k=2, スケールパラメータ θ=1/λ) の確率密度関数 (PDF) を計算
k = 2  # 形状パラメータ (k = 2)
theta = 1 / lambda_value  # スケールパラメータ (θ = 1/λ)

# ガンマ分布に基づくx軸の範囲
x = np.linspace(0, max(U), 1000)

# ガンマ分布の確率密度関数 (PDF) を計算
gamma_pdf = gamma.pdf(x, a=k, scale=theta)

# ガンマ分布の折れ線グラフをプロット
plt.plot(x, gamma_pdf, 'r-', label=f'ガンマ分布 (k={k}, θ={theta:.2f})')

# グラフの装飾
plt.axvline(expected_value, color='g', linestyle='--', label=f'理論値 = {expected_value:.2f}')
plt.xlabel('U = X1 + X2')
plt.ylabel('確率密度')
plt.title('Uの分布とガンマ分布の比較')
plt.legend()
plt.show()

# 実際のサンプル平均値と理論的期待値の比較
print(f'理論的な期待値: {expected_value}')
print(f'シミュレーションからの平均値: {np.mean(U)}')

理論的な期待値: 1.0
シミュレーションからの平均値: 0.9902440726761544