# 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