신뢰띠와 예측띠¶

import numpy as np

def generate_data(n, sigma, seed=0):
    """자료를 만든다. 참 모형은 y = 1 + 2x + 잡음 이다."""
    np.random.seed(seed)
    x = np.random.randn(n, 1)
    y = 1 + 2 * x + sigma * np.random.randn(n, 1)
    return x, y

# 전역 씨앗을 쓰는 함수의 재현성을 확인한다.
x_a, y_a = generate_data(100, 3)
x_b, y_b = generate_data(100, 3)
print(f"같은 인수로 두 번 부르면 같은가:  x {np.array_equal(x_a, x_b)},  y {np.array_equal(y_a, y_b)}")

x_s, y_s = generate_data(50, 3)
print(f"n=50 의 x 가 n=100 의 앞 50 과 같은가:  {np.allclose(x_s, x_a[:50])}")
print(f"n=50 의 y 가 n=100 의 앞 50 과 같은가:  {np.allclose(y_s, y_a[:50])}")
print(f"  n=50 의 잡음은 난수열의 50번째부터, n=100 의 잡음은 100번째부터 뽑힌다")

np.random.seed(0)
first, second = np.random.randn(3), np.random.randn(3)
print(f"\n씨앗을 한 번 고정한 뒤 연속 두 호출")
print(f"  첫 번째 {np.round(first, 4)}")
print(f"  두 번째 {np.round(second, 4)}   같은가? {np.allclose(first, second)}")

같은 인수로 두 번 부르면 같은가:  x True,  y True
n=50 의 x 가 n=100 의 앞 50 과 같은가:  True
n=50 의 y 가 n=100 의 앞 50 과 같은가:  False
  n=50 의 잡음은 난수열의 50번째부터, n=100 의 잡음은 100번째부터 뽑힌다

씨앗을 한 번 고정한 뒤 연속 두 호출
  첫 번째 [1.7641 0.4002 0.9787]
  두 번째 [ 2.2409  1.8676 -0.9773]   같은가? False

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

    beta_hat = r * (s_y / s_x) 이다. 최소제곱해와 정확히 같은 값이지만,
    회귀계수가 상관계수를 척도만 바꿔 옮긴 것임이 이 꼴에서 드러난다.
    """
    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

# 상관행렬의 두 비대각 원소, 그리고 척도를 바꿀 때 무엇이 변하는가.
corr_matrix = np.corrcoef(np.concatenate([x_a, y_a], axis=1), rowvar=False)
print("상관행렬")
print(np.round(corr_matrix, 6))
print(f"  [1,0] 과 [0,1] 의 차이 = {corr_matrix[1, 0] - corr_matrix[0, 1]:.2e}")
print()
y_hat_a, beta_a, y_bar_a, x_bar_a = estimate_regression_line(x_a, y_a)
print(f"기울기 {beta_a:.8f},  절편 = y-bar - b*x-bar = {y_bar_a - beta_a * x_bar_a:.8f}")
print()
print("척도를 바꾸면 무엇이 변하는가")
for label, xx, yy in [("x -> 10x", 10 * x_a, y_a), ("y -> 10y", x_a, 10 * y_a)]:
    r_new = np.corrcoef(np.concatenate([xx, yy], axis=1), rowvar=False)[1, 0]
    b_new = r_new * yy.std(ddof=1) / xx.std(ddof=1)
    print(f"  {label:10s} r = {r_new:.8f}   beta = {b_new:.8f}")
print(f"  원래      r = {beta_a / (y_a.std(ddof=1) / x_a.std(ddof=1)):.8f}   beta = {beta_a:.8f}")

상관행렬
[[1.       0.608065]
 [0.608065 1.      ]]
  [1,0] 과 [0,1] 의 차이 = -1.11e-16

기울기 2.34409530,  절편 = y-bar - b*x-bar = 1.22545926

척도를 바꾸면 무엇이 변하는가
  x -> 10x   r = 0.60806485   beta = 0.23440953
  y -> 10y   r = 0.60806485   beta = 23.44095301
  원래      r = 0.60806485   beta = 2.34409530

def calculate_residual_variance(y, y_hat, n):
    """잔차분산 s^2 을 구한다.

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

# s^2 을 r^2 으로 적는 세 가지 길이 같은지 본다.
from scipy.special import gammaln

n_a = 100
s_square_a, s_a = calculate_residual_variance(y_a, y_hat_a, n_a)
r_a = np.corrcoef(np.concatenate([x_a, y_a], axis=1), rowvar=False)[1, 0]
s_y_a = y_a.std(ddof=1)
S_yy_a = ((y_a - y_bar_a) ** 2).sum()
print(f"RSS                       = {((y_a - y_hat_a) ** 2).sum():.6f}")
print(f"S_yy * (1 - r^2)          = {S_yy_a * (1 - r_a ** 2):.6f}")
print(f"s^2 = RSS/(n-2)           = {s_square_a:.6f}")
print(f"(n-1)/(n-2) * s_y^2 (1-r^2) = {(n_a - 1) / (n_a - 2) * s_y_a ** 2 * (1 - r_a ** 2):.6f}")
print()
# s 는 sigma 의 불편추정량이 아니다
def unbias_factor(nu):
    return np.sqrt(2 / nu) * np.exp(gammaln((nu + 1) / 2) - gammaln(nu / 2))
print("E[s]/sigma = sqrt(2/nu) * Gamma((nu+1)/2) / Gamma(nu/2)")
for nu in (8, 18, 98):
    print(f"  자유도 {nu:3d}:  {unbias_factor(nu):.6f}   치우침 {unbias_factor(nu) - 1:+.6f}"
          f"   근사 -1/(4nu) = {-1 / (4 * nu):+.6f}")
print()
print(f"이 표본:  s = {s_a:.6f},  참 sigma = 3")
print(f"  치우침 보정한 추정값 s/{unbias_factor(n_a - 2):.6f} = {s_a / unbias_factor(n_a - 2):.6f}")

RSS                       = 951.454669
S_yy * (1 - r^2)          = 951.454669
s^2 = RSS/(n-2)           = 9.708721
(n-1)/(n-2) * s_y^2 (1-r^2) = 9.708721

E[s]/sigma = sqrt(2/nu) * Gamma((nu+1)/2) / Gamma(nu/2)
  자유도   8:  0.969311   치우침 -0.030689   근사 -1/(4nu) = -0.031250
  자유도  18:  0.986214   치우침 -0.013786   근사 -1/(4nu) = -0.013889
  자유도  98:  0.997452   치우침 -0.002548   근사 -1/(4nu) = -0.002551

이 표본:  s = 3.115882,  참 sigma = 3
  치우침 보정한 추정값 s/0.997452 = 3.123841

from scipy import stats

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

    두 식의 차이는 근호 안의 1 뿐이다. 평균을 맞히는 데는 추정오차만 들지만,
    개별 관측을 맞히려면 잡음 자체의 분산이 더 얹힌다.
    """
    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)
    ss_x = np.sum((x - x_bar) ** 2)

    # 평균반응의 신뢰구간. x 의 평균에서 멀어질수록 넓어진다.
    margin = t_val * s * np.sqrt((1 / n) + (x0 - x_bar) ** 2 / ss_x)
    lower = y0_hat - margin
    upper = y0_hat + margin

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

    return x0, lower, upper, lower2, upper2

# 코드가 그리는 범위와, 그 밖으로 나갔을 때 띠가 얼마나 넓어지는가.
x0_a, lo_a, up_a, lo2_a, up2_a = confidence_intervals(
    x_a, y_hat_a, beta_a, x_bar_a, y_bar_a, n_a, s_a
)
s_x_a = x_a.std(ddof=1)
S_xx_a = ((x_a - x_bar_a) ** 2).sum()
t_a = stats.t(n_a - 2).ppf(0.975)
print(f"코드가 그리는 x 범위 = [{x0_a.min():.4f}, {x0_a.max():.4f}]"
      f"   = x-bar + [{(x0_a.min() - x_bar_a) / s_x_a:.3f}, "
      f"{(x0_a.max() - x_bar_a) / s_x_a:.3f}] * s_x")
print()
print(f"{'k':>4s}{'x0 = x-bar + k*s_x':>20s}{'h':>12s}{'CI 반폭':>10s}{'배수':>9s}"
      f"{'sqrt(1+k^2 n/(n-1))':>22s}{'PI 반폭':>10s}{'PI 배수':>9s}")
half_ci_center = t_a * s_a * np.sqrt(1 / n_a)
half_pi_center = t_a * s_a * np.sqrt(1 + 1 / n_a)
for k in (0, 1, 2, 3, 5, 10):
    x0v = x_bar_a + k * s_x_a
    h = 1 / n_a + (x0v - x_bar_a) ** 2 / S_xx_a
    hc = t_a * s_a * np.sqrt(h)
    hp = t_a * s_a * np.sqrt(1 + h)
    print(f"{k:4d}{x0v:20.4f}{h:12.6f}{hc:10.4f}{hc / half_ci_center:9.4f}"
          f"{np.sqrt(1 + k ** 2 * n_a / (n_a - 1)):22.4f}{hp:10.4f}{hp / half_pi_center:9.4f}")
print()
print(f"격자 끝까지(k = 2.18)만 그리면 CI 반폭은 "
      f"{t_a * s_a * np.sqrt(1 / n_a + (2.18 * s_x_a) ** 2 / S_xx_a):.4f} 에서 멈춘다")
print(f"  k = 10 까지 가면 {t_a * s_a * np.sqrt(1 / n_a + (10 * s_x_a) ** 2 / S_xx_a):.4f}"
      f" 로 {np.sqrt(1 + 100 * n_a / (n_a - 1)) / np.sqrt(1 + 2.18 ** 2 * n_a / (n_a - 1)):.2f} 배가 된다")

코드가 그리는 x 범위 = [-2.5530, 2.2698]   = x-bar + [-2.579, 2.182] * s_x

   k  x0 = x-bar + k*s_x           h     CI 반폭       배수   sqrt(1+k^2 n/(n-1))     PI 반폭    PI 배수
   0              0.0598    0.010000    0.6183   1.0000                1.0000    6.2142   1.0000
   1              1.0728    0.020101    0.8767   1.4178                1.4178    6.2452   1.0050
   2              2.0857    0.050404    1.3882   2.2451                2.2451    6.3373   1.0198
   3              3.0987    0.100909    1.9642   3.1766                3.1766    6.4878   1.0440
   5              5.1246    0.262525    3.1682   5.1237                5.1237    6.9478   1.1180
  10             10.1894    1.020101    6.2452  10.1000               10.1000    8.7884   1.4142

격자 끝까지(k = 2.18)만 그리면 CI 반폭은 1.4892 에서 멈춘다
  k = 10 까지 가면 6.2452 로 4.19 배가 된다

import numpy as np
from scipy import stats

n, s = 100, 3.0
t_90 = stats.t(n - 2).ppf(0.95)
t_95 = stats.t(n - 2).ppf(0.975)

margin_90 = t_90 * s * np.sqrt(1 / n)
margin_95 = t_95 * s * np.sqrt(1 / n)
print(f"t_90 = {t_90:.4f}, margin_90 = {margin_90:.4f}")
print(f"t_95 = {t_95:.4f}, margin_95 = {margin_95:.4f}")
print(f"ratio = {margin_90 / margin_95:.4f}")

t_90 = 1.6606, margin_90 = 0.4982
t_95 = 1.9845, margin_95 = 0.5953
ratio = 0.8368

x30, y30 = generate_data(30, 1)
y_hat30, beta30, ybar30, xbar30 = estimate_regression_line(x30, y30)
_, s30 = calculate_residual_variance(y30, y_hat30, 30)