KaggleのHouse Pricesに挑戦(0.16781)

KaggleのHousePriceに挑戦した話

前回、Titanicに挑戦し、今回はHousePriceに挑戦しただけです。
ネットで見た感じラッソ回帰が主流?のようですが、アンサンブル学習をしたことがないので挑戦します。
今回も前回同様、まとめにプログラムの全容が載っているので結果だけ知りたい方はそちらに。
この記事をかなり参考にしました。初めのやっていることは同じ...

使用環境

• GoogleColaboratory
• Python

使用モデル

• Lasso
• ElasticNet
• Kernel Ridge Regression
• Gradient Boosting Regression
• XGBoost
• LightGBM
• Random Forest
• Support Vector
  上記のモデルを使用してアンサンブル学習を行ってみます。

相関関係を探る

test.py

import pandas as pd
import numpy as np
import seaborn as sns
from scipy import stats
from scipy.stats import norm, skew
from scipy.special import boxcox1p
import matplotlib.pyplot as plt

df = pd.read_csv('train-2.csv')

予測したい値のSalePriceと他のカラム全ての相関関係が知りたい。0.5以下の相関関係は学習には向かなそうなので...
数値以外のカラムもあるので、ダミー変数化して相関関係を出してみる。

HP_correlation.py

import pandas as pd

df = pd.read_csv('train-2.csv')

df_encoded = pd.get_dummies(df)
correlations = df_encoded.corr()['SalePrice'].sort_values(ascending=False)
strong_correlations = correlations[correlations >= 0.5]
print("SalePriceとの相関係数が0.5以上のカラム:")
print(strong_correlations)

上記の実行結果

カラム名	相関数
SalePrice	1.000000
OverallQual	0.790982
GrLivArea	0.708624
GarageCars	0.640409
GarageArea	0.623431
TotalBsmtSF	0.613581
1stFlrSF	0.605852
FullBath	0.560664
BsmtQual_Ex	0.553105
TotRmsAbvGrd	0.533723
YearBuilt	0.522897
YearRemodAdd	0.507101
KitchenQual_Ex	0.504094

0.6以上でも6カラムもある...一旦0.6以上だけで学習させるか!!

データの前処理

外れ値の削除

敷地面積と相関が強かったので、値段と敷地面積の関係をプロットしてみると

test.py

fig, ax = plt.subplots()
ax.scatter(x = df['GrLivArea'], y = df['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

となっており2点外れ値がいる...ということで削除！して再度プロットすると

test.py

df = df.drop(df[(df['GrLivArea']>4000) & (df['SalePrice']<300000)].index)
fig, ax = plt.subplots()
ax.scatter(df['GrLivArea'], df['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

うむ。良さそう。

分布を正規分布にする

とりあえず、今の状態を確認

test.py

sns.distplot(df['SalePrice'] , fit=norm)
fig = plt.figure()
res = stats.probplot(df['SalePrice'], plot=plt)
plt.show()

うむ。左にシフトしているので対数を取って正規化します。(右の場合は指数ですね。これは理系大学生なのでなんとなく...)

test.py

df['SalePrice'] = np.log1p(df['SalePrice'])
sns.distplot(df['SalePrice'] , fit=norm)
fig = plt.figure()
res = stats.probplot(df['SalePrice'], plot=plt)
plt.show()

正規分布にうまくいった!!

歪度を見てみる

正規化をうまくいった他に使用と思っている値の歪度を確認してみる。

test.py

numeric_feats = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
skewed_feats = df[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
skewness = pd.DataFrame({'Skew' :skewed_feats})
print("特定のカラムの歪度: ")
print(skewness)

結果↓

特定のカラムの歪度: 
                 Skew
GrLivArea    1.009951
1stFlrSF     0.886723
TotalBsmtSF  0.511177
OverallQual  0.200579
GarageArea   0.131612
GarageCars  -0.342025

歪んでそうなので、正規分布に変換

歪んでいるので、正規分布に変換

test.py

numeric_feats = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
skewed_feats = df[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)

# Box-Cox変換を適用
for feat in numeric_feats:
    df[feat] = boxcox1p(df[feat], 0.15)

# 歪度の再計算
skewed_feats_after = df[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)

skewness_after = pd.DataFrame({'Skew (After Box-Cox)' :skewed_feats_after})
print("Box-Cox変換後の歪度: ")
print(skewness_after)

Box-Cox変換後の歪度: 
             Skew (After Box-Cox)
1stFlrSF                 0.139021
GrLivArea                0.082854
OverallQual             -0.446390
GarageCars              -1.204495
GarageArea              -3.040319
TotalBsmtSF             -4.161816

これだけじゃわからないので、プロットしてみる。正規化なのでQ-Qを使用!

test.py

features = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
num_plots = len(features)

fig, axes = plt.subplots(nrows=1, ncols=num_plots, figsize=(15, 5))

for i, feature in enumerate(features):
    ax = axes[i]
    stats.probplot(df[feature], plot=ax)
    ax.set_title('Q-Q Plot for {}'.format(feature))
    ax.set_xlabel('Theoretical quantiles')
    ax.set_ylabel('Sample quantiles')

plt.tight_layout()
plt.show()

正規化できてそう!!

学習

待ちに待った学習を行う。

使用ライブラリ

HP.py

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
from sklearn.svm import SVR

データの読み込み

HP.py

df = pd.read_csv('train-2.csv')

予測データと学習データを宣言

HP.py

y = df.SalePrice
features = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
X = df[features]

テストデータと学主データに分ける。

HP.py

train_x, test_x, train_y, test_y = train_test_split(X, y, random_state=1)

ハイパーパラメータ

モデルのハイパーパラメータを先に宣言しておいて変更しやすいようにしておく。この辺の変数はネットで見たときの値からいい感じに...(適当だって！？)

HP.py

lasso_alpha = 0.0005 # 正則化項の強度

enet_alpha = 0.0005
enet_ratio = 0.9 # 正則の割合1に近いとL1正則化を使用

krr_alpha = 0.6
krr_degree = 2 # 多項式カーネルの次数を指定
krr_coef0 = 2.5 # カーネル関数の定数項の数

gbr_n_estimators = 3000 # 木の本数
gbr_learning_rate = 0.05 # 値が小さいほど学習が安定
gbr_depth = 4 # 深さ
min_samples_leaf = 15 # 葉ノードに必要な最小サンプル数
min_samples_split = 10 # ノードを分割するために必要な最小サンプル数

colsample_bytree = 0.4603 # 各木を構築する際に使用する特徴量の割合
gamma = 0.0468 # ツリーの成長を停止するために必要な最小のロス減少
xgb_learning_rate = 0.05
Xgb_n_estimators = 2200
min_child_weight = 1.7817 # 葉ノードに必要な最小の重み合計
reg_alpha = 0.4640 # L1正則化の強度
reg_lambda = 0.8571 # L2正則化の強度
subsample = 0.5213 # 各木を構築する際に使用するトレーニングデータの割合

num_leaves = 5 # 決定木の葉ノードの最大数
lgb_learning_rate = 0.05
lgb_n_estimators = 720
max_bin = 55 # 特徴量のバケット化の粒度
bagging_fraction = 0.8 # ブースティング中のトレーニングデータのサブサンプリング率
bagging_freq = 5 # バギングの頻度
feature_fraction = 0.2319 # 各木を構築する際に使用する特徴量の割合
feature_fraction_seed = 1
bagging_seed = 1 
min_data_in_leaf = 6 # 葉ノードに必要な最小サンプル数
min_sum_hessian_in_leaf = 11 # 各葉の損失関数の二階微分の合計の最小値

rf_n_estimators = 100

モデルの宣言

HP.py

lasso = make_pipeline(RobustScaler(), Lasso(alpha =lasso_alpha, random_state=1))
enet_model = ElasticNet(alpha=enet_alpha, l1_ratio=enet_ratio, random_state=1)
krr_model = KernelRidge(alpha=krr_alpha, kernel='polynomial', degree=krr_degree, coef0=krr_coef0)
gbr_model = GradientBoostingRegressor(n_estimators=gbr_n_estimators, learning_rate=gbr_learning_rate,
                                   max_depth=gbr_depth, max_features='sqrt',
                                   min_samples_leaf=min_samples_leaf, min_samples_split=min_samples_split, 
                                   loss='huber', random_state =1)
xgb_model = xgb.XGBRegressor(colsample_bytree=colsample_bytree, gamma=gamma, 
                             learning_rate=xgb_learning_rate, max_depth=3, 
                             min_child_weight=min_child_weight, n_estimators=Xgb_n_estimators,
                             reg_alpha=reg_alpha, reg_lambda=reg_lambda,
                             subsample=subsample, verbosity=1,
                             random_state =1, nthread = -1)
lgb_model = lgb.LGBMRegressor(objective='regression',num_leaves=num_leaves,
                              learning_rate=lgb_learning_rate, n_estimators=lgb_n_estimators,
                              max_bin = max_bin, bagging_fraction = bagging_fraction,
                              bagging_freq = bagging_freq, feature_fraction = feature_fraction,
                              feature_fraction_seed=feature_fraction_seed, bagging_seed=bagging_seed,
                              min_data_in_leaf =min_data_in_leaf, min_sum_hessian_in_leaf = min_sum_hessian_in_leaf)
rf_model = RandomForestRegressor(n_estimators=rf_n_estimators, random_state=1)
svm_model = SVR(kernel='rbf')

学習用の関数を汎用化

平均二乗誤差でモデルを評価する

HP.py

def train_and_score(model, X_train, y_train, X_test, y_test):
    model.fit(X_train, y_train)
    y_pred = model.predict(X_test)
    mse = mean_squared_error(y_test, y_pred)
    return np.sqrt(mse)

学習のために呼び出し

HP.py

lasso_score = train_and_score(lasso, train_x, train_y, test_x, test_y)
enet_score = train_and_score(enet_model, train_x, train_y, test_x, test_y)
krr_score = train_and_score(krr_model, train_x, train_y, test_x, test_y)
gbr_score = train_and_score(gbr_model, train_x, train_y, test_x, test_y)
xgb_score = train_and_score(xgb_model, train_x, train_y, test_x, test_y)
lgb_score = train_and_score(lgb_model, train_x, train_y, test_x, test_y)
rf_score = train_and_score(rf_model, train_x, train_y, test_x, test_y)
svm_score = train_and_score(svm_model, train_x, train_y, test_x, test_y)

スコアの表示をする。

HP.py

print("Lasso score: {:.4f}".format(lasso_score))
print("ElasticNet score: {:.4f}".format(enet_score))
print("Kernel Ridge Regression score: {:.4f}".format(krr_score))
print("Gradient Boosting Regression score: {:.4f}".format(gbr_score))
print("XGBoost score: {:.4f}".format(xgb_score))
print("LightGBM score: {:.4f}".format(lgb_score))
print("Random Forest score: {:.4f}".format(rf_score))
print("Support Vector Machine score: {:.4f}".format(svm_score))

その結果がこちら！(見やすいように数字の位置を揃えています。)

Lasso score:                        0.1757
ElasticNet score:                   0.1756
Kernel Ridge Regression score:      0.1669
Gradient Boosting Regression score: 0.1914
XGBoost score:                      0.1696
LightGBM score:                     0.1705
Random Forest score:                0.1780
Support Vector Machine score:       0.1718

平均化

ベースモデルの学習ができたので、スタッキングしていきます。
追加ライブラリ

HP.py

from sklearn.model_selection import KFold

メタモデルとしてRandomForestを使ってみる。

HP.py

def stack_models(base_models, meta_model, X_train, y_train, X_test, kfolds):
    meta_train = np.zeros((X_train.shape[0], len(base_models)))
    meta_test = np.zeros((X_test.shape[0], len(base_models)))
    for i, model in enumerate(base_models):
        fold_errors = []
        for train_idx, val_idx in kfolds.split(X_train):
            X_tr, y_tr = X_train.iloc[train_idx], y_train.iloc[train_idx]
            X_val, y_val = X_train.iloc[val_idx], y_train.iloc[val_idx]
            model.fit(X_tr, y_tr)
            y_pred = model.predict(X_val)
            fold_errors.append(mean_squared_error(y_val, y_pred))
            meta_train[val_idx, i] = y_pred
            meta_test[:, i] += model.predict(X_test)
        print(f"Model {i+1} MSE: {np.mean(fold_errors)}")
        meta_test[:, i] /= kfolds.n_splits
    meta_model.fit(meta_train, y_train)
    return meta_model.predict(meta_test)
    
X_train, y_train = train_x, train_y
X_test, y_test = test_x, test_y
base_models = [lasso, enet_model, krr_model, gbr_model, xgb_model, lgb_model, rf_model, svm_model]
meta_model = RandomForestRegressor()  #　ここを変更
kfolds = KFold(n_splits=5, shuffle=True, random_state=42)
stacked_pred = stack_models(base_models, meta_model, X_train, y_train, X_test, kfolds)
stacked_score = mean_squared_error(y_test, stacked_pred)
print(f"Stacked model MSE: {stacked_score}")

色々Warningは出ますが、実行に影響はないので無視して
結果(ファイル名はメタモデル名です。)

RandomForest.py

Model 1 MSE: 0.02926193853799905
Model 2 MSE: 0.029271270207278155
Model 3 MSE: 0.028641278640404593
Model 4 MSE: 0.028683352377638433
Model 5 MSE: 0.02733761558695501
Model 6 MSE: 0.03017459126827753
Model 7 MSE: 0.028084919454466856
Model 8 MSE: 0.031085991891324362
model MSE:   0.029865955985898612

メタモデルを変更してみる。比較で最も良かったモデルも載せる。

other.py

Model 5 MSE:           0.02733761558695501
Lasso model MSE:       0.0284434507513744
XGB model MSE:         0.03086616831795726
LGB model MSE:         0.029649113329076986
GBR model MSE:         0.03004591253068424
KRR model MSE:         0.028494784717844786
ENet model MSE:        0.15790200930837153
SVR model MSE:         0.028155102434682805

スタッキング無しで提出してみる。

スタッキングしない方が、性能が良いので、色々利点はありますが、スタッキングなしで一度提出してみます。

HP.py

df = pd.read_csv('train-2.csv')
y = df.SalePrice
features = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
X = df[features]
train_x, test_x, train_y, test_y = train_test_split(X, y, random_state=1)
xgb_model = xgb.XGBRegressor(colsample_bytree=colsample_bytree, gamma=gamma, 
                             learning_rate=xgb_learning_rate, max_depth=3, 
                             min_child_weight=min_child_weight, n_estimators=Xgb_n_estimators,
                             reg_alpha=reg_alpha, reg_lambda=reg_lambda,
                             subsample=subsample, verbosity=1,
                             random_state =1, nthread = -1)
xgb_model.fit(X_train, y_train)
test_data = pd.read_csv('test-2.csv')
features = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
X_test = test_data[features]
xgb_pred = xgb_model.predict(X_test)

predictions_df = pd.DataFrame({'Id': test_data['Id'], 'SalePrice': np.exp(xgb_pred)})
predictions_df.to_csv('predictions3.csv', index=False, float_format='%.6f')

として6桁まで保存し、提出したところ"1.17267"でした。んー。精度低い

スタッキング有りでやってみる。

スタッキングの際は最も性能の良かったLassoをメタモデルにする。

HP.py

import pandas as pd
import numpy as np
from scipy.special import boxcox1p
from sklearn.model_selection import KFold
from sklearn.linear_model import Lasso, ElasticNet
from sklearn.kernel_ridge import KernelRidge
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
import xgboost as xgb
import lightgbm as lgb
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error

# データの読み込み
train_df = pd.read_csv('train-4.csv')
test_df = pd.read_csv('test-3.csv')

# 外れ値の削除
train_df = train_df.drop(train_df[(train_df['GrLivArea']>4000) & (train_df['SalePrice']<300000)].index)

# 目的変数の対数変換
train_df['SalePrice'] = np.log1p(train_df['SalePrice'])

# 特徴量の選択
numeric_feats = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']

# 欠損値の補完と特徴量の変換
for feat in numeric_feats:
    train_df[feat].fillna(train_df[feat].mean(), inplace=True)
    train_df[feat] = boxcox1p(train_df[feat], 0.15)
    test_df[feat].fillna(test_df[feat].mean(), inplace=True)
    test_df[feat] = boxcox1p(test_df[feat], 0.15)

# 特徴量と目的変数の準備
X = train_df[numeric_feats]
y = train_df['SalePrice']

# モデルの定義
lasso = Lasso(alpha=0.001)
enet_model = ElasticNet(alpha=0.001, l1_ratio=0.5)
krr_model = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)
gbr_model = 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=1)
xgb_model = 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, verbosity=1, random_state=1, nthread=-1)
lgb_model = 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=1, bagging_seed=1, min_data_in_leaf=6,
                               min_sum_hessian_in_leaf=11)
rf_model = RandomForestRegressor(n_estimators=100, random_state=1)
svm_model = SVR(kernel='linear')

# ベースモデルとメタモデルの定義
base_models = [lasso, enet_model, krr_model, gbr_model, xgb_model, lgb_model, rf_model, svm_model]
meta_model = Lasso(alpha=0.001)
kfolds = KFold(n_splits=5, shuffle=True, random_state=42)

# スタッキングで予測を行う関数を作成
def stack_predict(base_models, meta_model, X_train, y_train, X_test, kfolds):
    meta_train = np.zeros((X_train.shape[0], len(base_models)))
    meta_test = np.zeros((X_test.shape[0], len(base_models)))
    
    for i, model in enumerate(base_models):
        fold_errors = []
        for train_idx, val_idx in kfolds.split(X_train):
            X_tr, y_tr = X_train.iloc[train_idx], y_train.iloc[train_idx]
            X_val, y_val = X_train.iloc[val_idx], y_train.iloc[val_idx]
            
            model.fit(X_tr, y_tr)
            y_pred = model.predict(X_val)
            fold_errors.append(mean_squared_error(y_val, y_pred))
            
            meta_train[val_idx, i] = y_pred
            meta_test[:, i] += model.predict(X_test)
        
        print(f"Model {i+1} MSE: {np.mean(fold_errors)}")
        meta_test[:, i] /= kfolds.n_splits
    
    meta_model.fit(meta_train, y_train)
    return meta_model.predict(meta_test)

# テストデータに対して予測を行う
test_predictions = stack_predict(base_models, meta_model, X, y, test_df[numeric_feats], kfolds)

# 予測結果をデータフレームに格納
predictions_df = pd.DataFrame({'Id': test_df['Id'], 'SalePrice': np.exp(test_predictions)})  # 対数変換を元に戻す

# CSVファイルに保存
predictions_df.to_csv('predictions.csv', index=False)

上記のプログラムを使用したところ、スコアは"0.16781"となり、断然良くなりました。
ハイパーパラメータの最適化まで行えばもう少しモデルとしては良くなりそうですが、今回はここまでにします。
ありがとうございました。

まとめ

これらを使用してデータが学習しやすい形であるかどうかを評価しました。

HP_correlation.py

import pandas as pd
import numpy as np
import seaborn as sns
from scipy import stats
from scipy.stats import norm, skew
from scipy.special import boxcox1p
import matplotlib.pyplot as plt

df = pd.read_csv('train-2.csv')

df_encoded = pd.get_dummies(df)
correlations = df_encoded.corr()['SalePrice'].sort_values(ascending=False)
strong_correlations = correlations[correlations >= 0.5]
print("SalePriceとの相関係数が0.5以上のカラム:")
print(strong_correlations)

# 敷地面積との関連が強かったためデータを表示。
fig, ax = plt.subplots()
ax.scatter(x = df['GrLivArea'], y = df['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

# 敷地面積が大きいが値段が低い外れ値が2点あるので、削除し再度表示。
df = df.drop(df[(df['GrLivArea']>4000) & (df['SalePrice']<300000)].index)
fig, ax = plt.subplots()
ax.scatter(df['GrLivArea'], df['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()
sns.distplot(df['SalePrice'] , fit=norm)
fig = plt.figure()
res = stats.probplot(df['SalePrice'], plot=plt)
plt.show()

# データが左に偏っているので、対数変換して正規分布にする。右に偏っている場合は指数変換
df['SalePrice'] = np.log1p(df['SalePrice'])
sns.distplot(df['SalePrice'] , fit=norm)
fig = plt.figure()
res = stats.probplot(df['SalePrice'], plot=plt)
plt.show()

# いい感じに正規分布になったので、今回学習に用いようと思っているデータの歪度を計算し、出力
numeric_feats = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
skewed_feats = df[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)
skewness = pd.DataFrame({'Skew' :skewed_feats})
print("歪度: ")
print(skewness)

# 歪んでいるので、正規分布に変換
numeric_feats = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
skewed_feats = df[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)

# Box-Cox変換を適用
for feat in numeric_feats:
    df[feat] = boxcox1p(df[feat], 0.15)

# 歪度の再計算
skewed_feats_after = df[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)

skewness_after = pd.DataFrame({'Skew (After Box-Cox)' :skewed_feats_after})
print("Box-Cox変換後のカラムの歪度: ")
print(skewness_after)

# Q-Qを使用して正規分布になったかを確認する。だいたい良さそう?
features = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
num_plots = len(features)

fig, axes = plt.subplots(nrows=1, ncols=num_plots, figsize=(15, 5))

for i, feature in enumerate(features):
    ax = axes[i]
    stats.probplot(df[feature], plot=ax)
    ax.set_title('Q-Q Plot for {}'.format(feature))
    ax.set_xlabel('Theoretical quantiles')
    ax.set_ylabel('Sample quantiles')

plt.tight_layout()
plt.show()

１回目の提出時はスタッキングを使用せずにXGBモデルだけを使用したものを提出しました。

HP_use_xgb.py

import pandas as pd
import numpy as np
import xgboost as xgb
from sklearn.model_selection import train_test_split

colsample_bytree = 0.4603 # 各木を構築する際に使用する特徴量の割合
gamma = 0.0468 # ツリーの成長を停止するために必要な最小のロス減少
xgb_learning_rate = 0.05
Xgb_n_estimators = 2200
min_child_weight = 1.7817 # 葉ノードに必要な最小の重み合計
reg_alpha = 0.4640 # L1正則化の強度
reg_lambda = 0.8571 # L2正則化の強度
subsample = 0.5213 # 各木を構築する際に使用するトレーニングデータの割合

df = pd.read_csv('train-2.csv')
df = df.drop(df[(df['GrLivArea']>4000) & (df['SalePrice']<300000)].index)
df['SalePrice'] = np.log1p(df['SalePrice'])
y = df.SalePrice
features = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
X = df[features]
for feat in features:
    df[feat] = boxcox1p(df[feat], 0.15)
train_x, test_x, train_y, test_y = train_test_split(X, y, random_state=1)
xgb_model = xgb.XGBRegressor(colsample_bytree=colsample_bytree, gamma=gamma, 
                             learning_rate=xgb_learning_rate, max_depth=3, 
                             min_child_weight=min_child_weight, n_estimators=Xgb_n_estimators,
                             reg_alpha=reg_alpha, reg_lambda=reg_lambda,
                             subsample=subsample, verbosity=1,
                             random_state =1, nthread = -1)
xgb_model.fit(X_train, y_train)
test_data = pd.read_csv('test-2.csv')
features = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']
X_test = test_data[features]
xgb_pred = xgb_model.predict(X_test)

predictions_df = pd.DataFrame({'Id': test_data['Id'], 'SalePrice': np.exp(xgb_pred)})
predictions_df.to_csv('predictions3.csv', index=False, float_format='%.6f')

最も良い精度を発揮したスタッキングしたモデルです。

HP_use_stack.py

import pandas as pd
import numpy as np
from scipy.special import boxcox1p
from sklearn.model_selection import KFold
from sklearn.linear_model import Lasso, ElasticNet
from sklearn.kernel_ridge import KernelRidge
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
import xgboost as xgb
import lightgbm as lgb
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error

# データの読み込み
train_df = pd.read_csv('train-4.csv')
test_df = pd.read_csv('test-3.csv')

# 外れ値の削除
train_df = train_df.drop(train_df[(train_df['GrLivArea']>4000) & (train_df['SalePrice']<300000)].index)

# 目的変数の対数変換
train_df['SalePrice'] = np.log1p(train_df['SalePrice'])

# 特徴量の選択
numeric_feats = ['OverallQual', 'GrLivArea', 'GarageCars', 'GarageArea', 'TotalBsmtSF', '1stFlrSF']

# 欠損値の補完と特徴量の変換
for feat in numeric_feats:
    train_df[feat].fillna(train_df[feat].mean(), inplace=True)
    train_df[feat] = boxcox1p(train_df[feat], 0.15)
    test_df[feat].fillna(test_df[feat].mean(), inplace=True)
    test_df[feat] = boxcox1p(test_df[feat], 0.15)

# 特徴量と目的変数の準備
X = train_df[numeric_feats]
y = train_df['SalePrice']

# モデルの定義
lasso = Lasso(alpha=0.001)
enet_model = ElasticNet(alpha=0.001, l1_ratio=0.5)
krr_model = KernelRidge(alpha=0.6, kernel='polynomial', degree=2, coef0=2.5)
gbr_model = 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=1)
xgb_model = 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, verbosity=1, random_state=1, nthread=-1)
lgb_model = 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=1, bagging_seed=1, min_data_in_leaf=6,
                               min_sum_hessian_in_leaf=11)
rf_model = RandomForestRegressor(n_estimators=100, random_state=1)
svm_model = SVR(kernel='linear')

# ベースモデルとメタモデルの定義
base_models = [lasso, enet_model, krr_model, gbr_model, xgb_model, lgb_model, rf_model, svm_model]
meta_model = Lasso(alpha=0.001)
kfolds = KFold(n_splits=5, shuffle=True, random_state=42)

# スタッキングで予測を行う関数を作成
def stack_predict(base_models, meta_model, X_train, y_train, X_test, kfolds):
    meta_train = np.zeros((X_train.shape[0], len(base_models)))
    meta_test = np.zeros((X_test.shape[0], len(base_models)))
    
    for i, model in enumerate(base_models):
        fold_errors = []
        for train_idx, val_idx in kfolds.split(X_train):
            X_tr, y_tr = X_train.iloc[train_idx], y_train.iloc[train_idx]
            X_val, y_val = X_train.iloc[val_idx], y_train.iloc[val_idx]
            
            model.fit(X_tr, y_tr)
            y_pred = model.predict(X_val)
            fold_errors.append(mean_squared_error(y_val, y_pred))
            
            meta_train[val_idx, i] = y_pred
            meta_test[:, i] += model.predict(X_test)
        
        print(f"Model {i+1} MSE: {np.mean(fold_errors)}")
        meta_test[:, i] /= kfolds.n_splits
    
    meta_model.fit(meta_train, y_train)
    return meta_model.predict(meta_test)

# テストデータに対して予測を行う
test_predictions = stack_predict(base_models, meta_model, X, y, test_df[numeric_feats], kfolds)

# 予測結果をデータフレームに格納
predictions_df = pd.DataFrame({'Id': test_df['Id'], 'SalePrice': np.exp(test_predictions)})  # 対数変換を元に戻す

# CSVファイルに保存
predictions_df.to_csv('predictions.csv', index=False)