機械学習 はじめよう

第14回 ベイズ線形回帰を実装してみよう

この記事を読むのに必要な時間:およそ 5 分

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

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

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

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

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

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

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

画像

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

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

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

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

ciはガウス基底の中心です。このガウス基底の線形和で求めたい関数を表すわけですから,回帰を行いたい区間をカバーするように取る必要があります。基本は先ほどの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を求めるまでのコードはそのまま使っていきましょう。

著者プロフィール

中谷秀洋(なかたにしゅうよう)

サイボウズ・ラボ(株)にてWebアプリ連携や自然言語処理を中心に研究開発を行いながら,英単語タイピングゲーム iVocaをサービス提供している。

URLhttp://d.hatena.ne.jp/n_shuyo/
Twitterhttp://twitter.com/shuyo

コメント

  • 今後の連載を心待ちにしています

    今月から連載を読みはじめて、追い付くことが出来ました。確率関係は苦手意識を持っており、今までだましだましやってきましたが、理論編がわかりやすく、なんとかついてくることができました。
    回を追うごとに少しずつ機械学習っぽくなっており、今後がすごく楽しみです。

    これからも継続的に拝読させていただきます。

    Commented : #4  itu (2012/05/22, 20:25)

  • ご教示頂き有難うございます。

    ご教示頂き有難うございます。 定数項の意味合い分かりました。α、βについても後の方でちゃんと説明が書いてあるのですね。α+β=1と勝手に思っていましたが、それぞれ独立なのですね。途中で投げ出すんじゃんかった。
    予測分布の墨絵のような図を見て、ベイズ線形回帰が何となく分かりかけてきました。この図でαを色々変えてみれば事前確率分布の意味なんかがよく分かりそうですね。
    誰かgnuplotでのこの絵の描き方教えて!

    Commented : #3  kiyo (2012/04/17, 09:53)

  • 定数項について

    筆者の中谷です。早速試していただいてありがとうございます。

    定数項は一般には付加することが多いですが、絶対に必要というわけではありません。結局それは基底関数の選び方の話なので、「問題を解きたい人がどんな答えが欲しいか」という部分にゆだねられます。つまり、定数項を入れる(あるいは入れない)のが正解とか一概には言えないのです。

    今回の例ではt=0上の点で対称なために、そもそも定数項があっても大きい値になる見込みがありませんが、もっと全体に上に移動した点であれば、定数項の有無で目に見える影響が出るでしょう。
    #その意味では例をもうちょっと工夫しておいた方が良かったかも?

    紙数の関係で本文では省略しましたが、定数項はグラフの上下方向の平行移動量に当たることから、正則化の対象から除外するということもよく行われます。
    ご興味があればそのあたりも交えて試していただけるとまたいろいろ掴めるのではないかなと思います。

    Commented : #2  中谷秀洋 (2012/04/15, 01:39)

  • お教えください

    待ちに待ったベイズ線形回帰実装編が出たので早速実装してみました。といっても新しい言語はわからないのでc++とgnuplotを使ってです。自分で実装していつも気になるのは正しいのかどうかなのですが・・・
    まず最初にガウス線形分布で、読み飛ばして、定数項の1を付けずに行ってしまいましたが、同様のグラフ結果を得ました。定数項は必要だったのですか?
    定数項を付けないままベイズ線形回帰まで行きましたが、α=0.1、β=0.9では同じグラフが得られず、α=0.01、β=0.99にすると、同じ過学習気味のグラフが得られました。
    前回の説明ではα/βは正則化定数なのでプログラムの動きとしては合ってそうな気がするのですが。
    ご教示いただければ幸甚です。

    Commented : #1  kiyo (2012/04/09, 15:58)

コメントの記入