パラメトリックとノンパラメトリックの狭間

※この記事は、Machine Learning Advent Calendar 2012（http://qiita.com/advent-calendar/2012/machinelearning）の10日目用に書かれています。

はじめに

Machine Learning Advent Calendar 2012の10日目を担当します、@risuoku です。
今回は、数ある機械学習手法の中で、以下の２つに焦点を当て、いくつかのアプローチを紹介します。

  1. 回帰分析
  2. クラスタリング

特に、パラメトリックな手法とノンパラメトリックな手法の違いや、それぞれの特徴の理解を目指しています。
また、@risuokuはPRMLをよりどころに機械学習を勉強してきました。なので、この本の影響を強く受けていることを初めにお伝えしておきます。

回帰分析

「形」の発見

突然ですが、以下の画像から、どんな知識が得られるでしょうか？

なんとなく、横軸の左より右の方が、縦軸の値が大きい気がしませんか？
「比例関係にある」という知識が得られたことになるかもしれません。

では、次の画像はどうでしょう？

今回は、全体的に左から右にいくにつれて、縦軸の値が小さくなってる気がします。
しかし、よく目をこらすと、横軸2付近で一番大きくなって、6あたりで一番小さくなる曲線に見えるかもしれません。

3枚目、これはどうでしょう？

これは、ぱっと見、傾向を掴むのが難しいかもしれません。
一度上がって、下がっているような曲線には見えるでしょうか？

実は、それぞれの画像に示された点は、ある曲線に正規乱数によるノイズをのせて生成した人工データです。
データ元となった曲線と、実際に生成されたデータを重ねると以下のようになります。


回帰分析のセクションでは、初めの3枚のようなデータから、データの背後に潜む「形」（この場合は、データ元となった曲線）を発見ないし推定することを目標とします。
また、以後は入力をx（画像では横軸）、出力をy（画像では縦軸）として扱います。

データの「特徴」

データないしオブジェクトというものは、何らかの特徴を持っています。
例えば、ドラクエのキャラは、HP、MP、攻撃力、守備力、素早さなどの特徴をもっています。
そして、ドラクエのキャラの「強さ」は、これらの特徴の組み合わせで決まります。
あるxに対してyを予測するといった時、xの特徴が手がかりとして使えそうです。
前節の1つめの例だと、特徴=$(1, x)$とすると、$y=w_0+w_1x$といった、それぞれの特徴に重みをつけた和で表現できそうです。
1つめの画像が1次関数に見えるというのは、xの特徴として1次の項までを考えることに相当します。
また、2つめが3次関数に見えるとすると、特徴=$(1, x, x^2, x^3)$ を使うことになります。
このように、データの特徴に重みを加えて線形結合した関数を、データにうまくフィットさせる手法を線形回帰といい、非常によく使われます。
この場合、重み(=パラメータ)が決まれば関数が決定するので、結局は関数を推定するというよりパラメータを推定する問題になります。
また、もっと大きな枠組みで、パラメトリック推定といいます。（多分）

パラメトリック推定の問題点

前節では、関数の形を決めてからパラメータを推定する線形回帰を紹介しました。しかし、すぐにこんな疑問がおきそうです。

  • 特徴の設計は何を手がかりにするのか？
  • そもそも、人間が何かを仮定してしまったら本末転倒ではないのか？

まず、特徴の設計（今回の例だと、多項式の次元）ですが、基本的には試行錯誤です。
例えば、2つめと3つめの例を1次関数で表現しようとしても無理があります。なので、1次関数でやってみて、ダメだったら2次関数、次は3次関数、みたいなことを繰り返してそれっぽいものを採用することになります。（これは面倒くさい！）
もちろん、全部を人間の目で評価する必要はなく、自動評価の指標が色々と提案されています。*1しかし、どうせなら関数の形までデータに基づいて推定できると嬉しいはずです。
また、人間が何かを仮定する件ですが、およそ科学という限り、何らかの仮定は必要になってきます。回帰問題を扱う際、多くの場合、実際のデータは正規分布に従うノイズがのった状態で生成されるという仮定が前提になります。これは基本的でありながら、非常に重要です。世の中の統計学者の大半はこの仮定が妥当であると信じています。

次節以降は、推定したい関数をより一般化した（=関数の形をあらかじめ決めない）アプローチを紹介します。
また、このような方法は、ノンパラメトリック推定といいます。（多分）

データ間の「関係」

最初に示した画像の3枚目で、x=2.25の時、yを推定したいとします。図に示すと以下のようになります。x=2.25は、青の点線で示されています。

さて、今回は推定したい関数の形を何も決めていません。見知らぬ土地に手ぶらで迷い込みました。こんな時、まずは何をすべきでしょうか？
やはり、現状把握でしょうか。現状把握するためには、周囲を観察します。周囲を観察するうちに、局所的かもしれませんが、何らかの「ルール」があることに気付くかもしれません。
あとは、自分の立ち位置に応じて、適切な行動を考えます。
x=2.25に相当するyを推定するため、とりあえずx=2.25±1.0くらいの範囲を観察します。この範囲は、x=2.25と強い関連があると考えられます。
すると、以下の図のように、この範囲でyの値が増加傾向にあることに気付きます。

このように見ると、x=2.25にあたるyの大まかな値が予測できた気がしますが、こんなので本当にいいのでしょうか？結果からいうと、この予測はいい線いっています。
すると、あるxの範囲について同様のことを繰り返せば、所望の関数が得られる気がします。

さて、次に進む前に、言葉の定義をしておきます。データ間の関係が推定に使えることがわかりました。ここで、データ間の関係を数値化する関数を準備します。このような関数は、「カーネル関数」という術語で知られています。
例えば、$k(x1,x2)=exp(-\beta{||x1-x2||}^2)$のようなカーネル関数を定義すると、x1とx2は近いほど関係が強く、遠いほど関係が弱いことを定量的に扱えます。

推定手法と呼ぶにはアドホックな感じがしますか？しかし、データ間の情報を推定に使うこのアプローチは、直感的に理解しやすいだけでなく、先に挙げたパラメトリックな手法と深い関連があることが知られています。*2

ガウス過程

前節では、データ間には関係があることをみて、関係を数値化するためにカーネル関数を導入しました。
ここで、推定したい関数fについて考えます。重要なポイントからいうと、fというのは、ガウス過程に従う確率変数です。
関数fの形を明示的に与えない代わりに、関数の値に確率的な制約（結局のところ、ノイズが正規分布に従う仮定に基づきます）をつけています。
fがガウス過程に従うことを形式的に表現すると、$(f(x_1),f(x_2),...,f(x_n))=N(0,C)$ となり、fの値がn次元多変量正規分布に従うことになります。（わかりづらいですが、平均の0はn次元です）
ここで重要なポイントが、共分散Cの意味です。共分散行列の各要素は各変数の関連度を表していました。各変数というのは、推定したい関数の値ですから、たしかにデータ間の関係を考えています。
実は、共分散行列の各要素の値にカーネル関数を使うと、上手く関数fを推定できることが知られています。細かい話は省きますが、p(訓練データ)とp(訓練データ,新しい入力)が計算できることから、p(新しい入力|訓練データ)が得られます。

デモ

理屈ばかりだといまいちイメージがつかめないかもしれないので、動くものを見ていきます（オンラインで学習していく様子を観察します）。
3つめの複雑な曲線から生成されたデータを入力として使います。
ソースコードはgithubに公開しています。
https://github.com/risuoku/AdventCalendar_2012_risuoku


まず１つめ。これは、パラメトリック推定です。特徴ベクトルの次元は6に設定してあります（つまり、5次関数による近似を目指しています）。
動画を見ればわかると思いますが、元の曲線のジグザクな傾向を捉えるのには失敗してます。ただ、大まかに上がって下がる傾向は捉えられています。
こういう時は、次元を上げると多少マシになります。なお、実線による曲線が3つプロットされているのは、パラメータが確率的な分布をもっていて、ランダムに3つサンプリングしているからです。*3


次は、ノンパラメトリック推定です。水色の帯みたいなものは、縦方向に見ます。これが細ければ細いほど、そのあたりにデータが生成されることに高い確信をもっています。
データが少ないころは帯が太い箇所が多いですが、これは情報が少ないので推定が困難なことを意味します。データが増えるにつれ、帯が細い場所が増えていき、システムが自信をつけていく様子が観察できます。
実際、データが増えるにつれ良い近似が得られていることが見てわかると思います。

結局どっちがいいのか？

パラメトリック推定とノンパラメトリック推定の2つを議論してきましたが、この節では、結局どっちがいいのかについて所感を軽く述べたいと思います。
非常に月並みというか、何も言ってないと同じ気もするのですが、結局のところ、「状況や問題による」ということになるのではないかと感じています。
まず、前節のデモでは、データの背後にある複雑な構造を捉えるのに、ノンパラメトリック推定がかなり有利なことを示しました。しかし、パラメトリック推定でも、特徴ベクトルに手を加えれば複雑さには対応できます。
また、実はノンパラメトリック推定のデモでは、データに対するフィッティングをやや大きめに設定してあります。こうすることで、複雑な現象に対しては有利になります。
一方で、設定を変えないままだと、1次関数のような単純な現象に対しては、過剰フィッティング（=過学習）することで汎化性能が弱くなります。
このことからわかるように、いくらノンパラメトリック（=no+parameter）といっても、パラメータの調整次第で上手くいったりいかなかったりします。
これでは、考えるべきパラメータは結局プラマイゼロということになって、何のためにノンパラメトリックにしたのかわからない、ということになりかねません（と、言ってしまうのはさすがに大げさだと思いますが・・）
まだあります。「ノンパラなのにパラメータ調整は必要」という以上に、むしろこっちがきつい気がするのですが、ノンパラメトリック推定は新しいデータに対する予測時の計算コストが訓練データの数に従って大きくなります。
この事実は、ガウス過程の意味を考えれば明らかです。実際、今回のサンプルでいうと、データ数100程度ですでにはっきりわかるくらい描画までにかかる時間が長くなります。

まとめ
  • 回帰分析におけるパラメトリック推定とノンパラメトリック推定を紹介しました
  • 実際に動かしてみました
文献案内
  • カーネル多変量解析：赤穂昭太郎：岩波書店

カーネル法についてとても詳しく解説されています。

クラスタリング

回帰分析のセクションが想像以上に手間取って、正直ちょっと辛くなってきています（汗）なので、ここから先はかなり端折りながら進めていきます。
クラスタリングということなのですが、とりあえずこのセクションではディリクレ過程だけ説明し、その他の話題は文献を紹介するにとどめます。
※注意
このエントリ全般に言えることですが、特に次のディリクレ過程は理解が怪しいところです。おかしいところがあれば教えていただけるとありがたいです。

ディリクレ過程

ディリクレ過程というのは、雑な言い方をすると、$\sum_{k=1}^{n+1}\pi_k=1$を満たす組$(\pi_1,...,\pi_{n+1})$を生成する「母艦」のようなものです。ただし、nまでをクラスとして認め、n+1はまだクラスになっていないとします。
ディリクレ分布と似ています。しかし、ディリクレ過程はより一般的な分布を考えることができます。
ディリクレ過程は、基底分布というものを考え、基底分布の定義域の区切り方により、n個のクラスを定義します。
例えば、基底分布にガウス分布を使い、3つのクラスを考えると、以下のようになります。

注意点として、ディリクレ分布においては各次元において混合比が決まるのに対し、ディリクレ過程では各クラス（この場合、$S_1, S_2, S_3$）で混合比が決まります。
重要なポイントは、基底分布の定義域の区切り方を自由にすることで、「クラスタ」に柔軟性を持たせることができるところです。
このことが、いわゆるノンパラベイズなクラスタリングによる、クラスタ数自動推定といった高度な技術を可能にしています。
ここでは、クラスを3つに明示的に区切りましたが、普通はn個の決まったクラス+未知のクラスのn+1個をもたせることで、理論上無限に拡張可能なモデルにしています。

（言い訳）stick breaking processとかchinese restaurant processの未知ゾーンをどう扱えば良いかいまだにわかってないので、とりあえずn+1個目のクラスタにおしこんでると解釈しています。

まとめ
  • ディリクレ過程について説明しました
文献案内
  • 最近のベイズ理論の進展と応用 (III) ノンパラメトリックベイズ：持橋大地

http://chasen.org/~daiti-m/paper/ieice10npbayes.pdf
日本語で読めるノンパラベイズ資料です。素人目にはよくまとまっている素晴らしい資料に見えます。ただ、気になったことがあって、『無限次元の離散分布によって「すかすか」にする（2ページ目左上）』とはどういうことでしょうか？
図には棒がたっているように見えます。しかし、この棒はどこの領域を占めているのでしょうか？もし、実数のうちある1つだとすると、非加算無限個の棒を考えることになります。しかし、ここで議論されている「無限」というのは明らかに「加算無限」です。
そこで、ある非加算無限の集合（例えば、[tex:$0

  • "Infinite LDA"-Implementing the HDP with minimum code complexity：Gregor Heinrich

http://www.arbylon.net/publications/ilda.pdf
タイトルに"Implementing"とついているとおり、応用寄りの資料です。巻末にはなんと無限次元のLDA（いわゆるHDP-LDA）のJava実装がのっています。とはいえ、前半は適度に理論的な話もあってバランスの良い構成です。また、趣味で書いているのかわかりませんが、かみ砕いた説明が多くて良いです。Y.W.TehやM.Jordanの書いたものを読んだけどちんぷんかんぷんだったという人はまずこのへんから読むとよいです。

  • Hierarchical Bayesian Nonparametric Models with Applications：Y.W.Teh, M.Jordan

http://www.stat.berkeley.edu/tech-reports/770.pdf
長いし難しいですが、面白いです。

  • The Infinite Gaussian Mixture Model：C.Rasmussen

http://www.gatsby.ucl.ac.uk/~edward/pub/inf.mix.nips.99.pdf
混合ガウスの混合数を無限個に拡張したモデルが説明されています。フルベイズだし仰々しい数式がつらつらと並んでいます。ただ、必要な数式は全部出てる感じなので、実装してみるだけならできると思います。
パラメータ推定には大部分ギブスサンプリングを使っていて、一部Adaptive Rejection Sampling(ARS)が必要になるようです。

細々としたこと

ベイズ線形回帰

パラメトリック推定で回帰問題を解くには、パラメータを最適化する必要があります。その際、基本的には2つの方針があります。まず、1つはいわゆる最小二乗法。汎化性能を高めるために正則化項をつけたりします。
2つめは確率モデルで考えて尤度を高める方針。ベイズ的なモデル化だと、パラメータに事前知識をつけたりします。
今回パラメトリック推定と言ってきたのは、正確にはベイズ的な確率モデルの最適化です。データが増えるほどパラメータの事後分布が更新されて賢くなっていくイメージです。
ちなみに、正則化項付き最小二乗法とベイズ的な確率モデルの最適化は数学的に等価です。

最適化について

今回、パラメータをどうやって学習するかはあえて触れませんでした。というのも、回帰の場合、最適解が決定的に得られるので、特に語ることはないと思ったからです（式を見ればわかる、ということになります）。
ちなみに、分類は大変です。分類の最適化を語るだけで1つのエントリが終わってしまいます。

カーネル関数の応用例

せっかくカーネル関数を取り上げたので、今回の主旨からは随分ずれますが、アレを紹介しておきます。カーネル関数の応用例として、最も有名な気がするサポートベクトルマシン（SVM）です。
SVMでカーネル関数を使うことで、非線形な特徴空間を扱えるメリットがあることは広く知られている通りです。しかしここでは特に、SVMでは新しいデータの予測にかかる計算コストが小さくてすむ点を強調しておきます。
というのも、SVMにおいては、各訓練データについて、「サポートベクトルか。さもなくば重み0か」という非常に強い制約（KKT条件の一部）がかかっているからです。*4
こういった理由で、ノンパラメトリック回帰において問題になっていた、訓練データが増えるにつれて計算コストが大きくなる現象を回避しています。

おわりに

最近勉強してたことまとめのような記事を書きました。やはり、アウトプットしてみることで考えが整理されることがあるようです。こういうのはきっかけがないとなかなか手が出ないので、機会を提供してくださった@naoya_tさんに感謝します。

*1:AICやBICなどが有名です

*2:雑な言い方ですが、特徴ベクトルを使って構成した関数を双対変換するとカーネル関数が自然と現れ、カーネル関数は2つの入力それぞれの特徴ベクトルの内積となります。カーネル関数と特徴ベクトルに関連がある事実は、論理遊びによる無味乾燥な産物というわけではなく、意味を考えれば自然に理解できます。

*3:特に断りませんでしたが、ベイズ的なモデル化をしています

*4:PRML風にいうと、「疎な解」を持ちます