ランダムフォレスト回帰【機械学習】

ランダムフォレストの回帰タスクについてまとめました。

:::note info
全く同じ記事をZennでも公開しています。
:::
https://zenn.dev/0w0kaomoji/articles/7a1218c23ba357

## 1. ランダムフォレストとは

- 複数の**決定木**（というモデル）を組み合わせる**アンサンブル手法**を用いたモデルです。
    
    ```mermaid
    graph TD
        A[入力データ] --> B[ブートストラップ]
        B --> C1[決定木 1]
        B --> C2[決定木 2]
        B --> Cn[決定木 3]
    
        C1 --> D[予測値]
        C2 --> D
        Cn --> D
    
        D --> E[予測値の平均]
        E --> F[最終的な出力]
    
    ```
    
- **クラス分類**や**回帰**に用いられます。
    
    
    :::note info
    **決定木(decision tree)**
    
    - 質問を重ねることで、**データを分割していく**モデル
    - 図１：予測値は区分線形関数として表される
        
        ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3428921/23a8c4a4-1112-f184-b471-1e496e1c2f92.png)

        
    - 図２：MSEが最小となる変数の閾値を設定していく（LSTATは変数名）
        
        ![image 1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3428921/885c234a-3541-a975-afb6-f20791d2dc4a.png)
        
:::

## 2. メリット

- アンサンブル手法なので、決定木よりも**外れ値**に対して頑健です。また**バギング**により、モデルのバイアスは変わらないまま**分散**が小さくなります。
    
    ![image 2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3428921/cecb9d7e-ad3a-b0a1-c452-c996a1561668.png)

    
    :::note info
    **バギング**(Bootstrap aggregating/Bagging)
    
    - **ブートストラップ**（後述）によりm個の訓練データセットを用意し、それぞれを学習させたあと、**m個の予測値を平均**して最終的な出力とするアンサンブル手法
    - クラスタリングでは、m個の予測値から多数決によって最終的な出力が決まる
:::

- 決定木を用いるので、**変数変換**に対しても頑健です。
- **欠測値**を推定してくれるので、埋める必要がありません。
- 分岐の繰り返しによって説明変数の**相互作用**を考慮するので、明示的に与える必要がありません。

## 3. アルゴリズム

1. 各決定木で、データを**ブートストラップ**します。
    
    :::note info
    **ブートストラップ(Bootstrap)**
    
    - サンプル$A$から、元のサンプルサイズと**同じサイズ**のサンプル$B$を**復元抽出**すること
    

2. 各ノードで、以下の操作を情報利得が**変化しなくなるまで**、繰り返します。
    1. 分岐に用いる**説明変数の候補**を選択します。
    2. 各ノードで、**分岐条件**を決定します。
        - 情報利得が**最大になる**ような分岐条件を求めます。
        - つまり、目的は子ノードの**不純度の総和を最小にする**ことです。
        :::note info
        
        **情報利得(information gain, IG)**
        
        - **親ノードの不純度**と、**子ノードの不純度の総和**の差
        - 各変数を以下のようにおく。
            
            $D_p$：親ノードに含まれるサンプル
            
            $N_p$：親ノードのサンプル数
            
            $D_j$：$j$番目の子ノードに含まれるサンプル
            
            $N_j$：$j$番目の子ノードのサンプル数
            
            $I(\cdot)$：不純度
            
            $f$：分割を行う説明変数
            
            このとき、**情報利得**は
            
            $$
            IG(D_p,f)=I(D_p)-\sum_{j=1}^{M}\frac{N_j}{N_p}I(D_j)
            $$
            
            と表される。      
    
        :::note info
        
        **不純度(impurity)**
        
        - 回帰では、**平均二乗誤差(MSE)**を用いる
            - 各変数を、以下のようにおく。
                
                $\hat{y_t}$：予測値
                
                $y^{(i)}$：真の目的変数
                
                このとき、**ノード$t$におけるMSE**は
                
                $$
                I(t)=MSE(t)=\frac{1}{N_t}\sum_{i\in D_t}(\hat{y_t}-y^{(i)})^2
                $$
                
                と表される。
                
                ここで、予測値は$D_t$の目的変数の平均とする。
                
                すなわち、
                
                $$
                \hat{y_t}=\frac{1}{N_t}\sum_{i\in  D_t}y^{(i)}
                $$
                  
       
        
3. 全ての木における**予測値を平均**します。
    - 最終的な予測値$\hat{y}^{(i)}$は以下のように表されます。
        
        $$
        \hat{y}^{(i)}=\frac{1}{M}\sum_{j=1}^{M}\hat{y}_j^{(i)}
        $$
        
        $M$：木の本数
        
        $\hat{y}_j^{(i)}$：木$j$における、$y^{(i)}$の予測値
        

## 4. 他のモデルとの違い

- **GBDT（勾配ブースティング決定木）との違い**
    - GBDTでは決定木を直列に並べるのに対して、ランダムフォレストでは**並列**に並べます。そのため、決定木の本数を増やすことで精度が悪くなることはありません。
    - アンサンブル手法としてGBDTではブースティングを用いるのに対し、ランダムフォレストでは**バギング**を用いています。
- **ランダムフォレストクラスタリング**との違い
    - 不純度の指標として、クラスタリングではジニ不純度やエントロピーを用いるのに対して、回帰では**平均二乗誤差**を使用します。
    - クラスタリングでは予測値は多数決により決定されますが、回帰では**平均**で計算されます。

## 5. サンプルコード

- ライブラリは**Scikit-learn**を使用します
    - Section 3 アルゴリズムの2-aにあるように、説明変数の候補を選択するステップがありますが、Scikit-learnのデフォルトの候補は**すべての説明変数** となっています。つまり、全ての説明変数が候補として選択される確率は１です。
    - 図２からもわかるように、Scikit-learnでは決定木は２分木で実装されています。
- コードは以下のようになります。
    
    ```python
    # dfについては補足にある図１のコードを参照
    X = df.iloc[:, :-1].values  # 特徴量データ（最後の列を除く）
    y = df['MEDV'].values       # 目的変数（MEDV）
    from sklearn.model_selection import train_test_split
    
    # データ分割（40%をテストデータに）
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=1)
    
    from sklearn.ensemble import RandomForestRegressor as RFR
    
    # ランダムフォレストモデルの設定
    forest = RFR(criterion='squared_error', # 不純度基準
                 random_state=1,            # 乱数シード
                 n_jobs=-1)                 # 並列処理
    
    forest.fit(X_train, y_train)  # モデルの学習
    
    # 予測値の取得
    y_train_pred = forest.predict(X_train)
    y_test_pred = forest.predict(X_test)
    
    from sklearn.metrics import mean_squared_error, r2_score
    
    # MSEとR^2スコアの表示
    print('MSE train: %.3f, test: %.3f' % (mean_squared_error(y_train, y_train_pred),
                                           mean_squared_error(y_test, y_test_pred)))
    print('R^2 train: %.3f, test: %.3f' % (r2_score(y_train, y_train_pred),
                                           r2_score(y_test, y_test_pred)))
    
    import matplotlib.pyplot as plt
    
    # 残差プロット
    plt.scatter(y_train, y_train_pred - y_train, 
                c='steelblue', edgecolor='white', marker='o', 
                s=35, alpha=0.9, label='training_data')
    plt.scatter(y_test, y_test_pred - y_test, 
                c='palevioletred', edgecolor='white', marker='s', 
                s=15, alpha=0.9, label='test_data')
    plt.xlabel('True values')
    plt.ylabel('Residuals')
    plt.legend(loc='upper right')
    plt.title('Residuals')
    plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='black')
    plt.tight_layout()
    plt.show()
    
    # 特徴量の重要度
    feature_importances = forest.feature_importances_
    
    import numpy as np
    
    # 特徴量重要度の棒グラフ
    plt.figure(figsize=(10, 5))
    y = feature_importances
    x = np.arange(len(y))
    plt.bar(x, y, align="center")
    plt.xticks(x, df.columns[:-1])
    plt.xlabel('Features')
    plt.ylabel('Importance')
    plt.title('Feature Importance')
    plt.tight_layout()
    plt.show()
    ```
    
    - 出力
      
        MSE train: 5.544, test: 15.618
        R^2 train: 0.927, test: 0.844
        
        ![image 3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3428921/954385ce-a52d-1852-7e9a-da5f1c65f4a2.png)

        
        ![image 4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3428921/e3c41039-e8df-57c0-c76f-a51c9d92542e.png)


## 6. 結果の解釈

- MSE
    - **過学習**する傾向にあることがわかります。
- 決定係数
    - モデルは**かなり良く**当てはまっていることがわかります。
- 残差プロット
    - 明らかに負の相関関係が見えることから、モデルが**データの情報をとらえきれていない**ことがわかります。
- 各説明変数の重要度
    - 特に“RM”と"LSTAT"という説明変数が、モデルに**大きな影響**を与えていることがわかります。
## 7. 感想

ランダムフォレストは、その構造がシンプルであるにもかかわらず、予想以上に高い精度を示すことに感心しました。また、モデルの構築プロセスが比較的容易であることも、実務上大きな利点だと思いました。

一方で、モデルの改善方法が明確でないという課題があり、その点では、統計的モデルのほうが扱いやすそうです。
## 参考文献

- Raschka, Sebastian, Vahid Mirjalili, 株式会社クイープ, and 福島真太朗. 2018. *［第2版］Python機械学習プログラミング 達人データサイエンティストによる理論と実践*. インプレス.
- 門脇大輔, 阪田隆司, 保坂桂佑, and 平松雄司. 2019. *Kaggleで勝つデータ分析の技術*.

https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html

https://scikit-learn.org/stable/auto_examples/ensemble/plot_forest_hist_grad_boosting_comparison.html#sphx-glr-auto-examples-ensemble-plot-forest-hist-grad-boosting-comparison-py

https://ja.wikipedia.org/wiki/%E3%83%90%E3%82%AE%E3%83%B3%E3%82%B0

https://bellcurve.jp/statistics/course/9706.html?srsltid=AfmBOopgESK3d1ZiHsbUDri-qsqmjQt-QqsOj4eOQLeWCv7Q2LhShub1

https://funatsu-lab.github.io/open-course-ware/machine-learning/random-forest/

https://dropout009.hatenablog.com/entry/2021/07/26/193907

- アンサンブルとバイアス、分散との関係を説明している

https://amalog.hateblo.jp/entry/decision-tree-scaling
- 変数変換に対する頑健性を検証している

https://qiita.com/g-k/items/e7dcc4d2b057dada405c
- ランダムフォレストの欠損値補完を説明している

https://blog.data-hacker.net/2020/07/randomforest.html
- 相互作用に関する検証をしている

## 補足

- 図１のコード
    
    ```python
    import pandas as pd
    import numpy as np
    import matplotlib.pyplot as plt
    from sklearn.tree import DecisionTreeRegressor
    def lin_regplot(X_test,X, y, model):
        plt.scatter(X, y, c='steelblue', edgecolor='white', s=70,label="Data")
        y_pred=model.predict(X_test)
        plt.plot(X_test, y_pred, color='orange', lw=2, label="Prediction")
        for i in range(1, len(y_pred)):
          if y_pred[i] != y_pred[i - 1]:  
              plt.axvline(x=X_test[i], color='green', linestyle='--', linewidth=1)
        
        return 
    
    df = pd.read_csv('https://raw.githubusercontent.com/rasbt/'
                     'python-machine-learning-book-2nd-edition'
                     '/master/code/ch10/housing.data.txt',
                     header=None,
                     sep='\s+')
    
    df.columns = ['CRIM', 'ZN', 'INDUS', 'CHAS', 
                  'NOX', 'RM', 'AGE', 'DIS', 'RAD', 
                  'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV']
    
    df=df.sample(100)
    
    X = df[['LSTAT']].values
    y = df['MEDV'].values
    
    tree = DecisionTreeRegressor(max_depth=2)
    tree.fit(X, y)
    
    X_test = np.arange(0, 40, 0.01)[:, np.newaxis]
    
    lin_regplot(X_test,X,y, tree)
    plt.xlabel('% lower status of the population [LSTAT]')
    plt.ylabel('Price in $1000s [MEDV]')
    plt.legend()
    plt.show()
    ```
    
- 図２のコード
    
    ```python
    from sklearn.tree import export_graphviz
    from graphviz import Source
    from sklearn.tree import plot_tree
    dot_data = export_graphviz(tree, 
                               filled=True, 
                               rounded=True,
                               feature_names=['LSTAT'], 
                               out_file=None)
    graph = Source(dot_data) 
    graph.format = 'png'
    graph.render('tree')
    ```