Kaggle チャレンジ 4日目 住宅価格問題を解いていく

シェアする

  • このエントリーをはてなブックマークに追加

スポンサーリンク

はじめに

 今回は、住宅の価格を色々な特徴量から予測していきたいと思います。今回もGoogleColaboratoryを使って進めていくので、はじめ方などは前回の記事を参考にしてください。

前回の記事:

Kaggle チャレンジ 1日目 タイタニックの問題からデータを読み解いてみる


やること

 GoogleColaboratoryとKernelsを使って住宅の価格を色々な特徴量から予測していきたいと思います。

・今回参考にするKernels: Stacked Regressions : Top 4% on LeaderBoard

GoogleColaboratoryの使い方はこちら

対象者

機械学習をKaggleを使って学びたい方、Kaggleに興味があるけどやり方がわからない方。

 

House Prices: Advanced Regression Techniques

Kernelsに沿って進めていきます。ところどころ日本語がおかしいところがあるかもしれませんが、よろしくお願いします。また、全てを翻訳せず、抜粋して取り上げていきます。

このコンペのKernelsで私はいくつかの素晴らしい記事を読みました。それが、次に挙げているものになります。

  1. Comprehensive data exploration with Python by Pedro Marcelino : すごく意欲的なデータ分析
  2. A study on Regression applied to the Ames dataset by Julien Cohen-Solal : 徹底した解析と線形回帰分析によって深くまで掘り下げているので、初心者はついていきやすいでしょう。
  3. Regularized Linear Models by Alexandru Papiu : モデリングとクロスバリデーションに関する素晴らしい入門記事

  1. 欠損値の補完
  2. 変換
  3. ラベルエンコード
  4. Box Cox 変換
  5. ダミー変数の取得

次に、多くのベースモデル(ほとんどSklearnベースのモデル+ DMLCのXGBoostとMicrosoftのLightGBMのSklearn API)を選択し、それらをスタッキング/アンサンブルする前にデータでクロスバリデーションする。ここでの重要な点は、(線形)モデルを外れ値に強くすることです。これにより、LBとクロスバリデーションの両方で結果が改善しました。

ライブラリの準備

 必要になるライブラリを用意しましょう。

In[1]:

#import some necessary librairies

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
%matplotlib inline
import matplotlib.pyplot as plt  # Matlab-style plotting
import seaborn as sns
color = sns.color_palette()
sns.set_style('darkgrid')
import warnings
def ignore_warn(*args, **kwargs):
    pass
warnings.warn = ignore_warn #ignore annoying warning (from sklearn and seaborn)


from scipy import stats
from scipy.stats import norm, skew #for some statistics


pd.set_option('display.float_format', lambda x: '{:.3f}'.format(x)) #Limiting floats output to 3 decimal points


from subprocess import check_output
print(check_output(["ls", "../input"]).decode("utf8")) #check the files available in the directory

データの準備

前回同様、ドライブからデータを引っ張ってきます。わからない方は、前回の記事を参考にしてください。

In[2]:

from pydrive.auth import GoogleAuth
from pydrive.drive import GoogleDrive
from google.colab import auth
from oauth2client.client import GoogleCredentials
auth.authenticate_user()
gauth = GoogleAuth()
gauth.credentials = GoogleCredentials.get_application_default()
drive = GoogleDrive(gauth)
id = '*************************' # 共有リンクで取得した id= より後の部分を*の部分に入力
downloaded = drive.CreateFile({'id': id})
downloaded.GetContentFile('train.csv') #ファイルの名前
id = '*************************' 
downloaded = drive.CreateFile({'id': id})
downloaded.GetContentFile('test.csv')
id = '*************************'
downloaded = drive.CreateFile({'id': id})
downloaded.GetContentFile('sample_submission.csv')

In[3]:

train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')

 データを用意できたら前回同様、最初に出力させてみてみましょう。このページでは、全て表示すると見切れてしまうので中間のカラムを省略しています。

In[4]:

##display the first five rows of the train dataset.
train.head(5)

Out[4]:

Id MSSubClass MSZoning LotFrontage SaleType SaleCondition SalePrice
0 1 60 RL 65.000 WD Normal 208500
1 2 20 RL 80.000 WD Normal 181500
2 3 60 RL 68.000 WD Normal 223500
3 4 70 RL 60.000 WD Abnorml 140000
4 5 60 RL 84.000 WD Normal 250000

5 rows × 81 columns

In[5]:

##display the first five rows of the test dataset.
test.head(5)
Out[5]:
Id MSSubClass MSZoning LotFrontage YrSold SaleType SaleCondition
0 1461 20 RH 80.000 2010 WD Normal
1 1462 20 RL 81.000 2010 WD Normal
2 1463 60 RL 74.000 2010 WD Normal
3 1464 60 RL 78.000 2010 WD Normal
4 1465 120 RL 43.000 2010 WD Normal

5 rows × 80 columns

カラムを見て気づくと思いますが、trainの方にのみ、SalePrice があります。これを予測していきます。

IDは予測に必要ないので削除します。削除する前と後のそれぞれデータセットのサイズを出力させています。

In[6]:

#check the numbers of samples and features
print("The train data size before dropping Id feature is : {} ".format(train.shape))
print("The test data size before dropping Id feature is : {} ".format(test.shape))

#Save the 'Id' column
train_ID = train['Id']
test_ID = test['Id']

#Now drop the  'Id' colum since it's unnecessary for  the prediction process.
train.drop("Id", axis = 1, inplace = True)
test.drop("Id", axis = 1, inplace = True)

#check again the data size after dropping the 'Id' variable
print("\nThe train data size after dropping Id feature is : {} ".format(train.shape)) 
print("The test data size after dropping Id feature is : {} ".format(test.shape))
Out[6]:
The train data size before dropping Id feature is : (1460, 81) 
The test data size before dropping Id feature is : (1459, 80) 

The train data size after dropping Id feature is : (1460, 80) 
The test data size after dropping Id feature is : (1459, 79)

外れ値

 エイムズの住宅データの文書は、トレーニングデータに外れ値が存在することを示しているので、これらの外れ値を調べましょう。

In[7]:

fig, ax = plt.subplots()
ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()
Out[7]:
 右下にはGrLivAreaが大きく、低価格の2つの大きな外れ値があります。 よってこの外れ値を削除します。
In[8]:
#Deleting outliers
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)

#Check the graphic again
fig, ax = plt.subplots()
ax.scatter(train['GrLivArea'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

Out[8]:

注意:
今回は非常に悪いデータだった(非常に低価格で非常に広い領域)ので、これらの2つを削除することに決めました。トレーニングデータには異常値が存在する可能性があります。 しかし、テストデータに外れ値があった場合、外れ値をすべて削除するとモデルが汎用性を失う可能性があります。

ターゲット変数

今回のターゲット変数はSalePriceです。なので、この変数について分析していきましょう。

In[9]:

sns.distplot(train['SalePrice'] , fit=norm);

# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

#Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('SalePrice distribution')

#Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)
plt.show()
Out[9]:
 mu = 180932.92 and sigma = 79467.79

 目標変数は右に偏っていますね。(線形)モデルは、正規分布したデータを愛用するので、この変数を変換してより正規分布にしましょう。

ターゲット変数のログ変換

In[10]:

#We use the numpy fuction log1p which  applies log(1+x) to all elements of the column
train["SalePrice"] = np.log1p(train["SalePrice"])

#Check the new distribution 
sns.distplot(train['SalePrice'] , fit=norm);

# Get the fitted parameters used by the function
(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

#Now plot the distribution
plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],
            loc='best')
plt.ylabel('Frequency')
plt.title('SalePrice distribution')

#Get also the QQ-plot
fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)
plt.show()
Out[10]:
 mu = 12.02 and sigma = 0.40

 データはより正規的に配置されているように見えます。

engineering 特徴

 まずトレーニングとテストのデータを同じデータフレームに連結しましょう。

In[11]:

ntrain = train.shape[0]
ntest = test.shape[0]
y_train = train.SalePrice.values
all_data = pd.concat((train, test)).reset_index(drop=True)
all_data.drop(['SalePrice'], axis=1, inplace=True)
print("all_data size is : {}".format(all_data.shape))
Out[11]:
all_data size is : (2917, 79)

欠損値

In[12]:

all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(20)
Out[12]:
Missing Ratio
PoolQC 99.691
MiscFeature 96.400
Alley 93.212
Fence 80.425
FireplaceQu 48.680
LotFrontage 16.661
GarageQual 5.451
GarageCond 5.451
GarageFinish 5.451
GarageYrBlt 5.451
GarageType 5.382
BsmtExposure 2.811
BsmtCond 2.811
BsmtQual 2.777
BsmtFinType2 2.743
BsmtFinType1 2.708
MasVnrType 0.823
MasVnrArea 0.788
MSZoning 0.137
BsmtFullBath 0.069

 表ではわかりにくいので、可視化してみましょう。

In[13]:

f, ax = plt.subplots(figsize=(15, 12))
plt.xticks(rotation='90')
sns.barplot(x=all_data_na.index, y=all_data_na)
plt.xlabel('Features', fontsize=15)
plt.ylabel('Percent of missing values', fontsize=15)
plt.title('Percent missing data by feature', fontsize=15)
Out[13]:
Text(0.5,1,'Percent missing data by feature')

データの相関

In[14]:

#Correlation map to see how features are correlated with SalePrice
corrmat = train.corr()
plt.subplots(figsize=(12,9))
sns.heatmap(corrmat, vmax=0.9, square=True)
Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x7efd7b454898>

欠損値の入力

  • PoolQC:NAは「プールなし」を意味します。 欠落値(+ 99%)は、一般的に普通の家はプールを持たないことを考えると、理にかなっています。

In[15]:

all_data["PoolQC"] = all_data["PoolQC"].fillna("None")

  • MiscFeature:NAと表示されます。「その他の機能はありません。」

In[16]:

all_data["MiscFeature"] = all_data["MiscFeature"].fillna("None")

  • Alley:NAは “路地なし”を意味します。

In[17]:

all_data["Alley"] = all_data["Alley"].fillna("None")
  • Fence:NAは「フェンスなし」を意味します。

In[18]:

all_data["Fence"] = all_data["Fence"].fillna("None")

  • FireplaceQu:NAが “暖炉なし”を意味します。

In[19]:

all_data["FireplaceQu"] = all_data["FireplaceQu"].fillna("None")

  • LotFrontage:家に接続されている道のエリアは、近隣の他の家と同様のエリアを持つ可能性が最も高いので、不足している値は近隣のLotFrontageの中央値で埋めることができます。

In[20]:

#Group by neighborhood and fill in missing value by the median LotFrontage of all the neighborhood
all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(
    lambda x: x.fillna(x.median()))

  • GarageYrBlt、GarageArea、GarageCars:不足しているデータを0に置き換える(ガレージなし=このようなガレージには車がないため)

In[21]:

for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
    all_data[col] = all_data[col].fillna(0)

  • BsmtFinSF1、BsmtFinSF2、BsmtUnfSF、TotalBsmtSF、BsmtFullBath、BsmtHalfBath:欠損値は、地下室がないから、ゼロである可能性が高い。

In[22]:

for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
    all_data[col] = all_data[col].fillna(0)

  • BsmtQual、BsmtCond、BsmtExposure、BsmtFinType1、BsmtFinType2:これらのカテゴリの全ての地下室関連の特徴量のNaNは、地下室がないことを意味します。

In[23]:

for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
    all_data[col] = all_data[col].fillna('None')

  • MasVnrAreaとMasVnrType: NA は、これらの家屋のための石工単板を意味しない可能性が最も高いです。 ここには0を、型にはNoneを指定できます。

In[24]:

all_data["MasVnrType"] = all_data["MasVnrType"].fillna("None")
all_data["MasVnrArea"] = all_data["MasVnrArea"].fillna(0)

  • MSZoning(一般的なゾーニングの分類): ‘RL’は最頻値です。 そこで、欠損値を ‘RL’に置き換えることができます。

In[25]:

all_data['MSZoning'] = all_data['MSZoning'].fillna(all_data['MSZoning'].mode()[0])

  • Utilities:このカテゴリ機能の場合、「NoSeWa」と「NA」を除き、すべてのレコードは「AllPub」です。 「NoSewa」の家がトレーニングセットに含まれているので、この機能は予測モデリングに役立ちません。 

In[26]:

all_data = all_data.drop(['Utilities'], axis=1)

  • Functional :NAは典型的であることを示します。

In[27]:

all_data["Functional"] = all_data["Functional"].fillna("Typ")

  • Electrical :NAが1つあります。 この特徴はの最頻値は’SBrkr’なので、欠損値を最頻値に置き換えましょう。

In[28]:

all_data['Electrical'] = all_data['Electrical'].fillna(all_data['Electrical'].mode()[0])

  • KitchenQual:KitchenQualの欠損値のために、1つのNA値とElectricalと同じだけ、私たちは ‘TA’(最頻値)に置き換えます。

In[29]:

all_data['KitchenQual'] = all_data['KitchenQual'].fillna(all_data['KitchenQual'].mode()[0])

  • Exterior1stとExterior2nd:再びExterior 1と2の両方に欠損値が1つしかありません。 私たちは、最頻値に置き換えます。

In[30]:

all_data['Exterior1st'] = all_data['Exterior1st'].fillna(all_data['Exterior1st'].mode()[0])
all_data['Exterior2nd'] = all_data['Exterior2nd'].fillna(all_data['Exterior2nd'].mode()[0])

  • SaleType:最頻値で欠損値を置き換えましょう。

In[31]:

all_data['SaleType'] = all_data['SaleType'].fillna(all_data['SaleType'].mode()[0])

  • MSSubClass:NAは建物クラスなしを意味します。 欠損値をNoneに置き換えることができます。

In[32]:

all_data['MSSubClass'] = all_data['MSSubClass'].fillna("None")

  • 欠損値が他にあるか確かめてみましょう。

In[33]:

#Check remaining missing values if any 
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head()
Missing Ratio

欠損値はありません。

  • 型の変換

In[34]:

#MSSubClass=The building class
all_data['MSSubClass'] = all_data['MSSubClass'].apply(str)


#Changing OverallCond into a categorical variable
all_data['OverallCond'] = all_data['OverallCond'].astype(str)


#Year and month sold are transformed into categorical features.
all_data['YrSold'] = all_data['YrSold'].astype(str)
all_data['MoSold'] = all_data['MoSold'].astype(str)

  • Label Encoding

In[35]:

from sklearn.preprocessing import LabelEncoder
cols = ('FireplaceQu', 'BsmtQual', 'BsmtCond', 'GarageQual', 'GarageCond', 
        'ExterQual', 'ExterCond','HeatingQC', 'PoolQC', 'KitchenQual', 'BsmtFinType1', 
        'BsmtFinType2', 'Functional', 'Fence', 'BsmtExposure', 'GarageFinish', 'LandSlope',
        'LotShape', 'PavedDrive', 'Street', 'Alley', 'CentralAir', 'MSSubClass', 'OverallCond', 
        'YrSold', 'MoSold')
# process columns, apply LabelEncoder to categorical features
for c in cols:
    lbl = LabelEncoder() 
    lbl.fit(list(all_data[c].values)) 
    all_data[c] = lbl.transform(list(all_data[c].values))

# shape        
print('Shape all_data: {}'.format(all_data.shape))
Shape all_data: (2917, 78)

  • 重要な特徴を1つ追加する

面積関連の特徴は住宅価格を決定するために非常に重要であるため、各住宅の地下1階と2階の合計面積を追加する

In[36]:

# Adding total sqfootage feature 
all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']

  • Skewed features

In[37]:

numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index

# Check the skew of all numerical features
skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
print("\nSkew in numerical features: \n")
skewness = pd.DataFrame({'Skew' :skewed_feats})
skewness.head(10)
Skew in numerical features:

Skew
MiscVal 21.940
PoolArea 17.689
LotArea 13.109
LowQualFinSF 12.085
3SsnPorch 11.372
LandSlope 4.973
KitchenAbvGr 4.301
BsmtFinSF2 4.145
EnclosedPorch 4.002
ScreenPorch 3.945

  • Box Cox Skewed featuresの変換

In[38]:

skewness = skewness[abs(skewness) > 0.75]
print("There are {} skewed numerical features to Box Cox transform".format(skewness.shape[0]))

from scipy.special import boxcox1p
skewed_features = skewness.index
lam = 0.15
for feat in skewed_features:
    #all_data[feat] += 1
    all_data[feat] = boxcox1p(all_data[feat], lam)
    
#all_data[skewed_features] = np.log1p(all_data[skewed_features])
There are 59 skewed numerical features to Box Cox transform

  • ダミー変数の取得

In[39]:

all_data = pd.get_dummies(all_data)
print(all_data.shape)
(2917, 220)

  • 新しいトレーニングとテストのデータを取得する

In[40]:

train = all_data[:ntrain]
test = all_data[ntrain:]

モデリング

ライブラリの準備

In[41]:

from sklearn.linear_model import ElasticNet, Lasso,  BayesianRidge, LassoLarsIC
from sklearn.ensemble import RandomForestRegressor,  GradientBoostingRegressor
from sklearn.kernel_ridge import KernelRidge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import RobustScaler
from sklearn.base import BaseEstimator, TransformerMixin, RegressorMixin, clone
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.metrics import mean_squared_error
import xgboost as xgb
import lightgbm as lgb

クロスバリデーションを定義する

 Sklearnのcross_val_score関数を使用します。クロスバリデーションの前にデータセットをシャッフルするために、1行のコードを追加します。

In[42]:

#Validation function
n_folds = 5

def rmsle_cv(model):
    kf = KFold(n_folds, shuffle=True, random_state=42).get_n_splits(train.values)
    rmse= np.sqrt(-cross_val_score(model, train.values, y_train, scoring="neg_mean_squared_error", cv = kf))
    return(rmse)

Base models and scores

  • LASSO Regression 

In[43]:

lasso = make_pipeline(RobustScaler(), Lasso(alpha =0.0005, random_state=1))
score = rmsle_cv(lasso)
print("\nLasso score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
Lasso score: 0.1115 (0.0074)

  • Elastic Net Regression 

In[44]:

ENet = make_pipeline(RobustScaler(), ElasticNet(alpha=0.0005, l1_ratio=.9, random_state=3))
score = rmsle_cv(ENet)
print("ElasticNet score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
ElasticNet score: 0.1116 (0.0074)

  • Kernel Ridge Regression 

In[45]:

KRR = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)
score = rmsle_cv(KRR)
print("Kernel Ridge score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
Kernel Ridge score: 0.1153 (0.0075)

  • Gradient Boosting Regression 

In[46]:

GBoost = GradientBoostingRegressor(n_estimators=3000, learning_rate=0.05,
max_depth=4, max_features='sqrt',
min_samples_leaf=15, min_samples_split=10, 
loss='huber', random_state =5)
score = rmsle_cv(GBoost)
print("Gradient Boosting score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
Gradient Boosting score: 0.1177 (0.0080)

  • XGBoost 

In[47]:

model_xgb = xgb.XGBRegressor(colsample_bytree=0.4603, gamma=0.0468, 
learning_rate=0.05, max_depth=3, 
min_child_weight=1.7817, n_estimators=2200,
reg_alpha=0.4640, reg_lambda=0.8571,
subsample=0.5213, silent=1,
random_state =7, nthread = -1)
score = rmsle_cv(model_xgb)
print("Xgboost score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
Xgboost score: 0.1164 (0.0069)

  • LightGBM

In[48]:

model_lgb = lgb.LGBMRegressor(objective='regression',num_leaves=5,
learning_rate=0.05, n_estimators=720,
max_bin = 55, bagging_fraction = 0.8,
bagging_freq = 5, feature_fraction = 0.2319,
feature_fraction_seed=9, bagging_seed=9,
min_data_in_leaf =6, min_sum_hessian_in_leaf = 11)
score = rmsle_cv(model_lgb)
print("LGBM score: {:.4f} ({:.4f})\n" .format(score.mean(), score.std()))
LGBM score: 0.1167 (0.0072)

Stacking models

 最も単純なstackのアプローチとして基本モデルの平均化があります。
ベースモデルを平均化するこの単純なアプローチから始めます。 、

  • 平均化されたベースモデルクラス

In[49]:

class AveragingModels(BaseEstimator, RegressorMixin, TransformerMixin):
    def __init__(self, models):
        self.models = models
        
    # we define clones of the original models to fit the data in
    def fit(self, X, y):
        self.models_ = [clone(x) for x in self.models]
        
        # Train cloned base models
        for model in self.models_:
            model.fit(X, y)

        return self
    
    #Now we do the predictions for cloned models and average them
    def predict(self, X):
        predictions = np.column_stack([
            model.predict(X) for model in self.models_
        ])
        return np.mean(predictions, axis=1)

  • 平均ベースモデルスコア

ENet、GBoost、KRR、そしてlassoの4つのモデルをここで平均します。 もちろん、もっと多くのモデルを簡単に追加することができます。

In[50]:

averaged_models = AveragingModels(models = (ENet, GBoost, KRR, lasso))

score = rmsle_cv(averaged_models)
print(" Averaged base models score: {:.4f} ({:.4f})\n".format(score.mean(), score.std()))
 Averaged base models score: 0.1091 (0.0075)

最も単純なスタッキングアプローチでさえ、本当にスコアを向上させるようです。他のものも試してみましょう。

メタモデルの追加

平均ベースモデルにメタモデルを追加し、これらのベースモデルの out-of-folds 予測を使用して、メタモデルをトレーニングします。

トレーニング部分の手順は、次のように記述できます。

  1. トータルトレーニングセットを2つの互いに素なセットに分割します(ここではtrainと holdout )
  2. 最初の部分(train)にいくつかのベースモデルを学習する
  3. これらのベースモデルを2番目の部分(ホールドアウト)でテストします。
  4. 3)からの予測を入力として使用し、メタモデルと呼ばれるより高いレベルの学習モデルを学習するための出力としての正しい応答(目標変数)を使用します。

最初の3つのステップは繰り返し実行されます。例えば、5つのスタッキングを取る場合、最初にトレーニングデータを5つに分割します。データを5個に分割した上で、そのうちの1つをテストデータとし、残る4個を学習用のデータとするアプローチです。

 予測部分では、テストデータ上のすべてのベースモデルの予測を平均し、それらをメタモデルとして使用して、最終的な予測をメタモデルで行います。

Stacking averaged Models Class

In[50]:

class StackingAveragedModels(BaseEstimator, RegressorMixin, TransformerMixin):
    def __init__(self, base_models, meta_model, n_folds=5):
        self.base_models = base_models
        self.meta_model = meta_model
        self.n_folds = n_folds
   
    # We again fit the data on clones of the original models
    def fit(self, X, y):
        self.base_models_ = [list() for x in self.base_models]
        self.meta_model_ = clone(self.meta_model)
        kfold = KFold(n_splits=self.n_folds, shuffle=True, random_state=156)
        
        # Train cloned base models then create out-of-fold predictions
        # that are needed to train the cloned meta-model
        out_of_fold_predictions = np.zeros((X.shape[0], len(self.base_models)))
        for i, model in enumerate(self.base_models):
            for train_index, holdout_index in kfold.split(X, y):
                instance = clone(model)
                self.base_models_[i].append(instance)
                instance.fit(X[train_index], y[train_index])
                y_pred = instance.predict(X[holdout_index])
                out_of_fold_predictions[holdout_index, i] = y_pred
                
        # Now train the cloned  meta-model using the out-of-fold predictions as new feature
        self.meta_model_.fit(out_of_fold_predictions, y)
        return self
   
    #Do the predictions of all base models on the test data and use the averaged predictions as 
    #meta-features for the final prediction which is done by the meta-model
    def predict(self, X):
        meta_features = np.column_stack([
            np.column_stack([model.predict(X) for model in base_models]).mean(axis=1)
            for base_models in self.base_models_ ])
        return self.meta_model_.predict(meta_features)

Stacking Averaged models Score

 同じ数のモデルを使用して2つのアプローチを比較可能にするために、我々はEnet KRRとGboostの平均値を取った後、lassoをメタモデルとして追加します。
In[51]:
stacked_averaged_models = StackingAveragedModels(base_models = (ENet, GBoost, KRR),
                                                 meta_model = lasso)

score = rmsle_cv(stacked_averaged_models)
print("Stacking Averaged models score: {:.4f} ({:.4f})".format(score.mean(), score.std()))
Stacking Averaged models score: 0.1085 (0.0074)

Ensembling StackedRegressor, XGBoost and LightGBM

 前に定義した StackedRegressor に XGBoost と LightGBM を追加します。

In[52]:

def rmsle(y, y_pred):
    return np.sqrt(mean_squared_error(y, y_pred))

Final Training and Prediction

  • StackedRegressor

In[53]:

stacked_averaged_models.fit(train.values, y_train)
stacked_train_pred = stacked_averaged_models.predict(train.values)
stacked_pred = np.expm1(stacked_averaged_models.predict(test.values))
print(rmsle(y_train, stacked_train_pred))
0.0781571937916
  • XGBoost

In[54]:

model_xgb.fit(train, y_train)
xgb_train_pred = model_xgb.predict(train)
xgb_pred = np.expm1(model_xgb.predict(test))
print(rmsle(y_train, xgb_train_pred))
0.0785165142425
  • LightGBM

In[55]:

model_lgb.fit(train, y_train)
lgb_train_pred = model_lgb.predict(train)
lgb_pred = np.expm1(model_lgb.predict(test.values))
print(rmsle(y_train, lgb_train_pred))
0.0716757468834
In[56]:
'''RMSE on the entire Train data when averaging'''

print('RMSLE score on train data:')
print(rmsle(y_train,stacked_train_pred*0.70 +
               xgb_train_pred*0.15 + lgb_train_pred*0.15 ))
RMSLE score on train data:
0.0752190464543

  • Ensemble prediction

In[57]:

ensemble = stacked_pred*0.70 + xgb_pred*0.15 + lgb_pred*0.15
  • Submission

In[58]:

sub = pd.DataFrame()
sub['Id'] = test_ID
sub['SalePrice'] = ensemble
sub.to_csv('submission.csv',index=False)

最後に、出力させたcsvファイルをローカルに落として提出しましょう。

from google.colab import files
files.download('submission.csv')

今回は、住宅の価格を色々なモデルで予測し、最後にスタッキングについても触れました。スタッキングはスコアを伸ばすには必要になってくるので色々変えて試してみてください。

最後まで読んでいただきありがとうございました。よろしければこの記事をシェアしていただけると励みになります。よろしくお願いします。

スポンサーリンク
レクタングル広告(大)
レクタングル広告(大)

シェアする

  • このエントリーをはてなブックマークに追加

フォローする

スポンサーリンク
レクタングル広告(大)