【House Prices】Kaggleの住宅価格の予測を40代サラリーマンが解く！初心者にもわかりやすく

# 各ライブラリのインポート
import numpy as np                # 数値計算ライブラリのnumpyをインポート
import pandas as pd               # データ操作ライブラリのpandasをインポート
import matplotlib.pyplot as plt   # データ可視化ライブラリのmatplotlibをインポート
import seaborn as sns             # データ可視化ライブラリのseabornをインポート
%matplotlib inline               # グラフをノートブック内に表示するためのマジックコマンド

# 訓練データとテストデータの読み込み
df_train = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/train.csv')  # 訓練データのCSVファイルを読み込み
df_test = pd.read_csv('/kaggle/input/house-prices-advanced-regression-techniques/test.csv')   # テストデータのCSVファイルを読み込み

# 訓練データのサイズを表示し、最初の3行を表示
print(f'df_train.shape: {df_train.shape}')      # 訓練データの行数と列数を表示
display(df_train.head(3))                      # 訓練データの最初の3行を表示

# テストデータのサイズを表示し、最初の3行を表示
print(f'df_test.shape: {df_test.shape}')        # テストデータの行数と列数を表示
display(df_test.head(3))                       # テストデータの最初の3行を表示

from scipy.stats import norm       # SciPyライブラリから確率分布用のnormモジュールをインポート
from scipy import stats           # SciPyライブラリから統計関連の機能を提供するstatsモジュールをインポート

# 分析対象のデータを抽出して変数"data"に格納
# data変数を用いることにより、変数の入れ替えを楽にできます
data = df_train['SalePrice']       

# 目的変数（SalePrice）の基本統計量を表示
print('-'*10, 'describe', '-'*10)
print(round(data.describe(), 1))

# 目的変数の歪度と尖度を計算して表示
skewness = data.skew()
kurtosis = data.kurtosis()
print(f'歪度 (Skewness): {skewness:.2f}')
print(f'尖度 (Kurtosis): {kurtosis:.2f}')

# 2つのサブプロットを持つ図を作成
fig, ax = plt.subplots(1, 2, figsize=(10, 4))

# ヒストグラムと正規分布を表示するサブプロット
sns.histplot(data, stat='density', kde=True, ax=ax[0])
ax[0].set_title('ヒストグラムと正規分布')

# x軸のラベルの文字サイズを調整
ax[0].tick_params(axis='x', labelsize=8)

xmin, xmax = ax[0].get_xlim()  # ax[0]の最大値と最小値を取得
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, np.mean(data), np.std(data))
ax[0].plot(x, p, 'k', linewidth=1)

# 正規確率プロットを表示するサブプロット
res = stats.probplot(data, plot=ax[1])
ax[1].set_title('正規確率プロット')

plt.tight_layout()  # グラフのレイアウトを調整

# グラフを表示
plt.show()

# 散布図グラフの描画
ax = sns.scatterplot(data=df_train, x='GrLivArea', y='SalePrice')

# 以下の部分は解析には不要 (説明のための矢印をつけています)
# 注釈を追加
for i in range(2):
    idx = df_train['GrLivArea'].nlargest(2).index[i]

    # x軸の最大値を見つける
    max_x_values = df_train.loc[idx,'GrLivArea']
    max_y_value = df_train.loc[idx,'SalePrice']

    # 矢印を追加
    plt.annotate('', xy=(max_x_values, max_y_value), xytext=(max_x_values + 250, max_y_value + 100),
                 arrowprops=dict(color='red', arrowstyle='->'))

# さらに注釈を追加
plt.annotate('', xy=(2000, 250000), xytext=(1000, 100000),
             arrowprops=dict(color='yellow', arrowstyle='->', alpha=0.8, linewidth=5))

plt.show()

# Seabornのboxplotを作成
plt.figure(figsize=(16, 5))  # プロットのサイズを設定
sns.boxplot(data=df_train, x='YearBuilt', y='SalePrice')

# x軸ラベルの回転とプロットの表示
plt.xticks(rotation=90)

# グラフを表示
plt.show()

# Overall Qualityに基づく箱ひげ図を描画
sns.boxplot(data=df_train, x='OverallQual', y='SalePrice', fliersize=1)

# グラフを表示
plt.show()

# ヒートマップの作成だけなら下記だけで描画可能
sns.heatmap(df_train.corr())
plt.show()

# ----------------- 以下は見やすくする & 説明のための点線枠を追加した場合のコード ----------------
# 図（Figure）と軸（Axes）を作成
fig, ax = plt.subplots(figsize=(10, 8))

# ヒートマップを作成
sns.heatmap(df_train.corr(), cmap='viridis', annot=True, fmt='.1f', annot_kws={'fontsize': 5})

# 赤い点線の四角を追加
rect = plt.Rectangle((0.1, 36.9), 37.8, 0.9, edgecolor='red', linewidth=2, fill=False, linestyle='dashed')
ax.add_patch(rect)

# オレンジの点線の四角を追加
xy = [11.9, 25.9]
hw = [2.1, 2.1]

for i, v in enumerate(xy):
    x, y = xy[i], xy[i]
    height, width = hw[i], hw[i]
    rect = plt.Rectangle((x, y), height, width, edgecolor='orange', linewidth=2, fill=False, linestyle='dashed')
    ax.add_patch(rect)

# グラフを表示
plt.show()

threshold = 0.4  # 相関係数のしきい値

# 相関係数が絶対値0.4以上の特徴量を抜き出し、スコアの高い順にソート
high_corr_cols = df_train.corr()['SalePrice'][abs(df_train.corr()['SalePrice']) >= threshold].sort_values(ascending=False).index

# ヒートマップを作成
sns.heatmap(df_train[high_corr_cols].corr(), cmap='viridis', annot=True, fmt='.1f', annot_kws={'fontsize': 10})
plt.show()

# DataFrameを結合
all_data = pd.concat([df_train, df_test])

# 欠損値の数を計算
missing_values_count = df.isnull().sum()

# 欠損値の数とパーセンテージを表形式で表示するためのDataFrameを作成
missing_values_table = pd.DataFrame({
    'Missing_total': missing_values_count,
    'Percent (%)': round((missing_values_count / len(df)) * 100, 1)
})

# 欠損値がある列の情報を抽出して降順でソート
df_missing = missing_values_table[missing_values_table['Missing_total'] > 0].sort_values(by='Missing_total', ascending=False)
    
# 欠損値がある列の情報を表示
display(df_missing)

# 図（Figure）のサイズを設定
plt.figure(figsize=(20, 4))

# Seabornを使用して箱ひげ図を作成
# x軸には "Neighborhood"、y軸には "LotFrontage" を設定
sns.boxplot(data=all_data, x='Neighborhood', y='LotFrontage')

# グラフを表示
plt.show()

# 各地域（Neighborhood）ごとの "LotFrontage" の中央値を計算し、新しいデータフレームに保存
df_group = all_data.groupby(by='Neighborhood')['LotFrontage'].agg('median')

# 欠損値を埋めるための関数を定義
def fillnaLot(row):
    # もし "LotFrontage" が欠損値（NaN）または空文字列の場合
    if pd.isna(row['LotFrontage']) or row['LotFrontage'] == '':
        # 対応する地域（Neighborhood）の中央値で欠損値を埋める
        return df_group[row['Neighborhood']]
    else:
        # それ以外の場合、元の値を保持
        return row['LotFrontage']

# 訓練データとテストデータに対して "LotFrontage" の欠損値を埋める
datasets = [df_train, df_test]
for df in datasets:
    df['LotFrontage'] = df.apply(fillnaLot, axis=1)

# 訓練データとテストデータのデータフレームをリストにまとめる
datasets = [df_train, df_test]

# 欠損値を 'None' で置き換える列のリスト
cols = ['PoolQC', 'MiscFeature', 'Alley', 'Fence', 'FireplaceQu', 'GarageCond', 'GarageQual', 'GarageFinish', 'GarageType',
       'BsmtCond', 'BsmtExposure', 'BsmtQual', 'BsmtFinType1', 'BsmtFinType2']

# 各データセットに対して、指定した列の欠損値を 'None' で置き換える
for df in datasets:
    for col in cols:
        df[col] = df[col].fillna('None')

# 訓練データとテストデータのデータフレームをリストにまとめる
datasets = [df_train, df_test]

# 欠損値を 0 で置き換える列のリスト
cols = ['GarageYrBlt', 'MasVnrArea', 'BsmtFullBath', 'BsmtHalfBath', 'GarageCars', 'GarageArea', 'TotalBsmtSF', 'BsmtUnfSF', 'BsmtFinSF2', 'BsmtFinSF1']

# 各データセットに対して、指定した列の欠損値を 0 で置き換える
for df in datasets:
    for col in cols:
        df[col] = df[col].fillna(0)

# 訓練データとテストデータのデータフレームをリストにまとめる
datasets = [df_train, df_test]

# 欠損値を最頻値で置き換える列のリスト
cols = ['MasVnrType', 'MSZoning', 'Functional', 'Electrical', 'KitchenQual', 'Exterior2nd', 'Exterior1st', 'SaleType']

# 各データセットに対して、指定した列の欠損値を最頻値で置き換える
for df in datasets:
    for col in cols:
        # 最頻値を取得し、欠損値をそれで置き換える
        df[col] = df[col].fillna(df[col].mode()[0])

# 訓練データとテストデータのデータフレームをリストにまとめる
datasets = [df_train, df_test]

# 'Utilities' 列を削除する
for df in datasets:
    df = df.drop('Utilities', axis=1)

# DataFrameを結合
all_data = pd.concat([df_train, df_test])

# 欠損値の数を計算
missing_values_count = df.isnull().sum()

# 欠損値の数とパーセンテージを表形式で表示するためのDataFrameを作成
missing_values_table = pd.DataFrame({
    'Missing_total': missing_values_count,
    'Percent (%)': round((missing_values_count / len(df)) * 100, 1)
})

# 欠損値がある列の情報を抽出して降順でソート
df_missing = missing_values_table[missing_values_table['Missing_total'] > 0].sort_values(by='Missing_total', ascending=False)
    
# 欠損値がある列の情報を表示
display(df_missing)

# GrLivArea（リビングエリアの面積）列で降順にソートし、最大値の2つの行のインデックスを取得
del_idx = df_train.nlargest(2, 'GrLivArea').index

# 削除前のデータフレームの形状を表示
print(f'before deleting: {df_train.shape}')

# del_idxで指定された行を削除
df_train = df_train.drop(del_idx, axis=0)

# 削除後のデータフレームの形状を表示
print(f'after deleting: {df_train.shape}')

# データセットリストを作成
datasets = [df_train, df_test]

# ラベル付けの対応を定義
label_mapping = {'Ex': 5, 'Gd': 4, 'TA': 3, 'Fa': 2, 'Po': 1, 'None': 0}  # ラベル対応の辞書

# ラベリングを行う列のリスト
cols = ['ExterQual', 'ExterCond', 'BsmtQual', 'BsmtCond', 'HeatingQC', 'KitchenQual',
       'FireplaceQu', 'GarageQual', 'GarageCond', 'PoolQC']

# データセットごとにラベリングを適用
for df in datasets:
    for col in cols:
        # 指定した辞書に基づいてラベルをマッピングし、整数型に変換
        df[col] = df[col].map(label_mapping).astype(int)

# データセットリストを作成
datasets = [df_train, df_test]

# ラベル付けの対応を定義
label_mapping = {'Gd': 4, 'Av': 3, 'Mn': 2, 'No': 1, 'None': 0}  # ラベル対応の辞書

# ラベリングを行う列のリスト
cols = ['BsmtExposure']

# データセットごとにラベリングを適用
for df in datasets:
    for col in cols:
        # 指定した辞書に基づいてラベルをマッピングし、整数型に変換
        df[col] = df[col].map(label_mapping).astype(int)

# データセットリストを作成
datasets = [df_train, df_test]

# ラベル付けの対応を定義
label_mapping = {'Fin': 3, 'RFn': 2, 'Unf': 1, 'None': 0}  # ラベル対応の辞書

# ラベリングを行う列のリスト
cols = ['GarageFinish']

# データセットごとにラベリングを適用
for df in datasets:
    for col in cols:
        # 指定した辞書に基づいてラベルをマッピングし、整数型に変換
        df[col] = df[col].map(label_mapping).astype(int)

# データセットリストを作成
datasets = [df_train, df_test]

# ラベル付けの対応を定義
label_mapping = {'GLQ': 6, 'ALQ': 5, 'BLQ': 4, 'Rec': 3, 'LwQ': 2, 'Unf': 1, 'None': 0}  # ラベル対応の辞書

# ラベリングを行う列のリスト
cols = ['BsmtFinType1', 'BsmtFinType2']

# データセットごとにラベリングを適用
for df in datasets:
    for col in cols:
        # 指定した辞書に基づいてラベルをマッピングし、整数型に変換
        df[col] = df[col].map(label_mapping).astype(int)

# SalePriceの対数変換
df_train['SalePrice'] = np.log(df_train['SalePrice'])

# 分析対象のデータを抽出して変数"data"に格納
data = df_train['SalePrice']

# サブプロットを作成
fig, ax = plt.subplots(1, 2, figsize=(10, 4))

# 歪度と尖度を計算して表示
skewness = data.skew()
kurtosis = data.kurtosis()
print(f'skew: {skewness:.2f}')  # 歪度を表示
print(f'kurt: {kurtosis:.2f}')  # 尖度を表示

# ヒストグラムと正規分布曲線の表示
sns.histplot(data, stat='density', kde=True, ax=ax[0])
ax[0].set_title('Histogram and Normal Distribution')

xmin, xmax = ax[0].get_xlim()  # ax[0]の最大値と最小値を取得
x = np.linspace(xmin, xmax, 100)
p = stats.norm.pdf(x, np.mean(data), np.std(data))
ax[0].plot(x, p, 'k', linewidth=1)

# 正規確率プロットの表示
res = stats.probplot(data, plot=ax[1])
ax[1].set_title('Normal Probability Plot')

plt.tight_layout()  # グラフのレイアウトを調整
plt.show()

# データセットのリストを作成
datasets = [df_train, df_test]

# 各データフレームに対して処理を行うループ
for df in datasets:
    # 変更前のデータフレームの形状を表示
    print(f"Before making: {df.shape}")
    
    # 新しい特徴量 'TotalSF' を追加し、'TotalBsmtSF'、'1stFlrSF'、'2ndFlrSF' を合計して代入
    df['TotalSF'] = df['TotalBsmtSF'] + df['1stFlrSF'] + df['2ndFlrSF']
    
    # 変更後のデータフレームの形状を表示
    print(f"After making: {df.shape}")

from scipy.stats import skew

# データの結合
ntrain = len(df_train)
all_data = pd.concat([df_train, df_test])

# 数値変数を特定
numeric_feats = all_data.dtypes[all_data.dtypes != "object"].index

# 各数値変数の歪度を計算し、降順にソート
skewed_feats = all_data[numeric_feats].apply(lambda x: skew(x.dropna())).sort_values(ascending=False)

# 歪度を含むデータフレームの作成
skewness = pd.DataFrame({'Skew' :skewed_feats})

# 最も歪んでいる数値変数の上位10件を表示
display(skewness.head(10))

# 歪度の絶対値が0.75より大きい特徴量を選択
skewness = skewness[abs(skewness) > 0.75]

# 歪度が0.75より小さい特徴量を削除
skewness = skewness.dropna(axis=0)

# Box-Cox変換が必要な歪んだ数値特徴量の数を表示
print(f"There are {skewness.shape[0]} skewed numerical features to Box Cox transform")

from scipy.special import boxcox1p

# Box-Cox変換を適用する数値変数の列を特定
indices = skewness.index

# 新しい列を格納するための空のリストを作成
skewness_after_boxcox = []

# Box-Cox変換のパラメータ（ラムダ）を設定
lam = 0.15

# 各数値変数に対してBox-Cox変換を適用し、新しい列として追加
for idx in indices:
    transformed_values = boxcox1p(all_data[idx], lam)
    all_data[idx] = transformed_values
    # Box-Cox変換後の歪度を計算し、リストに追加
    skewness_after_boxcox.append(skew(transformed_values))

# skewnessデータフレームにBox-Cox変換後の歪度を追加
skewness['Skew_BoxCox'] = skewness_after_boxcox

# Box-Cox変換後の歪度を表示
display(skewness.head(10))

# トレーニングデータの行数を取得
ntrain = len(df_train)

# ワンホットエンコーディングを実行
all_data = pd.get_dummies(all_data)

# トレーニングデータとテストデータに分割
df_train = all_data[:ntrain]  # 最初のntrain行はトレーニングデータ
df_test = all_data[ntrain:]   # 残りはテストデータ

# テストデータから 'SalePrice' 列を削除
df_test = df_test.drop(['SalePrice'], axis=1)

# トレーニングデータの形状を表示
print(f'df_train.shape: {df_train.shape}')  # テストデータの行数と列数を表示
display(df_train.head(3))  # テストデータの最初の3行を表示

# テストデータの形状を表示
print(f'df_test.shape: {df_test.shape}')  # テストデータの行数と列数を表示
display(df_test.head(3))  # テストデータの最初の3行を表示

# 必要なライブラリをインポート
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, Ridge, Lasso
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
from lightgbm import LGBMRegressor
from sklearn.metrics import mean_squared_error

# 説明変数と目的変数を選択
X = df_train.drop(['Id', 'SalePrice'], axis=1)
y = df_train['SalePrice']

# データの分割
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# モデルの構築と学習
models = [
    LinearRegression(),
    Ridge(alpha=1.0),
    Lasso(alpha=1.0),
    DecisionTreeRegressor(),
    RandomForestRegressor(),
    XGBRegressor(),
    LGBMRegressor()
]

# 結果を保存するためのデータフレームを作成
results = pd.DataFrame(columns=['Model', 'MSE'])

# 各モデルに対して学習と評価を行うループ
for model in models:
    model_name = model.__class__.__name__  # モデルの名前を取得
    model.fit(X_train, y_train)  # モデルの学習
    y_pred = model.predict(X_test)  # 予測
    
    mse = mean_squared_error(y_test, y_pred)  # 平均二乗誤差を計算
    
    # 結果をデータフレームに登録
    results = pd.concat([results, pd.DataFrame({'Model': [model_name], 'MSE': [mse]})], ignore_index=True)

# MSEによって結果をソート
results = results.sort_values(by='MSE')

# 結果を出力
display(results)

from sklearn.model_selection import GridSearchCV

# 開始値（最小値）と終了値（最大値）を指定
start_alpha = 0.01
end_alpha = 100

# 生成する値の個数を指定
num_values = 20

# 対数スケールでalphaの値を生成
alphas = np.logspace(np.log10(start_alpha), np.log10(end_alpha), num=num_values)

# Ridge回帰モデルを作成
ridge = Ridge()

# グリッドサーチのパラメータグリッドを定義
param_grid = {'alpha': alphas}

# クロスバリデーションとグリッドサーチを組み合わせて最適なalphaを探索
grid_search = GridSearchCV(ridge, param_grid, cv=5, scoring='neg_mean_squared_error')
grid_search.fit(X_train, y_train)

# 最適なalphaを表示
best_alpha = grid_search.best_params_['alpha']
print(f"Best alpha: {best_alpha}")

# 最適なalphaでRidge回帰モデルを再訓練
model = Ridge(alpha=best_alpha).fit(X, y)

# テストデータに対して予測を行い、exp(x) - 1を適用して元のスケールに戻す
sub_pred = np.expm1(model.predict(df_test.drop('Id', axis=1)))

# 提出用のデータフレームを作成
submission = pd.DataFrame({'Id': df_test['Id'], 'SalePrice': sub_pred})

# 予測結果をCSVファイルに保存
submission.to_csv("submission.csv", index=False)