2015 Q3(4)
重回帰モデルの重みの推定量の分散と、説明変数を減らした時のそれとの差を求め、説明変数間の相関係数によるどのような関数になるのか調べました。

コード
の変化が
と
に与える影響を、グラフで描画して確認します。
# 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()

が0に近いほど
は最小となり、相関が高まると分散が増加する一方で、
は一定であることが確認されました。
次に、横軸をとして、
と
の変化をプロットし、両者の差も併せて確認します。
# 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)
重回帰モデルの重みの推定量の分散を求め、説明変数間の相関がそれにどのように影響するか確認しました。

コード
の変化が
と
に与える影響を、グラフで描画して確認します。
# 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()

が0に近いほど
と
は最小となり、相関が高まると分散が増加することが確認されました。
2015 Q3(2)
重回帰モデルの重みβ1を3つの相関係数を使って表し、それが負になる必要十分条件を求めました。

コード
が大きな値をとる場合、
が負の値をとりやすくなることを確認するため、
と
は、正の数に固定した上で、グラフで可視化します。
# 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()

が1に近い値をとると、
が負の値をとる傾向があることが確認できました。
2015 Q3(1)
重回帰モデルにおいて正規方程式を用い、各重みの最小二乗推定量を求めました。

コード
重回帰モデルのシミュレーションを行い、パラメータβ0,β1,β2の推定量を計算し、真の値と比較します。
# 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
シミュレーションの結果、パラメータβ0,β1,β2の推定量はわずかな誤差はあるものの、真の値に近い値を示しました。
2015 Q2(5)
前問の検定がネイマン-ピアソンの基本定理に基づき一様最強力検定であることを示しました。

コード
有意水準α=0.05のもと、検出力を理論値とシミュレーションで描画し比較します。
# 2015 Q2(5) 2024.12.11
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# シミュレーション設定
alpha = 0.05 # 有意水準
z_alpha = norm.ppf(1 - alpha) # 臨界値 z_alpha
num_simulations = 10000 # シミュレーション回数
n = 10 # 標本サイズ
mu_values = np.linspace(0, 2, 10) # 真の平均 μ の範囲
# 臨界値の計算
critical_value = z_alpha / np.sqrt(n) # 棄却域の閾値
# 検出力を計算する関数
def compute_power_for_mu(mu, n, z_alpha):
return 1 - norm.cdf(z_alpha - mu * np.sqrt(n))
# シミュレーションで検出力を計算
empirical_powers = []
for mu in mu_values:
samples = np.random.normal(mu, 1 / np.sqrt(n), num_simulations) # 標本生成
rejection_rate = np.mean(samples > critical_value) # 棄却割合
empirical_powers.append(rejection_rate)
# 理論的な検出力を計算
theoretical_powers = [compute_power_for_mu(mu, n, z_alpha) for mu in mu_values]
# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(mu_values, theoretical_powers, label="理論的検出力", linestyle="--", color="blue")
plt.scatter(mu_values, empirical_powers, label="シミュレーションによる検出力", color="red", zorder=5)
plt.axhline(y=0.05, color="green", linestyle="--", label="有意水準 $\\alpha = 0.05$")
# グラフ設定
plt.xlabel("真の平均 $\\mu$", fontsize=12)
plt.ylabel("検出力 (Power)", fontsize=12)
plt.title("理論値とシミュレーションによる検出力の比較", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

検出力は理論値とシミュレーションでよく一致し、検定が正しく機能していることを確認しました。また、真の平均μ>0に対して検出力が単調増加することが分かり、この検定が一様最強力検定の特徴を備えていることが確認できました。