データの「分類」(classification)と「回帰」(regression)が機械学習の2大テーマであるが、実験観測データの回帰分析がもっとも使う頻度が高いであろうと考え、実践的なデータサイエンスとして、回帰問題に対する機械学習の手法を中心に事例、プログラム例をもとに講義、演習を行った。
線形単回帰
説明変数$x$と目的変数$y$の関係を$ y= w_0 + w_1 x$として、係数$w_0$, $w_1$を最小2乗法で推定
線形重回帰
ボストンの住宅価格と各地点の特徴データをもとに、説明変数が多数ある場合(説明変数$\{x_1, x_2, \cdots , x_M\}$について
$$y = w_0 +w_1 x_1 + w_2 x_2 + \cdots + w_M x_M$$
なる関係があると想定して係数を推定
非線形回帰
$x$と$y$が非線形な関係にある場合、べき級数近似
$$y = w_0 + w_1x + w_2x^2 + \cdots + w_M x^M$$
として係数を求める問題は説明変数を$\{x, x^2,\cdots , x^M\}$とすれば線形重回帰と同様である。
データ例 :$y=f(x)$の関係を求めるために、$x$から各べきの値を説明変数として並べる(numpyのvanderを用いて作成)
(追加) sckit-learnにはPolynomialFeaturesというメソッドもあるのでそれを利用してもよい。 https://cpp-learning.com/polynomial-regression/
過学習
多項式近似では、$M$が大きい(データ数と同程度)の場合は、過学習がおきる。
ridge回帰
$$ \tilde{E}({\bf w}) = \frac{1}{2} \sum_{n=1}^N \{y(x_n,{\bf w}) - t_n\}^2 + \frac{\lambda}{2}||{\bf w}||^2 $$
を最小にする$\{w_i,\ i=1,\cdots, M\}$を求める。
べき多項式以外のフィッティングへの一般化
非線形な関係のフィッテングは、$\{\phi_1(x) + \phi_2(x) + \cdots + \phi_M(x)\}$を適当に選び
$$ y(x, {\bf w}) = w_1 \phi_1(x) + w_2 \phi_2(x) + \cdots + w_M \phi_M(x) $$のように一般化できる。
特に、区間$[a,b]$の間のフィッティングを、区間を$M+1$個で割り、そこをピークするガウス関数
$$ \phi_i (x) = \exp \left( - \frac{(x - \mu_i)^2}{2\sigma^2}\right),\quad \mu_i = \frac{b-a}{M}i,\quad i=1,2,\cdots M $$を使うことが多い。これをガウス基底と呼ぶ。
パラメータの事前分布を仮定し、データが追加されるごとに逐次学習してパラメータの確率分布をベイズの定理に従って計算する手法。
(PRMLの図1.16):
ガウシアン基底を用いたベイズ統計回帰モデルの表式から、学習は、基底関数$\phi$の前の係数の確率分布ではなく、基底関数と係数をカーネル関数という形にまとめられることをみた。
$$ y(\boldsymbol{x}_i)=\sum_{k=1} w_k \phi_k (\boldsymbol{x}_i) +b $$例えば、RBF(ガウシアン)カーネル
$$ k(\boldsymbol{x},\boldsymbol{x'}) = \exp ( -\gamma ||\boldsymbol{x}- \boldsymbol{x'}||^2) $$非線形なサポートベクトル識別
SVMの回帰への応用: SVR
いくつかの改良
ニューラルネット(パーセプトロン)
決定木のアンサンブル学習(ランダムフォレスト)
データの正規化(reguralization)
from sklearn.preprocessing import StandardScaler
ハイパーパラメータの最適値探索 交差検定
scikitlearnのGridSearchの利用
import sklearn.model_selection as ms
est = ms.GridSearchCV(svm.LinearSVC(),
{'C': np.logspace(-3., 3., 10)})
自然科学(Natural Science)
⇒ 原理に基づくモデリングと数値実験(Numerical Simulation)
通底するのは数学 (その意味で、数学は自然科学の範疇には入らない)
とはいえ、創薬、大規模加速器実験や天文データ分析、複雑な機能性化合物の探索など自然科学ベースの対象物の探索にも用いられるようになり、一言では言い表せない状況も出現している。
目的変数、説明変数がそれぞれ一つのという単純な場合についてではあるが、本講義で取り上げた手法を並べて比較してみる。
データとしては、甲府の年平均気温の推移を用いる。
世界の年平均気温については、 https://toyoki-lab.ee.yamanashi.ac.jp/~toyoki/lectures/PracDataSci/Examples.html
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
plt.rcParams['font.family'] = 'IPAPGothic' # 全体のフォントを設定
# データの読み込み(月ごとの平均気温)
df_temperature = pd.read_csv("data/temperature_data_kofu.csv")
df_temperature['date'] = pd.to_datetime(df_temperature['年月'])
df_temperature['year'] = df_temperature['date'].dt.year
# 1920以降のデータを採用
origin_year=1920
df_temperature = df_temperature[df_temperature['year']>=origin_year]
# 年平均の算出
df_annual = df_temperature[["year","平均気温(℃)"]].groupby('year').mean().reset_index()
# 1945を除く
df_annual = df_annual[df_annual["year"]!=1945]
df_annual.plot(kind="line", x="year", y="平均気温(℃)")
<AxesSubplot: xlabel='year'>
scikit-learnの各種モデル(学習機械)のメソッドには統一感があり、プログラムをあまり書き換えることなく様々なモデルを試せます。
しかし、それぞれのモデルの特徴から、前処理に工夫が必要な場合があります。上のような単一の$x$、単一の$y$が非線形に関係している例をつかって、まとめます。
import numpy as np
# M次多項式近似の場合は、べき乗の値を列とするM+1列の行列(計画行列)
deg = 5
poly_feature = np.vander((df_annual['year'] - origin_year).to_list(), deg+1)
# SVR, RVRはカーネルを指定するだけなので、1列だけでよい
single_feature = np.array([(df_annual['year'] - origin_year).to_list()]).T
# 説明変数の正規化 (Neural Netについては、正規化データを与えないとよい結果が得られない
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler().fit(poly_feature) #正規化
scaled_poly_feature = scaler.transform(poly_feature)
#single_feature
target = df_annual["平均気温(℃)"].to_list()
# テストデータ (ここでは、推定曲線を出力させるだけなので、featureをそのまま使ってもよい)
n=100
test_single = np.linspace(0, 2022-origin_year, n)
test_single = test_single[:, np.newaxis]
test_poly = np.vander(np.linspace(0, 2022-origin_year, n)[:], deg+1)
scaled_test_poly = scaler.transform(test_poly)
# scikit-learnの各種モデル
# {model, title(name), 説明変数, 目的変数, test(predict)用データ}の配列を用意する
models = {}
# 非線形回帰:多項式近似 (vander行列をデータとして入力 ⇒ 線形重回帰へ入力)
from sklearn import linear_model as lm
models['LM'] = {"model": lm.LinearRegression(),
"name": "Simple Polynomial Model " + str(deg) +"次",
"feature": poly_feature, "target": target, "test_data": test_poly}
# SVR
from sklearn import svm
models["SVR"] = {"model": svm.SVR(kernel='rbf'),
"name": "Support Vector Regression",
"feature": single_feature, "target": target, "test_data": test_single}
# RVR (pipでインストールする必要あり)
from sklearn_rvm import EMRVR
models["RVR"] = {"model": EMRVR(kernel='rbf', gamma="scale"),
"name": "Relevant Vector Regression",
"feature": single_feature, "target": target, "test_data": test_single}
# Neural Network
from sklearn.neural_network import MLPRegressor as MLPR
models["NN"] = {"model": MLPR((100,), activation="identity", solver="lbfgs",
max_iter=2000, tol=1e-8),
"name": "Neural Network",
"feature": scaled_poly_feature, "target": target, "test_data": scaled_test_poly}
# Random Forest
from sklearn.ensemble import RandomForestRegressor as RFR
models["RFR"] = {"model": RFR(max_depth=3, n_estimators=100), # 深さをいろいろ変えてみる
"name": "Random Forest Regression",
"feature": poly_feature, "target": target, "test_data": test_poly}
# Mixture model
# 各モデルの学習 ⇒ 推定 ⇒ 値をmodelsの1要素に格納
for k in models.keys():
models[k]['model'].fit(models[k]['feature'], models[k]['target'])
models[k]['predict'] = models[k]['model'].predict(models[k]['test_data'])
# 結果の描画
fig, ax = plt.subplots()
# 表示するものをインタラクティブに選択したかったのだが。。。
models2show = ["LM", "SVR", "RVR", "NN", "RFR"]
for m in models2show:
ax.plot(test_single+origin_year, models[m]['predict'], label=models[m]['name'])
ax.legend()
df_annual.plot(kind='scatter',x="year", y="平均気温(℃)", ax=ax)
<AxesSubplot: xlabel='year', ylabel='平均気温(℃)'>
mixture modelは「手作り」のものなので上の一覧には入れられない。
参考に結果の図だけ貼っておく。
%%html
<link rel="stylesheet" type="text/css" href="custom.css">