# 2015 Q3(4)  2024.12.16

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
r12_squared_vals = np.linspace(0, 0.99, 50)  # r12^2の値（0から0.99まで）
beta2_squared_vals = np.linspace(0, 1.5, 50)  # β2^2の値（0から1.5まで）
sigma2 = 1  # 誤差分散 σ^2
S11 = 10  # S11の仮定値
S12 = 5  # S12の仮定値
S22 = 15  # S22の仮定値

# グリッド生成
R12_squared, Beta2_squared = np.meshgrid(r12_squared_vals, beta2_squared_vals)

# MSEの差の計算
Var_beta2_hat = (sigma2 / S22) * (1 / (1 - R12_squared))  # Var(β2)
MSE_diff = (S12 / S11)**2 * (Var_beta2_hat - Beta2_squared)  # MSEの差

# 3Dプロットの描画
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')

# 3Dサーフェスプロット
surf = ax.plot_surface(R12_squared, Beta2_squared, MSE_diff, cmap='coolwarm', edgecolor='none')

# 軸ラベルとタイトル
ax.set_title(r'3Dプロット: $\mathrm{MSE}(\hat{\beta}_1) - \mathrm{MSE}(\tilde{\beta}_1)$ の符号反転', fontsize=14)
ax.set_xlabel(r'$r_{12}^2$ (相関係数の二乗)', fontsize=12)
ax.set_ylabel(r'$\beta_2^2$ (回帰係数の二乗)', fontsize=12)
ax.set_zlabel(r'$\mathrm{MSE}(\hat{\beta}_1) - \mathrm{MSE}(\tilde{\beta}_1)$', fontsize=12)

# カラーバーの追加
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label='MSEの差')

plt.tight_layout()
plt.show()

# 2015 Q3(4)  2024.12.15

import numpy as np
import matplotlib.pyplot as plt

# 1. パラメータ設定
r12_vals = np.linspace(-0.99, 0.99, 200)  # r12の値（-0.99から0.99まで）
sigma2 = 1  # 分散 σ^2
S11 = 10  # S11の仮定値

# 重回帰と単回帰の分散計算
Var_beta1_hat = (sigma2 / S11) * (1 / (1 - r12_vals**2))  # 重回帰（r12に依存）
Var_beta1_tilde = sigma2 / S11  # 単回帰（r12=0の分散）

# 2. グラフのプロット
plt.figure(figsize=(10, 6))

plt.plot(r12_vals, Var_beta1_hat, label=r'$\mathrm{Var}(\hat{\beta}_1)$ (重回帰)', color='blue', linewidth=2)
plt.axhline(y=Var_beta1_tilde, color='orange', linestyle='--', label=r'$\mathrm{Var}(\tilde{\beta}_1)$ (単回帰, $r_{12}=0$)')

# 軸ラベルとタイトル
plt.title(r'$\mathrm{Var}(\tilde{\beta}_1)$ と $\mathrm{Var}(\hat{\beta}_1)$ の比較', fontsize=14, fontweight='bold')
plt.xlabel(r'$r_{12}$（$x_1$と$x_2$の相関）', fontsize=12)
plt.ylabel(r'分散', fontsize=12)

# 凡例とグリッド
plt.legend(fontsize=10, loc='upper right')
plt.grid(alpha=0.5)

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

# 2015 Q3(4)  2024.12.15

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
r12_vals = np.linspace(-0.99, 0.99, 200)  # r12の値（-0.99から0.99まで）
sigma2 = 1  # 分散 σ^2
S11 = 10  # S11の仮定値

# r12^2 の計算
r12_squared_vals = r12_vals**2  # r12^2 の値

# 重回帰と単回帰の分散計算
Var_beta1_hat = (sigma2 / S11) * (1 / (1 - r12_squared_vals))  # 重回帰（r12に依存）
Var_beta1_tilde = sigma2 / S11  # 単回帰（r12=0の分散）

# 分散の差の計算
Var_diff = Var_beta1_hat - Var_beta1_tilde  # 分散の差

# グラフのプロット
plt.figure(figsize=(10, 6))

# 分散のプロット
plt.plot(r12_squared_vals, Var_beta1_hat, label=r'$\mathrm{Var}(\hat{\beta}_1)$ (重回帰)', color='blue', linewidth=2)
plt.axhline(y=Var_beta1_tilde, color='orange', linestyle='--', label=r'$\mathrm{Var}(\tilde{\beta}_1)$ (単回帰, $r_{12}=0$)')

# 分散の差をプロット
plt.plot(r12_squared_vals, Var_diff, label=r'$\mathrm{Var}(\hat{\beta}_1) - \mathrm{Var}(\tilde{\beta}_1)$', color='green', linewidth=2, linestyle=':')

# 軸ラベルとタイトル
plt.title(r'$\mathrm{Var}(\tilde{\beta}_1)$, $\mathrm{Var}(\hat{\beta}_1)$ とその差の比較 ($r_{12}^2$)', fontsize=14, fontweight='bold')
plt.xlabel(r'$r_{12}^2$（相関係数の二乗）', fontsize=12)
plt.ylabel(r'分散 / 差', fontsize=12)

# 凡例とグリッド
plt.legend(fontsize=10, loc='upper right')
plt.grid(alpha=0.5)

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

# 2015 Q3(3)  2024.12.14

import numpy as np
import matplotlib.pyplot as plt

# 1. パラメータ設定
r12_vals = np.linspace(-0.99, 0.99, 200)  # r12の値（-0.99から0.99まで）
sigma2 = 1  # 分散 σ^2
S11 = 10  # S11の仮定値
S22 = 15  # S22の仮定値

# 2. 分散の計算
Var_beta1 = (sigma2 / S11) * (1 / (1 - r12_vals**2))  # Var(β1)
Var_beta2 = (sigma2 / S22) * (1 / (1 - r12_vals**2))  # Var(β2)

# 3. グラフのプロット
plt.figure(figsize=(10, 6))

plt.plot(r12_vals, Var_beta1, label=r'$\mathrm{Var}(\hat{\beta}_1)$', color='blue', linewidth=2)
plt.plot(r12_vals, Var_beta2, label=r'$\mathrm{Var}(\hat{\beta}_2)$', color='orange', linewidth=2)
plt.axhline(y=sigma2 / S11, color='blue', linestyle='--', label=r'$r_{12}=0$ の $\mathrm{Var}(\hat{\beta}_1)$')
plt.axhline(y=sigma2 / S22, color='orange', linestyle='--', label=r'$r_{12}=0$ の $\mathrm{Var}(\hat{\beta}_2)$')

# 軸ラベルとタイトル
plt.title(r'$\mathrm{Var}(\hat{\beta}_1)$ と $\mathrm{Var}(\hat{\beta}_2)$ の $r_{12}$ に対する変化', fontsize=14, fontweight='bold')
plt.xlabel(r'$r_{12}$（$x_1$と$x_2$の相関）', fontsize=12)
plt.ylabel(r'分散', fontsize=12)

# 凡例とグリッド
plt.legend(fontsize=10, loc='upper right')
plt.grid(alpha=0.5)

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

# 2015 Q3(2)  2024.12.13

import numpy as np
import matplotlib.pyplot as plt

# 1. パラメータの設定
n = 100  # サンプルサイズ
r12_vals = np.linspace(-0.99, 0.99, 100)  # r12 を -0.99 から 0.99 まで変化させる
r1y = 0.5  # r1y を固定
r2y = 0.6  # r2y を固定

# 2. Beta1 の計算
beta1_vals = []
for r12 in r12_vals:
    if abs(r12) >= 1:
        beta1_vals.append(np.nan)  # r12 = ±1 の場合は計算不能
        continue
    beta1 = (r1y - r12 * r2y) / (1 - r12**2)  # 簡略化された Beta1 の符号条件
    beta1_vals.append(beta1)

# 3. 結果の可視化
plt.figure(figsize=(8, 6))
plt.plot(r12_vals, beta1_vals, label=r'$\hat{\beta}_1$', color='blue')
plt.axhline(0, color='red', linestyle='--', label='ゼロライン (基準)')
plt.title(r'$\hat{\beta}_1$ と $r_{12}$ の関係', fontsize=14, fontweight='bold')
plt.xlabel(r'$r_{12}$（$x_1$と$x_2$の相関）', fontsize=12)
plt.ylabel(r'$\hat{\beta}_1$', fontsize=12)
plt.legend(fontsize=10, loc='upper left')
plt.grid(alpha=0.7)
plt.tight_layout()
plt.show()

# 2015 Q3(1)  2024.12.12

import numpy as np

# 1. パラメータの設定 (再現性のため同じ設定を使用)
n = 100
beta_0, beta_1, beta_2 = 2.0, 1.0, -0.5
sigma = 1.0

# 2. 説明変数と誤差の生成
x1 = np.random.randn(n)
x2 = np.random.randn(n)
epsilon = np.random.randn(n) * sigma

# 3. 応答変数の生成
y = beta_0 + beta_1 * x1 + beta_2 * x2 + epsilon

# 4. S11, S22, S12, S1y, S2y を計算
S11 = np.sum(x1**2)
S22 = np.sum(x2**2)
S12 = np.sum(x1 * x2)
S1y = np.sum(x1 * y)
S2y = np.sum(x2 * y)

# 5. 推定値の計算 (導出した式に基づく)
denominator = S11 * S22 - S12**2
beta1_hat = (S22 * S1y - S12 * S2y) / denominator
beta2_hat = (S11 * S2y - S12 * S1y) / denominator
beta0_hat = np.mean(y) - beta1_hat * np.mean(x1) - beta2_hat * np.mean(x2)

# 推定値の表示
print(f"推定値 (導出した式):")
print(f"  β0 = {beta0_hat:.3f}")
print(f"  β1 = {beta1_hat:.3f}")
print(f"  β2 = {beta2_hat:.3f}")

# 真の値と比較
print("\n真の値:")
print(f"  β0 = {beta_0:.3f}")
print(f"  β1 = {beta_1:.3f}")
print(f"  β2 = {beta_2:.3f}")

推定値 (導出した式):
  β0 = 1.909
  β1 = 0.863
  β2 = -0.505

真の値:
  β0 = 2.000
  β1 = 1.000
  β2 = -0.500

# 2016 Q3(1)  2024.11.20

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 50  # サンプルサイズ
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)

# プロット
plt.figure(figsize=(10, 6))
plt.scatter(xi, Yi, label='観測データ ($Y_i$)', color='blue', marker='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.xlabel('$x_i$')
plt.ylabel('$Y_i$')
plt.title('線形モデルの推定')
plt.legend()
plt.grid(True)
plt.show()

# 2016 Q3(1)  2024.11.20

import numpy as np
import matplotlib.pyplot as plt

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

# サンプルサイズを変化させて b0, b1 の期待値を検証
sample_sizes = np.arange(5, 101, 5)  # n を 5 から 100 まで 5 刻みで変化
true_beta = 2  # 真のβ
simulations = 1000  # 各 n に対してのシミュレーション回数

b0_means = []
b1_means = []

# 各サンプルサイズ n に対してシミュレーションを実行
for n in sample_sizes:
    b0_samples = []
    b1_samples = []
    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))
    b0_means.append(np.mean(b0_samples))
    b1_means.append(np.mean(b1_samples))

# 結果をプロット
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.axhline(y=true_beta, color='green', linestyle='solid', label='真の値 $\\beta=2$')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('推定量の期待値')
plt.title('推定量 $b_0, b_1$ の期待値とサンプルサイズの関係')
plt.legend()
plt.grid(True)
plt.show()

# 2016 Q3(1)  2024.11.20

import numpy as np
import matplotlib.pyplot as plt

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

b0_samples = []
b1_samples = []

# シミュレーションを実施
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))

# b0 と b1 の分布を重ねてプロット
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.axvline(x=true_beta, color='green', linestyle='dashed', label='真の値 $\\beta=2$')

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