2.3. 機械学習ライブラリ scikit-learn

Scikit-learn は Pythonで最も有名な機械学習ライブラリでしょう。 主な特徴は以下のようなものです。

  • 一般的な機械学習アルゴリズムが網羅されている
  • ドキュメントがよく整備されており、素人でも簡単に利用できる
  • 簡単に理解できる平易な実装となっている
  • 計算負荷の高いところはC++で実装されているため比較的高速

機械学習というと、IT関係、人工知能関係のごく限られた内容を想定されるかもしれませんが、 そんなことはありません。 機械学習とは統計と言い換えてもほとんど問題ないでしょう。 実験データや計算データなど、データ解析に関わる問題であれば関わる可能性も高いと思います。

主に

  • 回帰問題
  • 教師有り分類問題
  • 教師無し分類問題
  • 次元圧縮・特徴抽出
  • モデル選択
  • データ前処理

に関する様々なアルゴリズムが実装されています。

2.3.1. scikit-learn のインストール

Anaconda distribution を用いて Python 開発環境を整えた人は

conda install scikit-learn

としてください。Native の Python を用いている人は

pip install scikit-learn

とするとよいでしょう。

2.3.2. scikit-learn の使用例

scikit-learn では様々な手法が利用できますが、ここでは線形回帰を紹介します。 線形回帰とは、計測データ y を基底関数 \(\phi(x)\) の線形和で近似するものです。

\[y_i \approx \sum_j w_j \phi_j(x_i)\]

ここで、 \(w_j\) は線形結合の係数です。 \(\phi(x)\) にどういった式を用いるかで、 多項式近似や、スプライン近似などを表すことができる汎用的なモデルです。 scikit-learn では linear_model モジュールとして提供されています。 ここでは、前節で紹介したトムソン散乱による電子温度分布を複数のガウス関数の和で近似することにします。

最も簡単な近似法は、最小二乗法です。最小二乗法では、

\[\sum_i \left(y_i - \sum_j w_j \phi_j(x_i)\right)^2\]

を最小化する \(w_j\) を求めます。 このような最小二乗法は、ノイズがガウス分布に従っていると仮定したものです。 しかし、一般にノイズが綺麗なガウス分布に従っているとは限りません。 例えば 前節に示したトムソン散乱による電子温度分布にも、 外れ値と言えるようなデータ点も見受けられます。

このような外れ値を含むようなデータの近似アルゴリズムも様々存在します。 scikit-learn でも複数種類実装されていますが、 ここではHuber回帰という手法を用いたものを紹介します。

Huber回帰では、通常の最小二乗法と異なり、以下のコスト関数を最小化します。

\[\sum_i \mathcal{H}\left(\frac{y_i - \sum_j w_j \phi_j(x_i)}{\sigma}\right)\]

ここで、 \(\sigma\) は外れ値を除いたノイズの分散で、これもデータから推定します。 \(\mathcal{H}(z)\) はHuberのロス関数で、以下のように定義されるものです。

\[\begin{split}\mathcal{H}(z) = \begin{cases} \; z^2 & |z| \le 1 \\ \; 2|z| - 1 & |z| > 1 \end{cases}\end{split}\]

簡単には、近似値の近くにある点については通常の二乗誤差を考え、 \(\sigma\) より遠く離れている点については一乗誤差に緩和してロス関数に加えるものです。 遠く離れている点についてのコストが小さくなるため、外れ値に引っ張られることが少なくなります。

ここでは、scikit-learn を用いて Huber回帰を行ってみます。

In [1]: from sklearn import linear_model  # linear_model モジュールを用います
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-c379e4b67c23> in <module>()
----> 1 from sklearn import linear_model  # linear_model モジュールを用います

ImportError: No module named 'sklearn'

# data ここでは 6800 msに得られた Te の分布を解析します
In [2]: Te = thomson.sel(Time=6800, method='nearest')['Te'].values

In [3]: R = thomson['R'].values

# basis R:2500--5000 を10分割した点を中心とするガウス関数の和で近似しましょう
In [4]: centers = np.linspace(2500, 5000, 10)

In [5]: phi = np.exp(-((R.reshape(-1, 1) - centers) / 200)**2)

# 最小二乗法
In [6]: lin = linear_model.LinearRegression(fit_intercept=False)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-5ab8d10dc023> in <module>()
----> 1 lin = linear_model.LinearRegression(fit_intercept=False)

NameError: name 'linear_model' is not defined

# フィッティング
In [7]: lin.fit(phi, Te)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-189a96eb5d16> in <module>()
----> 1 lin.fit(phi, Te)

NameError: name 'lin' is not defined

# 求めたフィッティング係数を用いた予測
In [8]: Te_lin_fit = lin.predict(phi)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-31a114ce840b> in <module>()
----> 1 Te_lin_fit = lin.predict(phi)

NameError: name 'lin' is not defined

# ロバスト最小二乗法
In [9]: rob = linear_model.HuberRegressor(fit_intercept=False)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-9-7ce2dde34e7a> in <module>()
----> 1 rob = linear_model.HuberRegressor(fit_intercept=False)

NameError: name 'linear_model' is not defined

# フィッティング
In [10]: rob.fit(phi, Te)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-6e5aa8defcd8> in <module>()
----> 1 rob.fit(phi, Te)

NameError: name 'rob' is not defined

# 求めたフィッティング係数を用いた予測
In [11]: Te_rob_fit = rob.predict(phi)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-11-053e520b255a> in <module>()
----> 1 Te_rob_fit = rob.predict(phi)

NameError: name 'rob' is not defined

In [12]: plt.plot(R, Te, '--o', ms=3, label='data')
Out[12]: [<matplotlib.lines.Line2D at 0x7fce70bba358>]

In [13]: plt.plot(R, Te_lin_fit, label='linear regression', lw=2)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-13-b30e1a51b311> in <module>()
----> 1 plt.plot(R, Te_lin_fit, label='linear regression', lw=2)

NameError: name 'Te_lin_fit' is not defined

In [14]: plt.plot(R, Te_rob_fit, label='huber regression', lw=2)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-14-e4e2b9b54c2d> in <module>()
----> 1 plt.plot(R, Te_rob_fit, label='huber regression', lw=2)

NameError: name 'Te_rob_fit' is not defined

In [15]: plt.legend(loc='best')  # 凡例を表示する
Out[15]: <matplotlib.legend.Legend at 0x7fce70bbaac8>

In [16]: plt.xlabel('$R$ (mm)')
Out[16]: Text(0.5,0,'$R$ (mm)')

In [17]: plt.ylabel('$T_\mathrm{e}$ (eV)')