import numpy as np
import statsmodels.api as sm
from scipy import stats

np.random.seed(42)

# データ生成
n = 100
x1 = np.random.normal(0, 1, n)
x2 = np.random.normal(0, 1, n)
x3 = np.random.normal(0, 1, n)  # 影響のない変数

# 真のモデル: y = 2 + 3*x1 - 1.5*x2 + 0*x3 + ε
y = 2 + 3 * x1 - 1.5 * x2 + 0 * x3 + np.random.normal(0, 2, n)

# 計画行列（切片を追加）
X = np.column_stack([x1, x2, x3])
X_with_const = sm.add_constant(X)

# OLS回帰
model = sm.OLS(y, X_with_const).fit()
print(model.summary())

print("\n=== 回帰係数のt検定 ===")
for i, name in enumerate(['const', 'x1', 'x2', 'x3']):
    print(f"{name}: β̂ = {model.params[i]:.4f}, SE = {model.bse[i]:.4f}, "
          f"t = {model.tvalues[i]:.4f}, p = {model.pvalues[i]:.4f}")

print(f"\n=== モデル全体のF検定 ===")
print(f"F統計量 = {model.fvalue:.4f}")
print(f"p値 = {model.f_pvalue:.6f}")
print(f"R² = {model.rsquared:.4f}")
print(f"Adjusted R² = {model.rsquared_adj:.4f}")

# R²からFを算出して確認
p = 3
R2 = model.rsquared
F_from_R2 = (R2 / p) / ((1 - R2) / (n - p - 1))
print(f"\nR²から算出したF = {F_from_R2:.4f}（モデルのF = {model.fvalue:.4f}）")

import numpy as np
import statsmodels.api as sm

np.random.seed(42)

# データ生成
n = 100
x1 = np.random.normal(0, 1, n)
x2 = np.random.normal(0, 1, n)
y = 2 + 3 * x1 - 1.5 * x2 + np.random.normal(0, 2, n)

X = sm.add_constant(np.column_stack([x1, x2]))
model = sm.OLS(y, X).fit()

# 平方和の分解
y_hat = model.fittedvalues
y_bar = np.mean(y)
e = model.resid

TSS = np.sum((y - y_bar)**2)
SSR = np.sum((y_hat - y_bar)**2)
RSS = np.sum(e**2)

print("=== 平方和の分解 ===")
print(f"TSS (全変動)   = {TSS:.4f}")
print(f"SSR (回帰変動) = {SSR:.4f}")
print(f"RSS (残差変動) = {RSS:.4f}")
print(f"SSR + RSS      = {SSR + RSS:.4f}")
print(f"TSS = SSR + RSS? {np.isclose(TSS, SSR + RSS)}")

print(f"\n=== 検定統計量 ===")
p = 2  # 説明変数の数
MSR = SSR / p
MSE = RSS / (n - p - 1)
F = MSR / MSE

print(f"MSR = {MSR:.4f}")
print(f"MSE = {MSE:.4f}")
print(f"F統計量 = {F:.4f}")
print(f"R² = {SSR/TSS:.4f}")

# 交差項がゼロであることの確認
cross_term = np.sum(e * (y_hat - y_bar))
print(f"\n交差項 Σe_i(ŷ_i - ȳ) = {cross_term:.10f}（≈0）")

import numpy as np
from scipy import stats

def regression_t_test(X, y):
    """
    回帰係数のt検定をスクラッチで実装

    Parameters
    ----------
    X : ndarray
        計画行列（切片列を含む）
    y : ndarray
        応答変数

    Returns
    -------
    beta_hat : ndarray
        回帰係数の推定値
    se : ndarray
        標準誤差
    t_stats : ndarray
        t統計量
    p_values : ndarray
        p値
    """
    n, k = X.shape  # k = p + 1（切片を含む）

    # OLS推定
    XtX_inv = np.linalg.inv(X.T @ X)
    beta_hat = XtX_inv @ X.T @ y

    # 残差
    y_hat = X @ beta_hat
    e = y - y_hat

    # 誤差分散の推定（自由度 n - k）
    RSS = e.T @ e
    sigma2_hat = RSS / (n - k)

    # 回帰係数の共分散行列
    cov_beta = sigma2_hat * XtX_inv

    # 標準誤差
    se = np.sqrt(np.diag(cov_beta))

    # t統計量
    t_stats = beta_hat / se

    # p値（両側検定）
    df = n - k
    p_values = 2 * (1 - stats.t.cdf(np.abs(t_stats), df=df))

    return beta_hat, se, t_stats, p_values


# --- 具体例 ---
np.random.seed(42)

n = 100
x1 = np.random.normal(0, 1, n)
x2 = np.random.normal(0, 1, n)
x3 = np.random.normal(0, 1, n)

# 真のモデル: y = 2 + 3*x1 - 1.5*x2 + 0*x3 + ε
y = 2 + 3 * x1 - 1.5 * x2 + 0 * x3 + np.random.normal(0, 2, n)

# 計画行列（切片を追加）
X = np.column_stack([np.ones(n), x1, x2, x3])
names = ['const', 'x1', 'x2', 'x3']

beta_hat, se, t_stats, p_values = regression_t_test(X, y)

print("=== 回帰係数のt検定（スクラッチ実装） ===")
print(f"{'変数':<8} {'β̂':>10} {'SE':>10} {'t値':>10} {'p値':>10} {'有意?':>6}")
print("-" * 56)
for name, b, s, t, p in zip(names, beta_hat, se, t_stats, p_values):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    print(f"{name:<8} {b:>10.4f} {s:>10.4f} {t:>10.4f} {p:>10.4f} {sig:>6}")

print(f"\n真の値: β₀=2, β₁=3, β₂=-1.5, β₃=0")

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats

np.random.seed(42)

# データ生成
n = 150
x1 = np.random.normal(0, 1, n)
x2 = np.random.normal(0, 1, n)
y = 2 + 3 * x1 - 1.5 * x2 + np.random.normal(0, 2, n)

X = sm.add_constant(np.column_stack([x1, x2]))
model = sm.OLS(y, X).fit()

# 残差
residuals = model.resid
fitted = model.fittedvalues

# ハット行列の対角成分（てこ比）
H = X @ np.linalg.inv(X.T @ X) @ X.T
leverage = np.diag(H)

# 標準化残差
sigma_hat = np.sqrt(np.sum(residuals**2) / (n - 3))
standardized_resid = residuals / (sigma_hat * np.sqrt(1 - leverage))

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. 残差 vs 予測値
axes[0, 0].scatter(fitted, residuals, alpha=0.6, edgecolor='k', linewidth=0.5)
axes[0, 0].axhline(y=0, color='red', linestyle='--', linewidth=1.5)
axes[0, 0].set_xlabel('Fitted Values', fontsize=12)
axes[0, 0].set_ylabel('Residuals', fontsize=12)
axes[0, 0].set_title('Residuals vs Fitted Values', fontsize=13)
axes[0, 0].grid(True, alpha=0.3)

# 2. Q-Qプロット
stats.probplot(standardized_resid, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Normal Q-Q Plot', fontsize=13)
axes[0, 1].grid(True, alpha=0.3)

# 3. 標準化残差のヒストグラム
axes[1, 0].hist(standardized_resid, bins=25, density=True, alpha=0.7,
                color='steelblue', edgecolor='black')
x_range = np.linspace(-4, 4, 200)
axes[1, 0].plot(x_range, stats.norm.pdf(x_range), 'r-', linewidth=2,
                label='$N(0, 1)$')
axes[1, 0].set_xlabel('Standardized Residuals', fontsize=12)
axes[1, 0].set_ylabel('Density', fontsize=12)
axes[1, 0].set_title('Distribution of Standardized Residuals', fontsize=13)
axes[1, 0].legend(fontsize=11)
axes[1, 0].grid(True, alpha=0.3)

# 4. 標準化残差 vs てこ比
axes[1, 1].scatter(leverage, np.abs(standardized_resid), alpha=0.6,
                    edgecolor='k', linewidth=0.5)
axes[1, 1].axhline(y=2, color='red', linestyle='--', linewidth=1.5,
                    label='|Std. Residual| = 2')
axes[1, 1].axvline(x=2*(2+1)/n, color='orange', linestyle='--', linewidth=1.5,
                    label=f'Leverage = 2(p+1)/n = {2*3/n:.3f}')
axes[1, 1].set_xlabel('Leverage ($h_{ii}$)', fontsize=12)
axes[1, 1].set_ylabel('|Standardized Residuals|', fontsize=12)
axes[1, 1].set_title('Leverage vs |Standardized Residuals|', fontsize=13)
axes[1, 1].legend(fontsize=10)
axes[1, 1].grid(True, alpha=0.3)

plt.suptitle('Residual Diagnostics', fontsize=15, y=1.02)
plt.tight_layout()
plt.show()

# Shapiro-Wilk検定で正規性を確認
stat_sw, p_sw = stats.shapiro(residuals)
print(f"Shapiro-Wilk検定: W = {stat_sw:.4f}, p = {p_sw:.4f}")
print(f"正規性の帰無仮説を{'棄却しない' if p_sw > 0.05 else '棄却する'} (α=0.05)")

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats

np.random.seed(42)

# 単回帰のシミュレーション
n_simulations = 5000
n = 50
t_squared_vals = []
f_vals = []

for _ in range(n_simulations):
    x = np.random.normal(0, 1, n)
    # H1のもと: 真のβ = 1.0
    y = 1.0 + 1.0 * x + np.random.normal(0, 2, n)

    X = sm.add_constant(x)
    model = sm.OLS(y, X).fit()

    t_val = model.tvalues[1]  # x の t値
    f_val = model.fvalue       # F値

    t_squared_vals.append(t_val**2)
    f_vals.append(f_val)

t_squared_vals = np.array(t_squared_vals)
f_vals = np.array(f_vals)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# t² vs F の散布図
ax1.scatter(t_squared_vals, f_vals, alpha=0.3, s=10, c='steelblue')
ax1.plot([0, 30], [0, 30], 'r--', linewidth=2, label='$t^2 = F$ (identity)')
ax1.set_xlabel('$t^2$ (t-statistic squared)', fontsize=12)
ax1.set_ylabel('$F$ (F-statistic)', fontsize=12)
ax1.set_title('Simple Regression: $t^2$ vs $F$', fontsize=13)
ax1.legend(fontsize=11)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 30)
ax1.set_ylim(0, 30)

# 差のヒストグラム
diff = t_squared_vals - f_vals
ax2.hist(diff, bins=50, density=True, alpha=0.7, color='steelblue',
         edgecolor='black')
ax2.axvline(x=0, color='red', linestyle='--', linewidth=2)
ax2.set_xlabel('$t^2 - F$', fontsize=12)
ax2.set_ylabel('Density', fontsize=12)
ax2.set_title(f'$t^2 - F$ (mean = {np.mean(diff):.2e})', fontsize=13)
ax2.grid(True, alpha=0.3)

plt.suptitle('Equivalence of t-test and F-test in Simple Regression',
             fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# R²からF統計量を計算するデモ
n_values = [30, 50, 100, 200]
p = 3  # 説明変数の数
R2_range = np.linspace(0.01, 0.99, 100)

fig, ax = plt.subplots(figsize=(10, 6))

for n in n_values:
    F_values = (R2_range / p) / ((1 - R2_range) / (n - p - 1))
    ax.plot(R2_range, F_values, linewidth=2, label=f'n = {n}')

# 有意水準の臨界値を表示
for n in n_values:
    f_crit = stats.f.ppf(0.95, p, n - p - 1)
    R2_crit = (f_crit * p) / (f_crit * p + (n - p - 1))
    ax.plot(R2_crit, f_crit, 'rx', markersize=10, markeredgewidth=2)

ax.set_xlabel('$R^2$', fontsize=13)
ax.set_ylabel('F-statistic', fontsize=13)
ax.set_title(f'Relationship between $R^2$ and F-statistic (p = {p})', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 50)

plt.tight_layout()
plt.show()

# 各nでの臨界R²値を表示
print("=== 有意水準5%での臨界R²値 ===")
for n in n_values:
    f_crit = stats.f.ppf(0.95, p, n - p - 1)
    R2_crit = (f_crit * p) / (f_crit * p + (n - p - 1))
    print(f"n = {n}: F_crit = {f_crit:.4f}, R²_crit = {R2_crit:.4f}")