前回までに紹介したベイズ線形回帰を実装してみます。

ベイジアンという言葉に難しい印象を持たれている方もいるかもしれませんが、実装が劇的に難しくなったりはしませんから、ご安心ください。

ベイジアンに難しいところがあるとすれば、増えたパラメータをどう決めるかという点と、確率分布として求まる解をどう扱うかという点でしょうか。今回はそのあたりも含めて、見ていくことにしましょう。

環境はこれまでと同じPython&numpy&matplotlibを使用します。インストールなどがまだの方は連載第6回を参照ください。

普通の線形回帰のコードを復習

それでは、ベイズ線形回帰を解くコードを実際に書いていくのですが、第11回で書いた普通の線形回帰のコードに必要な部分を書き足す形で進めましょう。ただし、特徴関数φにはガウス基底を使うことにします。

ガウス基底は、次のような正規分布と同じ釣り鐘型をした関数です。ただし分布ではないので正規化は必要ありません。

ガウス基底の良いところは、データ点の近くはその情報を強く使って、離れるにしたがって影響を弱めるという自然なモデルを構成できるところです。

一方の多項式基底は見慣れた多項式という形式で解を得られるというメリットの反面、データが離れた点での推定に影響を強く及ぼすという特徴があります。解の表記に対する制約が強すぎるという言い方もできるでしょう。その点、ガウス基底は解の表記に対する事前の仮定を弱くできます。

ガウス基底を使う場合、上の式の中のc_iとsをどのように取るか決めなければなりません。

sはデータが影響を及ぼす距離をコントロールするパラメータです。大きくするほど遠くまで影響が届くようになります。そのため小さめの値が望ましいですが、あまり小さくすると、データの密度が低い箇所でどのデータ点も使えず、正しく推論ができないということが起きたりします。今回はs=0.1としておきましょう。

c_iはガウス基底の中心です。このガウス基底の線形和で求めたい関数を表すわけですから、回帰を行いたい区間をカバーするように取る必要があります。基本は先ほどのsに設定した値の間隔で取っていけば大丈夫でしょう。というわけで今回は[0,1]区間からs=0.1間隔で取った11個のガウス基底関数を使いましょう。

以下に第11回のコードについて、関数phiをガウス基底に変更したものと、その場合に推定された関数のグラフを掲げます。phiは11個のガウス基底関数に加えて、定数項を表すための1を加えて12次元のベクトルを返します。

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])

# 定数項＋ガウス基底(12次元) 
def phi(x): 
    s = 0.1 # ガウス基底の「幅」
    return np.append(1, np.exp(-(x - np.arange(0, 1 + s, s)) ** 2 / (2 * s * s)))

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()

12次元のガウス基底は4次元の多項式基底よりずっと自由度と表現力が高いため、推定された関数はかなり「おかしい」状態になってしまっているようです。いわゆる過学習という状態ですね。

これがベイジアンになることでどう変わるかも見て欲しいポイントです。

ちなみに今回、特徴関数をガウス基底に変更したわけですが、コードの変更点がその特徴関数phiの定義を書き換えるだけで済んでいます。このように特徴関数の作り方によらず、パラメータwを特徴行列PHIに対する行列計算だけで決めることができてしまうのが線形回帰の良いところでした。

その点はベイズ線形回帰も同様ですから、PHIを求めるまでのコードはそのまま使っていきましょう。

事後分布の平均と共分散を求める

先ほどのコードでのxlist=～で始まる行から下はグラフを描く部分です。

ということは、普通の線形回帰そのもののコードは次の一行だけだったことがわかります。

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

この一行をベイズ線形回帰にあわせて書き換えるのが第一ミッションです。

より具体的に言うと、前回紹介した事後分布N(w|μ_N,Σ_N)の平均μ_Nと共分散行列Σ_Nを求めるコードを書くことになります。

この式を実装するには、文字αとβをどうにかしなくてはなりません。これらの考察は後ほど見ることにして、ここではひとまずα=0.1, β=9.0とおいて進めることにします。

あとは単なる行列の計算ですから、ほぼ逐語訳でnumpyのコードに落とすことができます。

alpha = 0.1 
beta = 9.0 

Sigma_N = np.linalg.inv(alpha * np.identity(PHI.shape[1]) + beta * np.dot(PHI.T, PHI))
mu_N = beta * np.dot(Sigma_N, np.dot(PHI.T, t))

こうして求めたmu_Nを使ってグラフを描けばめでたしめでたしですが、せっかくですからベイズ線形回帰（青線）と普通の線形回帰（緑線）の両方をグラフに重ねて見比べましょう。

matplotlibのplot関数は続けて呼び出すことでグラフを重ねて描くことができます。色を指定するには、第3引数に'b'（青）,'g'（緑色）のように色を表す文字を指定します。

xlist = np.arange(0, 1, 0.01) 
plt.plot(xlist, [np.dot(w, phi(x)) for x in xlist], 'g') 
plt.plot(xlist, [np.dot(mu_N, phi(x)) for x in xlist], 'b') 
plt.plot(X, t, 'o') 
plt.show()

具体的には、事後確率が最大となる値をwの推定値としているわけですが、これは前回の最後で確認したように、正則化付きの線形回帰の結果と一致します。

このように、パラメータに事後分布を最大化する値を選ぶ手法を「MAP推定」（⁠Maximum a Posterior）と呼びます。

共分散の眺め方

これで一応ベイズ線形回帰を実装して問題を解いた格好にはなりましたが、しかしこれだけではわざわざベイジアンで解いた甲斐がありません。例えば共分散Σ_Nは実質使ってませんよね。

そこでまず、共分散とはなんだろうという話をしましょう。

一次元の分布の「分散」はわかりやすいです。正規分布で具体的に説明すると、分散がσ²であるとは、平均を中心として-σから+σの範囲に約7割の点が、-3σから+3σの範囲に99.7%の点が入るくらいの「ばらけ具合」を表しています。

余談ですが、学校や入試のテストで使われる偏差値とは平均50で分散10²、つまり40から60の範囲に約7割の、20から80の範囲に99.7%の受験者が入るくらいに点数を正規化したものだったりします（実際には正規分布ではないのでズレがあります⁠）⁠。

しかし多次元の「共分散」は行列になりますから、そこまで単純な話ではすみません。その正体を知るために、先ほどのサンプルコードで求めた共分散行列Σ_Nを表示して眺めてみましょう。ただしΣ_Nは12x12の結構大きい行列なのでそのまま表示すると見にくいですから、少し整形して表示します。

# print "\n".join(' '.join("% .2f" % x for x in y) for y in Sigma_N)

 2.94 -2.03 -0.92 -1.13 -1.28 -1.10 -1.21 -1.14 -1.23 -1.06 -0.96 -2.08
-2.03  2.33 -0.70  1.95  0.13  1.02  0.85  0.65  0.98  0.70  0.65  1.44
-0.92 -0.70  2.52 -1.86  1.97 -0.29  0.42  0.59  0.13  0.40  0.32  0.63
-1.13  1.95 -1.86  3.02 -1.66  1.50  0.17  0.29  0.73  0.33  0.36  0.82
-1.28  0.13  1.97 -1.66  2.82 -1.11  1.39  0.22  0.55  0.49  0.40  0.92
-1.10  1.02 -0.29  1.50 -1.11  2.39 -1.35  1.72 -0.29  0.53  0.46  0.69
-1.21  0.85  0.42  0.17  1.39 -1.35  2.94 -2.06  2.39 -0.02  0.25  1.01
-1.14  0.65  0.59  0.29  0.22  1.72 -2.06  4.05 -2.72  1.43  0.37  0.67
-1.23  0.98  0.13  0.73  0.55 -0.29  2.39 -2.72  3.96 -1.41  1.23  0.59
-1.06  0.70  0.40  0.33  0.49  0.53 -0.02  1.43 -1.41  3.30 -2.27  2.05
-0.96  0.65  0.32  0.36  0.40  0.46  0.25  0.37  1.23 -2.27  3.14 -0.86
-2.08  1.44  0.63  0.82  0.92  0.69  1.01  0.67  0.59  2.05 -0.86  2.45

まずこれを見ると、共分散行列は対称行列であることがわかります。

そこで「対角成分⁠」⁠、つまり左上の(1,1)成分から右下の(4,4)成分まで順に斜めに取ったものと、上三角行列に分けて、見やすくしましょう。

# 対角成分

(1)   (2)   (3)   (4)   (5)   (6)   (7)   (8)   (9)   (10)  (11)  (12) 
 2.94  2.33  2.52  3.02  2.82  2.39  2.94  4.05  3.96  3.30  3.14  2.45 

# 上三角行列の成分

    (2)   (3)   (4)   (5)   (6)   (7)   (8)   (9)   (10)  (11)  (12) 
(1) -2.03 -0.92 -1.13 -1.28 -1.10 -1.21 -1.14 -1.23 -1.06 -0.96 -2.08 
(2)       -0.70  1.95  0.13  1.02  0.85  0.65  0.98  0.70  0.65  1.44 
(3)             -1.86  1.97 -0.29  0.42  0.59  0.13  0.40  0.32  0.63 
(4)                   -1.66  1.50  0.17  0.29  0.73  0.33  0.36  0.82 
(5)                         -1.11  1.39  0.22  0.55  0.49  0.40  0.92 
(6)                               -1.35  1.72 -0.29  0.53  0.46  0.69 
(7)                                     -2.06  2.39 -0.02  0.25  1.01 
(8)                                           -2.72  1.43  0.37  0.67 
(9)                                                 -1.41  1.23  0.59 
(10)                                                      -2.27  2.05 
(11)                                                            -0.86 

これが正味の共分散行列の情報量となります。

実はこの対角成分と上三角行列をそれぞれ解釈することで、共分散行列を読み解くことができます。

まず対角成分ですが、これは対応するパラメータを単独で見た場合の分散にあたります[1]⁠。したがって、対角成分は一次元の正規分布の分散と同じように読むことができるわけです。

この例では一番小さいものでも2.33と、どの対角成分もあまり小さくありません。おそらくデータが少ないため、パラメータ推定の精度があまり高くないのだろうということが読みとれます。

そして三角行列の成分は、対応するパラメータ同士の相関を表しています。0だと相関なし（互いに独立⁠）⁠、正の値だと片方が平均より大きいときはもう片方も平均より大きくなる傾向が、逆に負の値なら平均より小さくなる傾向がある、という具合です。

上の例をじっくり眺めてみると、ガウス基底を使った場合、定数項の係数（第1パラメータ）とそれ以外のパラメータ、および隣同士のガウス基底の係数は負の相関があり、離れたガウス基底同士はおおむね正の相関があるがどれもだいたい弱い（0に近い⁠）⁠、ということが読みとれます。これは線形回帰がどのように関数を作るか少し考えてみれば納得いくのではないでしょうか。

実際に問題を解くときは共分散行列をわざわざこのように眺めることはあまりしません。しかし、どういう意味を持つのかわかっている方がやはりいいですよね。例えば計算間違いやバグのせいでとんでもない値が求まってしまったときにも気づくことができるかもしれません。

予測分布を描く

共分散行列の読み方が少しわかったとはいえ、まだ物足りない感じがあります。

もっと事後分布を分布として活用する方法の一つとして、予測分布があります。

イメージ的には、事後分布p(w|X)とはパラメータwの分布なわけですが、これをp(y|w,x)に「代入」すると、xを与えるとyの分布が返ってくる「関数」p(y|X,x)が得られます。

いえ、分布に分布を単純に代入なんてできませんし、分布が返ってくるものは関数の定義にはあてはまりませんよ。あくまでもイメージです。

このようにターゲットとなる変数の分布に事後分布を反映させたものを予測分布と呼びます。

今回のベイズ線形回帰の予測分布は、導出過程は都合上省略しますが、確率の加法定理・乗法定理などを駆使すると次のように得られます。

これは「xごとに決まるyの分布⁠」⁠、つまり2次元上の関数になりますので、これを描くには2次元では足りず、分布（確率密度）の高さを表すためにもう1次元必要になります。3次元のグラフを描いてもいいですが、ここでは濃淡で表すことにしましょう。

matplotlibのcontourf関数を使えばそういう図を簡単に描くことができるのですが、2次元上の関数値を計算して渡してあげる必要があり、そこだけは少々面倒なコードを書いてあげる必要があります。

# 正規分布の確率密度関数
def normal_dist_pdf(x, mean, var): 
    return np.exp(-(x-mean) ** 2 / (2 * var)) / np.sqrt(2 * np.pi * var)

# 2次形式( x^T A x を計算)
def quad_form(A, x):
    return np.dot(x, np.dot(A, x))

xlist = np.arange(0, 1, 0.01)
tlist = np.arange(-1.5, 1.5, 0.01)
z = np.array([normal_dist_pdf(tlist, np.dot(mu_N, phi(x)),
        1 / beta + quad_form(Sigma_N, phi(x))) for x in xlist]).T

plt.contourf(xlist, tlist, z, 5, cmap=plt.cm.binary)
plt.plot(xlist, [np.dot(mu_N, phi(x)) for x in xlist], 'r')
plt.plot(X, t, 'go')
plt.show()

濃い部分は確率密度関数の値が高い、つまり推定される関数がそこを通る可能性が高い部分です。

こうして描いた予測分布を見ると、データの密度が高いところは推定に自信があり、薄いところはデータ点の間隔が広いため推定に自信がないといったことがわかります。

普通の線形回帰では目立って「過学習」っぽい振る舞いをしていたx=0.7のあたりなどは特に色が薄いので、推定した関数のこのあたりはあまり信用できないな、ということがわかってしまうのもベイジアンならではですね。

パラメータαとβの決め方

最後に、後回しにしていたパラメータαとβについて少し考えてみましょう。

そもそもαとβとは何だったでしょう？　αは事前分布p(w)を導入するときに初めて登場しました。こういった、ベイジアンの事前分布に入っているパラメータは特にハイパーパラメータとも呼ばれます。

この式を見ると、αが大きいほど分散が小さくなる、つまりwが0に近い値だという事前知識が強くなります。この状態でベイズ線形回帰をとくと、wを0に近い値に推定しようとする力が強いため、いわゆる過学習しているような結果になるのが抑えられやすい反面、真の解にたどり着くまでに多くのデータを必要としてしまうかもしれません。

逆にαが小さいとwを押さえつける力が弱くなります。特にαが0の時は普通の線形回帰と一致します。そこで0に近い小さめの値、例えば0.1や0.01あたりを設定してまずは解いてみることが多いです。

もう一つのパラメータβは「精度」と呼ばれ、条件付き確率p(y|w,x)の導入で出てきました。

この式でβは、観測値tが真の値w^Tφ(x)からどれくらいぶれてもいいか（ノイズが乗ってもいいか）を表すパラメータとして与えられていました。βが大きいとぶれは少なく、小さいとぶれが大きくても許されます。

つまり、このβは「データを得るまでにノイズ要因がどれくらいあったか」といった事情を考慮しつつ、1以上のできるだけ大きめの値が好まれるようです。

ここまでαとβを問題を解きたい人が選ぶ話をしましたが、「⁠交差検定」によってパラメータを自動的に決める方法や、αにも事前分布を入ることでパラメータを直接決定しなくてもいいようにする手法など、他にもいろいろなパラメータの決め方はあるのですが、いろいろあるということは裏返せば正解がないということ。

一番必要なのは「経験と勘」……と言ってしまうと身も蓋もないですね。ただ、それは機械学習全般に当てはまることなので、いろいろな技術を試しておくのが使いこなせるようになる一番の近道なのかもしれません。

次回予告

次回からは分類問題の話に入っていき、何回かかけてロジスティック回帰にたどり着きたいと思います。お楽しみに。