house price - 3 - 모델

나는 토끼가 좋은 사람

2편에 이어서..

사실 전처리 과정이 다 끝난 건 아닌데, 새로 특성도 추가했으므로 Skew된 정도를 파악하고 Boxcox 변환으로 표준화 해주는 과정을 거친다.

1. Skew?

=> '비스듬한' 이런 뜻인데, 비대칭적인 데이터일수록 Skewness가 높다. 이를 펴줘야 더 좋은 모델을 만들 수 있다.

https://dining-developer.tistory.com/18

비대칭(skewed) 데이터를 처리하는 3가지 방법 / Skewed Data

2. Boxcox?

=> 표준화할 때 쓰이는 변환 중 하나다. 식을 자세히 들여다보면 확률과 통계에서 볼 수 있는 표준화 식과 상당히 유사하다.

https://seeyapangpang.tistory.com/35

Box Cox Transformation 파이썬 박스 콕스 변환 Skew 조절 [빅공남! 통계 같이해요]

# skew 라이브러리 불러오기.
from scipy.stats import skew

# one_hot_encoding의 열(column) 중에서 숫자형 열을 모두 numerics2에 저장하기.
numeric_dtypes = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']
numerics2 = []
for i in one_hot_encoding.columns:
    if one_hot_encoding[i].dtype in numeric_dtypes: 
        numerics2.append(i)

# 숫자형 열들을 skew한 정도에 따라서 내림차순으로 정리하기.
skew_features = one_hot_encoding[numerics2].apply(lambda x: skew(x)).sort_values(ascending=False)

# pandas의 dataframe 형태로 정리하기.
skews = pd.DataFrame({'skew':skew_features})

# boxcox 라이브러리 불러오기.
from scipy.special import boxcox1p
from scipy.stats import boxcox_normmax

# 위에서 정리한 skew 정도 중에 0.5 초과인 피처들을 따로 추출하기.
high_skew = skew_features[skew_features > 0.5]

# 인덱스만 추출하기.
skew_index = high_skew.index

# 추출한 인덱스를 갖고 boxcox 적용하기. (Sigmoid랑 식이 똑같다.)
for i in skew_index:
    one_hot_encoding[i]= boxcox1p(one_hot_encoding[i], boxcox_normmax(one_hot_encoding[i]+1))

이렇게 정리하다보니 미숙한 게 많이 보이는게, 표준화 혹은 정규화한 데이터가 이미 이 전에 있다는 사실이다.

만약에 이렇게 진행될 것이었다면, 처음부터 전처리 과정에 집중한 다음에 분포를 펴주는 작업을 하는 것이 더 좋았을 것 같다. 한편으로는, StandardScaler 등 다른 방법을 사용하지 않았음이 조금 아쉽다. 전문서적에서는 제일 처음 배우는 것인데도...

아무튼 그 때 당시에는 장기프로젝트 여러 개를 동시에 진행하다보니 조금 바빴던 게, 이렇게 다시 들여다보니 디테일이나 퀄리티 적인 측면에서 들어나는 것 같아 조금 부끄럽다..

# boxcox 작업까지 마친 데이터들을 다시 train과 test로 나눠준다.
one_hot_encoding_train = one_hot_encoding[:len_train]
one_hot_encoding_test = one_hot_encoding[len_train:]

# 여기서 one_hot_encoding_train과 test는 위에 동일한 이름의 변수가 있지만
# boxcox가 적용된 점이 다르다.

# 'SalePrice' 피처 제거해주기.
one_hot_encoding_train = one_hot_encoding_train.drop(['SalePrice'], axis=1)
one_hot_encoding_test = one_hot_encoding_test.drop(['SalePrice'], axis=1)

다만, 한 가지 칭찬해주고 싶은 점은, 주석을 꽤 꼼꼼하게 적었다는 것.

 

이렇게까지 했을 때 오차는?

#중 간 체 크
subset=one_hot_encoding_train[['OverallQual', 'YearBuilt', 'YearRemodAdd', 'TotalBsmtSF', '1stFlrSF',
       'GrLivArea', 'FullBath', 'TotRmsAbvGrd', 'GarageCars', 'GarageArea',
       'ExterQual_Gd', 'ExterQual_TA', 'Foundation_PConc',
       'KitchenQual_TA', 'Total_sqr_footage', 'Total_Bathrooms', 'TotalHouse',
       'TotalArea', 'GrLivArea_OverallQual']]
       
y = y_reg
X = subset

reg=LinearRegression()
reg.fit(X,y)

y_train_error = reg.predict(X)
rmse=(np.sqrt(mean_squared_error(np.expm1(y), np.expm1(y_train_error))))

print("RMSE is {}".format(rmse))

$24845 수준이 되었다.

https://brunch.co.kr/@itschloe1/11

Ridge와 Lasso Regression의 쉬운 풀이

1. Ridge

Ridge 모델을 도입해서 문제를 해결해보자.

#ridge를 해보자

kfolds = KFold(n_splits=5, shuffle=True, random_state=42)

#cv_rmse 함수는 cross_val_score 를 함수화한 것
def cv_rmse(model):
    rmse = np.sqrt(-cross_val_score(model, one_hot_encoding_train, y_reg, 
                                   scoring="neg_mean_squared_error", 
                                   cv = kfolds))
    return(rmse)

#pipeline을 만들어서 좀 더 유연하고 빠르게, 그리고 동시에 여러 모델을 만듬
def ridge_selector(k):
    ridge_model = make_pipeline(RobustScaler(),
                                RidgeCV(alphas = [k],
                                        cv=kfolds)).fit(one_hot_encoding_train, y_reg)
    
    ridge_rmse = cv_rmse(ridge_model).mean()
    return(ridge_rmse)

r_alphas = [.0001, .0003, .0005, .0007, .0009, .01, 0.05, 0.1, 0.3, 1, 3, 5, 10, 15, 20, 30, 50, 60, 70, 80]

ridge_scores = []
for alpha in r_alphas:
    score = ridge_selector(alpha)
    ridge_scores.append(score)
    
plt.plot(r_alphas, ridge_scores, label='Ridge')
plt.legend('center')
plt.xlabel('alpha')
plt.ylabel('score')

본격적으로 모델을 도입해서 문제를 해결해보자. 

Pipeline으로 모델을 만들었는데, Ridge는 규제항이 매우 중요하게 작용하는데, 규제항의 계수를 쉽게 찾을 수 있도록 여러 모델을 동시에 만들고 검증해서 가장 좋은 계수를 찾고 도와준다.

RobustScaler는 이름부터 알 수 있듯이, 정규화 하는 함수이고

RidgeCV는 Ridge + Cross Validation의 약자로서, 주 모델은 Ridge, 검증 모델은 kfolds 기법을 사용하겠다는 의미.

Score가 가장 낮은 지점에 최적화된 규제항을 찾을 수 있다.

#ridge를 좀 더 구체적으로 

alphas_alt = [9,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8,9.9,10,10.1,10.2]

ridge_model2 = make_pipeline(RobustScaler(),
                            RidgeCV(alphas = alphas_alt,
                                    cv=kfolds)).fit(one_hot_encoding_train, y_reg)

print("Best of alpha in ridge model :" ,ridge_model2.steps[1][1].alpha_)


print("Ridge rmse : ",cv_rmse(ridge_model2))
print("Ridge rmse : ",cv_rmse(ridge_model2).mean())

#cv_rmse 결과 평균 0.1107 정도의 값을 얻었다. 달러로 환산시 대략 19000~20000 달러 정도 예상

 

2. Lasso

#Lasso에 대한 것

from sklearn.linear_model import LassoCV

alphas2 = [0.00005, 0.0001, 0.0002, 0.0003, 0.0004, 0.0005,
           0.0006, 0.0007, 0.0008]


lasso_model2 = make_pipeline(RobustScaler(),
                             LassoCV(max_iter=1e5,
                                    alphas = alphas2,
                                    random_state = 42)).fit(one_hot_encoding_train, y_reg)
scores = lasso_model2.steps[1][1].mse_path_

plt.plot(alphas2, scores, label='Lasso')
plt.legend(loc='center')
plt.xlabel('alpha')
plt.ylabel('RMSE')
plt.tight_layout()
plt.show()

Lasso는 Ridge와 아주 유사하다.

역시 비슷한 방법으로 pipeline을 도입했고, 다음과 같은 결과를 얻을 수 있었다.

print("Best of alpha in lasso model :",lasso_model2.steps[1][1].alpha_)

print("lasso rmse : ",cv_rmse(lasso_model2))
print("lasso rmse : ",cv_rmse(lasso_model2).mean())

#cv_rmse 결과 평균 0.11046 정도의 값을 얻었다. 달러로 환산시 대략 19000~20000 달러 정도 예상

Lasso 모델 내부에서 가장 좋은 값의 alpha값을 가져오고 마무리했다.

여기까지, 전체적인 결과는 다음과 같다. (달러화 포함)

#전체의 결과를 포괄적으로 얘기한 것, CV 결과는 아니니 너무 신뢰하지는 말자

y_train_error_ridge = ridge_model2.predict(one_hot_encoding_train)

rmse=(np.sqrt(mean_squared_error(np.expm1(y_reg), np.expm1(y_train_error_ridge))))

print("train_RMSE(달러화)_on_RIDGE is {}".format(rmse))

rmse=(np.sqrt(mean_squared_error(y_reg, y_train_error_ridge)))

print("train_RMSE_on_RIDGE is {}".format(rmse))

y_train_error_lasso = lasso_model2.predict(one_hot_encoding_train)

rmse=(np.sqrt(mean_squared_error(np.expm1(y_reg), np.expm1(y_train_error_lasso))))

print("train_RMSE(달러화)_on_LASSO is {}".format(rmse))

rmse=(np.sqrt(mean_squared_error(y_reg, y_train_error_lasso)))

print("train_RMSE_on_LASSO is {}".format(rmse))

Ridge : 17642$

Lasso : 18433$

3. ElasticNet

Ridg는 L1 Norm, Lasso L2 Norm이라면, 엘라스틱 넷은 두 형태를 흡수한 형태인데, 다음과 같은 형태를 가진다.

출처 https://teddylee777.github.io/scikit-learn/scikit-learn-linear-with-regularizations

#ElasticNet

from sklearn.linear_model import ElasticNetCV

e_alphas = [0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 0.0007]
e_l1ratio = [0.8, 0.85, 0.9, 0.95, 0.99, 1]

elastic_cv = make_pipeline(RobustScaler(), 
                           ElasticNetCV(max_iter=1e5, alphas=e_alphas, 
                                        cv=kfolds, l1_ratio=e_l1ratio))

elastic_model3 = elastic_cv.fit(one_hot_encoding_train, y_reg)


print("Best of alpha in elastic_net :",elastic_model3.steps[1][1].alpha_)

print("Best of ratio in elastic_net :",elastic_model3.steps[1][1].l1_ratio_)

y_train_error_elasticNET = elastic_model3.predict(one_hot_encoding_train)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_reg), np.expm1(y_train_error_elasticNET))))
print("train_RMSE_for_ElasticNET is {}".format(rmse))

print("lasso rmse : ",cv_rmse(elastic_model3))
print("lasso rmse : ",cv_rmse(elastic_model3).mean())

#cv_rmse 결과 평균 0.110675 정도의 값을 얻었다. 달러로 환산시 대략 19000~20000 달러 정도 예상

ElasticNet : 18327$

4. RandomForest

여기서 그치지 않고 익숙한 모델인 결정나무를 도입해본다.

# RandomForestRegressor 모델 만들어보기~
from sklearn.ensemble import RandomForestRegressor

# 객채 생성하기.
rfModel = RandomForestRegressor(n_estimators=100)

# 모델 만들기.
rfModel.fit(x_train, y_train)

y_train_rf = rfModel.predict(x_train)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_train), np.expm1(y_train_rf))))
print("train_RMSE_for_RandomForestRegressor is {}".format(rmse))

print("RandomForestRegressor rmse : ",cv_rmse(rfModel))
print("RandomForestRegressor rmse mean : ",cv_rmse(rfModel).mean())

# 객채 생성하기.
rfModel = RandomForestRegressor(n_estimators=300)

# 모델 만들기.
rfModel.fit(x_train, y_train)

y_train_rf = rfModel.predict(x_train)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_train), np.expm1(y_train_rf))))
print("train_RMSE_for_RandomForestRegressor is {}".format(rmse))

print("RandomForestRegressor rmse : ",cv_rmse(rfModel))
print("RandomForestRegressor rmse mean : ",cv_rmse(rfModel).mean())

생성 나무 갯수를 100/300개로 하니, 10000$ 이하의 결과를 출력한다. 다만, 나무 갯수가 많으므로 아무래도 Overfitting 된 결과라는 생각이 들었다.

 

5.XGBoost

여러 모델을 Greedy Algorithm으로 결합해 만드는 XGBoost 모델도 적용해봤다.

https://brunch.co.kr/@snobberys/137

XGBoost 사용하기

 

윗 글이 참 정리가 잘되어 있다고 생각한다. 모델의 에러를 줄이는 과정이 잘 나타나있다.

학습을 할수록 오차를 줄이는 방향으로, Greedy 하게 움직인다.

# XGBoost

import xgboost as xgb
from xgboost import XGBRegressor

xgb = XGBRegressor(learning_rate =0.01, n_estimators=1000, max_depth=3,
                     min_child_weight=0 ,gamma=0, subsample=0.7,
                     colsample_bytree=0.7,objective= 'reg:linear',
                     nthread=4,scale_pos_weight=1,seed=27, reg_alpha=0.00006)

xgb = xgb.fit(one_hot_encoding_train, y_reg)
y_train_xgb = rfModel.predict(one_hot_encoding_train)

rmse=(np.sqrt(mean_squared_error(np.expm1(y_reg), np.expm1(y_train_xgb))))
print("train_RMSE_for_RandomForestRegressor is {}".format(rmse))

print("Xgboost model rmse : ",cv_rmse(xgb))
print("Xgboost model rmse : ",cv_rmse(xgb).mean())

Warning은 무시하자.

train에 한하여 9507$가 출력되었다.

6. GradientBoostingRegressor

GradientBoostingRegressor는 새로운 예측기를 반복마다 추가하는데, 전 모델의 오차를 가져와서 학습한다.

https://biology-statistics-programming.tistory.com/49

[Concept] 회귀를 위한 다양한 모델 (feat. Regressor)

결과가 엄청 좋게나왔는데, 반복수를 500개로 엄청 크게 했다는 점.

그리고 해당 모델을 잘모르는 시점에서 함 써보는 느낌으로 사용했다.

5/6은 솔직히 모델에 관해서 깊게 알지 못하는데 다음에는, 링크로 고수분들의 글이 아니라 내 글을 올릴 수 있도록 해봐야겠다... ㅠ

7. 앙상블 모델

앙상블 기법을 도입해 모델을 깎아보자.

여러 모델을 합쳐서 투표 및 연산해서 Target을 예측한다.

https://flonelin.wordpress.com/2016/08/02/stacking%EC%9D%B4%EB%9E%80

Stacking이란

https://techblog-history-younghunjo1.tistory.com/103

[ML] Scikit-learn을 이용한 Stacking 구현하기

(아래가 비교적 더 정리가 잘되어 있다.)

앙상블 모델은 여러 모델을 결합해서 만드는 모델이다.

Regressor 모델이 아니라 Classifier 모델이라면 개별 모델의 '투표'(각 모델이 예측한 값) 혹은 확률을 합하여 최종 예측 값을 결정하는 과정을 거치지만, Regressor 모델의 경우 개별 모델의 결과를 합치고 적절히 나눠서 결과를 출력한다.

train_ridge_preds = ridge_model.predict(one_hot_encoding_train)
train_lasso_preds = lasso_model.predict(one_hot_encoding_train)
train_elastic_preds = elastic_model.predict(one_hot_encoding_train)

stackX = pd.DataFrame(np.transpose(np.array([train_ridge_preds,train_lasso_preds,train_elastic_preds])))

model = Ridge(alpha=14.2)

stack_fit = model.fit(stackX, y_reg)

test_ridge_preds = ridge_model.predict(one_hot_encoding_test)
test_lasso_preds = lasso_model.predict(one_hot_encoding_test)
test_elastic_preds = elastic_model.predict(one_hot_encoding_test)

test_stackX = pd.DataFrame(np.transpose(np.array([test_ridge_preds,test_lasso_preds,test_elastic_preds])))
stack_pred = np.expm1(stack_fit.predict(test_stackX))
submission = pd.DataFrame({'Id': Id, 'SalePrice': stack_pred})
submission.to_csv('Ridge_stacking.csv', index=False)

간단하게, 서로 다른 새 모델의 예측값을 DataFrame으로 새로 만들어서 합쳐주고 해당 값들을 다시 Ridge model로 학습시켜줬다.

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(stackX, y_reg, test_size = .30, random_state=0)

y_train_error = stack_fit.predict(x_train)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_train), np.expm1(y_train_error))))
print("stack... is {}".format(rmse))


rmse=(np.sqrt(mean_squared_error(y_train, y_train_error)))
print("stack... is {}".format(rmse))

y_test_error = stack_fit.predict(x_test)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_test), np.expm1(y_test_error))))
print("test stack... is {}".format(rmse))


rmse=(np.sqrt(mean_squared_error(y_test, y_test_error)))
print("test stack... is {}".format(rmse))

여기서 Test stack이 오차가 제공되는데, 윗사진을 보다시피 kaggle에서 예시로 추가한 Submission.csv 파일을 test 파일로 오해했기 때문에 벌어졌다.

해당 모델이 실제 Test.csv 를 얼마나 잘 예측하는 지에 대한 지표는 Kaggle 상위권의 고수 분들의 예측결과와 비교해볼 예정이다.

train_ridge_preds = ridge_model.predict(one_hot_encoding_train)
train_lasso_preds = lasso_model.predict(one_hot_encoding_train)
train_elastic_preds = elastic_model.predict(one_hot_encoding_train)

stackX = pd.DataFrame(np.transpose(np.array([train_ridge_preds,train_lasso_preds,train_elastic_preds])))

model = Lasso(alpha=0.00044)

stack_fit = model.fit(stackX, y_reg)

test_ridge_preds = ridge_model.predict(one_hot_encoding_test)
test_lasso_preds = lasso_model.predict(one_hot_encoding_test)
test_elastic_preds = elastic_model.predict(one_hot_encoding_test)

test_stackX = pd.DataFrame(np.transpose(np.array([test_ridge_preds,test_lasso_preds,test_elastic_preds])))
stack_pred = np.expm1(stack_fit.predict(test_stackX))
submission = pd.DataFrame({'Id': Id, 'SalePrice': stack_pred})
submission.to_csv('lasso_stacking.csv', index=False)

Lasso 모델로 학습시켜보자

y_train_error = stack_fit.predict(x_train)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_train), np.expm1(y_train_error))))
print("stack... is {}".format(rmse))


rmse=(np.sqrt(mean_squared_error(y_train, y_train_error)))
print("stack... is {}".format(rmse))

y_test_error = stack_fit.predict(x_test)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_test), np.expm1(y_test_error))))
print("test stack... is {}".format(rmse))


rmse=(np.sqrt(mean_squared_error(y_test, y_test_error)))
print("test stack... is {}".format(rmse))

Stacking을 하니 기존 모델에 비해 확실히 더 좋은 결과가 나온다.

 

똑같은 방법으로 ElasticNet까지 Stacking을 진행했다.

train_ridge_preds = ridge_model.predict(one_hot_encoding_train)
train_lasso_preds = lasso_model.predict(one_hot_encoding_train)
train_elastic_preds = elastic_model.predict(one_hot_encoding_train)

stackX = pd.DataFrame(np.transpose(np.array([train_ridge_preds,train_lasso_preds,train_elastic_preds])))

model = ElasticNet(alpha=0.0005, l1_ratio = 0.85)

stack_fit = model.fit(stackX, y_reg)

test_ridge_preds = ridge_model.predict(one_hot_encoding_test)
test_lasso_preds = lasso_model.predict(one_hot_encoding_test)
test_elastic_preds = elastic_model.predict(one_hot_encoding_test)

test_stackX = pd.DataFrame(np.transpose(np.array([test_ridge_preds,test_lasso_preds,test_elastic_preds])))
stack_pred = np.expm1(stack_fit.predict(test_stackX))
submission = pd.DataFrame({'Id': Id, 'SalePrice': stack_pred})
submission.to_csv('ElasticNet_stacking.csv', index=False)

y_train_error = stack_fit.predict(x_train)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_train), np.expm1(y_train_error))))
print("stack... is {}".format(rmse))


rmse=(np.sqrt(mean_squared_error(y_train, y_train_error)))
print("stack... is {}".format(rmse))

y_test_error = stack_fit.predict(x_test)
rmse=(np.sqrt(mean_squared_error(np.expm1(y_test), np.expm1(y_test_error))))
print("test stack... is {}".format(rmse))


rmse=(np.sqrt(mean_squared_error(y_test, y_test_error)))
print("test stack... is {}".format(rmse))

 

본격적으로 고수의 예측 모델과 비교해보자

3년 전 당시, 4000개 이상의 좋아요를 받은, 높은 순위의 고수분의 데이터와 비교했다. (일단 300개 까지의 비교지만 300까지가 아니어도 전체적으로 비슷한 양상을 보인다.. 그 양상은 아래를 확인해보자)

첫번째 비교는 Ridge...

RMSE between gosu is is 9761.704269580167

거의 유사한 결과를 얻을 수 있었다.

두번째 비교는? Lasso

RMSE between gosu is is 10720.443267414317

이것 역시 아주 유사하다!!

세번째는 ElasticNet

RMSE between gosu is is 10704.434230909543

네번째는 적절히 Stacking 시켜서 출력해줬다.

RMSE between gosu is is 10244.368460758933

 

결과적으로 스태킹한 모델이 제일 유사한 결과를 얻을 수 있었다.

 

PBL은 이렇게 마지막 비교를 보여주고 끝냈다. 개인적으로는 그때 당시에 제대로 진지하게 공부하지 못하게 급하게 진행한 티가 나서 조금 아쉬웠던 것 같다. 이 때를 반면교사 삼아 제대로 발전하길!