#2018 Q4(2) 2024.10.19

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n_simulations = 100000  # シミュレーション回数
rho = 0.7  # ρの値
x_t = 1.0  # X_t の初期値

# シミュレーション結果を格納するリスト
X_t1_samples = []

for _ in range(n_simulations):
    # 1. X_t -> Y_t のステップ
    y_t = np.random.normal(loc=rho * x_t, scale=np.sqrt(1 - rho**2))
    
    # 2. Y_t -> X_{t+1} のステップ
    x_t1 = np.random.normal(loc=rho * y_t, scale=np.sqrt(1 - rho**2))
    
    # 結果をリストに追加
    X_t1_samples.append(x_t1)

# シミュレーション結果を配列に変換
X_t1_samples = np.array(X_t1_samples)

# シミュレーション結果の期待値と分散を計算
simulated_mean = np.mean(X_t1_samples)
simulated_variance = np.var(X_t1_samples)

# 理論値の期待値と分散
expected_mean = rho**2 * x_t  # 理論上の期待値
expected_variance = 1 - rho**4  # 理論上の分散

# 結果の表示
print(f"理論上の期待値: {expected_mean}, シミュレーションによる期待値: {simulated_mean}")
print(f"理論上の分散: {expected_variance}, シミュレーションによる分散: {simulated_variance}")

# ヒストグラムの表示
plt.hist(X_t1_samples, bins=50, density=True, alpha=0.6, color='g', edgecolor='black', label='シミュレーション')
plt.axvline(expected_mean, color='r', linestyle='--', label=f'理論上の期待値: {expected_mean:.3f}')
plt.title(r"$X_{t+1}$ のシミュレーション (マルコフ性)")
plt.xlabel(r"$X_{t+1}$ の値")
plt.ylabel("頻度")
plt.legend()
plt.grid(True)
plt.show()

理論上の期待値: 0.48999999999999994, シミュレーションによる期待値: 0.4909842816275996
理論上の分散: 0.7599, シミュレーションによる分散: 0.7608618418358614