Kaggle チャレンジ

正月なまりを吹っ飛ばすべく、今回はKaggleチャレンジをします。

Kaggleとは

世界中のデータサイエンティスト達がデータ分析の腕を競い合うプラットフォームです。企業や研究者が提供するデータセットに対し、適切な予測モデルを構築しその予測精度を競います。上位入賞者には賞金が出ます。
様々なデータセットが手に入り、他の競技者の解説（カーネル）を見ることができるので、データ分析の勉強にも非常に有用です。

Kaggleを始める方法

1. Kaggleのサイトに接続してアカウントを作成します。
2. Jupyter Notebook（Webブラウザで動作するpythonの対話環境）をインストールします。
3. コンペに参加

これだけです！

住宅価格予測

今回は、Kaggleの回帰問題のチュートリアルである「House Prices: Advanced Regression Techniques」をやってみます。（cybozu02さんのカーネルを参考にしております）

www.kaggle.com

1. ライブラリの準備

今回使用するライブラリをインポートします。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import PowerTransformer

• numpy： 数値計算を効率的に行うための拡張モジュール
• pandas： データ解析を支援する機能を提供するライブラリ
• matplotlib： グラフ描画ライブラリ
• seaborn： 可視化ライブラリ
• PowerTransformer： 0や負値を含んだ変数も対数変換ができるライブラリ

2. データの読み込み

こちらのリンクから以下のデータセットをダウンロードして読み込みます。

• train.csv： 学習データ
• test.csv： テストデータ

#pandasを使用してcsvデータ読み込み
train = pd.read_csv("train.csv")
test_x = pd.read_csv("test.csv")

#学習データの特徴量を出力
train.columns

学習データの特徴量は79個あります。

• MSSubClass： 住居のタイプ
  
• MSZoning： 住居のゾーン
  
• LotFrontage： 隣接した道路の長さ
  
• LotArea： 敷地面積
  
• Street： 道路のタイプ
  
• Alley： 路地のタイプ
  
• LotShape： 土地の形状
  
• LandContour： 土地の平坦さ
  
• Utilities： ガス・電気・水の利用
  
• LotConfig： 土地の構成
  
• LandSlope：土地の傾斜
  
• Neighborhood： 近所
  

等です。
実践さながらのデータでわくわくしてきますね！

3. データ分析

まずは、予測対象の変数（目的変数）である「SalesPrice（住宅価格）」について確認します。

plt.hist(train_saleprice, bins=100)
plt.show()

左に偏っており正規分布になっていません。機械学習において目的変数が正規分布に従ったほうが精度が良くなる、ということが知られているので後ほど正規分布に変換します。

次に、特徴量を確認します。

• TotalBsmtSF： 地下室の広さ
• 1stFlrSF： １階の広さ
• GrLivArea： リビングの広さ
• GarageArea： 車庫の広さ
• LotArea： 土地の広さ

の広さが、住宅価格に影響するのでは、と考えます。
相関関係を見てみましょう。

#地下室の広さと住宅価格の散布図を作成
data = pd.concat([train["TotalBsmtSF"],train["SalePrice"]],axis=1)
plt.figure(figsize=(20, 10))
plt.scatter(train["TotalBsmtSF"],train["SalePrice"])
plt.xlabel("TotalBsmtSF")
plt.ylabel("SalePrice")

#１階の広さと住宅価格の散布図を作成
data = pd.concat([train["1stFlrSF"],train["SalePrice"]],axis=1)
plt.figure(figsize=(20, 10))
plt.scatter(train["1stFlrSF"],train["SalePrice"])
plt.xlabel("1stFlrSF")
plt.ylabel("SalePrice")

#リビングの広さと住宅価格の散布図を作成
data = pd.concat([train["GrLivArea"],train["SalePrice"]],axis=1)
plt.figure(figsize=(20, 10))
plt.scatter(train["GrLivArea"],train["SalePrice"])
plt.xlabel("GrLivArea")
plt.ylabel("SalePrice")

#車庫の広さと住宅価格の散布図を作成
data = pd.concat([train["GarageArea"],train["SalePrice"]],axis=1)
plt.figure(figsize=(20, 10))
plt.scatter(train["GarageArea"],train["SalePrice"])
plt.xlabel("GarageArea")
plt.ylabel("SalePrice")

#土地の広さと住宅価格の散布図を作成
data = pd.concat([train["LotArea"],train["SalePrice"]],axis=1)
plt.figure(figsize=(20, 10))
plt.scatter(train["LotArea"],train["SalePrice"])
plt.xlabel("LotArea")
plt.ylabel("SalePrice")

広くなるにつれて住宅価格が上昇する傾向が見られますね。外れ値がいくつかあるので除外しておきます。

#外れ値を除外
train = train.drop(train[(train['TotalBsmtSF']>6000) & (train['SalePrice']<200000)].index)
train = train.drop(train[(train['1stFlrSF']>4000) & (train['SalePrice']<200000)].index)
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<200000)].index)
train = train.drop(train[(train['GarageArea']>1200) & (train['SalePrice']<300000)].index)
train = train.drop(train[(train['LotArea']>80000) & (train['SalePrice']<400000)].index)

4. 特徴量の作成

より良い特徴量を作成する特徴量エンジニアリングを進めていきます。

目的変数（住宅価格）を学習データから切り出します。また、学習データとテストデータをまとめて処理するため一旦マージしておきます。

#学習データを目的変数（住宅価格）とそれ以外に分ける
train_x = train.drop("SalePrice",axis=1)
train_y = train["SalePrice"]

#学習データとテストデータをマージ
all_data = pd.concat([train_x,test_x],axis=0,sort=True)

#IDのカラムは学習に不要なので別の変数に格納
train_ID = train['Id']
test_ID = test_x['Id']

all_data.drop("Id", axis = 1, inplace = True)

5. 欠損値の処理

分析に利用するデータには多くの場合、何らかの理由により記録されなかった値、欠損値が含まれます。欠損値があると統計処理ができなくなるので、削除や0埋め等が必要です。まずは欠損値がどのくらいあるのかを確認します。

#データの欠損値
all_data_na = all_data.isnull().sum()[all_data.isnull().sum()>0].sort_values(ascending=False)

#欠損値の数をグラフ化
plt.figure(figsize=(20,10))
plt.xticks(rotation='90')
sns.barplot(x=all_data_na.index, y=all_data_na)

#欠損値があるカラムをリスト化
na_col_list = all_data.isnull().sum()[all_data.isnull().sum()>0].index.tolist()
all_data[na_col_list].dtypes.sort_values()

• PoolQC： プールの質を表す。プールがない場合にはNAとなる。
• MiscFeature： 備え付けられている特別な設備（エレベータ等）を表す。何もない場合はNAとなる。
• Alley： 物件にアクセスするための路地のタイプを表す。該当しない場合はNAとなる。
• Fence： 敷地などを仕切る囲いの質を表す。囲いがない場合はNAとなる。
• FireplaceQu： 暖炉の質を表す。暖炉がない場合はNAとなる。

が多く欠損しています。
対象物そのものが存在しない場合、欠損として扱われているようです。

• float型の場合は「0」
• object型の場合は「None」

で置換することにしましょう。

物件に隣接した道路の長さ（LotFrontage）については平均値で穴埋めします。

#物件に隣接した道路の長さ（LotFrontage）の欠損値を平均値で穴埋め
all_data['LotFrontage'] = all_data['LotFrontage'].fillna(all_data['LotFrontage'].median())

#欠損値が存在する、かつfloat型の場合は「0」で置換
float_list = all_data[na_col_list].dtypes[all_data[na_col_list].dtypes == "float64"].index.tolist()
all_data[float_list] = all_data[float_list].fillna(0)

#欠損値が存在する、かつobject型の場合は「None」で置換
obj_list = all_data[na_col_list].dtypes[all_data[na_col_list].dtypes == "object"].index.tolist()
all_data[obj_list] = all_data[obj_list].fillna("None")

6. 数値変数の処理

• MSSubClass： 住宅の種類
• YrSold： 販売年
• MoSold： 販売月

については、数値ですが、数や量で測れない変数なので、カテゴリ変数として扱うことにします。

#カテゴリ変数に変換する
all_data['MSSubClass'] = all_data['MSSubClass'].apply(str)
all_data['YrSold'] = all_data['YrSold'].astype(str)
all_data['MoSold'] = all_data['MoSold'].astype(str)

7. 目的変数を対数変換

目的変数（住宅価格）が正規分布になっていないため、対数変換をすることで正規分布にします。

#目的変数を対数変換する
train_y = np.log1p(train_y)

#可視化
plt.figure(figsize=(20, 10))
sns.distplot(train_y)

8. 説明変数を対数変換

説明変数についても正規分布になっていない（歪度が0.5よりも大きい）ものは対数変換していきます。0や負値を含んだ変数も対数変換ができる「Yeo-Johnson変換」を使用します。

#数値の説明変数のリスト
num_feats = all_data.dtypes[all_data.dtypes != "object" ].index

#各説明変数の歪度を計算
skewed_feats = all_data[num_feats].apply(lambda x: x.skew()).sort_values(ascending = False)

#歪度の絶対値が0.5より大きい変数だけに絞る
skewed_feats_over = skewed_feats[abs(skewed_feats) > 0.5].index

#Yeo-Johnson変換
pt = PowerTransformer()
pt.fit(all_data[skewed_feats_over])

#変換後のデータで各列を置換
all_data[skewed_feats_over] = pt.transform(all_data[skewed_feats_over])

9. 特徴量の追加

新しい特徴量を追加することで、予測モデルの性能が向上することがあります。以下の特徴量を追加します。

#プールの有無（広さが0より大きい場合は有）
all_data['haspool'] = all_data['PoolArea'].apply(lambda x: 1 if x > 0 else 0)

#２階の有無（広さが0より大きい場合は有）
all_data['has2ndfloor'] = all_data['2ndFlrSF'].apply(lambda x: 1 if x > 0 else 0)

#ガレージの有無（広さが0より大きい場合は有）
all_data['hasgarage'] = all_data['GarageArea'].apply(lambda x: 1 if x > 0 else 0)

#地下室の有無（広さが0より大きい場合は有）
all_data['hasbsmt'] = all_data['TotalBsmtSF'].apply(lambda x: 1 if x > 0 else 0)

#暖炉の有無
all_data['hasfireplace'] = all_data['Fireplaces'].apply(lambda x: 1 if x > 0 else 0)

10. カテゴリ変数の処理

カテゴリ変数を、One-Hotエンコーディングにより「0」と「1」だけの数列に変換します。

#カテゴリ変数となっているカラム
cal_list = all_data.dtypes[all_data.dtypes=="object"].index.tolist()

#pandasのget_dummies関数でOne-Hotエンコーディング
all_data = pd.get_dummies(all_data,columns=cal_list)

11. データの分割

マージした学習データとテストデータを再度分割します。

#学習データとテストデータに再分割
train_x = all_data.iloc[:train_x.shape[0],:].reset_index(drop=True)
test_x = all_data.iloc[train_x.shape[0]:,:].reset_index(drop=True)

12. 予測モデルの作成＆実行

勾配ブースティング（アンサンブル学習と決定木を組み合わせた手法）のライブラリ「XGBoost」でモデルを作成します。

• max_depth： 木の深さの最大値
• eta： 学習率を調整
• objective： 損失関数を指定

等のハイパーパラメータについて、設定値を色々試してみたいですね。

from sklearn.ensemble import RandomForestRegressor,  GradientBoostingRegressor
from sklearn.model_selection import KFold, cross_val_score, train_test_split
from sklearn.metrics import mean_squared_error
import xgboost as xgb
from sklearn.model_selection import GridSearchCV

#データの分割
train_x, valid_x, train_y, valid_y = train_test_split(
        train_x,
        train_y,
        test_size=0.3,
        random_state=0)

#特徴量と目的変数をxgboostのデータ構造に変換する
dtrain = xgb.DMatrix(train_x, label=train_y)
dvalid = xgb.DMatrix(valid_x,label=valid_y)

#パラメータを指定して勾配ブースティング
num_round = 5000
evallist = [(dvalid, 'eval'), (dtrain, 'train')]

evals_result = {}

#パラメータ
param = {
    'max_depth': 3,
    'eta': 0.01,
    'objective': 'reg:squarederror',
}

#学習の実行
bst = xgb.train(
    param, dtrain,
    num_round,
    evallist,
    evals_result=evals_result,
)    

#学習曲線を可視化する
plt.figure(figsize=(20, 10))
train_metric = evals_result['train']['rmse']
plt.plot(train_metric, label='train rmse')
eval_metric = evals_result['eval']['rmse']
plt.plot(eval_metric, label='eval rmse')
plt.grid()
plt.legend()
plt.xlabel('rounds')
plt.ylabel('rmse')
plt.ylim(0, 0.3)
plt.show()

13. 提出用csvファイルの作成

いよいよ最後です。予測結果をcsvファイルに出力します。予測結果のSalePrice（住宅価格）は対数をとった値なので、Exponentialをかけてあげる必要があります。

#予測結果を出力
dtest = xgb.DMatrix(test_x)
my_submission = pd.DataFrame()
my_submission["Id"] = test_ID
my_submission["SalePrice"] = np.exp(bst.predict(dtest))
my_submission.to_csv('submission.csv', index=False)

1998位/5343人でした！

以下は今回参考にさせていただいたcybozu02さんのカーネルです。ありがとうございました。 www.kaggle.com