선형회귀의 정규성 확인¶

import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(7)
n = 120

# X는 균등, Y는 X에 선형으로 의존하되 오차의 분산이 X와 함께 커진다.
# 이렇게 두면 선형성은 성립하고 등분산성만 깨져, 각 진단이 무엇을
# 잡아내고 무엇을 놓치는지 구분해 볼 수 있다.
X = rng.uniform(0, 10, n)
Y = 2.0 + 1.5 * X + rng.normal(0, 0.5 + 0.35 * X, n)

df = pd.DataFrame({"X": X, "Y": Y})
model = sm.OLS(Y, sm.add_constant(X)).fit()
residuals = model.resid
fitted = model.fittedvalues

print(f"beta_hat = {model.params.round(4)}")
print(f"R^2 = {model.rsquared:.4f}")

beta_hat = [1.5933 1.5317]
R^2 = 0.7813

# c_i = (x_i - xbar) / S_xx 가 가중치다. sum c_i = 0, sum c_i x_i = 1 을 먼저 본다.
Sxx = ((X - X.mean()) ** 2).sum()
c = (X - X.mean()) / Sxx
print(f"S_xx        = {Sxx:.4f}")
print(f"sum c_i     = {c.sum():+.2e}   (0 이어야 한다)")
print(f"sum c_i x_i = {(c * X).sum():.10f}   (1 이어야 한다)")

# 이 자료는 sigma_i 를 우리가 안다. 참 표준오차를 직접 계산할 수 있다.
sigma_i = 0.5 + 0.35 * X
se_true = np.sqrt((c ** 2 * sigma_i ** 2).sum())
se_ols = model.bse[1]
print(f"\n참 SE(b1)   = sqrt(sum c_i^2 sigma_i^2) = {se_true:.6f}")
print(f"OLS 보고 SE = s/sqrt(S_xx)              = {se_ols:.6f}")
print(f"보고값이 참값의 {se_ols / se_true:.1%}")

# 같은 X 를 고정한 채 오차만 다시 뽑아 흔들림을 직접 재 본다.
sim = np.random.default_rng(123)
b1 = 1.5 + (sim.normal(0, 1, (20000, n)) * sigma_i * c).sum(axis=1)
print(f"\n모의 20000회: 평균 {b1.mean():.4f}, 표준편차 {b1.std(ddof=1):.6f}")

S_xx        = 1004.3944
sum c_i     = -2.78e-17   (0 이어야 한다)
sum c_i x_i = 1.0000000000   (1 이어야 한다)

참 SE(b1)   = sqrt(sum c_i^2 sigma_i^2) = 0.082525
OLS 보고 SE = s/sqrt(S_xx)              = 0.074595
보고값이 참값의 90.4%

모의 20000회: 평균 1.5009, 표준편차 0.082505

import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm, skew, kurtosis

# 정규성은 자료가 아니라 잔차에 요구되는 가정이다. 게다가 계수 추정의
# 불편성에는 필요 없고, 작은 표본에서 t 검정과 신뢰구간을 쓰기 위해 필요하다.
residuals = model.resid

fig, ax = plt.subplots(figsize=(8, 5))
ax.hist(residuals, bins=30, density=True, edgecolor='black', alpha=0.7, label='Residuals')

# 잔차의 평균과 표준편차로 만든 정규곡선을 겹쳐 그린다.
x = np.linspace(residuals.min(), residuals.max(), 100)
ax.plot(x, norm.pdf(x, loc=residuals.mean(), scale=residuals.std()),
        'r-', linewidth=2, label='Normal PDF')

ax.set_xlabel('Residuals')
ax.set_ylabel('Density')
ax.set_title('Histogram of Residuals with Normal Overlay')
ax.legend()
plt.show()

# 그림이 말해 주지 않는 수치들
print(f"잔차 평균   = {residuals.mean():+.3e}   (정규방정식이 강제한다)")
print(f"왜도        = {skew(residuals):+.4f}")
print(f"첨도        = {kurtosis(residuals, fisher=False):.4f}")
print(f"겹친 곡선의 척도 residuals.std() = {residuals.std():.4f}")
print(f"회귀의      sigma_hat            = {np.sqrt(model.mse_resid):.4f}")

# 계급 수를 바꾸면 가장 높은 막대의 높이가 얼마나 달라지는가
for b in (10, 30, 60):
    h, _ = np.histogram(residuals, bins=b, density=True)
    print(f"bins={b:>3}: 최고 막대 밀도 {h.max():.4f}")

# 이 그림이 전혀 보지 못하는 것: 오차의 분산이 x 와 함께 커진다는 사실
lo, hi = residuals[X < 5], residuals[X >= 5]
print(f"\nX<5  쪽 잔차 표준편차 = {lo.std(ddof=1):.4f}  (n={lo.size})")
print(f"X>=5 쪽 잔차 표준편차 = {hi.std(ddof=1):.4f}  (n={hi.size})")

잔차 평균   = +3.523e-15   (정규방정식이 강제한다)
왜도        = -0.0508
첨도        = 3.5595
겹친 곡선의 척도 residuals.std() = 2.3443
회귀의      sigma_hat            = 2.3641
bins= 10: 최고 막대 밀도 0.2459
bins= 30: 최고 막대 밀도 0.3136
bins= 60: 최고 막대 밀도 0.3320

X<5  쪽 잔차 표준편차 = 1.3132  (n=57)
X>=5 쪽 잔차 표준편차 = 3.0126  (n=63)

import scipy.stats as stats
import matplotlib.pyplot as plt

residuals = model.resid

# Q-Q 그림이 히스토그램보다 낫다. 계급 수에 휘둘리지 않고, 정규에서
# 벗어나는 자리가 꼬리인지 가운데인지까지 알려 준다.
fig, ax = plt.subplots(figsize=(6, 6))
stats.probplot(residuals, dist="norm", plot=ax)
ax.set_title('Q-Q Plot of Residuals')
plt.show()

# probplot 은 그림만 그리는 것이 아니라 기준선의 계수도 돌려준다.
(osm, osr), (slope, intercept, r) = stats.probplot(residuals, dist="norm")
print(f"기준선 기울기 = {slope:.4f}   (sigma 의 추정값)")
print(f"기준선 절편   = {intercept:+.2e}   (0 이어야 한다)")
print(f"상관계수  r   = {r:.6f},   r^2 = {r ** 2:.6f}")

# 양 끝 두 점이 얼마나 바깥에 있는가
print(f"\n가장 작은 잔차 {osr[0]:+.3f}  (기준선 예측 {slope * osm[0] + intercept:+.3f})")
print(f"가장 큰   잔차 {osr[-1]:+.3f}  (기준선 예측 {slope * osm[-1] + intercept:+.3f})")

기준선 기울기 = 2.3627   (sigma 의 추정값)
기준선 절편   = +3.10e-15   (0 이어야 한다)
상관계수  r   = 0.990542,   r^2 = 0.981173

가장 작은 잔차 -7.036  (기준선 예측 -5.970)
가장 큰   잔차 +6.518  (기준선 예측 +5.970)

import numpy as np
from scipy.stats import shapiro, norm

residuals = model.resid

# Shapiro-Wilk 는 작은 표본에서 검정력이 좋다. 다만 표본이 수천을 넘으면
# 실질적으로 무시할 만한 이탈에도 유의하게 나오므로 그림과 함께 읽는다.
stat, p_value = shapiro(residuals)
print(f'Shapiro-Wilk statistic: {stat:.4f}')
print(f'Shapiro-Wilk p-value: {p_value:.4f}')

if p_value > 0.05:
    print('No significant evidence of non-normality (fail to reject H0)')
else:
    print('Significant evidence of non-normality (reject H0)')

# (1) 의 등호 조건을 직접 만들어 본다. 정렬된 표본을 정규점수 그 자체로 두면
# x_(i) - xbar 가 a_i 에 비례하는 상황에 가장 가까워진다.
i = np.arange(1, 121)
perfect = norm.ppf((i - 0.375) / (120 + 0.25))
print(f'\n정규점수 격자 자체의 W = {shapiro(perfect).statistic:.6f}')

# 반대쪽 극단: 치우친 분포에서는 W 가 1 에서 멀어진다.
expo = np.sort(np.random.default_rng(0).exponential(size=120))
print(f'지수분포 표본의    W = {shapiro(expo).statistic:.6f}, '
      f'p = {shapiro(expo).pvalue:.2e}')
print(f'실제 잔차의        W = {stat:.6f}')

Shapiro-Wilk statistic: 0.9836
Shapiro-Wilk p-value: 0.1524
No significant evidence of non-normality (fail to reject H0)

정규점수 격자 자체의 W = 0.999275
지수분포 표본의    W = 0.842123, p = 5.45e-10
실제 잔차의        W = 0.983570

from scipy.stats import anderson

residuals = model.resid

# Anderson-Darling 은 꼬리 쪽 이탈에 특히 민감하다. p-값 대신 유의수준별
# 기각값을 돌려주므로, 통계량과 기각값을 견주어 판단한다.
result = anderson(residuals, dist='norm')
print(f'Anderson-Darling statistic: {result.statistic:.4f}')
print()
for i in range(len(result.critical_values)):
    sig_level = result.significance_level[i]
    crit_value = result.critical_values[i]
    status = 'REJECT' if result.statistic > crit_value else 'Fail to reject'
    print(f'At {sig_level}% significance: Critical value = {crit_value:.4f} → {status}')

Anderson-Darling statistic: 0.9338

At 15.0% significance: Critical value = 0.5580 → REJECT
At 10.0% significance: Critical value = 0.6360 → REJECT
At 5.0% significance: Critical value = 0.7630 → REJECT
At 2.5% significance: Critical value = 0.8900 → REJECT
At 1.0% significance: Critical value = 1.0590 → Fail to reject

import numpy as np
from scipy.stats import norm, anderson, shapiro

residuals = model.resid
n = len(residuals)
xs = np.sort(residuals)

# 척도를 ddof=1 로 잡아야 scipy 와 맞는다. ddof=0 을 쓰면 값이 달라진다.
for ddof in (0, 1):
    F = norm.cdf((xs - xs.mean()) / xs.std(ddof=ddof))
    i = np.arange(1, n + 1)
    A2 = -n - np.sum((2 * i - 1) / n * (np.log(F) + np.log(1 - F[::-1])))
    print(f"정의대로 계산 (ddof={ddof}): A^2 = {A2:.12f}")
print(f"scipy.stats.anderson       : A^2 = {anderson(residuals).statistic:.12f}")

# 왜 꼬리에 민감한가: 가장 바깥 두 점이 기대보다 얼마나 더 바깥인가
F = norm.cdf((xs - xs.mean()) / xs.std(ddof=1))
print(f"\n가장 작은 잔차의 F     = {F[0]:.5f}  (기대 1/(n+1) = {1 / (n + 1):.5f})")
print(f"가장 큰   잔차의 1 - F = {1 - F[-1]:.5f}")
print(f"log F 의 끝항          = {np.log(F[0]):.3f}")

# 양 끝 두 개씩을 덜어 내면 두 통계량이 어떻게 움직이는가
trim = xs[2:-2]
print(f"\n{'':14}{'전체':>12}{'양끝 2개씩 제거':>18}")
print(f"{'A^2':14}{anderson(residuals).statistic:>12.4f}{anderson(trim).statistic:>18.4f}")
print(f"{'Shapiro-Wilk W':14}{shapiro(residuals).statistic:>12.4f}{shapiro(trim).statistic:>18.4f}")

정의대로 계산 (ddof=0): A^2 = 0.923067249972
정의대로 계산 (ddof=1): A^2 = 0.933839718132
scipy.stats.anderson       : A^2 = 0.933839718132

가장 작은 잔차의 F     = 0.00140  (기대 1/(n+1) = 0.00826)
가장 큰   잔차의 1 - F = 0.00281
log F 의 끝항          = -6.571

                        전체         양끝 2개씩 제거
A^2                 0.9338            0.7793
Shapiro-Wilk W      0.9836            0.9821

import numpy as np
from scipy.stats import chi2
from statsmodels.stats.stattools import jarque_bera

residuals = model.resid

# Jarque-Bera 는 왜도와 첨도만 본다. 정규분포는 왜도 0, 첨도 3 이므로
# 이 둘이 얼마나 벗어났는지를 하나의 통계량으로 묶은 것이다.
# 큰 표본에서만 믿을 만하다.
jb_stat, jb_pvalue, skew, kurtosis = jarque_bera(residuals)
print(f'Jarque-Bera statistic: {jb_stat:.4f}')
print(f'Jarque-Bera p-value: {jb_pvalue:.4f}')
print(f'Skewness: {skew:.4f}')
print(f'Kurtosis: {kurtosis:.4f}')

# 손으로 다시 계산한다.
n = len(residuals)
skew_part = n / 6 * skew ** 2
kurt_part = n / 6 * (kurtosis - 3) ** 2 / 4
print(f'\n왜도 몫  = (n/6) S^2        = {skew_part:.6f}   ({skew_part / jb_stat:5.1%})')
print(f'첨도 몫  = (n/6)(K-3)^2 / 4 = {kurt_part:.6f}   ({kurt_part / jb_stat:5.1%})')
print(f'합       = {skew_part + kurt_part:.6f}   (statsmodels: {jb_stat:.6f})')

# 자유도 2 의 카이제곱은 생존함수가 exp(-x/2) 로 닫힌 꼴이다.
print(f'\np = exp(-JB/2)      = {np.exp(-jb_stat / 2):.6f}')
print(f'p = chi2(2).sf(JB)  = {chi2.sf(jb_stat, 2):.6f}')

Jarque-Bera statistic: 1.6167
Jarque-Bera p-value: 0.4456
Skewness: -0.0508
Kurtosis: 3.5595

왜도 몫  = (n/6) S^2        = 0.051610   ( 3.2%)
첨도 몫  = (n/6)(K-3)^2 / 4 = 1.565108   (96.8%)
합       = 1.616719   (statsmodels: 1.616719)

p = exp(-JB/2)      = 0.445589
p = chi2(2).sf(JB)  = 0.445589