# 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