ホーム » Qchan の投稿 (ページ 11)

作者アーカイブ: Qchan

投稿一覧

2017 Q1(3)

この前と同じ問ですが、解き方が間違っていたので、再チャレンジです。

 

コード

Xの標本平均の尖度(過剰尖度)が\beta_2\left(\overline{X}\right) = \frac{\beta_2}{n}になるか確認するためにシミュレーションを行います。母集団分布はガンマ分布とします。比較のためシミュレーション回数は前問(歪度)と同じ回数とする。

# 2017 Q1(3)  2024.10.26

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

# シミュレーションパラメータの設定
num_simulations = 30000  # シミュレーションの繰り返し回数
sample_sizes = range(10, 110, 10)  # サンプルサイズの設定

# 尖度を持つ母集団分布としてガンマ分布を使用(shape=2で正の尖度を持つ)
shape_param = 2.0  # ガンマ分布のshapeパラメータ
scale_param = 2.0  # ガンマ分布のscaleパラメータ
population_dist = gamma(a=shape_param, scale=scale_param)

# 理論的な過剰尖度(ガンマ分布の尖度の計算)
beta_2_theoretical = 6 / shape_param  # ガンマ分布の過剰尖度

# シミュレーション結果を保存するリスト
sample_mean_excess_kurtosis = []

# 各サンプルサイズでのシミュレーションの実行
for n in sample_sizes:
    kurt_values = []
    
    # シミュレーションを繰り返し
    for _ in range(num_simulations):
        sample = population_dist.rvs(size=n)  # サンプルを抽出
        sample_mean = np.mean(sample)  # 標本平均
        kurt_values.append(sample_mean)
    
    # 標本平均の尖度を計算し、過剰尖度に変換(-3)
    kurtosis_of_sample_mean = kurtosis(kurt_values, fisher=False) - 3  # 尖度を計算し、過剰尖度に変換
    sample_mean_excess_kurtosis.append(kurtosis_of_sample_mean)

# 理論値の計算
theoretical_excess_kurtosis = [beta_2_theoretical / n for n in sample_sizes]

# グラフの作成
plt.figure(figsize=(10, 6))

# シミュレーション結果のプロット
plt.plot(sample_sizes, sample_mean_excess_kurtosis, marker='o', label='シミュレーション結果 (過剰尖度)')
plt.plot(sample_sizes, theoretical_excess_kurtosis, color='r', linestyle='--', label='理論値: $\\frac{\\beta_2}{n}$')

plt.title('$\\bar{X}$ の過剰尖度のシミュレーションと理論値の比較')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('$\\bar{X}$ の過剰尖度')
plt.legend()
plt.grid(True)

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

尖度は理論値に沿うものの、歪度(前問)より更に外れ値の影響を受けやすいため分散が大きい。

シミュレーション回数を10倍にしてみます。

尖度は理論値に近づきました。

尖度の計算でマイナス3がついている

尖度の計算でマイナス3がついていることを腹落ちさせるために計算しました。

2017 Q1(2)

標本平均の歪度を求めました。

 

コード

Xの標本平均の歪度が\beta_1\left(\overline{X}\right) = \frac{\beta_1}{\sqrt{n}}になるか確認するためにシミュレーションを行います。母集団分布はガンマ分布とします。

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

# シミュレーションパラメータの設定
num_simulations = 30000  # シミュレーションの繰り返し回数
sample_sizes = range(10, 110, 10)  # サンプルサイズの設定

# 歪度を持つ母集団分布としてガンマ分布を使用(shape=2で正の歪度を持つ)
shape_param = 2.0  # ガンマ分布のshapeパラメータ
scale_param = 2.0  # ガンマ分布のscaleパラメータ
population_dist = gamma(a=shape_param, scale=scale_param)

# 理論的な歪度(ガンマ分布の歪度の計算)
beta_1_theoretical = 2 / np.sqrt(shape_param)  # ガンマ分布の理論的な歪度

# シミュレーション結果を保存するリスト
sample_mean_skewness = []

# 各サンプルサイズでのシミュレーションの実行
for n in sample_sizes:
    skew_values = []
    
    # シミュレーションを繰り返し
    for _ in range(num_simulations):
        sample = population_dist.rvs(size=n)  # サンプルを抽出
        sample_mean = np.mean(sample)  # 標本平均
        skew_values.append(sample_mean)
    
    # 標本平均の歪度を計算
    skewness_of_sample_mean = skew(skew_values)
    sample_mean_skewness.append(skewness_of_sample_mean)

# 理論値の計算
theoretical_skewness = [beta_1_theoretical / np.sqrt(n) for n in sample_sizes]

# グラフの作成
plt.figure(figsize=(10, 6))

# シミュレーション結果のプロット
plt.plot(sample_sizes, sample_mean_skewness, marker='o', label='シミュレーション結果')
plt.plot(sample_sizes, theoretical_skewness, color='r', linestyle='--', label='理論値: $\\frac{\\beta_1}{\sqrt{n}}$')

plt.title('$\\bar{X}$ の歪度のシミュレーションと理論値の比較')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('$\\bar{X}$ の歪度')
plt.legend()
plt.grid(True)

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

歪度は外れ値の影響を受けやすいものの概ね理論値と一致しました。

2017 Q1(1)

不偏分散が母分散の不偏推定量であることを示しました。

 

コード

サンプルサイズnを変化させてE[\bar{X}],V[\bar{X}],E[T^2],E[S^2]をそれぞれシミュレーションして描画します。

import numpy as np
import matplotlib.pyplot as plt

# シミュレーションパラメータの設定
mu = 5.0             # 母平均
sigma = 2.0          # 母標準偏差
num_simulations = 1000  # シミュレーションの繰り返し回数

# n の範囲を設定 (100 から 10000 まで 100 刻み)
n_values = range(100, 10100, 100)

# 各 n に対する統計量を保存するリスト
E_X_bar_list = []
V_X_bar_list = []
E_T2_list = []
E_S2_list = []

# 各 n に対するシミュレーションの実行
for n in n_values:
    sample_means = []
    T2_values = []
    S2_values = []
    
    # シミュレーションの実行
    for _ in range(num_simulations):
        sample = np.random.normal(mu, sigma, n)
        sample_mean = np.mean(sample)
        sample_means.append(sample_mean)
        
        # T^2 計算 (母平均使用)
        T2 = np.mean((sample - mu) ** 2)
        T2_values.append(T2)
        
        # S^2 計算 (標本平均使用)
        S2 = np.var(sample, ddof=1)
        S2_values.append(S2)
    
    # 結果をリストに保存
    E_X_bar_list.append(np.mean(sample_means))
    V_X_bar_list.append(np.var(sample_means))
    E_T2_list.append(np.mean(T2_values))
    E_S2_list.append(np.mean(S2_values))

# グラフの作成(理論値を追加)
plt.figure(figsize=(14, 10))

# 理論値
theoretical_mean = mu       # E[X_bar] の理論値
theoretical_variance = sigma**2 / np.array(n_values)  # V[X_bar] の理論値
theoretical_variance_true = sigma**2  # 母分散

# E[X_bar] のプロット
plt.subplot(2, 2, 1)
plt.plot(n_values, E_X_bar_list, marker='o', label='シミュレーション')
plt.axhline(y=theoretical_mean, color='r', linestyle='--', label='理論値: $E[\overline{X}] = \mu$')
plt.title('$E[\overline{X}]$')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('$E[\overline{X}]$')
plt.legend()

# V[X_bar] のプロット
plt.subplot(2, 2, 2)
plt.plot(n_values, V_X_bar_list, marker='o', label='シミュレーション')
plt.plot(n_values, theoretical_variance, color='r', linestyle='--', label='理論値: $V[\overline{X}] = \\frac{\sigma^2}{n}$')
plt.title('$V[\overline{X}]$')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('$V[\overline{X}]$')
plt.legend()

# E[T^2] のプロット
plt.subplot(2, 2, 3)
plt.plot(n_values, E_T2_list, marker='o', label='シミュレーション')
plt.axhline(y=theoretical_variance_true, color='r', linestyle='--', label='理論値: $E[T^2] = \sigma^2$')
plt.title('$E[T^2]$')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('$E[T^2]$')
plt.legend()

# E[S^2] のプロット
plt.subplot(2, 2, 4)
plt.plot(n_values, E_S2_list, marker='o', label='シミュレーション')
plt.axhline(y=theoretical_variance_true, color='r', linestyle='--', label='理論値: $E[S^2] = \sigma^2$')
plt.title('$E[S^2]$')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('$E[S^2]$')
plt.legend()

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

サンプルサイズnが増加するほど、各推定量は理論値に収束していく傾向に見えます。E[T^2]E[S^2]は外れ値の影響を受けやすいため収束速度が緩やかです。

2018 Q5(3)

一様分布の3つの順序統計量のうち第3と第1の確率変数の差の期待値と分散を求めました。

 

コード

同時確率密度関数f13(y1,y3)を視覚化してみます。またzが一定になる線も描画してみます。

# 2018 Q5(3)  2024.10.23

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 定義域 (y1 と y3 は [0, 1] の範囲で変動)
y1_values = np.linspace(0, 1, 100)
y3_values = np.linspace(0, 1, 100)

# メッシュグリッドを作成
Y1, Y3 = np.meshgrid(y1_values, y3_values)

# 同時確率密度関数 f13(y1, y3) = 6(y3 - y1), ただし y3 >= y1
f13 = np.where(Y3 >= Y1, 6 * (Y3 - Y1), 0)

# 等高線のための Z の値を設定 (例として Z = 0.2, 0.5, 0.8 を使用)
Z_values = [0.2, 0.5, 0.8]

# 3D グラフの描画を再度行う
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 3Dプロット (同時確率密度関数 f13(y1, y3))
ax.plot_surface(Y1, Y3, f13, cmap='viridis', alpha=0.8)

# 各 Z = Y3 - Y1 の線をプロット
for z_value in Z_values:
    y1_line = np.linspace(0, 1 - z_value, 100)
    y3_line = y1_line + z_value
    ax.plot(y1_line, y3_line, zs=0, zdir='z', label=f'$Z = {z_value}$', lw=2)

# 軸ラベル
ax.set_title('同時確率密度関数 $f_{13}(y_1, y_3)$ と $Z = Y_3 - Y_1$', fontsize=14)
ax.set_xlabel('$y_1$', fontsize=12)
ax.set_ylabel('$y_3$', fontsize=12)
ax.set_zlabel('$f_{13}(y_1, y_3)$', fontsize=12)

# 凡例
ax.legend()

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

Z=z で垂直に切った断面の面積がf_Z(z)に相当します。

f_Z(z) = \int_0^{1-z} f_{13}(y_1, y_1 + z) \, dy_1を解いてZの確率密度関数をf_Z(z) = 6z(1 - z), \quad 0 \leq z \leq 1のように求めます。順序統計量Y1とY3をシミュレーションし、理論と一致するか確認します。

# 2018 Q5(3)  2024.10.23

import numpy as np
import matplotlib.pyplot as plt

# サンプルデータの生成 (Z = Y3 - Y1 としてシミュレーション)
n_samples = 10000  # サンプル数
samples = np.random.uniform(0, 1, (n_samples, 3))  # 一様分布からサンプリング
Y1_samples = np.min(samples, axis=1)  # 最小値 Y1
Y3_samples = np.max(samples, axis=1)  # 最大値 Y3

# Z = Y3 - Y1 を計算
Z_samples = Y3_samples - Y1_samples

# Zの期待値と分散を計算
expected_Z = np.mean(Z_samples)
variance_Z = np.var(Z_samples)

# 理論値との比較
theoretical_Z_mean = 0.5
theoretical_Z_var = 1/20

# 期待値と分散を表示
print(f"シミュレーションによる E[Z](期待値): {expected_Z:.4f}, 理論値: {theoretical_Z_mean:.4f}")
print(f"シミュレーションによる Var[Z](分散): {variance_Z:.4f}, 理論値: {theoretical_Z_var:.4f}")

# Zの理論的な確率密度関数の定義
z_values = np.linspace(0, 1, 500)
f_z_theoretical = 6 * z_values * (1 - z_values)

# ヒストグラムと理論的なPDFを重ねて描画
plt.figure(figsize=(8, 6))

# サンプルのヒストグラムを描画 (確率密度として正規化)
plt.hist(Z_samples, bins=50, density=True, alpha=0.5, color='purple', label='Z samples (Y3 - Y1)')

# 理論的な確率密度関数のプロット
plt.plot(z_values, f_z_theoretical, label='Theoretical $f_Z(z) = 6z(1-z)$', color='black', linewidth=2)

# グラフの設定
plt.title('Z = Y3 - Y1 の確率密度関数 (シミュレーションと理論)', fontsize=14)
plt.xlabel('Z', fontsize=12)
plt.ylabel('密度', fontsize=12)
plt.legend(loc='upper right', fontsize=10)
plt.grid(True)

# グラフの表示
plt.show()
シミュレーションによる E[Z](期待値): 0.4992, 理論値: 0.5000
シミュレーションによる Var[Z](分散): 0.0500, 理論値: 0.0500

シミュレーションによる順序統計量Y3とY1の差の分布と期待値と分散が理論値と一致しました。

また興味深いことに、ZとY2は同じ分布になっていることが分かりました。

2018 Q5(2)

一様分布の3つの順序統計量のうち真ん中の確率密度関数と確率の計算をしました。

 

コード

シミュレーションにより順序統計量Y2の分布と期待値が理論値と一致するか確認します。

# 2018 Q5(2)  2024.10.22

import numpy as np
import matplotlib.pyplot as plt

# 定義域
y_values = np.linspace(0, 1, 500)

# 確率密度関数 f2(y) の定義
f2_y = 6 * y_values * (1 - y_values)

# シミュレーションの設定
n_samples = 10000   # サンプル数

# 一様分布から3つの変数をサンプリング
samples = np.random.uniform(0, 1, (n_samples, 3))

# 順序統計量の第2要素を抽出 (Y2)
Y2_samples = np.median(samples, axis=1)  # 中央値

# 期待値の計算
expected_Y2 = np.mean(Y2_samples)

# 理論値との比較
theoretical_Y2 = 0.5  # 中央値

# 期待値を表示
print(f"シミュレーションによる E[Y2](中間値の期待値): {expected_Y2:.4f}, 理論値: {theoretical_Y2:.4f}")

# ヒストグラムの作成
plt.figure(figsize=(8, 6))

# 理論密度関数 f2(y) のプロット
plt.plot(y_values, f2_y, label='$f_2(y) = 6y(1-y)$', color='red')

# サンプルのヒストグラムを追加 (確率密度として正規化)
plt.hist(Y2_samples, bins=50, density=True, alpha=0.5, color='red', label='Y2 samples (median)')

# グラフの設定
plt.title('確率密度関数 $f_2(y)$ とヒストグラム (シミュレーションと理論)', fontsize=14)
plt.xlabel('y', fontsize=12)
plt.ylabel('密度', fontsize=12)
plt.legend(loc='upper right', fontsize=10)
plt.grid(True)

# グラフの表示
plt.show()
シミュレーションによる E[Y2](中間値の期待値): 0.5006, 理論値: 0.5000

シミュレーションによる順序統計量Y1とY3の分布と期待値が理論値と一致しました。

次に、Y1~Y3の確率密度関数と累積分布関数を重ねて見てみます。

# 2018 Q5(2)  2024.10.22

import numpy as np
import matplotlib.pyplot as plt

# 定義域
y_values = np.linspace(0, 1, 500)

# 確率密度関数の定義
f1_y = 3 * (1 - y_values) ** 2
f2_y = 6 * y_values * (1 - y_values)
f3_y = 3 * y_values ** 2

# CDF (累積分布関数) の計算
cdf1_y = np.cumsum(f1_y) * (y_values[1] - y_values[0])
cdf2_y = np.cumsum(f2_y) * (y_values[1] - y_values[0])
cdf3_y = np.cumsum(f3_y) * (y_values[1] - y_values[0])

# グラフの作成
plt.figure(figsize=(10, 8))

# 確率密度関数 (PDF) のプロット
plt.subplot(2, 1, 1)
plt.plot(y_values, f1_y, label='$f_1(y)$', color='blue')
plt.plot(y_values, f2_y, label='$f_2(y)$', color='red')
plt.plot(y_values, f3_y, label='$f_3(y)$', color='green')
plt.title('確率密度関数 (PDF)', fontsize=14)
plt.xlabel('y', fontsize=12)
plt.ylabel('密度', fontsize=12)
plt.legend(loc='upper right', fontsize=10)
plt.grid(True)

# 累積分布関数 (CDF) のプロット
plt.subplot(2, 1, 2)
plt.plot(y_values, cdf1_y, label='$F_1(y)$', color='blue')
plt.plot(y_values, cdf2_y, label='$F_2(y)$', color='red')
plt.plot(y_values, cdf3_y, label='$F_3(y)$', color='green')
plt.title('累積分布関数 (CDF)', fontsize=14)
plt.xlabel('y', fontsize=12)
plt.ylabel('累積確率', fontsize=12)
plt.legend(loc='lower right', fontsize=10)
plt.grid(True)

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

Y1は左寄り、Y2は中央寄り、Y3は右寄りになっていて、順序統計量の性質を直感的に表されています。Y2の確率密度関数は左右対称なので期待値が0.5になることも分かります。Y2の累積分布関数が直線でないことも面白いです。

2018 Q5(1)

一様分布の順序統計量の確率密度関数と期待値を求めました。

 

コード

シミュレーションにより順序統計量Y1とY3の分布と期待値が理論値と一致するか確認します。

# 2018 Q5(1)  2024.10.21

import numpy as np
import matplotlib.pyplot as plt

# 定義域
y_values = np.linspace(0, 1, 500)

# 確率密度関数の定義
f1_y = 3 * (1 - y_values) ** 2
f3_y = 3 * y_values ** 2

# シミュレーションの設定
n_samples = 10000   # サンプル数

# 一様分布から3つの変数をサンプリング
samples = np.random.uniform(0, 1, (n_samples, 3))

# 順序統計量を求める
Y1_samples = np.min(samples, axis=1)  # 最小値
Y3_samples = np.max(samples, axis=1)  # 最大値

# 期待値の計算
expected_Y1 = np.mean(Y1_samples)
expected_Y3 = np.mean(Y3_samples)

# 理論値との比較
theoretical_Y1 = 1/4
theoretical_Y3 = 3/4

# 期待値を表示
print(f"シミュレーションによる E[Y1](最小値の期待値): {expected_Y1:.4f}, 理論値: {theoretical_Y1:.4f}")
print(f"シミュレーションによる E[Y3](最大値の期待値): {expected_Y3:.4f}, 理論値: {theoretical_Y3:.4f}")

# ヒストグラムの作成
plt.figure(figsize=(8, 6))

# 理論密度関数のプロット
plt.plot(y_values, f1_y, label='$f_1(y) = 3(1-y)^2$', color='blue')
plt.plot(y_values, f3_y, label='$f_3(y) = 3y^2$', color='green')

# サンプルのヒストグラムを追加 (確率密度として正規化)
plt.hist(Y1_samples, bins=50, density=True, alpha=0.5, color='blue', label='Y1 samples (min)')
plt.hist(Y3_samples, bins=50, density=True, alpha=0.5, color='green', label='Y3 samples (max)')

# グラフの設定
plt.title('確率密度関数 $f_1(y)$ と $f_3(y)$ (シミュレーションと理論)', fontsize=14)
plt.xlabel('y', fontsize=12)
plt.ylabel('密度', fontsize=12)
plt.legend(loc='upper right', fontsize=10)
plt.grid(True)

# グラフの表示
plt.show()
シミュレーションによる E[Y1](最小値の期待値): 0.2499, 理論値: 0.2500
シミュレーションによる E[Y3](最大値の期待値): 0.7487, 理論値: 0.7500

シミュレーションによる順序統計量Y1とY3の分布と期待値が理論値と一致しました。

2018 Q4(3)

マルコフ連鎖する条件付き正規分布の任意のtに対する分布の証明をやりました。

 

コード

X_tの条件付き分布がN(\rho^{2t} x_0, 1 - \rho^{4t})に従うのかシミュレーションをしてみます。

# 2018 Q4(3)  2024.10.20

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
t_max = 50  # 時間 t の最大値
n_simulations = 10000  # シミュレーション回数
rho = 0.7  # ρの値
x_0 = 1.0  # X_0 の初期値

# シミュレーション結果を格納するリスト
X_t_samples = np.zeros((t_max, n_simulations))
X_t_samples[0, :] = x_0  # X_0 = x_0

# マルコフプロセスをシミュレーション
for t in range(1, t_max):
    # X_t を次のステップで生成
    mean = rho**2 * X_t_samples[t-1, :]
    variance = 1 - rho**4
    X_t_samples[t, :] = np.random.normal(loc=mean, scale=np.sqrt(variance), size=n_simulations)

# シミュレーション結果の平均と分散を計算
simulated_means = np.mean(X_t_samples, axis=1)
simulated_variances = np.var(X_t_samples, axis=1)

# 理論値の期待値と分散を計算
theoretical_means = rho**(2 * np.arange(t_max)) * x_0
theoretical_variances = 1 - rho**(4 * np.arange(t_max))

# グラフを描画
plt.figure(figsize=(10, 6))

# 期待値のプロット
plt.subplot(2, 1, 1)
plt.plot(np.arange(t_max), simulated_means, label="シミュレーション期待値", color='blue', linestyle='--')
plt.plot(np.arange(t_max), theoretical_means, label="理論期待値", color='red')
plt.title("期待値の推移")
plt.xlabel("時間 t")
plt.ylabel("期待値 E[X_t]")
plt.legend()
plt.grid(True)

# 分散のプロット
plt.subplot(2, 1, 2)
plt.plot(np.arange(t_max), simulated_variances, label="シミュレーション分散", color='blue', linestyle='--')
plt.plot(np.arange(t_max), theoretical_variances, label="理論分散", color='red')
plt.title("分散の推移")
plt.xlabel("時間 t")
plt.ylabel("分散 Var[X_t]")
plt.legend()
plt.grid(True)

# グラフを表示
plt.tight_layout()
plt.show()

シミュレーションと理論値は一致しました。

2018 Q4(2)-2

マルコフ連鎖する条件付き正規分布の分散の問題をやりました。

 

コード

V_{X_{t+1}}[X_{t+1} \mid X_t = x_t] = 1 - \rho^4をシミュレーションしてみます。

#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

シミュレーション結果は理論値に一致しました。

2018 Q4(2)-1

マルコフ連鎖する条件付き正規分布の期待値の問題をやりました。

 

コード

E_{X_{t+1}}[X_{t+1} \mid X_t = x_t] = \rho^2 x_tをシミュレーションしてみます。

#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

シミュレーション結果は理論値に一致しました。