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

線形回帰の復習

今回は連載第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]

推定した関数を図示する

こうしてwを求めることができましたが、数字だけ見てても、問題が解けた!!という実感は湧きませんよね。そこで、matplotlibで図を描いてみましょう。

そのためにまず（式1）または（式3）を使って、関数f(x)を用意する必要があります。

# (式 1) の場合
def f(w, x):
    return np.dot(w, phi(x))

# (式 3) の場合
def f(w, x):
    return w[0] + w[1] * x + w[2] * (x ** 2) + w[3] * (x ** 3)

どちらでもかまいませんが、（⁠式1）に基づくコードは、phi(x)の定義を書き変えてもそのまま使えます。

関数を曲線として描くにはpyplot.plot関数を用います。点列を与えて曲線を描くのがmatplotlibの流儀です。

xlist = np.arange(0, 1, 0.01)
ylist = [f(w, x) for x in xlist]

plt.plot(xlist, ylist)

xlistは、グラフを描く範囲 [0, 1] を等間隔に刻んだ配列（ベクトル）です。

np.arange()の3つ目の引数0.01は刻み幅です。これを0.1などにするとカクカクしたグラフが描かれてしまいますので、曲線っぽく見えるくらいまで十分細かくしましょう。

ylistは、xlistのそれぞれの入力xに対応する関数の値f(x)のリストです。ここでもリスト内包が活躍していますね。

あとはデータ点の分布を描いて、曲線と見比べたいですよね。先ほど曲線を描くのに使ったplot()の第三引数に'o'を与えると、点の分布を描くことができます。

plt.plot(X, t, 'o')

最後にplt.show()を呼ぶと、描いたグラフが表示されます。

plt.show()

様々な条件で試す

ここまでのコードをまとめましょう。

import numpy as np
import matplotlib.pyplot as plt

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.79, -0.89, -0.04])

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

PHI = np.array([phi(x) for x in X])
w = np.linalg.solve(np.dot(PHI.T, PHI), np.dot(PHI.T, t))

xlist = np.arange(0, 1, 0.01)
ylist = [np.dot(w, phi(x)) for x in xlist]

plt.plot(xlist, ylist)
plt.plot(X, t, 'o')
plt.show()

見通しの良い短いコードですから、いろいろ書き換えて楽しむことができます。

例えば、データを増やしたり減らしたりして、それに適した関数がきちんと得られるかを見てみる、というのが序の口でしょうが、もっと本質的な部分に触ってみたいという向きには、「⁠基底関数の変更」と「正則化の導入」あたりがおすすめです。

基底関数を変更する

上のコードではφ(x)の次数、つまり基底関数の個数はM=4でしたが、これを変えてみましょう。

例えばM=8にするには、phi(x)の返値を [1, x, x**2, x**3, x**4, x**5, x**6, x**7] にすればOKです。

あるいは、多項式以外の基底関数を試してみることもできます。例えばガウス基底は次の数式で表されます。

ここでは [0, 1] 区間をM等分する点を、sは正の定数をとります。sは0.1くらいから始めて、いろいろ変えてみるといいでしょう。

このように基底関数を取り替えても、phi(x)の定義を書き換えるだけでよいのが、線形回帰の良いところですね。

正則化項を入れる

「正則化」については連載の第9回をご覧ください。

実装に関係のある結果だけ説明すると、正則化を入れることでパラメータwを求める式が以下のように変わります。しかし、それ以外は今までの線形回帰そのままです。

Iは単位行列、λは正の定数です。これなら上のコードを一ケ所ちょこっと書き換えるだけで済みそうですよね。

λをいくつにするかが迷うところですが、0.01あたりから始めて、増やしたり減らしたりいろいろ試してみてください。ちなみに、基底関数を [1, x, ..., x**7] にしておくと、正則化の効果を実感しやすいですよ。

このように条件を少しずつ変えながら結果を確認しておくことで、線形回帰の「手触り⁠」⁠、つまり傾向や特性のようなものを感じることができるかもしれません。

そういった蓄積は、機械学習を実際に使いこなそうというときにきっと役にたちますから、是非いろいろ試してみてください。

次回は理論編。ベイズ線形回帰を紹介します。