# 2016 Q3(2)  2024.11.21

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 10  # サンプルサイズ
true_beta = 2  # 真のβ
sigma = 5  # 誤差の標準偏差

# 1. xi を一様分布から生成
xi = np.random.uniform(1, 10, n)

# 2. εi を正規分布から生成
epsilon_i = np.random.normal(0, sigma, n)

# 3. Yi を計算
Yi = true_beta * xi + epsilon_i

# 4. 推定量 b0, b1 を計算
b0 = np.mean(Yi / xi)
b1 = np.sum(Yi) / np.sum(xi)

# 5. 最小二乗推定量 b2 を計算
b2 = np.sum(xi * Yi) / np.sum(xi**2)

# プロット
plt.figure(figsize=(10, 6))
plt.scatter(xi, Yi, label='観測データ ($Y_i$)', color='blue', marker='x')  # マーカーを 'x' に変更
plt.plot(np.sort(xi), true_beta * np.sort(xi), label=f'真の直線 ($\\beta={true_beta}$)', color='green', linestyle='dashed')
plt.plot(np.sort(xi), b0 * np.sort(xi), label=f'推定直線 ($b_0={b0:.2f}$)', color='red')
plt.plot(np.sort(xi), b1 * np.sort(xi), label=f'推定直線 ($b_1={b1:.2f}$)', color='orange')
plt.plot(np.sort(xi), b2 * np.sort(xi), label=f'推定直線 ($b_2={b2:.2f}$)', color='purple')

plt.xlabel('$x_i$')
plt.ylabel('$Y_i$')
plt.title('線形モデルの推定')
plt.legend()
plt.grid(True)
plt.show()

# 2016 Q3(2)  2024.11.21

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
sigma = 5  # 誤差の標準偏差
true_beta = 2  # 真のβ
sample_sizes = np.arange(5, 101, 5)  # n を 5 から 100 まで 5 刻みで変化
simulations = 1000  # 各 n に対してのシミュレーション回数

b0_means = []
b1_means = []
b2_means = []  # b2 の期待値を記録

# 各サンプルサイズ n に対してシミュレーションを実行
for n in sample_sizes:
    b0_samples = []
    b1_samples = []
    b2_samples = []  # b2 を記録
    for _ in range(simulations):
        xi = np.random.uniform(1, 10, n)
        epsilon_i = np.random.normal(0, sigma, n)
        Yi = true_beta * xi + epsilon_i
        b0_samples.append(np.mean(Yi / xi))
        b1_samples.append(np.sum(Yi) / np.sum(xi))
        b2_samples.append(np.sum(xi * Yi) / np.sum(xi**2))  # b2 を計算
    b0_means.append(np.mean(b0_samples))
    b1_means.append(np.mean(b1_samples))
    b2_means.append(np.mean(b2_samples))  # b2 の期待値を記録

# 結果をプロット
plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, b0_means, label='$b_0$の期待値', marker='o', linestyle='dashed', color='red')
plt.plot(sample_sizes, b1_means, label='$b_1$の期待値', marker='x', linestyle='dashed', color='orange')
plt.plot(sample_sizes, b2_means, label='$b_2$の期待値', marker='s', linestyle='dashed', color='purple')  # b2 を追加
plt.axhline(y=true_beta, color='green', linestyle='solid', label='真の値 $\\beta=2$')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('推定量の期待値')
plt.title('推定量 $b_0, b_1, b_2$ の期待値とサンプルサイズの関係')
plt.legend()
plt.grid(True)
plt.show()

# 2016 Q3(2)  2024.11.21

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 100  # サンプルサイズ
true_beta = 2  # 真のβ
sigma = 1  # 誤差の標準偏差
simulations = 10000  # シミュレーション回数

b0_samples = []
b1_samples = []
b2_samples = []  # b2 の分布を記録

# シミュレーションを実施
for _ in range(simulations):
    xi = np.random.uniform(1, 10, n)
    epsilon_i = np.random.normal(0, sigma, n)
    Yi = true_beta * xi + epsilon_i
    b0_samples.append(np.mean(Yi / xi))
    b1_samples.append(np.sum(Yi) / np.sum(xi))
    b2_samples.append(np.sum(xi * Yi) / np.sum(xi**2))  # b2 を計算

# b0, b1, b2 の分布を重ねてプロット
plt.figure(figsize=(10, 6))
plt.hist(b0_samples, bins=30, density=True, color='red', alpha=0.5, label='$b_0$')
plt.hist(b1_samples, bins=30, density=True, color='orange', alpha=0.5, label='$b_1$')
plt.hist(b2_samples, bins=30, density=True, color='purple', alpha=0.5, label='$b_2$')  # b2 を追加
plt.axvline(x=true_beta, color='green', linestyle='dashed', label='真の値 $\\beta=2$')

plt.title('$b_0$, $b_1$, $b_2$ の分布 (n=100)')
plt.xlabel('推定値')
plt.ylabel('頻度')
plt.legend()
plt.grid(True)
plt.show()