★操作方法:

上から順に自分で実行したい場合は、

  1. メニューから Cell -> All Output -> Clear とする
  2. 上から順にコード部分をクリックして、メニューの実行ボタン(右向き三角ボタン)を押していく

逆に全部計算された状態のノートブックを得たい場合は、

  1. メニューから Cell -> Run All とする

です。

多項式近似とAICによるモデル選択の例

ここではデータ点を最小二乗法で多項式近似するに際して、多項式次数をAIC(赤池情報量基準)を用いて最適化してみます。

まず

  • matplotlib をipython内で使えるようにする (matplotlibなどで描写したグラフをノートの中にそのまま表示したりできる)
  • NumPy と matplotlib を import

を行います。

In [9]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

多項式近似される元データを作成します。$x$ として $0 \sim 3$ の区間を30等分して両端含め31点のベクトル値。 $y$ として $y \sim \sin x + \mathcal{N}(\mu=0, \sigma=0.1)$ というデータを作成します。

python の sin 関数はスカラー値(単一の数値)を引数にとってスカラー値を返しますが、NumPyで定義されたsin (np.sin) はベクトル値を引数にとって、各値に対する sin の値のベクトルを返します。

In [10]:
xs = np.linspace(0.0, 3.0, 31)
ys = np.sin(xs) + np.random.normal(0.0, 0.1, len(xs))

下記で計算している xp は曲線プロット用の変数で、線形回帰の計算とは関係ありません。 matplotlib では表示するデータ点の形状を指定できます。

In [11]:
xp = np.linspace(-0.1, 3.1, 100)
plt.plot(xs,ys,"o", xp, np.sin(xp), "-")
plt.ylim(-0.2, 1.5)
Out[11]:
(-0.2, 1.5)

最小自乗法は、誤差 $y_i - f(x_i)$ が正規分布 $\mathcal{N}(0, \sigma^2)$ に従うと仮定し、尤度

$$ \mathcal{L} = \prod_{i=1}^n \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \{- \frac{(y_i - f(x_i))^2}{2 \sigma^2}\} $$

を最大にするように関数 $f$ を決定します。簡単な計算により対数尤度 $\ln \mathcal{L}$ は、

$$ \ln \mathcal{L} = - \frac{n}{2}\ln(2\pi\sigma^2) - \frac{ns^2}{2\sigma^2}, \;\; s^2 = \frac{1}{n} \sum_i (y_i - f(x_i))^2 $$

となります。($s^2$ は標本分散)

多項式で最小自乗法を行う場合は、NumPy の多項式の係数を与える関数 polyfit と、多項式係数から多項式の値を計算する関数 poly1d を利用することができます。

次数が低すぎても、モデルの一致度は悪く、

In [12]:
xp = np.linspace(-0.1, 3.1, 1000)
plt.plot(xs,ys, "o", xp, np.poly1d(np.polyfit(xs,ys,1))(xp), "-", xp, np.poly1d(np.polyfit(xs,ys,2))(xp), "--")
plt.ylim(-0.2, 1.2)
Out[12]:
(-0.2, 1.2)

次数が大きすぎると過学習を起こします。(悪条件になっていてソルバーも精度が出ておらず、本来はすべての点を通るはずですが、通っていません。)

★注意:この下のスクリプトは実行すると警告が出ますが障害ではありません(警告が出るのは想定された結果です)。悪条件で計算精度が出ていないという警告です。

In [13]:
plt.plot(xs,ys, "o", xp, np.poly1d(np.polyfit(xs,ys,30))(xp), "-")
plt.xlim(-0.2, 3.1)
plt.ylim(-1, 2)
/Users/takashi_miyamoto/.pyenv/versions/anaconda-2.1.0/lib/python2.7/site-packages/numpy/lib/polynomial.py:588: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)
Out[13]:
(-1, 2)

そこで、AICを用いてモデル(今回は次数)を選択することにします。モデルのパラメータ数を$k$とすると、AICは

$$ \mathrm{AIC} = -2 \ln \mathcal{L} + 2 k $$

で定義されます。

$\sigma^2$ は現実の問題では未知の値であるため、測定値から求めた標本分散 $s^2$ で代用すると

$$ \mathrm{AIC} = n \ln(2\pi\sigma^2) + \frac{ns^2}{\sigma^2} + 2 k \sim n \ln (2\pi s^2) + n + 2 k = n \ln s^2 + 2k + (\text{モデルに依存しない項})$$

となります。

線形回帰で得た$i$次の多項式関数を返す関数を fit(i)、その関数と真の値 $y$ との残差平方和を返す関数 s2(i) を定義します。 多項式の次数を1次から15次まで変化させて、次数と残差平方和の対数のグラフを作成してみます。

自作した s2 のようなスカラー値をとってスカラー値を返す関数を、ベクトル値をとってベクトル値を返す関数に変換するためには、np.vectorize という関数を使います。

In [14]:
def fit(i):
    return np.poly1d(np.polyfit(xs,ys,i))
def s2(i):
    return ((ys - fit(i)(xs))**2).sum()/len(xs)
degs = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])
plt.plot(degs, np.vectorize(s2)(degs))
plt.yscale('log')

次に赤池情報量基準を返す関数 aic(i) を定義し、プロットしてみます。

In [15]:
def aic(i):
    return len(xs)*np.log(s2(i))+2*(i+1)
aics = np.vectorize(aic)(degs)
plt.plot(degs, aics)
Out[15]:
[<matplotlib.lines.Line2D at 0x10a7372d0>]

最初のデータの関数形や誤差の大きさ、データ点数など色々変化させて試してみましょう。

時間が余った人は下記も試してみてください。

自由課題1

また、今回は元の関数が $\sin x$ と有限項の多項式では書けない関数でしたが、 $$y \sim x - \frac{x^3}{6} + \mathcal{N}(\mu=0, \sigma=0.1)$$ のような3次式であったならば、AICで正しく次数を3次と推定できるでしょうか?

In [15]:
 

自由課題2

(先に boston.ipynb を終わらせてから。) boston.ipynb で学んだ Ridge, Lasso などを用いて多項式近似を試してみましょう。

多項式近似の係数行列は Vandermonde の行列と呼ばれ、numpy ライブラリで簡単に生成できます。ここでは2次式で試してみましょう。 通常の線形回帰は下記のようになります。

In [16]:
X2 = np.polynomial.polynomial.polyvander(xs, 2)
from sklearn.linear_model import LinearRegression
lr = LinearRegression(normalize=True)
lr.fit(X2, ys)
predicted = lr.predict(X2)
plt.plot(xs, ys, "o", xs, predicted, "-")
Out[16]:
[<matplotlib.lines.Line2D at 0x10ab53fd0>,
 <matplotlib.lines.Line2D at 0x10ac05e10>]

次数を変えるとか Ridge, Lasso なども自分で試してみてください。

In [16]: