「スモールデータ解析と機械学習」を寄り道写経～第４章「線形回帰モデルにおける入力変数選択」

2024年3月30日 11:47

第４章「線形回帰モデルにおける入力変数選択」

書籍の著者　藤原幸一先生

この記事は、テキスト「スモールデータ解析と機械学習」第４章「線形回帰モデルにおける入力変数選択」の通称「寄り道写経」を取り扱います。

テキスト4.7節「NIRスペクトルの検量線入力波長選択」の寄り道写経に取り組みます。
ではテキストを開いてスモールデータの旅に出発です🚀

このシリーズは書籍「スモールデータ解析と機械学習」（オーム社、「テキスト」と呼びます）の機械学習の理論・数式とPythonプログラムを参考にしながら、テキストにはプログラムの紹介が無いけれども気になったテーマ、または、テキストのプログラム以外の方法を試したいテーマを「実験的」にPythonコード化する寄り道写経ドキュメンタリーです。

はじめに

テキスト「スモールデータ解析と機械学習」のご紹介

テキストは、2022年2月に発売され、スモールデータと呼ばれる小さなサンプルサイズのデータを用いた機械学習の取り組み方について、理論（数式）とPythonサンプルプログラムを提案する素晴らしい実用書です。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「スモールデータ解析と機械学習」第１版第１刷、著者藤原幸一、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

第４章線形回帰モデルにおける入力変数選択

主に Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。

NIRスペクトルの検量線入力波長選択

4.7節「NIRスペクトルの検量線入力波長選択」の寄り道写経に取り組みます。
テキストは、前回利用した「ディーゼル燃料のNIRスペクトルデータ」を用いて、次のアルゴリズムによる説明変数選択を実行するシナリオになっています。
・Lasso
・PLS-beta（20%と50%）
・VIP法
・NCSC-VS

寄り道写経的には、①テキストオリジナルコードの誤植・ワーニング対策と、②図4.17「入力波長選択結果」の一括描画の実装に取り組みます！

なお、この記事のコードを動かすには、テキストの処理の流れに沿って、テキストオリジナルコードを実行しておく必要があります。
「テキストオリジナルコード」はオーム社HPからダウンロードできる書籍掲載のプログラムであり、書籍購入者のみ利用できるとのことです。

ひとまず共通的に利用するライブラリをインポートします。

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

① 誤植・ワーニング対策

プログラム4.8「Lassoによる検量線入力波長選択」からプログラム4.10「NCSC-VSによる検量線入力波長選択」で発生するワーニングの対策案を考えてみます。

■ ワーニング発生例
前回記事でも取り上げた「RMSEおよび相関係数を算出する関数linear.pred_eval」で発生するワーニングの対策です。
簡単にワーニング内容、対策を書きます。

【ワーニングメッセージ（部分）】
RuntimeWarning: Degrees of freedom <= 0 for slice
RuntimeWarning: divide by zero encountered in divide
RuntimeWarning: invalid value encountered in multiply

【対策の概要】
RMSEおよび相関係数を算出する関数「linear.pred_eval」に与える目的変数の正解値 yval と予測値 y_hat を ravel() flatten() で平坦化します。

【コード変更（案）】

rmse, r = linear.pred_eval(yval.ravel(), y_hat.ravel(), plotflg=True) # ★ravel

■ 誤植例１
プログラム4.9「PLS-betaとVIP法による検量線入力波長選択」のPLS-beta(20%)にかかる誤植です。
【該当箇所】

# PLS-beta(20%)
beta, _, _, _, _, _, _ = simpls(X_, y_, R=7) # 前章のクロスバリデーションの結果より
sel_var_beta20 = pls_beta(beta, share=0.50)

【コード変更（案）】
３行目の関数 pls_beta の引数 shareを 0.20 に変更します。

# PLS-beta(20%)
beta, _, _, _, _, _, _ = simpls(X_, y_, R=7) # 前章のクロスバリデーションの結果より
sel_var_beta20 = pls_beta(beta, share=0.20)  # ★shareは0.50ではなく0.20

■ 誤植例２
プログラム4.9「PLS-betaとVIP法による検量線入力波長選択」のVIP法にかかる誤植です。
【該当箇所】

# VIP
beta, W, P, Q, T, U, cont = simpls(X_, y_, R=7)
sel_var_vip, vip = pls_vip(W, D, T, threshold=1)

【コード変更（案）】
この変更内容が正解かどうかは分かりません。
・２行目で利用する関数を simpls から nipals_pls1 へ変更します。
・３行目の関数 pls_vip の引数 threshold を削除します

# VIP
beta, W, _, D, T = nipals_pls1(X_, y_, R=7)  # ★nipals_pls1に変更
sel_var_vip, vip = pls_vip(W, D, T)          # ★threshold引数を削除

② おまけ

プログラム4.8「Lassoによる検量線入力波長選択」からプログラム4.10「NCSC-VSによる検量線入力波長選択」で算出する「潜在変数の数 optLV」と、pred_eval 関数実行時に「plotflg=True」にして描画したグラフを共有します。
なお、算出される値は実行のたびに変化するのでご留意ください。

■ Lasso
相関係数 r はテキストの 0.78 より大きく、RMSEはテキストの 2.21 より小さいです。

■ PLS-beta(20%)
潜在変数の数はテキストの 12 と同じです。
相関係数 r はテキストの 0.81 と同じ、RMSEはテキストの 2.11 と同じです。

■ PLS-beta(50%)
潜在変数の数はテキストの 7 より小さいです。
相関係数 r はテキストの 0.83 より小さく、RMSEはテキストの 2.00 より大きいです。

■ VIP
潜在変数の数はテキストの 7 と同じです。
相関係数 r はテキストの 0.83 と同じ、RMSEはテキストの 2.00 と同じです。

■ NCSC-VS
潜在変数の数はテキストの 8 より小さいです。
相関係数 r はテキストの 0.84 より小さく、RMSEはテキストの 1.95 より大きいです。

③ 図4.17 入力波長選択結果の一括描画

プログラム4.11「選択された入力波長の可視化」を引用して、PLS-beta(50%)、VIP法、Lasso、NCSC-VSを一括で可視化するコード改造を行います。

## 描画用データの作成 ★NCSC-VSの選択変数数がテキストよりもかなり多い
# Lassoの選択された(とみなした)変数のインデックス取得
sel_var_lasso = np.where(beta_lasso.ravel() > 0.001)[0]
# 4つのsel_varをリスト化
sel_var = [sel_var_beta50, sel_var_vip, sel_var_lasso, sel_var_ncsc]
# 4つのタイトルをリスト化
title = ['PLS-beta(50%)', 'VIP', 'Lasso', 'NCSC-VS']
# x軸の波長の値の生成
wave = np.arange(750, 1551, 2)

## 描画処理
fig, ax = plt.subplots(2, 2, figsize=(10, 7))
for i in range(len(sel_var)):
    # axesの位置の算出
    pos = divmod(i, 2)
    # X_trainの値の描画（青線）
    for n in range(X_train.shape[0]):
        ax[pos].plot(wave, X_train[n, :], color='royalblue', alpha=0.05)
    # 選択された変数の描画（薄赤色）
    xmin, xmax = ax[pos].get_xlim()
    ymin, ymax = ax[pos].get_ylim()
    ax[pos].vlines(wave[sel_var[i]], ymin, ymax, lw=0.5, color='tomato',
                   alpha=0.2)
    # 選択された変数の数を表示
    num_text = f'選択された変数の数: {len(sel_var[i])}'
    posx = (xmax - xmin)*0.1 + xmin
    posy = (ymax - ymin)*0.85 + ymin
    ax[pos].text(posx, posy, num_text)
    # 修飾
    ax[pos].set(xlabel='波長[nm]', ylabel='スペクトル強度', title=title[i])

plt.tight_layout();
plt.show()

【実行結果】
テキストの可視化結果と比べると、今ひとつ、観測値にピークを捉えられていない印象です。
NCSC-VSの選択変数数は多すぎますね・・・（汗）

第４章の寄り道写経は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルはＲ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

この記事が参加している募集

新生活をたのしく

39,047件

この記事が気に入ったらサポートをしてみませんか？