# 2016 Q3(3)  2024.11.22

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
sigma = 1  # 誤差の標準偏差
true_beta = 2  # 真のβ
sample_sizes = np.arange(10, 101, 10)  # サンプルサイズ n を 10 から 100 まで 10 刻みで変化

# 理論的な分散を格納するリスト
variances_b0 = []
variances_b1 = []
variances_b2 = []

# 各サンプルサイズ n に対して分散を計算
for n in sample_sizes:
    xi = np.arange(1, n + 1)  # 一様に分布した x_i を生成 (例: 1, 2, ..., n)
    sum_xi = np.sum(xi)
    sum_xi_squared = np.sum(xi**2)
    sum_1_over_xi_squared = np.sum(1 / xi**2)
    
    # 分散を理論式で計算
    V_b0 = (sigma**2 / n**2) * sum_1_over_xi_squared
    V_b1 = (sigma**2 * n) / (sum_xi**2)
    V_b2 = sigma**2 / sum_xi_squared
    
    variances_b0.append(V_b0)
    variances_b1.append(V_b1)
    variances_b2.append(V_b2)

# 理論値からの標準偏差として各分散の平方根を計算
std_b0 = np.sqrt(variances_b0)
std_b1 = np.sqrt(variances_b1)
std_b2 = np.sqrt(variances_b2)

# グラフ作成
plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, std_b0, label='$|b_0 - \\beta|$ (標準偏差)', marker='o', linestyle='dashed', color='red')
plt.plot(sample_sizes, std_b1, label='$|b_1 - \\beta|$ (標準偏差)', marker='x', linestyle='dashed', color='orange')
plt.plot(sample_sizes, std_b2, label='$|b_2 - \\beta|$ (標準偏差)', marker='s', linestyle='dashed', color='purple')

plt.axhline(y=0, 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(3)  2024.11.22

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
sigma = 1  # 誤差の標準偏差
sample_sizes = np.arange(10, 101, 10)  # サンプルサイズ n を 10 から 100 まで 10 刻みで変化

# 理論的な分散を格納するリスト
variances_b0 = []
variances_b1 = []
variances_b2 = []

# 各サンプルサイズ n に対して分散を計算
for n in sample_sizes:
    xi = np.arange(1, n + 1)  # 一様に分布した x_i を生成 (例: 1, 2, ..., n)
    sum_xi = np.sum(xi)
    sum_xi_squared = np.sum(xi**2)
    sum_1_over_xi_squared = np.sum(1 / xi**2)
    
    # 分散を理論式で計算
    V_b0 = (sigma**2 / n**2) * sum_1_over_xi_squared
    V_b1 = (sigma**2 * n) / (sum_xi**2)
    V_b2 = sigma**2 / sum_xi_squared
    
    variances_b0.append(V_b0)
    variances_b1.append(V_b1)
    variances_b2.append(V_b2)

# グラフ作成
plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, variances_b0, label='$V[b_0]$', marker='o', linestyle='dashed', color='red')
plt.plot(sample_sizes, variances_b1, label='$V[b_1]$', marker='x', linestyle='dashed', color='orange')
plt.plot(sample_sizes, variances_b2, label='$V[b_2]$', marker='s', linestyle='dashed', color='purple')

plt.xlabel('サンプルサイズ $n$')
plt.ylabel('分散')
plt.title('推定量 $b_0$, $b_1$, $b_2$ の理論的分散')
plt.legend()
plt.grid(True)
plt.show()

# 2016 Q3(3)  2024.11.22

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
sigma = 1  # 誤差の標準偏差
true_beta = 2  # 真のβ
sample_sizes = np.arange(10, 101, 10)  # サンプルサイズ n を 10 から 100 まで 10 刻みで変化

# 理論的な分散を格納するリスト
variances_b0 = []
variances_b1 = []
variances_b2 = []

# 各サンプルサイズ n に対して分散を計算
for n in sample_sizes:
    xi = np.arange(1, n + 1)  # 一様に分布した x_i を生成 (例: 1, 2, ..., n)
    sum_xi = np.sum(xi)
    sum_xi_squared = np.sum(xi**2)
    sum_1_over_xi_squared = np.sum(1 / xi**2)
    
    # 分散を理論式で計算
    V_b0 = (sigma**2 / n**2) * sum_1_over_xi_squared
    V_b1 = (sigma**2 * n) / (sum_xi**2)
    V_b2 = sigma**2 / sum_xi_squared
    
    variances_b0.append(V_b0)
    variances_b1.append(V_b1)
    variances_b2.append(V_b2)

# 理論値からの標準偏差として各分散の平方根を計算
std_b0 = np.sqrt(variances_b0)
std_b1 = np.sqrt(variances_b1)
std_b2 = np.sqrt(variances_b2)

# グラフ作成
plt.figure(figsize=(10, 6))
plt.plot(sample_sizes, std_b0, label='$|b_0 - \\beta|$ (標準偏差)', marker='o', linestyle='dashed', color='red')
plt.plot(sample_sizes, std_b1, label='$|b_1 - \\beta|$ (標準偏差)', marker='x', linestyle='dashed', color='orange')
plt.plot(sample_sizes, std_b2, label='$|b_2 - \\beta|$ (標準偏差)', marker='s', linestyle='dashed', color='purple')

plt.axhline(y=0, 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 = 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()