前回の掲載からしばらく間が空いてしまいましたが，今後は中谷の方で連載を進めていくことになりました。理論編と実践編を交互に進めていくスタイルは継続していきますので，引き続きよろしくお願いします。

線形回帰の復習

今回は連載第8回と第9回で紹介した線形回帰を実装してみる実践編です。

まずは簡単に復習しましょう。

回帰とは，与えられたデータに適した（データを上手く説明できる）関数を求める手法です。点の近くを通る曲線を見つけるときにも用いられます。中でも，基本形として選んだ関数の線形和（式1）から関数を探すのが線形回帰です。

ここでφ(x)=(φ_m(x))を基底関数といい，何か適切な関数を選んで固定します。その選び方によってモデルの性能や得られる関数の形などが決まるので，基底関数は解きたい問題にあわせて選ぶ必要があります。

しかし今はわかりやすさを優先して，シンプルな多項式基底（式2）を採用しましょう。基底関数に何を選んでも線形回帰の解き方自体は同じになりますので，心配しなくても大丈夫です。

この多項式基底を使った場合，f(x)は（式3）のような形になり，(M-1)次式を求めるモデルになることがわかります。

基底関数φ(x)はこのように選んで固定しつつ，パラメータwは自由に動かします。そうするとf(x)が変化しますから，その中で「与えられたデータ点に一番適したw」を見つけるのが線形回帰のアプローチです。

f(x)がデータ点 (x_n, t_n) (n=1, ..., N) の「一番近く」を通るようなwは，次の連立方程式の解として与えられます。

これを導くには二乗誤差やその偏微分などの理解が必要ですが，そのあたりのことは連載の第8回と9回を参照してください。

行列を使ってこの連立方程式を整理すると，解wは（式4）で求められます。

ただし，Φは（式5）のように定義される行列とします。

こうして求めたパラメータwを（式1）に代入すると，データ点の近くを通る関数f(x)が得られます。

ちなみにM=2のときは，データに一番近い直線を求める「最小二乗法」に一致します。

実装のための準備

それでは線形回帰を実装していきましょう。

環境には，これまでと同じくPythonとnumpy/matplotlibを利用します。インストール手順などは連載の第6回および第7回を参照してください。

機械学習に限らず，数式で説明されたアルゴリズムを実装するときに一番最初にやるべきことは，使われている記号（変数や関数）がどういう形式をしていて，どのように初期化したり定義したりしなければいけないかを確認することです。

そこで，今回の数式の中で使われている変数と関数を改めて確認しておきましょう。

注

x_nがD次元空間の点の場合，XはN×D次元の行列として与えられます。今回は簡単のためD=1としていますが，D≧2の場合も線形回帰のアルゴリズムは同じです。

当たり前のことをわざわざ繰り返しているだけに見えるかもしれませんが，これがなかなか大切なことです。いざ実装してみようとしたら上手くいかない，という場合，実は変数や関数の形を勘違いしていた，ということも結構多いのです。

また，記号それぞれがスカラー値なのかベクトルなのか行列なのか，ベクトルや行列ならその次元はいくつか，最初に確認しておくことで無駄なつまずきを減らせます。

こういったことは紙の上で丁寧に確認するのがおすすめです。すべてコンピュータの上で済ませられると思ったら大間違いですよ！

線形回帰の実装

変数の形がわかったところで，順にコードを書いていきましょう。

まずはnumpyとmatplotlibを利用するための宣言を書きます。

import numpy as np
import matplotlib.pyplot as plt

Xとtはデータでしたから，その値は計算するものではなく，何か与えられるものになります。

今回は簡単のために，今回の記事のために生成したデータをコードの中に直書きします。

X = np.array([0.02, 0.12, 0.19, 0.27, 0.42, 0.51, 0.64, 0.84, 0.88, 0.99])
t = np.array([0.05, 0.87, 0.94, 0.92, 0.54, -0.11, -0.78, -0.89, -0.79, -0.04])

Xとtを外部のファイルから読みこむようにすれば，いろいろなデータで試すことができるようになりますが，それは各自の課題として取り組んでみてください。

次は行列Φの計算です。

まずは基底関数φ(x)にあたるものとして， [φ₁(x), φ₂(x), ……, φ_M(x)] を返す関数phi(x)を定義しましょう。

そろそろMも決めなくてはいけませんね。ここではひとまずM=4とします。

今回は多項式基底を使うのでしたから，（式2）の通り，φ_m(x)=x^m-1をm=1から4まで考えて，次のような関数を定義しましょう。

def phi(x):
    return [1, x, x**2, x**3]

このphi(x)とXから（式5）にしたがって行列Φを作るには，次のように書きます。

PHI = np.array([phi(x) for x in X])

[～ for ～ in ～] というのは「リスト内包」と呼ばれるPythonの機能で，ループを回しながら配列を作ってくれます。そしてnp.arrayは配列をベクトルや行列に変換します。

phi(x)は [φ₁(x), φ₂(x), ……, φ_M(x)] を返す関数でしたから， [phi(x) for x in X] は [φ₁(x_n), φ₂(x_n), ……, φ_M(x_n)] が並んだ配列となります。これをnp.arrayで変換することにより，求めたい行列Φが手に入ります。

# print PHI で PHI の中身を確認すると：
[[  1.00000000e+00   2.00000000e-02   4.00000000e-04   8.00000000e-06]
 [  1.00000000e+00   1.20000000e-01   1.44000000e-02   1.72800000e-03]
 [  1.00000000e+00   1.90000000e-01   3.61000000e-02   6.85900000e-03]
 [  1.00000000e+00   2.70000000e-01   7.29000000e-02   1.96830000e-02]
 [  1.00000000e+00   4.20000000e-01   1.76400000e-01   7.40880000e-02]
 [  1.00000000e+00   5.10000000e-01   2.60100000e-01   1.32651000e-01]
 [  1.00000000e+00   6.40000000e-01   4.09600000e-01   2.62144000e-01]
 [  1.00000000e+00   8.40000000e-01   7.05600000e-01   5.92704000e-01]
 [  1.00000000e+00   8.80000000e-01   7.74400000e-01   6.81472000e-01]
 [  1.00000000e+00   9.90000000e-01   9.80100000e-01   9.70299000e-01]]

最後に，（式4）にしたがってパラメータwを推定する数式を実装します。numpyには逆行列を求めるnp.linalg.inv関数がありますから，それを使うと次のように書けます。

w = np.dot(np.linalg.inv(np.dot(PHI.T, PHI)), np.dot(PHI.T, t))

また，連立方程式の解を求めるnp.linalg.solve関数もあります。これを使うと次のようにも書けます。

# np.linalg.solve(A, b) は A^{-1}b を返す
w = np.linalg.solve(np.dot(PHI.T, PHI), np.dot(PHI.T, t)) 

np.linalg.solveの方がコードが短く，実行速度も速いのでこちらを使うことにしましょう。

# print w で w の中身を確認すると：
[ -0.14051447  11.51076413 -33.6161329   22.33919877]