단순 OLS 추정량의 신뢰구간¶

from scipy import stats

def main():
    """회귀 출력표의 기울기와 표준오차만으로 신뢰구간을 만든다.

    표준오차는 이미 주어진 값이라고 보고, t 임계값만 자유도 n-2 로 구한다.
    """
    beta_1_hat = 0.164
    n = 20
    df = n - 2
    confidence_level = 0.95
    alpha = 1 - confidence_level
    t_star = stats.t(df).ppf(1 - alpha / 2)
    standard_error = 0.057
    margin_of_error = t_star * standard_error
    print(f"{confidence_level:.0%} confidence interval of the slope")
    print(f"{beta_1_hat:.4f} ± {margin_of_error:.4f}")

if __name__ == "__main__":
    main()

95% confidence interval of the slope
0.1640 ± 0.1198

from scipy import stats

beta_1_hat, se_reported, n = 0.164, 0.057, 20
df = n - 2
t_star, z_star = stats.t(df).ppf(0.975), stats.norm.ppf(0.975)
print(f"t_{df}(0.975) = {t_star:.6f},   z(0.975) = {z_star:.6f},   비 = {t_star / z_star:.4f}")
print(f"정규로 대신하면 오차한계가 {z_star * se_reported:.4f} 로 "
      f"{t_star * se_reported:.4f} 보다 {1 - z_star / t_star:.1%} 짧아진다")
print()
print(f"표의 T = 2.862,   표의 Coef/SE = {beta_1_hat / se_reported:.4f}")
print(f"  T = 2.862 를 믿으면 SE = 0.164/2.862 = {beta_1_hat / 2.862:.6f}")
print(f"  p-값:  T=2.862 -> {2 * stats.t(df).sf(2.862):.4f},"
      f"   T={beta_1_hat / se_reported:.4f} -> {2 * stats.t(df).sf(beta_1_hat / se_reported):.4f}")
print(f"  R^2 = T^2/(T^2+df) = {2.862 ** 2 / (2.862 ** 2 + df):.4f}")
print()
lo, hi = beta_1_hat - t_star * se_reported, beta_1_hat + t_star * se_reported
print(f"95% CI = ({lo:.4f}, {hi:.4f})")
print(f"  0 을 담는가? {lo <= 0 <= hi}")
print(f"  |T| = {beta_1_hat / se_reported:.4f} > t_crit = {t_star:.4f} 인가? "
      f"{beta_1_hat / se_reported > t_star}")

t_18(0.975) = 2.100922,   z(0.975) = 1.959964,   비 = 1.0719
정규로 대신하면 오차한계가 0.1117 로 0.1198 보다 6.7% 짧아진다

표의 T = 2.862,   표의 Coef/SE = 2.8772
  T = 2.862 를 믿으면 SE = 0.164/2.862 = 0.057303
  p-값:  T=2.862 -> 0.0104,   T=2.8772 -> 0.0100
  R^2 = T^2/(T^2+df) = 0.3127

95% CI = (0.0442, 0.2838)
  0 을 담는가? False
  |T| = 2.8772 > t_crit = 2.1009 인가? True

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

# 이 쪽의 계산에서 scipy 가 꼭 필요한 곳은 t 분위수 하나다.
z = stats.norm.ppf(0.975)
print(f"{'자유도':>8s}{'t(0.975)':>12s}{'정규 z':>10s}{'z + (z^3+z)/(4v)':>20s}{'보정의 오차':>14s}")
for nu in (8, 18, 48, 98, 998):
    t_exact = stats.t(nu).ppf(0.975)
    t_approx = z + (z ** 3 + z) / (4 * nu)
    print(f"{nu:8d}{t_exact:12.6f}{z:10.6f}{t_approx:20.6f}{t_approx - t_exact:+14.6f}")
print()
print("numpy 만으로 할 수 있는 일과 못 하는 일")
print(f"  평균·표준편차·상관: numpy 로 충분  (예: corrcoef)")
print(f"  t 분위수: 역함수가 초기함수가 아니므로 수치해법이 필요하다")
print(f"  그림: matplotlib 는 통계와 무관하고 그리기만 맡는다")

     자유도    t(0.975)      정규 z    z + (z^3+z)/(4v)        보정의 오차
       8    2.306004  1.959964            2.256498     -0.049506
      18    2.100922  1.959964            2.091757     -0.009165
      48    2.010635  1.959964            2.009386     -0.001248
      98    1.984467  1.959964            1.984171     -0.000297
     998    1.962344  1.959964            1.962341     -0.000003

numpy 만으로 할 수 있는 일과 못 하는 일
  평균·표준편차·상관: numpy 로 충분  (예: corrcoef)
  t 분위수: 역함수가 초기함수가 아니므로 수치해법이 필요하다
  그림: matplotlib 는 통계와 무관하고 그리기만 맡는다

def generate_data(n, sigma, seed=0):
    """모의 회귀자료를 만든다. 참 모형은 y = 1 + 2x + 잡음 이다.

    매개변수
    --------
    n : 관측 수
    sigma : 잡음의 표준편차
    seed : 난수 씨앗

    돌려주는 값
    ----------
    x, y : 모양 (n, 1) 인 배열
    """
    np.random.seed(seed)
    x = np.random.randn(n, 1)
    y = 1 + 2 * x + sigma * np.random.randn(n, 1)
    return x, y

# 참 모형이 y = 1 + 2x + 3*eps,  x ~ N(0,1) 이므로 이론값을 먼저 적을 수 있다.
n_demo, sigma_demo = 100, 3
print(f"이론:  Var(y) = 2^2*Var(x) + sigma^2 = 4 + {sigma_demo ** 2} = {4 + sigma_demo ** 2}")
print(f"       rho^2  = 4/(4+{sigma_demo ** 2}) = {4 / (4 + sigma_demo ** 2):.6f}")
print(f"       E[S_xx] = (n-1)*Var(x) = {n_demo - 1}")
print(f"       se(beta1) = sigma/sqrt(S_xx) = {sigma_demo}/sqrt({n_demo - 1}) "
      f"= {sigma_demo / np.sqrt(n_demo - 1):.6f}")
print()
x_demo, y_demo = generate_data(n_demo, sigma_demo)
S_xx_demo = ((x_demo - x_demo.mean()) ** 2).sum()
r_demo = np.corrcoef(np.concatenate([x_demo, y_demo], axis=1), rowvar=False)[1, 0]
print(f"표본:  Var(y) = {y_demo.var(ddof=1):.4f}")
print(f"       r^2    = {r_demo ** 2:.6f}")
print(f"       S_xx   = {S_xx_demo:.4f}")
print(f"       x-bar  = {x_demo.mean():+.6f}  (이론 0),"
      f"   s_x = {x_demo.std(ddof=1):.6f}  (이론 1)")
print()
print(f"상대 어긋남:  Var(y) {y_demo.var(ddof=1) / 13 - 1:+.2%},"
      f"  r^2 {r_demo ** 2 / (4 / 13) - 1:+.2%},  S_xx {S_xx_demo / 99 - 1:+.2%}")
print(f"분산 추정의 상대 표준오차 ~ sqrt(2/(n-1)) = {np.sqrt(2 / (n_demo - 1)):.4f}")

이론:  Var(y) = 2^2*Var(x) + sigma^2 = 4 + 9 = 13
       rho^2  = 4/(4+9) = 0.307692
       E[S_xx] = (n-1)*Var(x) = 99
       se(beta1) = sigma/sqrt(S_xx) = 3/sqrt(99) = 0.301511

표본:  Var(y) = 15.2488
       r^2    = 0.369743
       S_xx   = 101.5827
       x-bar  = +0.059808  (이론 0),   s_x = 1.012960  (이론 1)

상대 어긋남:  Var(y) +17.30%,  r^2 +20.17%,  S_xx +2.61%
분산 추정의 상대 표준오차 ~ sqrt(2/(n-1)) = 0.1421

def estimate_regression_line(x, y):
    """상관계수 공식으로 기울기와 절편을 추정한다.

    beta_hat = r * (s_y / s_x) 이다. 최소제곱해와 같은 값이지만, 회귀계수가
    상관계수를 두 변수의 척도로 되돌린 것임이 이 꼴에서 드러난다.

    돌려주는 값
    ----------
    y_hat : 적합값
    beta_hat : 추정된 기울기
    y_bar, x_bar : 표본평균
    """
    x_bar = x.mean()
    y_bar = y.mean()
    s_x = x.std(ddof=1)
    s_y = y.std(ddof=1)
    r = np.corrcoef(np.concatenate([x, y], axis=1), rowvar=False)[1, 0]
    beta_hat = r * s_y / s_x
    y_hat = beta_hat * (x - x_bar) + y_bar
    return y_hat, beta_hat, y_bar, x_bar

# 세 가지 꼴이 같은 수를 주는지 확인한다.
y_hat_demo, beta_hat_demo, y_bar_demo, x_bar_demo = estimate_regression_line(x_demo, y_demo)
S_xy_demo = ((x_demo - x_bar_demo) * (y_demo - y_bar_demo)).sum()
print(f"r * s_y / s_x   = {beta_hat_demo:.10f}")
print(f"S_xy / S_xx     = {S_xy_demo / S_xx_demo:.10f}")
print(f"두 값의 차이     = {abs(beta_hat_demo - S_xy_demo / S_xx_demo):.2e}")
print()
beta_0_demo = y_bar_demo - beta_hat_demo * x_bar_demo
print(f"절편 = y-bar - b*x-bar = {beta_0_demo:.6f}   (참값 1)")
print(f"  x = 0 에서의 적합값  = {float(beta_hat_demo * (0 - x_bar_demo) + y_bar_demo):.6f}")
print(f"  x-bar 에서의 적합값  = {float(beta_hat_demo * (x_bar_demo - x_bar_demo) + y_bar_demo):.6f}"
      f"   (= y-bar = {y_bar_demo:.6f})")
print()
s2_demo = ((y_demo - y_hat_demo) ** 2).sum() / (n_demo - 2)
se_beta = np.sqrt(s2_demo / S_xx_demo)
t_c = stats.t(n_demo - 2).ppf(0.975)
print(f"기울기 추정값 {beta_hat_demo:.6f}   (참값 2)")
print(f"  se = s/sqrt(S_xx) = {se_beta:.6f}   (보기 3 의 이론값 0.301511)")
print(f"  95% CI = ({beta_hat_demo - t_c * se_beta:.6f}, {beta_hat_demo + t_c * se_beta:.6f})"
      f"   참값 2 를 담는가? {beta_hat_demo - t_c * se_beta <= 2 <= beta_hat_demo + t_c * se_beta}")
print(f"  (추정값 - 참값)/se = {(beta_hat_demo - 2) / se_beta:+.4f}")

r * s_y / s_x   = 2.3440953006
S_xy / S_xx     = 2.3440953006
두 값의 차이     = 1.33e-15

절편 = y-bar - b*x-bar = 1.225459   (참값 1)
  x = 0 에서의 적합값  = 1.225459
  x-bar 에서의 적합값  = 1.365655   (= y-bar = 1.365655)

기울기 추정값 2.344095   (참값 2)
  se = s/sqrt(S_xx) = 0.309151   (보기 3 의 이론값 0.301511)
  95% CI = (1.730594, 2.957596)   참값 2 를 담는가? True
  (추정값 - 참값)/se = +1.1130

def calculate_residual_variance(y, y_hat, n):
    """잔차분산 s^2 과 그 제곱근 s 를 구한다.

    n 이 아니라 n-2 로 나눈다. 절편과 기울기 둘을 자료에서 추정하느라
    자유도를 둘 잃었기 때문이다. 그래야 s^2 이 sigma^2 의 불편추정량이 된다.
    """
    s_square = np.sum((y - y_hat) ** 2) / (n - 2)
    s = np.sqrt(s_square)
    return s_square, s

# 이 표본의 s^2, 그리고 n 으로 나눈 값
s_square_demo, s_demo = calculate_residual_variance(y_demo, y_hat_demo, n_demo)
rss_demo = ((y_demo - y_hat_demo) ** 2).sum()
print(f"RSS = {rss_demo:.4f}")
print(f"  n-2 로 나눈 s^2 = {s_square_demo:.6f}   (참 sigma^2 = 9)")
print(f"  n   로 나눈 값  = {rss_demo / n_demo:.6f}")
print(f"  이론 치우침 -2*sigma^2/n = {-2 * 9 / n_demo:.4f}")
print(f"  s^2 의 상대 표준오차 sqrt(2/(n-2)) = {np.sqrt(2 / (n_demo - 2)):.4f}")
print(f"  이 표본의 어긋남 {s_square_demo / 9 - 1:+.4f} = "
      f"{(s_square_demo / 9 - 1) / np.sqrt(2 / (n_demo - 2)):+.3f} 표준오차")
print()
# 2000 번 되풀이해 두 추정량의 평균을 본다. 몬테카를로 오차를 함께 적는다.
reps = 2000
rng = np.random.default_rng(1)
s2_list = np.empty(reps)
for b in range(reps):
    xb_ = rng.standard_normal((n_demo, 1))
    yb_ = 1 + 2 * xb_ + sigma_demo * rng.standard_normal((n_demo, 1))
    yh_, _, _, _ = estimate_regression_line(xb_, yb_)
    s2_list[b] = ((yb_ - yh_) ** 2).sum() / (n_demo - 2)
mc_se = s2_list.std(ddof=1) / np.sqrt(reps)
print(f"모의실험 {reps} 회")
print(f"  mean(s^2)        = {s2_list.mean():.4f} ± {mc_se:.4f}   (참값 9)")
print(f"  mean(RSS/n)      = {s2_list.mean() * (n_demo - 2) / n_demo:.4f}"
      f" ± {mc_se * (n_demo - 2) / n_demo:.4f}   (이론 {9 * (n_demo - 2) / n_demo:.2f})")
print(f"  sd(s^2)          = {s2_list.std(ddof=1):.4f}"
      f"   (이론 sigma^2*sqrt(2/(n-2)) = {9 * np.sqrt(2 / (n_demo - 2)):.4f})")

RSS = 951.4547
  n-2 로 나눈 s^2 = 9.708721   (참 sigma^2 = 9)
  n   로 나눈 값  = 9.514547
  이론 치우침 -2*sigma^2/n = -0.1800
  s^2 의 상대 표준오차 sqrt(2/(n-2)) = 0.1429
  이 표본의 어긋남 +0.0787 = +0.551 표준오차

모의실험 2000 회
  mean(s^2)        = 8.9534 ± 0.0286   (참값 9)
  mean(RSS/n)      = 8.7744 ± 0.0280   (이론 8.82)
  sd(s^2)          = 1.2769   (이론 sigma^2*sqrt(2/(n-2)) = 1.2857)

def confidence_intervals(x, y_hat, beta_hat, x_bar, y_bar, n, s):
    """평균반응 E[y|x] 의 신뢰구간과 개별관측 y|x 의 예측구간을 구한다.

    두 구간의 차이는 근호 안의 1 하나뿐이다. 평균을 맞히는 일에는 추정의
    오차만 들어가지만, 개별 관측값을 맞히려면 그 위에 잡음 자체의 분산이
    더 얹히기 때문이다. 그래서 예측구간이 언제나 더 넓다.

    돌려주는 값
    ----------
    x0 : 그림을 그릴 x 격자
    lower, upper : 평균반응 신뢰구간의 위아래
    lower2, upper2 : 예측구간의 위아래
    """
    x0 = np.linspace(x.min(), x.max(), 20)
    y0_hat = beta_hat * (x0 - x_bar) + y_bar
    t_val = stats.t(n - 2).ppf(0.975)

    # 평균반응의 신뢰구간. (x0 - x_bar)^2 항 때문에 x 의 평균에서 멀어질수록
    # 구간이 넓어진다. 그래서 띠가 가운데가 잘록한 모래시계 모양이 된다.
    margin = t_val * s * np.sqrt(
        (1 / n) + (x0 - x_bar) ** 2 / np.sum((x - x_bar) ** 2)
    )
    lower = y0_hat - margin
    upper = y0_hat + margin

    # 예측구간. 근호 안에 1 이 더 있다. 이 1 이 잡음의 분산 몫이며,
    # n 을 아무리 키워도 사라지지 않는다.
    margin2 = t_val * s * np.sqrt(
        1 + (1 / n) + (x0 - x_bar) ** 2 / np.sum((x - x_bar) ** 2)
    )
    lower2 = y0_hat - margin2
    upper2 = y0_hat + margin2

    return x0, lower, upper, lower2, upper2

# 지렛값 h 와 두 폭을 자리마다 재어 본다.
x0_demo, lo_demo, up_demo, lo2_demo, up2_demo = confidence_intervals(
    x_demo, y_hat_demo, beta_hat_demo, x_bar_demo, y_bar_demo, n_demo, s_demo
)
t_val_demo = stats.t(n_demo - 2).ppf(0.975)
print(f"t_{n_demo - 2}(0.975) = {t_val_demo:.6f},  s = {s_demo:.6f},  S_xx = {S_xx_demo:.4f}")
print()
print(f"{'x0':>10s}{'h':>12s}{'CI 폭':>11s}{'PI 폭':>11s}{'PI/CI':>10s}{'sqrt((1+h)/h)':>16s}")
for x0v, label in [(x_bar_demo, "x-bar"), (float(x_demo.min()), "최소"), (float(x_demo.max()), "최대")]:
    h = 1 / n_demo + (x0v - x_bar_demo) ** 2 / S_xx_demo
    w_ci = 2 * t_val_demo * s_demo * np.sqrt(h)
    w_pi = 2 * t_val_demo * s_demo * np.sqrt(1 + h)
    print(f"{x0v:10.4f}{h:12.6f}{w_ci:11.4f}{w_pi:11.4f}{w_pi / w_ci:10.4f}"
          f"{np.sqrt((1 + h) / h):16.4f}   ({label})")
print()
h_min = 1 / n_demo
print(f"x0 = x-bar 에서 h = 1/n = {h_min:.6f}")
print(f"  그때 PI/CI = sqrt(n+1) = {np.sqrt(n_demo + 1):.4f}")
print(f"  띠가 가장 좁은 자리가 x-bar 인가? "
      f"{np.argmin(up_demo - lo_demo) == np.argmin(np.abs(x0_demo - x_bar_demo))}")
print()
print(f"n -> 무한일 때:  CI 폭 -> 0,  PI 폭 -> 2*t*s = {2 * t_val_demo * s_demo:.4f}")
print(f"  참값으로는 2*z*sigma = {2 * stats.norm.ppf(0.975) * 3:.4f}")
print(f"  실제 x-bar 에서의 PI 폭 {2 * t_val_demo * s_demo * np.sqrt(1 + h_min):.4f}"
      f" = 극한값의 {np.sqrt(1 + h_min):.4f} 배")

t_98(0.975) = 1.984467,  s = 3.115882,  S_xx = 101.5827

        x0           h       CI 폭       PI 폭     PI/CI   sqrt((1+h)/h)
    0.0598    0.010000     1.2367    12.4284   10.0499         10.0499   (x-bar)
   -2.5530    0.077204     3.4362    12.8352    3.7353          3.7353   (최소)
    2.2698    0.058078     2.9803    12.7208    4.2683          4.2683   (최대)

x0 = x-bar 에서 h = 1/n = 0.010000
  그때 PI/CI = sqrt(n+1) = 10.0499
  띠가 가장 좁은 자리가 x-bar 인가? True

n -> 무한일 때:  CI 폭 -> 0,  PI 폭 -> 2*t*s = 12.3667
  참값으로는 2*z*sigma = 11.7598
  실제 x-bar 에서의 PI 폭 12.4284 = 극한값의 1.0050 배

def plot_intervals(x, y, y_hat, x0, lower, upper, lower2, upper2):
    """두 구간을 나란히 그린다. 오른쪽 띠가 훨씬 넓은 것이 요점이다."""
    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 4))

    # 왼쪽: 평균반응의 신뢰구간
    ax0.plot(x, y, 'o', alpha=0.5)
    ax0.plot(x, y_hat, '--b', label='Fitted line')
    ax0.plot(x0, upper, '--r', label='95% CI bounds')
    ax0.plot(x0, lower, '--r')
    ax0.set_title('95% Confidence Interval for $E[y]$')
    ax0.legend()

    # 오른쪽: 개별 관측의 예측구간
    ax1.plot(x, y, 'o', alpha=0.5)
    ax1.plot(x, y_hat, '--b', label='Fitted line')
    ax1.plot(x0, upper2, '--r', label='95% PI bounds')
    ax1.plot(x0, lower2, '--r')
    ax1.set_title('95% Prediction Interval for $y$')
    ax1.legend()

    plt.tight_layout()
    plt.show()

# 두 패널이 각자 y 축 범위를 잡는다. 그것이 비교를 어떻게 왜곡하는가.
span_left = max(float(y_demo.max()), up_demo.max()) - min(float(y_demo.min()), lo_demo.min())
span_right = max(float(y_demo.max()), up2_demo.max()) - min(float(y_demo.min()), lo2_demo.min())
print(f"왼쪽 패널에 담기는 세로 범위  = {span_left:.4f}")
print(f"오른쪽 패널에 담기는 세로 범위 = {span_right:.4f}")
print(f"  두 범위의 비 = {span_right / span_left:.4f}")
print()
w_ci = (up_demo - lo_demo).mean()
w_pi = (up2_demo - lo2_demo).mean()
print(f"띠의 평균 실제 폭:  CI {w_ci:.4f},  PI {w_pi:.4f},  비 {w_pi / w_ci:.4f}")
print(f"화면에서 보이는 비(각 패널 높이로 나눈 뒤):"
      f"  {(w_pi / span_right) / (w_ci / span_left):.4f}")
print(f"  곧 그림은 실제 {w_pi / w_ci:.2f} 배 차이를 "
      f"{(w_pi / span_right) / (w_ci / span_left):.2f} 배로 보이게 한다")
print()
print(f"x 가 정렬되어 있는가? {bool(np.all(np.diff(x_demo.ravel()) >= 0))}")
print(f"  적합값이 x 의 일차함수이므로 정렬하지 않아도 같은 직선 위를 왕복할 뿐이다")
print(f"  x0 격자는 linspace 로 만들었으므로 정렬되어 있다: "
      f"{bool(np.all(np.diff(x0_demo) >= 0))}")
print(f"  격자점 {len(x0_demo)}개 — 띠가 꺾은선으로 그려진다")

왼쪽 패널에 담기는 세로 범위  = 18.9203
오른쪽 패널에 담기는 세로 범위 = 24.0830
  두 범위의 비 = 1.2729

띠의 평균 실제 폭:  CI 2.0895,  PI 12.5594,  비 6.0106
화면에서 보이는 비(각 패널 높이로 나눈 뒤):  4.7221
  곧 그림은 실제 6.01 배 차이를 4.72 배로 보이게 한다

x 가 정렬되어 있는가? False
  적합값이 x 의 일차함수이므로 정렬하지 않아도 같은 직선 위를 왕복할 뿐이다
  x0 격자는 linspace 로 만들었으므로 정렬되어 있다: True
  격자점 20개 — 띠가 꺾은선으로 그려진다

n = 100
sigma = 3

# 참 모형은 y = 1 + 2x + 잡음 이다. 참 sigma 를 알고 있으므로 아래에서
# 추정값 s^2 이 sigma^2 = 9 근처로 나오는지 확인할 수 있다.
x, y = generate_data(n, sigma)

y_hat, beta_hat, y_bar, x_bar = estimate_regression_line(x, y)

s_square, s = calculate_residual_variance(y, y_hat, n)
print(f"True σ²: {sigma**2}")
print(f"Estimated s²: {s_square:.4f}")

# 구간 계산
x0, lower, upper, lower2, upper2 = confidence_intervals(
    x, y_hat, beta_hat, x_bar, y_bar, n, s
)

# 그림으로 확인
plot_intervals(x, y, y_hat, x0, lower, upper, lower2, upper2)

True σ²: 9
Estimated s²: 9.7087

# 이 표본의 s^2 이 chi-square 표집분포에서 어디에 놓이는가
dfree = n - 2
q = stats.chi2(dfree).cdf(dfree * s_square / sigma ** 2)
print(f"s^2 = {s_square:.4f},  참 sigma^2 = {sigma ** 2}")
print(f"  (n-2)s^2/sigma^2 = {dfree * s_square / sigma ** 2:.4f}  ~ chi2_{dfree}")
print(f"  그 분포의 {q:.4f} 분위수  ->  양쪽 꼬리로 보면 흔한 값")
print(f"  sigma^2 의 95% 신뢰구간 = "
      f"({dfree * s_square / stats.chi2(dfree).ppf(0.975):.4f}, "
      f"{dfree * s_square / stats.chi2(dfree).ppf(0.025):.4f})   9 를 담는가? "
      f"{dfree * s_square / stats.chi2(dfree).ppf(0.975) <= 9 <= dfree * s_square / stats.chi2(dfree).ppf(0.025)}")
print()
# 두 구간의 포함률을 모의실험으로 확인한다 (2000 회, 몬테카를로 오차를 함께 적는다)
reps = 2000
rng = np.random.default_rng(7)
x0_fixed = 1.0                      # 이 자리에서 두 구간을 재어 본다
mean_true = 1 + 2 * x0_fixed        # 참 평균반응
hit_ci = np.zeros(reps, dtype=bool)
hit_pi = np.zeros(reps, dtype=bool)
for b in range(reps):
    xb_ = rng.standard_normal((n, 1))
    yb_ = 1 + 2 * xb_ + sigma * rng.standard_normal((n, 1))
    yh_, bh_, yb_bar, xb_bar = estimate_regression_line(xb_, yb_)
    s2_, s_ = calculate_residual_variance(yb_, yh_, n)
    Sxx_ = ((xb_ - xb_bar) ** 2).sum()
    h_ = 1 / n + (x0_fixed - xb_bar) ** 2 / Sxx_
    t_ = stats.t(n - 2).ppf(0.975)
    center = bh_ * (x0_fixed - xb_bar) + yb_bar
    hit_ci[b] = abs(center - mean_true) <= t_ * s_ * np.sqrt(h_)
    y_new = mean_true + sigma * rng.standard_normal()      # 새 관측값 하나
    hit_pi[b] = abs(center - y_new) <= t_ * s_ * np.sqrt(1 + h_)
for name, hits in [("평균반응 신뢰구간", hit_ci), ("개별관측 예측구간", hit_pi)]:
    p_hat = hits.mean()
    mc = np.sqrt(p_hat * (1 - p_hat) / reps)
    print(f"{name}  포함률 {p_hat:.4f} ± {mc:.4f}   (목표 0.95,"
          f"  z = {(p_hat - 0.95) / mc:+.2f})")

s^2 = 9.7087,  참 sigma^2 = 9
  (n-2)s^2/sigma^2 = 105.7172  ~ chi2_98
  그 분포의 0.7206 분위수  ->  양쪽 꼬리로 보면 흔한 값
  sigma^2 의 95% 신뢰구간 = (7.4752, 13.1233)   9 를 담는가? True

평균반응 신뢰구간  포함률 0.9490 ± 0.0049   (목표 0.95,  z = -0.20)
개별관측 예측구간  포함률 0.9495 ± 0.0049   (목표 0.95,  z = -0.10)