회귀문제에서 Random Forest, SVM, XGBoost 튜닝하는법 (in python)

Tuning practice (SVM / RF / XGBoost) - Regression (PRSA data)¶

오늘은 튜닝이 필요한 머신러닝 모형들 중에 state of the art의 예측력으로 알려진 SVM(Support Vector Machine), RF(Random Forest), XGBoost 모델의 회귀문제 튜닝 예시를 가져와보았다. 먼저 데이터는 UCI Repository에 있는 PASA_data 이며, 중국의 미세먼지 데이터이다. 데이터의 속성은 아래와 같다.

Attribute Information: (PASA_data)¶

No: row number

year: year of data in this row

month: month of data in this row

day: day of data in this row

hour: hour of data in this row

pm2.5: PM2.5 concentration (ug/m^3)

DEWP: Dew Point (℃)

TEMP: Temperature (℃)

PRES: Pressure (hPa)

cbwd: Combined wind direction

Iws: Cumulated wind speed (m/s)

Is: Cumulated hours of snow

Ir: Cumulated hours of rain

EDA and Data preprocessing¶

In [1]:

import pandas as pd
import numpy as np
import seaborn as sns
pd.set_option('display.max_rows',None)

데이터를 불러온다.

In [2]:

df=pd.read_csv("PRSA_data_2010.1.1-2014.12.31.csv")

In [3]:

df.head(10)

Out[3]:

	No	year	month	day	hour	pm2.5	DEWP	TEMP	PRES	cbwd	Iws	Is	Ir
0	1	2010	1	1	0	NaN	-21	-11.0	1021.0	NW	1.79	0	0
1	2	2010	1	1	1	NaN	-21	-12.0	1020.0	NW	4.92	0	0
2	3	2010	1	1	2	NaN	-21	-11.0	1019.0	NW	6.71	0	0
3	4	2010	1	1	3	NaN	-21	-14.0	1019.0	NW	9.84	0	0
4	5	2010	1	1	4	NaN	-20	-12.0	1018.0	NW	12.97	0	0
5	6	2010	1	1	5	NaN	-19	-10.0	1017.0	NW	16.10	0	0
6	7	2010	1	1	6	NaN	-19	-9.0	1017.0	NW	19.23	0	0
7	8	2010	1	1	7	NaN	-19	-9.0	1017.0	NW	21.02	0	0
8	9	2010	1	1	8	NaN	-19	-9.0	1017.0	NW	24.15	0	0
9	10	2010	1	1	9	NaN	-20	-8.0	1017.0	NW	27.28	0	0

In [4]:

df.describe()

Out[4]:

	No	year	month	day	hour	pm2.5	DEWP	TEMP	PRES	Iws	Is	Ir
count	43824.000000	43824.000000	43824.000000	43824.000000	43824.000000	41757.000000	43824.000000	43824.000000	43824.000000	43824.000000	43824.000000	43824.000000
mean	21912.500000	2012.000000	6.523549	15.727820	11.500000	98.613215	1.817246	12.448521	1016.447654	23.889140	0.052734	0.194916
std	12651.043435	1.413842	3.448572	8.799425	6.922266	92.050387	14.433440	12.198613	10.268698	50.010635	0.760375	1.415867
min	1.000000	2010.000000	1.000000	1.000000	0.000000	0.000000	-40.000000	-19.000000	991.000000	0.450000	0.000000	0.000000
25%	10956.750000	2011.000000	4.000000	8.000000	5.750000	29.000000	-10.000000	2.000000	1008.000000	1.790000	0.000000	0.000000
50%	21912.500000	2012.000000	7.000000	16.000000	11.500000	72.000000	2.000000	14.000000	1016.000000	5.370000	0.000000	0.000000
75%	32868.250000	2013.000000	10.000000	23.000000	17.250000	137.000000	15.000000	23.000000	1025.000000	21.910000	0.000000	0.000000
max	43824.000000	2014.000000	12.000000	31.000000	23.000000	994.000000	28.000000	42.000000	1046.000000	585.600000	27.000000	36.000000

In [5]:

df2=df.dropna() #remove the rows contain any NA.

In [6]:

df2.describe()

Out[6]:

	No	year	month	day	hour	pm2.5	DEWP	TEMP	PRES	Iws	Is	Ir
count	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000	41757.000000
mean	22279.380104	2012.042771	6.513758	15.685514	11.502311	98.613215	1.750174	12.401561	1016.442896	23.866747	0.055344	0.194866
std	12658.168415	1.415311	3.454199	8.785539	6.924848	92.050387	14.433658	12.175215	10.300733	49.617495	0.778875	1.418165
min	25.000000	2010.000000	1.000000	1.000000	0.000000	0.000000	-40.000000	-19.000000	991.000000	0.450000	0.000000	0.000000
25%	11464.000000	2011.000000	4.000000	8.000000	5.000000	29.000000	-10.000000	2.000000	1008.000000	1.790000	0.000000	0.000000
50%	22435.000000	2012.000000	7.000000	16.000000	12.000000	72.000000	2.000000	14.000000	1016.000000	5.370000	0.000000	0.000000
75%	33262.000000	2013.000000	10.000000	23.000000	18.000000	137.000000	15.000000	23.000000	1025.000000	21.910000	0.000000	0.000000
max	43824.000000	2014.000000	12.000000	31.000000	23.000000	994.000000	28.000000	42.000000	1046.000000	565.490000	27.000000	36.000000

In [7]:

df2=df2.reset_index()
df2.head()

Out[7]:

	index	No	year	month	day	hour	pm2.5	DEWP	TEMP	PRES	cbwd	Iws	Is	Ir
0	24	25	2010	1	2	0	129.0	-16	-4.0	1020.0	SE	1.79	0	0
1	25	26	2010	1	2	1	148.0	-15	-4.0	1020.0	SE	2.68	0	0
2	26	27	2010	1	2	2	159.0	-11	-5.0	1021.0	SE	3.57	0	0
3	27	28	2010	1	2	3	181.0	-7	-5.0	1022.0	SE	5.36	1	0
4	28	29	2010	1	2	4	138.0	-7	-5.0	1022.0	SE	6.25	2	0

In [8]:

type(df2['cbwd'][0])

Out[8]:

str

In [9]:

df2.info() #one categorical variable

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 41757 entries, 0 to 41756
Data columns (total 14 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   index   41757 non-null  int64  
 1   No      41757 non-null  int64  
 2   year    41757 non-null  int64  
 3   month   41757 non-null  int64  
 4   day     41757 non-null  int64  
 5   hour    41757 non-null  int64  
 6   pm2.5   41757 non-null  float64
 7   DEWP    41757 non-null  int64  
 8   TEMP    41757 non-null  float64
 9   PRES    41757 non-null  float64
 10  cbwd    41757 non-null  object 
 11  Iws     41757 non-null  float64
 12  Is      41757 non-null  int64  
 13  Ir      41757 non-null  int64  
dtypes: float64(4), int64(9), object(1)
memory usage: 4.5+ MB

모델에 풍향변수를 설명변수로 사용하기 위해, 더미변수화 해준다.

In [10]:

cbwd_cat = pd.get_dummies(df2['cbwd'])
cbwd_cat.head() #change the categorical variable to dummy variables

Out[10]:

	NE	NW	SE	cv
0	0	0	1	0
1	0	0	1	0
2	0	0	1	0
3	0	0	1	0
4	0	0	1	0

In [11]:

df2 = pd.concat([df2,cbwd_cat],axis=1) #concatenate
df2.head()

Out[11]:

	index	No	year	month	day	hour	pm2.5	DEWP	TEMP	PRES	cbwd	Iws	Is	Ir	NE	NW	SE	cv
0	24	25	2010	1	2	0	129.0	-16	-4.0	1020.0	SE	1.79	0	0	0	0	1	0
1	25	26	2010	1	2	1	148.0	-15	-4.0	1020.0	SE	2.68	0	0	0	0	1	0
2	26	27	2010	1	2	2	159.0	-11	-5.0	1021.0	SE	3.57	0	0	0	0	1	0
3	27	28	2010	1	2	3	181.0	-7	-5.0	1022.0	SE	5.36	1	0	0	0	1	0
4	28	29	2010	1	2	4	138.0	-7	-5.0	1022.0	SE	6.25	2	0	0	0	1	0

In [12]:

quant_cols=['DEWP','TEMP','PRES','Iws','Is','Ir','pm2.5'] #select quantitative variables
qual_cols=['month','hour','NE','NW','SE','cv'] #select categorical data

In [13]:

corr_df=pd.DataFrame(np.corrcoef(np.array(df2[quant_cols]).T))

In [14]:

corr_df.columns=quant_cols

In [15]:

corr_df.index = quant_cols

변수별 상관계수는 아래와 같다.

In [16]:

import matplotlib.pyplot as plt
import seaborn as sns
ax = sns.heatmap(corr_df, annot=True, annot_kws=dict(color='g'), cmap='Greys')
plt.show()

풍속과 pm2.5 간에는 약간의 음의 상관관계가 있어보이고 이슬점과는 양의 상관관계가 있어보인다.

In [17]:

corr_df['pm2.5'] #correlation with pm2.5

Out[17]:

DEWP     0.171423
TEMP    -0.090534
PRES    -0.047282
Iws     -0.247784
Is       0.019266
Ir      -0.051369
pm2.5    1.000000
Name: pm2.5, dtype: float64

언뜻보기에는 4~9월에는 미세먼지 농도가 낮고 10월과 1월, 2월 미세먼지 농도가 높은편인것으로 보인다.

In [18]:

sns.barplot(x='month', y='pm2.5', data=df2)

Out[18]:

<AxesSubplot: xlabel='month', ylabel='pm2.5'>

하루중에는 새벽시간대에 미세먼지 농도가 높은것으로 보인다.

In [19]:

sns.barplot(x='hour', y='pm2.5', data=df2)

Out[19]:

<AxesSubplot: xlabel='hour', ylabel='pm2.5'>

풍속이 CV일때 미세먼지 농도가 가장 높은것으로 보인다. CV가 뭐의 줄임말인지는 모르겠지만..

In [20]:

sns.barplot(x='cbwd', y='pm2.5', data=df2)

Out[20]:

<AxesSubplot: xlabel='cbwd', ylabel='pm2.5'>

Modeling¶

이제 본격적으로 모델링에 들어가겠다.추후 객관적인 예측성능 비교를 위해 test set을 미리 분리해둔다.

In [21]:

cols = quant_cols+qual_cols
cols.remove('pm2.5')
X,y = np.array(df2[cols]), np.array(df2["pm2.5"])

In [22]:

X[0:2,:]

Out[22]:

array([[-1.60e+01, -4.00e+00,  1.02e+03,  1.79e+00,  0.00e+00,  0.00e+00,
         1.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  1.00e+00,  0.00e+00],
       [-1.50e+01, -4.00e+00,  1.02e+03,  2.68e+00,  0.00e+00,  0.00e+00,
         1.00e+00,  1.00e+00,  0.00e+00,  0.00e+00,  1.00e+00,  0.00e+00]])

In [23]:

X.shape

Out[23]:

(41757, 12)

In [24]:

np.random.seed(1)
train_idx=np.random.choice(41757,31757,replace=False)
test_idx = np.array(list(set(range(41757))-set(train_idx))) #Take 10000 samples as the test set.

추후에 공정한 예측력 평가를 위해 test set을 따로 떼어둔다.

In [25]:

X_train = X[train_idx,]
y_train = y[train_idx,]

X_test = X[test_idx,]
y_test = y[test_idx,]

In [26]:

print(X_train.shape)
print(X_test.shape)
print(y_train.shape)
print(y_test.shape)

(31757, 12)
(10000, 12)
(31757,)
(10000,)

0. Base (mean of y_train)¶

Baseline 모델로써, 설명변수가 전혀 없는 모델인 y의 평균을 예측값으로 사용해보았다. RMSE가 대략 93정도 나온다.

In [27]:

np.sqrt(np.mean((y_test-np.mean(y_train))**2)) # RMSE of the base (y_mean)

Out[27]:

93.37689990778144

1. OLS¶

비선형 머신러닝 모형에 들어가기에 앞서 단순 선형 회귀 모형을 적합해보았다.

In [28]:

from sklearn.linear_model import LinearRegression

In [29]:

ols_regr = LinearRegression()
ols_regr.fit(X_train,y_train)

Out[29]:

LinearRegression()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

In [30]:

ols_pred = ols_regr.predict(X_test)

pm2.5값이 0보다 작을일은 없으므로, OLS 선형회귀모형의 예측값이 0보다 작은 경우는 0으로 처리해주자.

In [31]:

ols_pred[ols_pred<0,]=0
ols_pred

Out[31]:

array([ 43.44065755, 178.57921986,  67.07205565, ..., 131.9691329 ,
        56.35336049,  40.32194833])

단순 선형 회귀 모형의 RMSE는 대략 79 정도 나온다.

In [32]:

np.sqrt(np.mean((y_test-ols_pred)**2)) #RMSE of OLS

Out[32]:

79.23636933768623

In [33]:

ols_regr.coef_

Out[33]:

array([  4.40215676,  -6.58870573,  -1.61106348,  -0.20029049,
        -4.12796092,  -6.32214068,  -1.15246671,   1.32612634,
       -11.72228683, -14.11758984,  13.0342684 ,  12.80560827])

2. RandomForest¶

이제 대표적인 머신러닝 모형 중 하나인 랜덤포레스트를 사용해보자.

In [34]:

from sklearn.ensemble import RandomForestRegressor
import time

먼저, default 값으로써 n_estimator 즉, 나무의 개수는 500개, max_features 즉, 트리가 가지를 칠 때, 참고하는 랜덤한 변수의 최대 개수를 4개로 설정하여 진행해보자.

In [35]:

rf_regr = RandomForestRegressor(n_estimators=500, criterion='squared_error', max_features=4, n_jobs=-1, random_state=0) #criterion is mse

In [37]:

rf_pred = rf_regr.predict(X_test)

pm2.5값이 0보다 작을일은 없으므로, 랜덤포레스트의 예측값이 0보다 작은 경우는 0으로 처리해주자. 물론 트리기반의 랜덤포레스트가 학습과정에서 0보다 작은 반응변수를 만난적이 없으므로, 예측값이 0보다 작을일은 없다. 따라서 그냥 이 처리를 해주지 않아도 무방하다.

In [38]:

rf_pred[rf_pred<0,]=0
rf_pred

Out[38]:

array([ 32.594, 180.072,  18.088, ..., 175.024,  36.154,  25.918])

튜닝 전 랜덤포레스트의 RMSE는 54.918이다.

In [39]:

np.sqrt(np.mean((y_test-rf_pred)**2)) #RMSE of default RF

Out[39]:

54.91833675909555

이제 좀 더 정교하게 튜닝파라미터들을 튜닝해보자. 여기서는 그리드 서치(grid search)와 5 fold cross validation을 사용하였다.

In [40]:

# 튜닝

param_grid = {'n_estimators':[500,1000,1500,2000],
             'max_features':[2,4,6,8,10,12]}

In [42]:

from sklearn.model_selection import GridSearchCV

grid_search = GridSearchCV(rf_regr,param_grid,cv=5,scoring='neg_root_mean_squared_error') 

In [43]:

grid_search.fit(X_train,y_train) 

Out[43]:

GridSearchCV(cv=5,
             estimator=RandomForestRegressor(max_features=4, n_estimators=500,
                                             n_jobs=-1, random_state=0),
             param_grid={'max_features': [2, 4, 6, 8, 10, 12],
                         'n_estimators': [500, 1000, 1500, 2000]},
             scoring='neg_root_mean_squared_error')

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

그리드 서치 결과 최적의 파라미터는 max_features 10, n_estimators 2000으로 밝혀졌다.

In [44]:

grid_search.best_params_

Out[44]:

{'max_features': 10, 'n_estimators': 2000}

그리드 서치 결과 최적의 파라미터에서의 validation score(RMSE)는 54.358이다.

In [45]:

grid_search.best_score_

Out[45]:

-54.367712944460585

다음은 그리드서치 과정을 상세하게 보여주는 테이블이다.

In [46]:

import pandas as pd
results = pd.DataFrame(grid_search.cv_results_)
display(results.head(5))

	mean_fit_time	std_fit_time	mean_score_time	std_score_time	param_max_features	param_n_estimators	params	split0_test_score	split1_test_score	split2_test_score	split3_test_score	split4_test_score	mean_test_score	std_test_score	rank_test_score
0	2.286065	0.228809	0.163676	0.007771	2	500	{'max_features': 2, 'n_estimators': 500}	-59.259871	-58.038878	-60.189908	-59.472687	-58.940318	-59.180333	0.703074	24
1	4.368623	0.122684	0.437073	0.056706	2	1000	{'max_features': 2, 'n_estimators': 1000}	-59.127964	-58.012705	-60.117883	-59.408901	-58.844824	-59.102455	0.689867	23
2	7.429689	0.593906	0.790882	0.081900	2	1500	{'max_features': 2, 'n_estimators': 1500}	-59.102795	-57.969256	-60.130077	-59.341624	-58.863323	-59.081415	0.700331	22
3	9.928616	0.827491	1.218432	0.017750	2	2000	{'max_features': 2, 'n_estimators': 2000}	-59.073619	-57.982308	-60.122460	-59.322955	-58.881232	-59.076515	0.691399	21
4	3.411477	0.214035	0.220291	0.116477	4	500	{'max_features': 4, 'n_estimators': 500}	-56.417897	-55.196058	-57.473388	-56.882318	-56.328091	-56.459550	0.751270	20

다음은 그리드서치 결과를 시각화로 보여준 것이다. max_features와 n_estimators의 값에 따라 validation score(RMSE)가 어떻게 변하는지 보여준다.

In [47]:

import mglearn

scores=np.array(results.mean_test_score).reshape(6,4)

mglearn.tools.heatmap(scores,xlabel='n_estimators', xticklabels=param_grid['n_estimators'], ylabel='max_features',yticklabels=param_grid['max_features'],cmap='viridis')

Out[47]:

<matplotlib.collections.PolyCollection at 0x13c5878e0>

그리드서치로 찾은 최적의 랜덤포레스트로 test set에 대한 예측값을 구해낸다.

In [48]:

rft_pred = grid_search.predict(X_test)

In [49]:

rft_pred[rft_pred<0,]=0
rft_pred

Out[49]:

array([ 41.407 , 214.9665,  18.566 , ..., 206.4885,  32.243 ,  25.8585])

튜닝이 완료된 랜덤포레스트의 RMSE는 52.942가 나왔다. 튜닝을 하지 않았을때에 비해 RMSE가 2정도 낮춰졌다.

In [50]:

np.sqrt(np.mean((y_test-rft_pred)**2)) #RMSE of tuned RF

Out[50]:

52.94212291306083

3. SVM (RBF kernel)¶

다음은 SVM을 사용해보겠다. SVM에는 결정경계를 정해주는 함수인 커널의 종류를 설정해주어야 하는데, 일반적인 비선형 관계를 잘 잡아내주는 RBF 커널을 사용해보도록 하겠다.

In [51]:

from sklearn.svm import SVR

SVM은 커널을 사용하기 때문에 변수들의 스케일을 맞춰주는것이 좋다. 따라서 먼저 0과 1 사이의 값으로 Normalization을 진행해주도록 하자.

In [52]:

#Implement normalization for svm
Xn_train = X_train-np.min(X_train,axis=0)
Xn_train = Xn_train/(np.max(X_train,axis=0)-np.min(X_train,axis=0))

Xn_test = X_test-np.min(X_train,axis=0)
Xn_test = Xn_test/(np.max(X_train,axis=0)-np.min(X_train,axis=0))

SVM은 rbf커널을 사용할 경우, 정규화 파라미터인 C와 Gamma 값을 지정해주어야 한다.

In [53]:

svm_regr = SVR(kernel='rbf',C=2**-3, gamma=2**3)

SVM은 rbf커널을 사용할 경우, 시간복잡도가 O(n^3)으로 알려져있다. 지금 훈련데이터가 30000개가 넘는것을 감안했을때, 튜닝을 하기 전에 미리 걸리는 시간을 확인해보는게 좋다. 필자의 cpu상으로 svm이 한번 돌아갈때, 대략 20초정도 걸리는것으로 확인되었다.

In [54]:

start=time.time()
svm_regr.fit(Xn_train,y_train)
print(time.time()-start) #one implementation, about 20 seconds... 

19.887809991836548

SVM은 inference단계도 꽤 오래걸린다. 확인해보니 10초정도 소요된다. 따라서 하나의 파라미터 조합 당 대략 30초정도 걸릴것으로 예상할 수 있다.

In [55]:

start=time.time()
svm_pred = svm_regr.predict(Xn_test) 
print(time.time()-start) #inference is also demanding, about 10 seconds...

11.451613187789917

위에서 했던것처럼 pm2.5값이 0보다 작을일은 없으므로, SVM의 예측값이 0보다 작은 경우는 0으로 처리해주자.

In [56]:

svm_pred[svm_pred<0,]=0
svm_pred

Out[56]:

array([55.38548456, 85.88205656, 63.85062897, ..., 77.50121515,
       55.14311824, 53.71417643])

튜닝을 하지 않은 SVM의 RMSE는 92.930으로 확인되었다. 매우 안좋은 수치이다.

In [57]:

np.sqrt(np.mean((y_test-svm_pred)**2)) #RMSE of rough SVM (very bad rmse!)

Out[57]:

92.92957249195852

이제 튜닝을 진행하자. 튜닝은 coarse tuning -> fine tuning으로 진행되었다. 즉, 먼저 넓은 범위에서 튜닝을 진행하고 범위를 좁혀서 점차 세세하게 튜닝이 진행되었다. 처음 튜닝파라미터의 범위는 https://www.csie.ntu.edu.tw/~cjlin/papers/guide/guide.pdf 을 참고하였다. 최종적으로 좁혀진 그리드서치 범위는 아래와 같다.

In [58]:

# 튜닝

expc_list = [6,7,8,9] 
expg_list = [3,4,5,6] 

param_grid = {'C':[2**i for i in expc_list],
             'gamma':[2**i for i in expg_list]} #coarse tuning -> fine tuning

In [59]:

param_grid

Out[59]:

{'C': [64, 128, 256, 512], 'gamma': [8, 16, 32, 64]}

추후 visualization을 위해서 파라미터값을 소수점 3자리로 만든 수치들을 저장해두자.

In [60]:

param_grid_visual = {'C':list(np.round(np.array([2**i for i in expc_list]),3)),
             'gamma':list(np.round(np.array([2**i for i in expg_list]),3))} #save this for visualization later

SVM은 5 fold cross validation하기에 시간이 너무 오래걸리므로 cv를 2로 설정해주었다. 앞서 언급했던것처럼 coarse tuning -> fine tuning이 진행되었고 아래는 fine tuning 과정이다. (coarse tuning과정은 본 노트북에서 생략하였다.)

In [61]:

from sklearn.model_selection import GridSearchCV

grid_search = GridSearchCV(svm_regr,param_grid,cv=2,scoring='neg_root_mean_squared_error',n_jobs=-1)

fine tuning과정에서 대략 13분정도 소요되는 것을 확인할 수 있다.

In [62]:

start=time.time()
grid_search.fit(Xn_train,y_train) #grid search of SVM 
print(time.time()-start) #required time

758.729752779007

In [105]:

758/60 #mins

Out[105]:

12.633333333333333

fine tuning으로 찾아진 SVM의 최적의 파라미터는 아래와 같다.

In [64]:

grid_search.best_params_

Out[64]:

{'C': 512, 'gamma': 16}

In [65]:

import pandas as pd
results = pd.DataFrame(grid_search.cv_results_)
display(results.head(5))

	mean_fit_time	std_fit_time	mean_score_time	std_score_time	param_C	param_gamma	params	split0_test_score	split1_test_score	mean_test_score	std_test_score	rank_test_score
0	10.798589	0.184446	17.107800	0.020478	64	8	{'C': 64, 'gamma': 8}	-65.463382	-68.218947	-66.841165	1.377783	15
1	11.774513	0.142468	17.335311	0.058818	64	16	{'C': 64, 'gamma': 16}	-64.276424	-66.908181	-65.592303	1.315879	12
2	17.736599	0.157609	17.826401	0.057385	64	32	{'C': 64, 'gamma': 32}	-64.104630	-66.133478	-65.119054	1.014424	10
3	34.398991	0.254529	18.215857	0.007710	64	64	{'C': 64, 'gamma': 64}	-64.948414	-67.155527	-66.051971	1.103557	14
4	13.704728	0.015989	18.607772	0.007982	128	8	{'C': 128, 'gamma': 8}	-64.420055	-67.048551	-65.734303	1.314248	13

SVM fine tuning 결과를 아래와 같이 시각화할수 있다. 평가지표가 negative RMSE 이므로, 앞의 부호 minus를 지우면 2 fold cross validation 상의 RMSE라고 볼 수 있다.

In [66]:

import mglearn

%matplotlib inline
import matplotlib.pylab as plt

plt.rcParams["figure.figsize"] = (10,10)
plt.rcParams['lines.linewidth'] = 2
plt.rcParams['lines.color'] = 'r'
plt.rcParams['axes.grid'] = True 

scores=np.array(results.mean_test_score).reshape(4,4)

mglearn.tools.heatmap(scores,xlabel='gamma', xticklabels=param_grid_visual['gamma'], ylabel='C',yticklabels=param_grid_visual['C'],cmap='viridis')

Out[66]:

<matplotlib.collections.PolyCollection at 0x13f671d00>

이제 fine tuning된 SVM의 test set 예측값을 구한다.

In [67]:

svmt_pred = grid_search.predict(Xn_test)

In [68]:

svmt_pred[svmt_pred<0,]=0
svmt_pred

Out[68]:

array([ 30.96605221, 151.90373598,   8.89091142, ..., 144.55076265,
         7.83334625,  24.76907943])

튜닝이 완료된 SVM의 RMSE는 59.851이 나왔다. 튜닝을 하지 않았을때에 비해 RMSE가 30이상 낮춰졌다.

In [69]:

np.sqrt(np.mean((y_test-svmt_pred)**2)) #RMSE of tuned SVM

Out[69]:

59.851333148561515

4. XGBoost¶

In [70]:

#Import libraries:
import pandas as pd
import numpy as np
import xgboost as xgb
from xgboost.sklearn import XGBClassifier
from sklearn.model_selection import train_test_split
from sklearn import metrics   #Additional scklearn functions

import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 12, 4

XGBoost를 돌리기 위해서는 xgboost 패키지가 필요하다. 이번 분석에 사용된 xgboost 패키지의 버전은 1.6.2였다.

In [71]:

print(xgb.__version__)

1.6.2

XGBoost의 튜닝은 다른 모델들에 비해서 특히 까다로운데, 튜닝파라미터의 개수가 훨씬 많기 때문이다. 본 노트북에서 XGBoost의 튜닝 방법은 https://www.analyticsvidhya.com/blog/2016/03/complete-guide-parameter-tuning-xgboost-with-codes-python/ 을 참조하였다. 위 링크에는 분류문제에서 XGBoost를 어떻게 튜닝을 하는지 나와있는데 이를 회귀문제에서 적용할 수 있도록 살짝 변형하였다.

먼저 아래와 같은 함수를 정의해준다. 다른 파라미터들이 정해졌을 때, 최적의 n_estimator를 찾아주는 함수이다. n_estimator는 XGBoost의 성능에 특히 중요하므로, 이를 최적화하는 함수를 이렇게 따로 정의해둔다.

In [72]:

def modelfit(alg, dtrain, predictors,useTrainCV=True, cv_folds=5, early_stopping_rounds=50): #find the best n_estimators
    
    if useTrainCV:
        xgb_param = alg.get_xgb_params()
        xgtrain = xgb.DMatrix(dtrain[predictors].values, label=dtrain[target].values)
        cvresult = xgb.cv(xgb_param, xgtrain, num_boost_round=alg.get_params()['n_estimators'], nfold=cv_folds,
            metrics='rmse', early_stopping_rounds=early_stopping_rounds)
        alg.set_params(n_estimators=cvresult.shape[0])
    
    #print the best n_estimators
    #print(cvresult)
    print("n_estimators:", cvresult.shape[0])

이제 임의의 XGBoost 모델을 만들어보자. 튜닝파라미터는 learning_rate, n_estimators, max_depth, min_child_weight, gamma, subsample, colsample_bytree, reg_alpha(아래에서 튜닝), reg_lambda(아래에서 튜닝) 등이 있으며, 대략적 성능을 보기 위해 먼저 적당히 아무값이나 넣어보았다.

In [73]:

import xgboost as xgb 
from xgboost.sklearn import XGBRegressor
xgb_regr = xgb.XGBRegressor( #just arbitrary xgboost
    learning_rate=0.1,
    n_estimators=1000,
    max_depth=5,
    min_child_weight=1,
    gamma=0,
    subsample=0.8,
    colsample_bytree=0.8,
    objective='reg:squarederror',
    nthread=3,
    scale_pos_weight=1,
    seed=1
    )

In [74]:

#make new format of train set for xgboost
train=df2.iloc[train_idx,]
train.reset_index(drop=True,inplace=True)
print(train.head())

   index     No  year  month  day  hour  pm2.5  DEWP  TEMP    PRES cbwd  \
0    631    632  2010      1   27     7  140.0   -11  -4.0  1015.0   NW   
1  11516  11517  2011      4   25    20   93.0     7  14.0  1003.0   SE   
2  42226  42227  2014     10   26    10   10.0    -8  15.0  1027.0   NE   
3  29825  29826  2013      5   27    17   86.0    15  23.0  1005.0   NW   
4   7882   7883  2010     11   25    10   80.0   -12   0.0  1027.0   NW   

     Iws  Is  Ir  NE  NW  SE  cv  
0  18.79   0   0   0   1   0   0  
1  33.98   0   0   0   0   1   0  
2  28.16   0   0   1   0   0   0  
3   1.79   0   0   0   1   0   0  
4   5.37   0   0   0   1   0   0  

In [75]:

#make new format of test set for xgboost
test=df2.iloc[test_idx,]
test.reset_index(drop=True,inplace=True)
print(test.head())

   index     No  year  month  day  hour  pm2.5  DEWP  TEMP    PRES cbwd  \
0  34717  34718  2013     12   17    13   24.0   -19   1.0  1033.0   NE   
1     33     34  2010      1    2     9  132.0    -7  -5.0  1025.0   SE   
2  34726  34727  2013     12   17    22   13.0   -21  -4.0  1037.0   NE   
3  34729  34730  2013     12   18     1   12.0   -21  -6.0  1037.0   NE   
4     38     39  2010      1    2    14  158.0    -9  -5.0  1025.0   SE   

     Iws  Is  Ir  NE  NW  SE  cv  
0   3.13   0   0   1   0   0   0  
1  14.30   0   0   0   0   1   0  
2  33.08   0   0   1   0   0   0  
3  46.94   0   0   1   0   0   0  
4  31.73   0   0   0   0   1   0  

In [107]:

cols

Out[107]:

['DEWP',
 'TEMP',
 'PRES',
 'Iws',
 'Is',
 'Ir',
 'month',
 'hour',
 'NE',
 'NW',
 'SE',
 'cv']

추후 위에서 정의했던 modelfit 함수에 사용하기 위해 target과 predictors 객체에 반응변수와 설명변수 리스트를 저장해둔다.

In [76]:

target = 'pm2.5'
predictors = cols

In [77]:

xgb_regr.fit(train[predictors],train[target])
xgb_pred=xgb_regr.predict(test[predictors])

튜닝하지 않은 XGBoost의 RMSE는 56.916이다. 꽤 준수한 성능을 보임을 확인할 수 있다.

In [78]:

np.sqrt(np.mean((y_test-xgb_pred)**2)) #not tuned xgboost, already beat SVM.

Out[78]:

56.91627766681695

이제 본격적으로 XGBoost를 튜닝해보도록 하겠다. XGBoost의 튜닝파라미터는 총 10개로 이것 전체를 한꺼번에 그리드서치하기에는 시간이 비현실적으로 오래걸린다. 따라서 실제로 XGBoost를 튜닝할때는 중요한 파라미터부터 하나씩 차례대로 튜닝을 진행한다. XGBoost의 파라미터들 중 n_estimators와 learning_rate가 가장 중요하므로 이를 먼저 튜닝한다.

위에서 임의로 정한 XGBoost 파라미터들 상태에서 최적의 n_estimator를 modelfit 함수를 통해 구하자. 이때, 5-fold cross validation상 최적의 n_estimators 는 767로 밝혀졌다.

In [79]:

# find the best n_estimators given other parameters are fixed
modelfit(xgb_regr,train,predictors) 

n_estimators: 767

XGBoost의 n_estimators가 크면, model fitting이 오래걸릴수 있다. 따라서 learning_rate를 조금 증가시켜서(0.1->0.5) 최적의 n_estimators를 줄여보겠다.

In [80]:

import xgboost as xgb
from xgboost.sklearn import XGBRegressor
xgb1_regr = xgb.XGBRegressor(
    learning_rate=0.5,
    n_estimators=1000,
    max_depth=5,
    min_child_weight=1,
    gamma=0,
    subsample=0.8,
    colsample_bytree=0.8,
    objective='reg:squarederror',
    nthread=3,
    scale_pos_weight=1,
    seed=1
    )

learning_rate를 0.5로 하니 최적의 n_estimators는 100이 되었다. 이정도면 납득할만한 n_estimators라 생각하고 그대로 사용하도록 하겠다.

In [81]:

modelfit(xgb1_regr,train,predictors) #select this n_estimators.

n_estimators: 100

n_estimators가 정해졌으면, 이제 차례대로 파라미터들을 튜닝해주어야 한다. XGBoost를 구성하는 트리에 대한 파라미터인 max_depth와 min_child_weight를 정해보자.

In [82]:

from sklearn.model_selection import GridSearchCV
param_test1={
    'max_depth':range(3,15,2),
    'min_child_weight':range(1,8,2)
}

gsearch1 = GridSearchCV(estimator=XGBRegressor(learning_rate=0.1, n_estimators=100, max_depth=5, min_child_weight=1,
                                              gamma=0, subsample=0.8, colsample_bytree=0.8,
                                              objective='reg:squarederror', scale_pos_weight=1, seed=1),
                       param_grid = param_test1, scoring='neg_root_mean_squared_error', n_jobs=-1, cv=5)

5 fold cross validation상, 최적의 max_depth는 13, min_child_weight는 7로 정해졌다.

In [83]:

gsearch1.fit(train[predictors],train[target])
gsearch1.best_params_, gsearch1.best_score_ 

Out[83]:

({'max_depth': 13, 'min_child_weight': 7}, -56.12428533566238)

다음은 gamma에 대해 튜닝해보자.

In [84]:

param_test2={
    'gamma':[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7]
}

gsearch2 = GridSearchCV(estimator=XGBRegressor(learning_rate=0.1, n_estimators=100, max_depth=13, min_child_weight=7,
                                              gamma=0, subsample=0.8, colsample_bytree=0.8,
                                              objective='reg:squarederror', scale_pos_weight=1, seed=1),
                       param_grid = param_test2, scoring='neg_root_mean_squared_error', n_jobs=3, cv=5)

gamma는 원래 값인 0이 최적의 파라미터로 선정되었다. 따라서 바로 전 과정의 validation 스코어와 동일하다.

In [85]:

gsearch2.fit(train[predictors],train[target])
gsearch2.best_params_, gsearch2.best_score_

Out[85]:

({'gamma': 0}, -56.12428533566238)

지금까지 파라미터 max_depth, min_child_weight, gamma 를 튜닝하였다. 다른 파라미터들을 튜닝하기에 앞서, 지금 상태의 max_depth, min_child_weight, gamma 상태에서 최적의 n_estimators를 다시 찾는다. 그렇게 찾아진 최적의 n_estimators는 73이다.

In [86]:

xgb2_regr = xgb.XGBRegressor(
    learning_rate=0.1,
    n_estimators=1000,
    max_depth=13,
    min_child_weight=7,
    gamma=0,
    subsample=0.8,
    colsample_bytree=0.8,
    objective='reg:squarederror',
    nthread=3,
    scale_pos_weight=1,
    seed=1
    )

modelfit(xgb2_regr, train, predictors)

n_estimators: 73

위에서 정해진 n_estimators, max_depth, min_child_weight, gamma를 고정시킨뒤 subsample, colsample_bytree 파라미터를 튜닝해보자.

In [87]:

param_test3 = {
    'subsample': [0.5,0.6,0.7,0.8,0.9],  #[i/10.0 for i in range(6,10)],
    'colsample_bytree': [0.5,0.6,0.7,0.8,0.9]  #[1/10.0 for i in range(6,10)]
}

gsearch3 = GridSearchCV(estimator=XGBRegressor(learning_rate=0.1, n_estimators=73, max_depth=13, min_child_weight=7,
                                              gamma=0, subsample=0.8, colsample_bytree=0.8,
                                              objective='reg:squarederror', scale_pos_weight=1, seed=1),
                       param_grid = param_test3, scoring='neg_root_mean_squared_error', n_jobs=3, cv=5)

이때 찾아진 최적의 파라미터는 0.9로 밝혀졌다.

In [88]:

gsearch3.fit(train[predictors],train[target])
gsearch3.best_params_, gsearch3.best_score_

Out[88]:

({'colsample_bytree': 0.9, 'subsample': 0.9}, -55.46338744252314)

지금까지 구한 파라미터들(n_estimators, max_depth, min_child_weight, gamma, subsample, colsample_bytree)를 고정시킨뒤, 정규화 파라미터인 reg_alpha, reg_lambda를 튜닝한다. 처음에는 범위를 넓게 잡는 coarse tuning을 진행한다.

In [89]:

param_test4 = {
    'reg_alpha': [0,10**-5,10**-2,0.1,1,100],
    'reg_lambda': [0,10**-5,10**-2,0.1,1,100]  
}

gsearch4 = GridSearchCV(estimator=XGBRegressor(learning_rate=0.1, n_estimators=73, max_depth=13, min_child_weight=7,
                                              gamma=0, subsample=0.9, colsample_bytree=0.9,
                                              objective='reg:squarederror', scale_pos_weight=1, seed=1),
                       param_grid = param_test4, scoring='neg_root_mean_squared_error', n_jobs=3, cv=5)

그렇게 찾아진 reg_alpha는 100, reg_lambda는 0.1로 밝혀졌다.

In [90]:

gsearch4.fit(train[predictors],train[target])
gsearch4.best_params_, gsearch4.best_score_

Out[90]:

({'reg_alpha': 100, 'reg_lambda': 0.1}, -55.46229828534162)

이제 reg_alpha와 reg_lambda 값을 위에서 찾아진 값들 주변을 더 세세하게 찾아보자.

In [91]:

param_test4 = {
    'reg_alpha': [10,20,50,100],
    'reg_lambda': [0.05,0.1,0.2,0.5]  
}

gsearch4 = GridSearchCV(estimator=XGBRegressor(learning_rate=0.1, n_estimators=73, max_depth=13, min_child_weight=7,
                                              gamma=0, subsample=0.9, colsample_bytree=0.9,
                                              objective='reg:squarederror', scale_pos_weight=1, seed=1),
                       param_grid = param_test4, scoring='neg_root_mean_squared_error', n_jobs=3, cv=5)

그렇게 찾아진 reg_alpha는 50, reg_lambda는 0.05 이다.

In [92]:

gsearch4.fit(train[predictors],train[target])
gsearch4.best_params_, gsearch4.best_score_

Out[92]:

({'reg_alpha': 50, 'reg_lambda': 0.05}, -55.38576886555485)

이제 마지막으로 위에서 정해진 파라미터들을 고정시키고, n_estimators를 다시 튜닝하여준다. 이렇게 찾아진 n_estimators는 77로 밝혀졌다.

In [93]:

xgb3_regr = xgb.XGBRegressor(
    learning_rate=0.1,
    n_estimators=1000,
    max_depth=13,
    min_child_weight=7,
    gamma=0,
    subsample=0.9,
    colsample_bytree=0.9,
    objective='reg:squarederror',
    nthread=3,
    scale_pos_weight=1,
    reg_alpha=50,
    reg_lambda=0.05,
    seed=1
    )

modelfit(xgb3_regr, train, predictors)

n_estimators: 77

혹시 learning_rate를 낮추면 성능이 더 좋아질수도 있다. 따라서 learning_rate를 0.01로 줄여서, 그에 최적이 되는 n_estimators를 찾아보자.

In [94]:

xgb4_regr = xgb.XGBRegressor(
    learning_rate=0.01,
    n_estimators=1000,
    max_depth=13,
    min_child_weight=7,
    gamma=0,
    subsample=0.9,
    colsample_bytree=0.9,
    objective='reg:squarederror',
    nthread=3,
    scale_pos_weight=1,
    reg_alpha=50,
    reg_lambda=0.05,
    seed=1
    )

modelfit(xgb4_regr, train, predictors)

n_estimators: 890

위에서 정한 learning_rate가 0.1일때와 0.01일때중 어느것이 더 좋은지 5 fold cross validation을 통해 결정한다.

In [95]:

param_test5 = {
    'learning_rate': [0.1,0.01],  #[i/10.0 for i in range(6,10)],
    'n_estimators': [77,890]  #[1/10.0 for i in range(6,10)]
}

gsearch5 = GridSearchCV(estimator=XGBRegressor(learning_rate=0.1, n_estimators=5, max_depth=13, min_child_weight=7,
                                              gamma=0, subsample=0.9, colsample_bytree=0.9,
                                              objective='reg:squarederror', scale_pos_weight=1,reg_alpha=50,
                                               reg_lambda=0.05, seed=1),
                       param_grid = param_test5, scoring='neg_root_mean_squared_error', n_jobs=3, cv=5)

그렇게 5 fold cross validation을 진행한 결과 learning_rate값이 0.01, n_estimators값이 890일때, 최적의 모델이 됨을 확인할 수 있었다.

In [96]:

gsearch5.fit(train[predictors],train[target])
gsearch5.best_params_, gsearch5.best_score_

/Users/minsoo/torch_ground/lib/python3.9/site-packages/joblib/externals/loky/process_executor.py:702: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.
  warnings.warn(

Out[96]:

({'learning_rate': 0.01, 'n_estimators': 890}, -54.600899387228786)

지금까지 XGBoost를 튜닝한 값들을 대입하여 모델링한다.

In [97]:

xgbt_regr = xgb.XGBRegressor(
    learning_rate=0.01,
    n_estimators=890,
    max_depth=13,
    min_child_weight=7,
    gamma=0,
    subsample=0.9,
    colsample_bytree=0.9,
    objective='reg:squarederror',
    nthread=3,
    scale_pos_weight=1,
    reg_alpha=50,
    reg_lambda=0.05,
    seed=1
    )

튜닝된 XGBoost를 통해 test_set을 예측해본다.

In [98]:

xgbt_regr.fit(train[predictors],train[target])
xgbt_pred=xgbt_regr.predict(test[predictors])

In [99]:

xgbt_pred

Out[99]:

array([ 41.60661 , 198.37718 ,  13.011922, ..., 204.16077 ,  47.88999 ,
        26.341387], dtype=float32)

튜닝된 XGBoost의 RMSE는 52.887로 튜닝된 Random Forest와 SVM에 비해 준수하게 나왔다.

In [100]:

np.sqrt(np.mean((y_test-xgbt_pred)**2)) #tuned xgboost. it outperforms tuned RandomForest and SVM.

Out[100]:

52.88686343503121

test_set의 실제 PM2.5 농도값과 XGBoost로 예측된 PM2.5 농도값을 처음 5개 비교해보면 아래와 같다.

In [114]:

y_test[1:6]

Out[114]:

array([132.,  13.,  12., 158., 154.])

In [115]:

xgbt_pred[1:6]

Out[115]:

array([198.37718 ,  13.011922,  13.743099, 123.998924, 134.80574 ],
      dtype=float32)

오늘은 Random Forest, Support Vector Machine, XGBoost 모델의 회귀 문제에서 튜닝 방법에 대해 알아보았다. 머신러닝에서 튜닝은 정말로 중요하며, 튜닝을 통해서 모델의 예측 성능이 크게 개선될 수 있음을 확인할 수 있었다. 머신러닝에서 튜닝 방법은 상당히 휴리스틱한 면이 강하며, 본 노트북에 제시된 튜닝 방법이 최선이라고 말하진 못하겠다. 다만, 필자가 실무에서 실제로 사용하는 튜닝 방법에 대해서 포스팅해보았다. (더 좋은 방법이 있다면 공유해주시면 감사합니다.) 다음번엔 분류문제에서 Random Forest, Support Vector Machine, XGBoost 모델을 어떻게 튜닝할 수 있는지 포스팅해보도록 하겠다.