import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 1000  # サンプル数
mu, sigma = 0, 1  # 正規分布の平均と標準偏差
iterations = 10000  # シミュレーション回数

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

# シミュレーションを実行
for _ in range(iterations):
    X_i = np.random.normal(mu, sigma, size=n)
    X_j = np.random.normal(mu, sigma, size=n)
    X_k = np.random.normal(mu, sigma, size=n)
    
    Y_i = X_i - mu
    Y_j = X_j - mu
    Y_k = X_k - mu
    
    Y_i_cubed_values.append(np.mean(Y_i**3))
    Y_i_squared_Y_j_values.append(np.mean(Y_i**2 * Y_j))
    Y_i_Y_j_Y_k_values.append(np.mean(Y_i * Y_j * Y_k))

# グラフを横に並べて表示
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# (1) E[Y_i^3] の分布
axes[0].hist(Y_i_cubed_values, bins=50, alpha=0.75)
axes[0].axvline(np.mean(Y_i_cubed_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[0].set_title(r"正規分布における $E[Y_i^3]$ の分布")
axes[0].set_xlabel('値')
axes[0].set_ylabel('頻度')
axes[0].legend()

# (2) E[Y_i^2 Y_j] の分布
axes[1].hist(Y_i_squared_Y_j_values, bins=50, alpha=0.75)
axes[1].axvline(np.mean(Y_i_squared_Y_j_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[1].set_title(r"正規分布における $E[Y_i^2 Y_j]$ の分布")
axes[1].set_xlabel('値')
axes[1].set_ylabel('頻度')
axes[1].legend()

# (3) E[Y_i Y_j Y_k] の分布
axes[2].hist(Y_i_Y_j_Y_k_values, bins=50, alpha=0.75)
axes[2].axvline(np.mean(Y_i_Y_j_Y_k_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[2].set_title(r"正規分布における $E[Y_i Y_j Y_k]$ の分布")
axes[2].set_xlabel('値')
axes[2].set_ylabel('頻度')
axes[2].legend()

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

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 1000  # サンプル数
scale = 1.0  # 指数分布のスケールパラメータ（1/λ）
iterations = 10000  # シミュレーション回数

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

# シミュレーションを実行
for _ in range(iterations):
    X_i = np.random.exponential(scale=scale, size=n)
    X_j = np.random.exponential(scale=scale, size=n)
    X_k = np.random.exponential(scale=scale, size=n)
    
    Y_i = X_i - np.mean(X_i)
    Y_j = X_j - np.mean(X_j)
    Y_k = X_k - np.mean(X_k)
    
    Y_i_cubed_values.append(np.mean(Y_i**3))
    Y_i_squared_Y_j_values.append(np.mean(Y_i**2 * Y_j))
    Y_i_Y_j_Y_k_values.append(np.mean(Y_i * Y_j * Y_k))

# グラフを横に並べて表示
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# (1) E[Y_i^3] の分布
axes[0].hist(Y_i_cubed_values, bins=50, alpha=0.75)
axes[0].axvline(np.mean(Y_i_cubed_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[0].set_title(r"指数分布における $E[Y_i^3]$ の分布")
axes[0].set_xlabel('値')
axes[0].set_ylabel('頻度')
axes[0].legend()

# (2) E[Y_i^2 Y_j] の分布
axes[1].hist(Y_i_squared_Y_j_values, bins=50, alpha=0.75)
axes[1].axvline(np.mean(Y_i_squared_Y_j_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[1].set_title(r"指数分布における $E[Y_i^2 Y_j]$ の分布")
axes[1].set_xlabel('値')
axes[1].set_ylabel('頻度')
axes[1].legend()

# (3) E[Y_i Y_j Y_k] の分布
axes[2].hist(Y_i_Y_j_Y_k_values, bins=50, alpha=0.75)
axes[2].axvline(np.mean(Y_i_Y_j_Y_k_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[2].set_title(r"指数分布における $E[Y_i Y_j Y_k]$ の分布")
axes[2].set_xlabel('値')
axes[2].set_ylabel('頻度')
axes[2].legend()

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

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 1000  # サンプル数
shape, scale = 2.0, 2.0  # ガンマ分布の形状パラメータとスケールパラメータ
iterations = 10000  # シミュレーション回数

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

# シミュレーションを実行
for _ in range(iterations):
    X_i = np.random.gamma(shape, scale, size=n)
    X_j = np.random.gamma(shape, scale, size=n)
    X_k = np.random.gamma(shape, scale, size=n)
    
    Y_i = X_i - np.mean(X_i)
    Y_j = X_j - np.mean(X_j)
    Y_k = X_k - np.mean(X_k)
    
    Y_i_cubed_values.append(np.mean(Y_i**3))
    Y_i_squared_Y_j_values.append(np.mean(Y_i**2 * Y_j))
    Y_i_Y_j_Y_k_values.append(np.mean(Y_i * Y_j * Y_k))

# グラフを横に並べて表示
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# (1) E[Y_i^3] の分布
axes[0].hist(Y_i_cubed_values, bins=50, alpha=0.75)
axes[0].axvline(np.mean(Y_i_cubed_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[0].set_title(r"ガンマ分布における $E[Y_i^3]$ の分布")
axes[0].set_xlabel('値')
axes[0].set_ylabel('頻度')
axes[0].legend()

# (2) E[Y_i^2 Y_j] の分布
axes[1].hist(Y_i_squared_Y_j_values, bins=50, alpha=0.75)
axes[1].axvline(np.mean(Y_i_squared_Y_j_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[1].set_title(r"ガンマ分布における $E[Y_i^2 Y_j]$ の分布")
axes[1].set_xlabel('値')
axes[1].set_ylabel('頻度')
axes[1].legend()

# (3) E[Y_i Y_j Y_k] の分布
axes[2].hist(Y_i_Y_j_Y_k_values, bins=50, alpha=0.75)
axes[2].axvline(np.mean(Y_i_Y_j_Y_k_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[2].set_title(r"ガンマ分布における $E[Y_i Y_j Y_k]$ の分布")
axes[2].set_xlabel('値')
axes[2].set_ylabel('頻度')
axes[2].legend()

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

# 2021 Q4(1)  2024.8.31

import numpy as np

# パラメータの設定
n = 100  # サンプル数
mu = 0  # 平均
sigma = 1  # 標準偏差
iterations = 10000  # シミュレーション回数

# 各シミュレーションで得た標本平均を格納するリスト
sample_means = []

for _ in range(iterations):
    # 正規分布に従う乱数を生成
    X = np.random.normal(loc=mu, scale=sigma, size=n)
    
    # 標本平均を計算
    X_bar = np.mean(X)
    sample_means.append(X_bar)

# 標本平均の分散を計算
estimated_variance_of_sample_mean = np.var(sample_means)

# 理論的な分散
theoretical_variance = sigma**2 / n

print(f"理論的な分散: {theoretical_variance}")
print(f"推定された標本平均の分散: {estimated_variance_of_sample_mean}")

理論的な分散: 0.01
推定された標本平均の分散: 0.009838287336429174

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 100  # サンプル数
mu = 0  # 平均
sigma = 1  # 標準偏差
iterations = 10000  # シミュレーション回数

# 分散の推定値を格納するリスト
estimated_variances = []

for _ in range(iterations):
    # 正規分布に従う乱数を生成
    X = np.random.normal(loc=mu, scale=sigma, size=n)
    
    # 標本平均の分散を計算
    estimated_variance = np.var(X) / n
    estimated_variances.append(estimated_variance)

# 理論的な分散
theoretical_variance = sigma**2 / n

# 視覚化
plt.figure(figsize=(10, 6))
plt.hist(estimated_variances, bins=100, alpha=0.75, label='推定された分散', range=(0.005, 0.015))
plt.axvline(theoretical_variance, color='red', linestyle='dashed', linewidth=2, label=f'理論的な分散 = {theoretical_variance}')
plt.xlabel('分散')
plt.ylabel('頻度')
plt.title(f'{iterations}回のシミュレーションにおける推定された分散の分布\n(サンプルサイズ = {n})')
plt.legend()
plt.grid(True)
plt.show()

# 2021 Q4(1)  2024.8.31

import numpy as np

# パラメータの設定
n_simulations = 10000  # シミュレーション回数
mu = 0  # 平均
sigma = 1  # 標準偏差

# X1, X2を正規分布から生成
X1 = np.random.normal(loc=mu, scale=sigma, size=n_simulations)
X2 = np.random.normal(loc=mu, scale=sigma, size=n_simulations)

# (X1 - X2)^3を計算
expression = (X1 - X2)**3

# 期待値を計算
expected_value = np.mean(expression)

print(f"(X1 - X2)^3 の期待値: {expected_value}")

(X1 - X2)^3 の期待値: 0.044962030052379295

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 10000  # シミュレーション回数
mu = 0  # 平均
sigma = 1  # 標準偏差

# (X1 - X2)^3 の推定値を格納するリスト
expected_values = []

for _ in range(n):
    # 正規分布に従う X1, X2 を生成
    X1 = np.random.normal(loc=mu, scale=sigma)
    X2 = np.random.normal(loc=mu, scale=sigma)
    
    # (X1 - X2)^3 を計算
    value = (X1 - X2)**3
    expected_values.append(value)

# 理論的な期待値
theoretical_expectation = 0

# 視覚化
plt.figure(figsize=(10, 6))
plt.hist(expected_values, bins=100, alpha=0.75, label='(X1 - X2)^3 の推定値', range=(-30, 30))
plt.axvline(theoretical_expectation, color='red', linestyle='dashed', linewidth=2, label=f'理論的な期待値 = {theoretical_expectation}')
plt.xlabel('値')
plt.ylabel('頻度')
plt.title(f'{n}回のシミュレーションにおける (X1 - X2)^3 の期待値の分布')
plt.legend()
plt.grid(True)
plt.show()