Patsy를 이용한 스플라인¶

import pandas as pd

# "Practical Statistics for Data Scientists" 저장소의 자료. 탭으로 구분되어 있다.
url = ("https://raw.githubusercontent.com/gedeck/"
       "practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/house_sales.csv")
house = pd.read_csv(url, sep='\t')

print(f"{len(house)}건, 열 {house.shape[1]}개")
print(house[['AdjSalePrice', 'SqFtTotLiving', 'YrBuilt']].describe().round(1).to_string())

22687건, 열 22개
       AdjSalePrice  SqFtTotLiving  YrBuilt
count       22687.0        22687.0  22687.0
mean       565233.3         2080.2   1971.2
std        385402.9          913.7     30.3
min          3368.0          370.0   1900.0
25%        360563.0         1420.0   1950.0
50%        471315.0         1910.0   1977.0
75%        649411.0         2540.0   2000.0
max      11644855.0        10740.0   2015.0

import numpy as np
from scipy import stats

price = house['AdjSalePrice'].values
age = 2024 - house['YrBuilt'].values

print(f"평균/중앙값 = {price.mean() / np.median(price):.3f},  최댓값/중앙값 = {price.max() / np.median(price):.1f}")
print(f"가격 왜도 = {stats.skew(price):.2f}")
yrs = np.unique(house['YrBuilt'])
print(f"서로 다른 건축연도 {yrs.size}개, 빠진 해 {len(set(range(1900, 2016)) - set(yrs.tolist()))}개")
counts = house['YrBuilt'].value_counts()
print(f"연도별 건수 최소 {counts.min()}({counts.idxmin()}) ~ 최대 {counts.max()}({counts.idxmax()})")
print(f"나이 20년 이하 {(age <= 20).sum()}건,  100년 이상 {(age >= 100).sum()}건")

평균/중앙값 = 1.199,  최댓값/중앙값 = 24.7
가격 왜도 = 5.85
서로 다른 건축연도 116개, 빠진 해 0개
연도별 건수 최소 20(2015) ~ 최대 1266(2006)
나이 20년 이하 4170건,  100년 이상 2256건

import numpy as np
import pandas as pd
from patsy import dmatrix, build_design_matrices
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error

# 건축연도 대신 집의 나이를 쓴다. 나이와 가격의 관계는 직선이 아니다.
age = 2024 - house['YrBuilt'].values
df = pd.DataFrame({'age': age, 'price': house['AdjSalePrice'].values})

# patsy 의 bs() 가 B-스플라인 기저를 만들어 준다. df=4 면 매듭을 자료의
# 분위수에 맞춰 알아서 놓는다. 끝의 -1 은 patsy 가 붙이는 절편을 빼는 것이고,
# 절편은 아래 LinearRegression 이 따로 넣는다.
bs_design = dmatrix("bs(age, df=4, degree=3, include_intercept=False) - 1",
                    {"age": df['age']}, return_type='dataframe')
bs_model = LinearRegression().fit(bs_design, df['price'])
bs_r2 = r2_score(df['price'], bs_model.predict(bs_design))

def patsy_knots(design):
    """patsy 가 design_info 에 기억해 둔 매듭벡터를 꺼낸다."""
    info = list(design.design_info.factor_infos.values())[0]
    return list(info.state['transforms'].values())[0]._all_knots

print("기저 열 개수 =", bs_design.shape[1])
print("매듭벡터     =", patsy_knots(bs_design), "  길이", len(patsy_knots(bs_design)))
print(f"내부 매듭 {patsy_knots(bs_design)[4]:.0f},  나이의 중앙값 {np.median(df['age']):.0f}")
print(f"R^2 = {bs_r2:.4f}")

기저 열 개수 = 4
매듭벡터     = [  9.   9.   9.   9.  47. 124. 124. 124. 124.]   길이 9
내부 매듭 47,  나이의 중앙값 47
R^2 = 0.0346

# 매듭 자리를 직접 정할 수도 있다. 관계가 꺾인다고 볼 만한 근거가 있으면
# 분위수에 맡기는 것보다 낫다.
knots_custom = [20, 40, 60]
bs_custom_design = dmatrix(
    f"bs(age, knots={knots_custom}, degree=3, include_intercept=False) - 1",
    {"age": df['age']}, return_type='dataframe'
)
bs_custom_model = LinearRegression().fit(bs_custom_design, df['price'])

print("기저 열 개수 =", bs_custom_design.shape[1])
print("매듭벡터     =", patsy_knots(bs_custom_design))
print(f"지정 매듭 [20, 40, 60]  R^2 = "
      f"{r2_score(df['price'], bs_custom_model.predict(bs_custom_design)):.4f}")

# 열 개수가 같은 자동 배치와 견준다.
auto6 = dmatrix("bs(age, df=6, degree=3, include_intercept=False) - 1",
                {"age": df['age']}, return_type='dataframe')
auto6_model = LinearRegression().fit(auto6, df['price'])
print(f"자동 매듭 {patsy_knots(auto6)[4:7]}  R^2 = "
      f"{r2_score(df['price'], auto6_model.predict(auto6)):.4f}")

기저 열 개수 = 6
매듭벡터     = [  9   9   9   9  20  40  60 124 124 124 124]
지정 매듭 [20, 40, 60]  R^2 = 0.0352
자동 매듭 [24. 47. 74.]  R^2 = 0.0349

# 자연 3차 스플라인. 양 끝에서 직선이 되도록 묶어 두어, 자료가 드문
# 바깥쪽에서 곡선이 크게 튀는 일을 막는다. B-스플라인의 약점을 고친 것이다.
cs_design = dmatrix("cr(age, df=4) - 1",
                    {"age": df['age']}, return_type='dataframe')
cs_model = LinearRegression().fit(cs_design, df['price'])
cs_r2 = r2_score(df['price'], cs_model.predict(cs_design))

print("기저 열 개수 =", cs_design.shape[1], " 매듭 =", patsy_knots(cs_design).round(2))
print(f"행별 합  최소 {cs_design.sum(axis=1).min():.12f}  최대 {cs_design.sum(axis=1).max():.12f}")

# 절편 열을 하나 더 붙이면 어떻게 되는가. 특이값이 바로 말해 준다.
X1 = np.column_stack([np.ones(len(cs_design)), cs_design.values])
print("절편을 붙인 설계행렬의 특이값 =", np.linalg.svd(X1, compute_uv=False).round(4))
print(f"열은 {X1.shape[1]}개인데 계수(rank)는 {np.linalg.matrix_rank(X1)}")
print(f"LinearRegression() 계수의 최대 절대값 = 10^{np.log10(np.abs(cs_model.coef_).max()):.1f}")

# 절편을 떼고 다시 적합한다.
cs_stable = LinearRegression(fit_intercept=False).fit(cs_design, df['price'])
Ck = pd.DataFrame(np.asarray(build_design_matrices(
    [cs_design.design_info], {"age": patsy_knots(cs_design)})[0]),
    columns=cs_design.columns)
B = Ck.values
print(f"매듭에서의 기저값: 대각 {np.diag(B).round(6)},  비대각 최대 "
      f"{np.abs(B - np.diag(np.diag(B))).max():.1e}")
print("계수             =", cs_stable.coef_.round(0))
print("매듭에서의 적합값 =", cs_stable.predict(Ck).round(0))
print(f"R^2 = {cs_r2:.4f} (절편 포함) 대 "
      f"{r2_score(df['price'], cs_stable.predict(cs_design)):.4f} (절편 없이)")
print(f"예측의 최대 차이 = "
      f"{np.abs(cs_model.predict(cs_design) - cs_stable.predict(cs_design)).max():.0f}")

기저 열 개수 = 4  매듭 = [  9.    47.33  85.67 124.  ]
행별 합  최소 1.000000000000  최대 1.000000000000
절편을 붙인 설계행렬의 특이값 = [175.5037  77.3597  52.5336  29.7136   0.    ]
열은 5개인데 계수(rank)는 4
LinearRegression() 계수의 최대 절대값 = 10^17.2
매듭에서의 기저값: 대각 [1. 1. 1. 1.],  비대각 최대 4.6e-17
계수             = [713887. 515523. 508845. 673888.]
매듭에서의 적합값 = [713887. 515523. 508845. 673888.]
R^2 = 0.0301 (절편 포함) 대 0.0301 (절편 없이)
예측의 최대 차이 = 5936

age_grid = np.linspace(df['age'].min(), df['age'].max(), 300)

# 예측할 때는 반드시 원래 기저를 그대로 써야 한다. 새 자료로 dmatrix 를
# 다시 부르면 매듭이 그 자료의 분위수로 다시 정해져, 훈련 때와 다른 기저가
# 만들어진다. design_info 를 넘겨 주는 것이 그 실수를 막는 방법이다.
bs_grid = build_design_matrices([bs_design.design_info], {"age": age_grid})[0]
bs_pred = bs_model.predict(np.asarray(bs_grid))

cs_grid = build_design_matrices([cs_design.design_info], {"age": age_grid})[0]
cs_pred = cs_model.predict(np.asarray(cs_grid))

# 해서는 안 되는 쪽. 식 문자열로 dmatrix 를 다시 부른다.
bs_wrong_design = dmatrix("bs(age, df=4, degree=3, include_intercept=False) - 1",
                          {"age": age_grid}, return_type='dataframe')
bs_wrong = bs_model.predict(np.asarray(bs_wrong_design))

print(f"올바른 기저의 내부 매듭   = {patsy_knots(bs_design)[4]:.1f}  (자료 22687건의 중앙값)")
print(f"다시 부른 기저의 내부 매듭 = {patsy_knots(bs_wrong_design)[4]:.1f}  (격자 300점의 중앙값)")
gap = np.abs(bs_wrong - bs_pred)
print(f"두 예측의 최대 차이 = {gap.max():,.0f}   (age = {age_grid[gap.argmax()]:.1f})")
print(f"양 끝 age = 9, 124 에서의 차이 = {gap[0]:.0f}, {gap[-1]:.0f}")
print(f"예측값 자체의 범위 = {bs_pred.min():,.0f} ~ {bs_pred.max():,.0f}")

올바른 기저의 내부 매듭   = 47.0  (자료 22687건의 중앙값)
다시 부른 기저의 내부 매듭 = 66.5  (격자 300점의 중앙값)
두 예측의 최대 차이 = 39,608   (age = 46.7)
양 끝 age = 9, 124 에서의 차이 = 0, 0
예측값 자체의 범위 = 462,491 ~ 656,578

for d in range(3, 9):
    design = dmatrix(f"bs(age, df={d}, degree=3, include_intercept=False) - 1",
                     {"age": df['age']}, return_type='dataframe')
    model = LinearRegression().fit(design, df['price'])
    grid_design = build_design_matrices([design.design_info],
                                       {"age": age_grid})[0]
    print(f"df={d}: R^2 = {model.score(design, df['price']):.4f}")

df=3: R^2 = 0.0292
df=4: R^2 = 0.0346
df=5: R^2 = 0.0351
df=6: R^2 = 0.0349
df=7: R^2 = 0.0353
df=8: R^2 = 0.0354

age_extrap = np.array([0.0, 5.0, 145.0, 150.0])

# 자연 스플라인: 선형으로 외삽한다.
cs_extrap = build_design_matrices([cs_design.design_info],
                                  {"age": age_extrap})[0]
print("Natural:", cs_model.predict(np.asarray(cs_extrap)).round(0))

# B-스플라인: 매듭 바깥에서는 예외를 던진다.
try:
    bs_extrap = build_design_matrices([bs_design.design_info],
                                      {"age": age_extrap})[0]
    print("B-spline:", bs_model.predict(np.asarray(bs_extrap)).round(0))
except Exception as exc:
    print(f"B-spline: {type(exc).__name__}")

Natural: [772512. 741152. 775648. 801248.]
B-spline: PatsyError

PatsyError: Error evaluating factor: NotImplementedError:
some data points fall outside the outermost knots,
and I'm not sure how to handle them.

from sklearn.linear_model import Ridge

design = dmatrix("bs(age, df=8, degree=3, include_intercept=False) - 1",
                 {"age": df['age']}, return_type='dataframe')
ridge_model = Ridge(alpha=1e4).fit(design, df['price'])