회귀 (대마초 가격과 인구통계)¶

high_q = (250 + 0.0009 * (income - 52_000) * 1000 / 1000
          - 1.2e-6 * pop + 18 * (pct_white - 0.72)
          + np.random.normal(0, 12, n_states))

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
import statsmodels.formula.api as smf

np.random.seed(42)


def build_dataset(n_states=50):
    """주 단위 인공자료를 만든다.

    실제 자료로 바꾸어도 아래 절차는 그대로다.
    HighQ(고품질 대마초 가격)를 인구, 1인당 소득, 백인 비율로 설명한다.
    """
    states = [f"state_{i:02d}" for i in range(n_states)]
    pop = np.random.lognormal(mean=15.0, sigma=0.9, size=n_states)
    income = np.random.normal(52_000, 9_000, n_states)
    pct_white = np.clip(np.random.normal(0.72, 0.13, n_states), 0.2, 0.95)
    pct_black = np.clip(np.random.normal(0.12, 0.08, n_states), 0.01, 0.40)
    pct_hispanic = np.clip(1 - pct_white - pct_black, 0.01, None)
    # 참 관계: 소득이 높을수록 비싸고, 인구가 많을수록 미세하게 싸다.
    high_q = (250 + 0.0009 * (income - 52_000) * 1000 / 1000
              - 1.2e-6 * pop + 18 * (pct_white - 0.72)
              + np.random.normal(0, 12, n_states))
    return pd.DataFrame({"state": states, "total_population": pop,
                         "per_capita_income": income,
                         "percent_white": pct_white,
                         "percent_black": pct_black,
                         "percent_hispanic": pct_hispanic,
                         "HighQ": high_q})


df = build_dataset()

# 검정용으로 떼어 둘 주를 지정한다. 실제 이름 대신 인덱스로 고른다.
TEST_IDX = {3, 8, 15, 22, 29, 34, 41, 47}
df["state"] = [f"state_{i:02d}" for i in range(len(df))]

test_states = {f"state_{i:02d}" for i in TEST_IDX}
train = df[~df['state'].isin(test_states)].copy()
test = df[df['state'].isin(test_states)].copy()

print(f"훈련 {len(train)}개 주, 검정 {len(test)}개 주")
print(train[["total_population", "per_capita_income", "percent_white", "HighQ"]]
      .describe().round(2).to_string())

훈련 42개 주, 검정 8개 주
       total_population  per_capita_income  percent_white   HighQ
count             42.00              42.00          42.00   42.00
mean         3502514.62           52041.15           0.71  246.38
std          3509386.94            7639.03           0.13   13.55
min           560338.80           28422.29           0.47  222.09
25%          1274815.60           47535.00           0.61  237.75
50%          2570441.13           51698.97           0.72  246.33
75%          4269259.39           56278.24           0.79  252.53
max         17314422.21           66081.79           0.95  297.57

# build_dataset 의 식에 적힌 참 모수를 꺼내 적는다.
TRUE_INTERCEPT = 250 - 0.0009 * 52_000 - 18 * 0.72
print(f"참 절편      = 250 - 0.0009*52000 - 18*0.72 = {TRUE_INTERCEPT:.4f}")
print(f"참 b_pop     = -1.2e-06")
print(f"참 b_income  =  0.0009")
print(f"참 b_white   =  18")
print(f"참 sigma     =  12")
print()
pop = train["total_population"]
print(f"인구  최솟값 {pop.min():,.0f}  최댓값 {pop.max():,.0f}  최댓값/최솟값 = {pop.max() / pop.min():.1f}")
print(f"      평균 {pop.mean():,.0f}  중앙값 {pop.median():,.0f}  (평균이 중앙값보다 크다)")
print(f"      왜도 {pop.skew():.4f}")
print(f"log(인구) 왜도 {np.log(pop).skew():.4f}")

참 절편      = 250 - 0.0009*52000 - 18*0.72 = 190.2400
참 b_pop     = -1.2e-06
참 b_income  =  0.0009
참 b_white   =  18
참 sigma     =  12

인구  최솟값 560,339  최댓값 17,314,422  최댓값/최솟값 = 30.9
      평균 3,502,515  중앙값 2,570,441  (평균이 중앙값보다 크다)
      왜도 2.4519
log(인구) 왜도 0.2292

# 먼저 인구 하나만 넣은 단변량 모형으로 기준선을 잡는다.
model1 = LinearRegression().fit(train[['total_population']], train['HighQ'])
pred1 = model1.predict(test[['total_population']])
rmse1 = np.sqrt(np.mean((test['HighQ'] - pred1) ** 2))
print(f"단변량 RMSE = {rmse1:.2f}")

단변량 RMSE = 12.66

# 참 모수와 훈련자료의 표준편차로 이론 예측오차를 계산한다.
sd_pop = train["total_population"].std()
sd_income = train["per_capita_income"].std()
sd_white = train["percent_white"].std()
print(f"훈련자료 표준편차:  인구 {sd_pop:,.0f}   소득 {sd_income:,.0f}   백인비율 {sd_white:.4f}")
print()
print("각 설명변수가 1 표준편차 오를 때 HighQ 의 참 변화량")
for name, beta, sd in [("인구", -1.2e-6, sd_pop), ("소득", 0.0009, sd_income),
                       ("백인비율", 18.0, sd_white)]:
    print(f"  {name:8s} {beta:+.2e} x {sd:12,.4f} = {beta * sd:+8.4f}")
print()
theory_total = np.sqrt((0.0009 * sd_income) ** 2 + (18 * sd_white) ** 2
                       + (1.2e-6 * sd_pop) ** 2 + 12.0 ** 2)
theory_uni = np.sqrt(theory_total ** 2 - (1.2e-6 * sd_pop) ** 2)
print(f"이론 sd(HighQ)                = {theory_total:.4f}   (관측 {train.HighQ.std():.4f})")
print(f"인구만 쓴 모형의 이론 예측오차 = {theory_uni:.4f}")
print(f"세 변수를 다 쓴 모형의 이론값  = 12.0000  (= 참 sigma)")
print()
print(f"실제 검정 RMSE                = {rmse1:.4f}")
print(f"검정자료가 {len(test)}개이므로 RMSE 의 상대 표준오차 ~ 1/sqrt(2n) = {1 / np.sqrt(2 * len(test)):.4f}")
print(f"이론값에서의 상대 차이        = {rmse1 / theory_uni - 1:+.4f}")

훈련자료 표준편차:  인구 3,509,387   소득 7,639   백인비율 0.1252

각 설명변수가 1 표준편차 오를 때 HighQ 의 참 변화량
  인구       -1.20e-06 x 3,509,386.9355 =  -4.2113
  소득       +9.00e-04 x   7,639.0328 =  +6.8751
  백인비율     +1.80e+01 x       0.1252 =  +2.2541

이론 sd(HighQ)                = 14.6316   (관측 13.5481)
인구만 쓴 모형의 이론 예측오차 = 14.0124
세 변수를 다 쓴 모형의 이론값  = 12.0000  (= 참 sigma)

실제 검정 RMSE                = 12.6631
검정자료가 8개이므로 RMSE 의 상대 표준오차 ~ 1/sqrt(2n) = 0.2500
이론값에서의 상대 차이        = -0.0963

# 변수를 셋으로 늘려 RMSE 가 실제로 줄어드는지 본다.
formula = "HighQ ~ total_population + per_capita_income + percent_white"
sm_model = smf.ols(formula=formula, data=train).fit()
# summary()는 실행 날짜와 시각을 함께 찍으므로 계수 표만 인쇄한다.
print(sm_model.summary().tables[1])

pred3 = sm_model.predict(test)
rmse3 = np.sqrt(np.mean((test['HighQ'] - pred3) ** 2))
print(f"다변량 RMSE = {rmse3:.2f}")

=====================================================================================
                        coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------------
Intercept           187.9016     22.531      8.340      0.000     142.290     233.513
total_population  -2.889e-08   5.85e-07     -0.049      0.961   -1.21e-06    1.15e-06
per_capita_income     0.0007      0.000      2.409      0.021       0.000       0.001
percent_white        32.0081     17.538      1.825      0.076      -3.496      67.512
=====================================================================================
다변량 RMSE = 10.32

# 참 모수가 95% 신뢰구간에 들어오는지 확인한다.
truth = {"Intercept": 190.24, "total_population": -1.2e-6,
         "per_capita_income": 0.0009, "percent_white": 18.0}
ci = sm_model.conf_int()
print(f"{'항':20s}{'추정값':>14s}{'95% CI 하한':>15s}{'95% CI 상한':>15s}{'참값':>14s}  판정")
for name, true_value in truth.items():
    low, high = ci.loc[name]
    verdict = "포함" if low <= true_value <= high else "벗어남"
    print(f"{name:20s}{sm_model.params[name]:14.6g}{low:15.6g}{high:15.6g}"
          f"{true_value:14.6g}  {verdict}")
print()
print(f"잔차 표준오차 s = {np.sqrt(sm_model.mse_resid):.4f}   (참 sigma = 12)")
print(f"훈련 R^2        = {sm_model.rsquared:.4f}")
print(f"  이론 R^2 = 1 - 12^2/14.6316^2 = {1 - 144 / 14.6316 ** 2:.4f}")
print()
print(f"단변량 검정 RMSE = {rmse1:.4f}   다변량 검정 RMSE = {rmse3:.4f}"
      f"   비 = {rmse3 / rmse1:.4f}")

항                              추정값      95% CI 하한      95% CI 상한            참값  판정
Intercept                  187.902         142.29        233.513        190.24  포함
total_population      -2.88857e-08   -1.21254e-06    1.15477e-06      -1.2e-06  포함
per_capita_income      0.000688352    0.000109966     0.00126674        0.0009  포함
percent_white              32.0081        -3.4962        67.5125            18  포함

잔차 표준오차 s = 12.9413   (참 sigma = 12)
훈련 R^2        = 0.1543
  이론 R^2 = 1 - 12^2/14.6316^2 = 0.3274

단변량 검정 RMSE = 12.6631   다변량 검정 RMSE = 10.3194   비 = 0.8149

# 주별로 실제값과 예측값을 나란히 놓아 어느 주에서 크게 빗나갔는지 본다.
# 평균 오차 하나만 보아서는 놓치는 것이 있다.
result = pd.DataFrame({
    'state': test['state'].values,
    'actual': test['HighQ'].values,
    'predicted': np.round(pred3.values, 2),
})
result['error'] = result['actual'] - result['predicted']
print(result.round(2).to_string(index=False))

   state  actual  predicted  error
state_03  250.27     246.82   3.45
state_08  256.90     249.80   7.10
state_15  268.43     256.34  12.09
state_22  257.06     252.26   4.80
state_29  237.86     232.26   5.60
state_34  258.80     237.71  21.09
state_41  254.95     258.07  -3.12
state_47  231.03     242.62 -11.59

# 여덟 오차로 RMSE 를 손으로 계산하고, 치우침과 극단값의 몫을 본다.
errors = result['error'].values
print(f"오차 = {np.round(errors, 2)}")
print(f"제곱합 = {np.sum(errors ** 2):.4f},  나누기 {len(errors)} = {np.mean(errors ** 2):.4f}")
print(f"제곱근 = {np.sqrt(np.mean(errors ** 2)):.4f}   (보기 3 의 RMSE = {rmse3:.4f})")
print()
print(f"평균 오차(치우침) = {errors.mean():+.4f}")
print(f"오차의 표준편차   = {errors.std(ddof=1):.4f}")
print(f"평균 오차의 표준오차 = {errors.std(ddof=1) / np.sqrt(len(errors)):.4f}"
      f"   t = {errors.mean() / (errors.std(ddof=1) / np.sqrt(len(errors))):.4f}")
print(f"RMSE^2 = 치우침^2 + (n-1)/n * 표준편차^2 = "
      f"{errors.mean() ** 2:.4f} + {(len(errors) - 1) / len(errors) * errors.var(ddof=1):.4f}"
      f" = {errors.mean() ** 2 + (len(errors) - 1) / len(errors) * errors.var(ddof=1):.4f}")
print()
worst = np.argmax(np.abs(errors))
print(f"가장 큰 오차: {result['state'].iloc[worst]}  오차 {errors[worst]:+.2f}")
print(f"  제곱합에서 차지하는 몫 = {errors[worst] ** 2 / np.sum(errors ** 2):.4f}")
print(f"  이 주를 빼면 RMSE = {np.sqrt(np.mean(np.delete(errors, worst) ** 2)):.4f}")
print(f"오차가 양수인 주 {np.sum(errors > 0)}개, 음수인 주 {np.sum(errors < 0)}개")

오차 = [  3.45   7.1   12.09   4.8    5.6   21.09  -3.12 -11.59]
제곱합 = 851.9593,  나누기 8 = 106.4949
제곱근 = 10.3196   (보기 3 의 RMSE = 10.3194)

평균 오차(치우침) = +4.9279
오차의 표준편차   = 9.6930
평균 오차의 표준오차 = 3.4270   t = 1.4380
RMSE^2 = 치우침^2 + (n-1)/n * 표준편차^2 = 24.2847 + 82.2102 = 106.4949

가장 큰 오차: state_34  오차 +21.09
  제곱합에서 차지하는 몫 = 0.5222
  이 주를 빼면 RMSE = 7.6257
오차가 양수인 주 6개, 음수인 주 2개

features = ['total_population', 'per_capita_income',
            'percent_white', 'percent_black', 'percent_hispanic']
corr = train[['HighQ'] + features].corr()
print(corr['HighQ'].sort_values(ascending=False))

HighQ                1.000000
per_capita_income    0.280779
percent_white        0.154687
total_population     0.054487
percent_hispanic    -0.072988
percent_black       -0.087616
Name: HighQ, dtype: float64

formula_full = ("HighQ ~ total_population + per_capita_income + "
                "percent_white + percent_black + percent_hispanic")
sm_full = smf.ols(formula=formula_full, data=train).fit()
pred_full = sm_full.predict(test)
rmse_full = np.sqrt(np.mean((test['HighQ'] - pred_full) ** 2))

from sklearn.model_selection import train_test_split

X_all = df[['total_population', 'per_capita_income', 'percent_white']]
y_all = df['HighQ']
X_tr, X_te, y_tr, y_te = train_test_split(X_all, y_all,
                                          test_size=0.2, random_state=42)
model_split = LinearRegression().fit(X_tr, y_tr)
rmse_split = np.sqrt(np.mean((y_te - model_split.predict(X_te)) ** 2))

sm_reduced = smf.ols("HighQ ~ total_population", data=train).fit()
f_stat = ((sm_reduced.ssr - sm_model.ssr) / 2) / sm_model.mse_resid
from scipy import stats
p_val = 1 - stats.f.cdf(f_stat, 2, sm_model.df_resid)