2019-04-19

ベイズ統計の理論と方法：ノート2.5

本読み機械学習

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。おかしい点がありましたらコメント等でご指摘いただけますと幸いです。

ベイズ統計の理論と方法

作者: 渡辺澄夫
出版社/メーカー: コロナ社
発売日: 2012/03/01
メディア: 単行本
購入: 1人クリック: 4回
この商品を含むブログ (8件) を見る

前回：ノート2 ／次回：まだ

f:id:cookie-box:20190101155733p:plain:w60

前回は1章の 1.2.2 節以降を読みましたが、そこで学んだ主なことは以下でした。

ベイズ推測は確率モデルが指数型分布である場合は、事前分布から事後分布への更新が「ハイパーパラメータの更新」の形になる。
サンプル $X^n$ から逆温度 $\beta=1$ でベイズ推測を行うとき、
点 $x$ の真の分布と予測分布 $p(x|X^n)$ の交差エントロピー（汎化誤差） $G_n$ の期待値は、
サンプル $X^{n+1}$ に対して $\varphi(w) , \, p(x|w)$ を仮定したときの自由エネルギー $F_{n+1}(1)$ の期待値から、
サンプル $X^n$ に対して $\varphi(w) , \, p(x|w)$ を仮定したときの自由エネルギー $F_n(1)$ の期待値を差し引いたものに等しい。

後者は、「ベイズ推測をするということはどんな予測をするということなのか」に踏み入ってきていますね。つまり、「ベイズ推測をするということは、誤差が平均的にこれくらいになる予測をするということなのだ」ということですが…しかし、私たちの手元にあるのはある1通りのサンプルの出方のみなので、自由エネルギーの「期待値」を知ることはできません。

f:id:cookie-box:20190101160814p:plain:w60

第2章ではもう少し詳しく、そもそも確率モデルによって真の分布が実現可能なのか、パラメータ集合 $W$ の中に最適な $w$ があるのかなどを考えて、どんな条件下で何が成り立つのかを議論するのかな。

勝手に先をぱらぱら読まないでください…しかし、1.6 節の「本書の概略」をみておくと、2章で汎化と推測の間に何が成り立つかを調べて、3章では事後分布が正規分布で近似できる理想的な場合、4章ではそうでない場合により詳しく何がいえるかをみていくのでしょうか。5章は現実の計算方法のようですね。6章はベイズ推測によって現実の目的を達成したいときここまでやるべきということ、といった感じがしますが…。7章は一歩下がって、なぜベイズ推測をするのか、ベイズ推測に限らず統計的推測をするときの心構えのようなものがかかれている感じがしますね。なぜこれが最終章なのかわかりませんが…学生からの質問が多い内容だから付け足したとかでしょうか…。

その理由のベイズ推測は難しそうだね。早速2章に進む？

いえ、少し気になる点が。17ページの「『パラメータ集合がコンパクトであり事前分布が一定値ならば、ベイズ推測は最尤推測と同じである』という言明は、本書の定義を採用するならば、正しくない」というくだりです。この本の定義においては「『ベイズ推測は最尤推測と同じである』ということはない」ということなので別にいいといえばいいんですが…でも、ベイズ推測の定義によってはその条件下でベイズ推測と最尤推測が同じになりうることを示唆していますよね。それってどういうことでしょうか？この本の定義ならベイズ推測と最尤推測は以下ですよね。

$\displaystyle p^{\ast}(x) = \int_W p(x|w)p(w|X^n)dw \, , \quad p(w|X^n) = \frac{1}{Z_n(\beta)} \varphi(w) \prod_{i=1}^n p(X_i | w)^{\beta}$
$\displaystyle p_{ML}(x) = p(x|w_{ML}) \, , \quad \displaystyle w_{ML} = \underset{w \in W}{\rm arg max} \prod_{i=1}^n p(X_i | w)$

これらが同じにはあまりみえませんが…。

…同じにみえるかどうかはおいておいて、ベイズ推測と最尤推測の例を考えてみようか。事前分布は一定値で…データが離散値だからあまりよくない例かもしれないんだけど、コインを投げて表が出る確率を予測したいとして、確率モデルを $p( {\rm head} |w) = w,$ $p( {\rm tail} |w) = 1-w$ とし、 $w \in [0, 1]$ とする。事前分布は $\varphi(w)$ はこの区間内だけ1をとり、この区間外では0をとる分布ね。例えば観測されたデータが「表、裏、裏、裏」だったとすると、最尤推測の結果は尤度関数 $w(1-w)^3$ を最大にするのが $w_{ML} = 1/4$ だから $p_{ML}( {\rm head}) = 1/4$ だけど、 $\beta=1$ でベイズ推測すると、

$Z_n(1) = \int_{0}^{1} w (1-w)^3 dw = 1 / 20$
$p(w|X^n) = 20 \cdot \varphi(w) w (1-w)^3$
$p^{\ast}( {\rm head} ) = 20 \cdot \int_{0}^{1} w^2 (1-w)^3 dw = 20 \cdot 1/60 = 1/3$
$p^{\ast}( {\rm tail} ) = 20 \cdot \int_{0}^{1} w (1-w)^4 dw = 20 \cdot 1/30 = 2/3$

$p_{ML}( {\rm head}) \neq p^{\ast}({\rm head})$ ですね…。最尤推測もベイズ推測もさすがに「表の出る確率の方が低いだろう」となるのは一緒ですが、ベイズ推測は「3分の1は表が出るんじゃないか」ということですから、最尤推測の「4分の1は表が出るんじゃないか」より表が出る確率を大きく見積もっていますが…。

ちなみに $\beta=2$ にするとこうなるね。

$Z_n(2) = \int_{0}^{1} w^2 (1-w)^6 dw = 1 / 252$
$p(w|X^n) = 252 \cdot \varphi(w) w^2 (1-w)^6$
$p^{\ast}( {\rm head} ) = 252 \cdot \int_{0}^{1} w^3 (1-w)^6 dw = 252 \cdot 1/840 = 3/10$
$p^{\ast}( {\rm tail} ) = 252 \cdot \int_{0}^{1} w^2 (1-w)^7 dw = 252 \cdot 1/360 = 7/10$

さっきより微妙に最尤推測に近づきましたね。微妙にですが。

ただコインを投げた回数が4回は少ないかもね。「表1回、裏3回」じゃなく「表1回、裏99回」だったらどうだろう。最尤推測の結果は $p_{ML}({\rm head}) = 1/100$ になるはずだ。でもベイズ推測は $\beta=1$ だったら、

$Z_n(1) = \int_{0}^{1} w (1-w)^{99} dw = 1 / 10100$
$p(w|X^n) = 10100 \cdot \varphi(w) w (1-w)^{99}$
$p^{\ast}( {\rm head} ) = 10100 \cdot \int_{0}^{1} w^2 (1-w)^{99} dw = 10100 \cdot 1/515100 = 1/51$
$p^{\ast}( {\rm tail} ) = 10100 \cdot \int_{0}^{1} w (1-w)^{100} dw = 10100 \cdot 1/10302 = 50/51$

やっぱりベイズ推測の方が表が出る確率を大きく見積もっているようですね。…しかし副部長、定積分の計算速くないですか？

以下のサイトに以下のコマンドを打ち込んだだけだからね。
Wolfram|Alpha: Computational Intelligence

integrate w (1-w)^3 dw from w = 0 to 1
integrate w^2 (1-w)^3 dw from w = 0 to 1
integrate w (1-w)^4 dw from w = 0 to 1

integrate w^2 (1-w)^6 dw from w = 0 to 1
integrate w^3 (1-w)^6 dw from w = 0 to 1
integrate w^2 (1-w)^7 dw from w = 0 to 1

integrate w (1-w)^99 dw from w = 0 to 1
integrate w^2 (1-w)^99 dw from w = 0 to 1
integrate w (1-w)^100 dw from w = 0 to 1

手抜きだった。

最初は手計算してたけど途中でめんどくさくなってきたからね。

…あの、 $\beta=1$ の場合ですが、「表1回、裏3回」の場合も「表1回、裏99回」の場合も、ベイズ推測の結果は最尤推測の結果と一致していませんが、「表1回、裏1回」を付け足した「表2回、裏4回」「表2回、裏100回」に対する最尤推測の結果とは一致していませんか？？

…「表1回、裏n-1回」のベイズ推測の結果が「表2回、裏n回」に対する最尤推測の結果と一致するね。最尤推測だと尤度関数 $w(1-w)^{n-1}$ の最大点 $w_{ML}$ をとるわけだけど、ベイズ推測だとこの $w(1-w)^{n-1}$ を事後分布としてこれで確率モデルを平均するから、表と裏が出る確率の予測値の比は $w^2(1-w)^{n-1}$ と $w(1-w)^n$ をそれぞれ $[0, 1]$ で定積分した比になる。この比は $2:n$ になる。前者を部分積分すれば後者の $2/n$ 倍になることがただちにわかるからね。

最尤推測では $w(1-w)^{n-1}$ の最大点をとってしまうのに対して、ベイズ推測ではこれを事後分布として平均することが表1回と裏1回分の差につながっている…。

こうしてみると「事前分布が一定値である」と「パラメータに事前分布などない（頻度論）」は明確に違うんだね。パラメータ集合のどこでも一定値をとる事前分布は一見「データを観測するまではパラメータはどの点なのかわからない」といった差し障りがなさそうなものにみえるけど、その実「パラメータ集合のどの点である可能性も等しくある」という確固たる信念だってわけだね。だから、一定値の事前分布でのベイズ推測はコインの表が出る確率を大きめに見積もる。コインを4回投げて表が出た回数が1回だけだったにもかかわらず、「表が出る確率は 0～1 のどの可能性も等しくあるはずだ」という事前の信念が強いからね。…もっとも $w(1-w)^{n-1}$ を事後分布として考えるとしても、17ページの最後の行にあるように、それで平均するんじゃなくMAP推定するなら最尤推測の結果と同じになるよ。だってMAP推定は「 $w(1-w)^{n-1}$ の最大点をとってしまう」だからね。最尤推測と同じだ。

確かに。じゃあ「パラメータ集合がコンパクトであり事前分布が一定値ならば、ベイズ推測は最尤推測と同じである」といいうのはベイズ推測＝MAP推測と定義した場合のことだったんでしょうか？

あるいはこの本の定義のベイズ推測でも逆温度 $\beta$ を無限大にすれば最尤推測になるけど…そんなベイズ推測の定義は考えにくいか。他にどんなベイズ推測の定義を想定してそのくだりがかかれたかはよくわからないな。そもそも上で考えた例も適当かわからないし。

いまいちすっきりし切りませんね…まあいいです。そういえば、この本は章末問題がありますよね。1章の章末問題をやってみたんですが、1問目から解けなかったんです。だいたい、左辺では $\beta$ に関する下限だったのが右辺では $w$ に関する下限になっているんですよ…。

自由エネルギーの逆温度に関する下限が、確率モデルの負の対数尤度のパラメータに関する下限に等しい、か。逆温度 $\beta$ が大きいほど分配関数は大きくなるから自由エネルギーを小さくできそうだけど、自由エネルギーには逆温度の逆数がかかっているから… $\beta \to +\infty$ でどこかに収束するのかな？まあその収束先が右辺なんだろうけど…。

どうすればこの右辺になるのかかなり考えてもわからなかったんですが、この本よくみたら章末問題の解答がついているんですよ。もうこの本に対する好感度が非常に上がりましたね。

そ、そっか…。ていうか先に解答ついてるか確認しようよ…。

つまり、分配関数は確率モデルの負の対数尤度でかけますから、そこを確率モデルの負の対数尤度を最小にする（分配関数を最大にする） $\hat{w}$ における負の対数尤度に置き換えれば自由エネルギーを下から評価できます。自由エネルギーはこれより小さくならない、とできます。しかし $F_n(\beta) \geqq n L_n(\hat{w})$ とかけたところでこれは $F_n(\beta)$ の下限が $n L_n(\hat{w})$ であることを意味しませんね（というネタがここ数日ツイッターで流行っているような）。なので、あるときに $F_n(\beta)$ は $n L_n(\hat{w})$ に近くなることを示さなければなりません。分配関数はパラメータの積分範囲を狭くすれば小さくなりますから、確率モデルの負の対数尤度が最小値＋ ε 未満となる領域だけに積分範囲を制限することで自由エネルギーを上から抑えることはできます。任意の ε > 0 について ε よりも小さくできます。のでこの問題の左辺と右辺の下限は等しかったんです。

解答の前半は違和感ない気がするかな…後半が、「分配関数の対数に逆温度の逆数をかけたものはこれより小さくならない」ってことだけど、分配関数の積分の領域をすごく狭くしてるイメージがあるから、そんなに狭くして下限が抑えられるのがちょっと不思議…だけど、逆温度が大きいときはこの狭い領域に分配関数は局在してるわけだし、逆温度を小さくすれば分配関数は局在しないけど「逆温度の逆数」が大きいから結局「分配関数の対数に逆温度の逆数をかけたもの」は大きくなってしまうって感じなのかな…。

（ノート3があれば）つづく

2019-04-17

雑記

【下のノートについて】多変量正規分布があったとき、 $x$ から $y$ への線形回帰と $y$ から $x$ への線形回帰は一致しないらしいんですが、あまり納得できなかったのでどう説明したら納得できるか考えてみたんですが、やっぱりあまり納得できなかったし、あとグラフに確率密度の等高線を描いた方が一目瞭然だったんじゃないかと思いました。「こう考えれば自明では」というのがあったら教えてください。よろしくお願いします。
gistaeb8a4a5dd292f783c69c28907da96ee

2019-04-10

ベイズ統計の理論と方法：ノート2

本読み機械学習

ベイズ統計の理論と方法

作者: 渡辺澄夫
出版社/メーカー: コロナ社
発売日: 2012/03/01
メディア: 単行本
購入: 1人クリック: 4回
この商品を含むブログ (8件) を見る

前回：ノート1 ／次回：まだ

前回のノートでは確率と確率密度がごっちゃになっているところがあります。事象が起きる確率に対しては事象の選択情報量ですが、確率変数の確率密度関数に対しては選択情報量とはよばないと思います。よばないと思いますがよび方に困るので、以下「選択情報量のような量」とよんでいます（自由エネルギーのことなんですが、自由エネルギーって何っていうことばとして）。
- そんな変な呼び方をしなくても「負の対数尤度」などとよべばさしつかえはないですが、個人的にそうよぶとそれって大きいほどどういう意味だっけというのがなんかすぐわかんないので採用しません。

前回読んだ1～10ページの内容をふりかえります。この本はベイズ推測について解説する本であり、ベイズ推測とは何かというのは 1.1 節に定義がありました。つまり、「ベイズ推測する」とは

サンプル $x^n$ 内の各点が独立に従っている分布は、確率モデル $p(x|w)$ の事後分布 $p(w|X^n)$ による平均 $\mathbb{E}_w \bigl[p(x|w) \bigr] \equiv p(x|X^n)$ だろうと考える（★）

ことに他なりません。ここで事後分布 $p(w|X^n)$ というのは私たちが何となく知る事後分布ではなく、逆温度 $\beta$ を用いて定義される (1.5) 式であることに注意してください（といっても、 $\beta = 1$ の場合は $p(X^n|w)$ と $\varphi(w)$ からベイズの定理によって導かれる、私たちがよく知る事後分布と同じですが）。

しかし、「～だろうと考える（★）」というだけでは、「じゃあそう考えれば？」という話です。この本で私たちが学ぶのは、「『～だろうと考える（★）』というのは結局どういうことなのか」でしょう。

6ページの1.～3.の問題を考えていくことでそれに迫るんだね。きっと「『～だろうと考える（★）』ことで、こんな条件下で、こんな風によい推測を達成できる」という感じの出口になるのかな。そういう出口があることがわかっているから、私たちはベイズ推測を利用することができる。そういう出口を知らずにベイズ推測をするなら、「理由はないけどとりあえず（★）だろうと考えました」ということにしかならない。

しかし「よい推測」とは何かという問題があります。私たちは現実に推測を行う場面で真の分布を知り得ないので、真の分布にどれだけ近づけたかを確認することはできません。しかしそれでも、 $\bigl( \, q(x), \, p(x|w), \, \varphi(w) \, \bigr)$ に拠らない数学的な法則が存在して、推測の「限界」を議論することができるということですが…。1章の続きではその足掛かりとして、「自由エネルギー $F_n(\beta)$ 」と「汎化損失 $G_n$ 」が定義されました。

自由エネルギー $F_n(\beta)$ は事前分布 $\varphi(w)$ と確率モデル $p(x|w)$ と逆温度 $\beta$ とサンプル $X^n$ による量で、分配関数 $Z_n(\beta)$ の対数の温度倍のマイナス1倍です。 $\beta=1$ のとき $F_n(1)$ は周辺尤度（その確率モデルでその事前分布にしたがうあらゆるパラメータの下でサンプル $X^n$ が観測される確率密度） $Z_n(1)$ の対数のマイナス1倍に等しいです。確率密度を確率のようなものと思えば、 $F_n(1)$ は「その確率モデルでその事前分布にしたがうあらゆるパラメータの下でサンプル $X^n$ が観測されるという事象の選択情報量のような量」です。これは、 $F_n(1)$ をあらゆるサンプル $X^n$ の現れ方について平均すると（＝ $q(x^n)$ で平均すると）、 $q(x^n)$ と $Z_n(1)$ の交差エントロピーになります。つまり、「真の分布」と「事前分布で平均した確率モデル」がどれだけ似ているかを測る指標になりそうです。もっとも、知り得ない $q(x^n)$ で平均するということは不可能なので議論が必要ですが…。
汎化損失 $G_n$ は真の分布 $q(x)$ と予測分布 $p^\ast (x)$ の交差エントロピーです。しかし、やはり知り得ない $q(x)$ で平均するということは不可能なので、経験損失 $T_n$ から見積るということですが…。 $T_n$ は「予測分布 $p^\ast (x)$ の下で各サンプル $X_i$ が観測される事象の選択情報量のような量の全サンプル平均」です。 $T_n$ は全サンプルからつくった櫛形の（ $n$ 本のデルタ関数が立った）経験分布と予測分布の交差エントロピーととらえることもできるかもしれません。

そして、 $\beta=1$ のとき $p^\ast (X_{n+1}) = Z_{n+1}(1) /Z_n(1)$ が成り立つので、この両辺の逆数の対数より、「予測分布（ただし $\beta=1$ ）の下で未知データ $X_{n+1}$ が観測される事象の選択情報量のような量」が「事前分布で平均した確率モデルの下でサンプル $X^{n+1}$ が観測される事象の選択情報量のような量（ $F_{n+1}(1)$ ）」から「事前分布で平均した確率モデルの下でサンプル $X^{n}$ が観測される事象の選択情報量のような量（ $F_{n}(1)$ ）」を差し引いたものに等しいことがわかります。前回はここまで読みました。…ここから $X^{n+1}$ で平均してその下の式になりますか？？右辺第1項は $X^{n+1}$ の関数ですが、右辺第2項は $X^{n}$ のみの関数であるようにみえますが…それに左辺は…左辺はこれ何の関数でしたっけ？？

右辺第1項と右辺第2項は丁寧にかくとこうかな。

$\displaystyle \int q(x_1) \cdots q(x_n) q(x_{n+1}) F_{n+1} (1) dx_1 \cdots dx_n dx_{n+1} = \mathbb{E} \bigl[ F_{n+1}(1) \bigr]$
$\displaystyle \int q(x_1) \cdots q(x_n) q(x_{n+1}) F_n (1) dx_1 \cdots dx_n dx_{n+1}$
$\displaystyle = \int q(x_1) \cdots q(x_n) F_n (1) dx_1 \cdots dx_n \cdot \int q(x_{n+1}) dx_{n+1}$
$= \mathbb{E} \bigl[ F_n(1) \bigr]$

左辺は先に $x_{n+1}$ について積分しようか。

$\displaystyle - \int q(x_1) \cdots q(x_n) q(x_{n+1}) \log p^\ast (x_{n+1}) dx_1 \cdots dx_n dx_{n+1}$
$\displaystyle = \int q(x_1) \cdots q(x_n) G_n dx_1 \cdots dx_n$
$\displaystyle = \mathbb{E} \bigl[ G_n \bigr]$

$\mathbb{E} \bigl[ \cdot \bigr]$ はあらゆるサンプル $x^{n+1}$ の出方に対する平均だから、10ページの一番下の式は確率的に変動する項を含まないね。任意の $\bigl( \, q(x), \, p(x|w), \, \varphi(w) \,\bigr)$ に対して、 $\beta=1$ のベイズ推測を実施したときの汎化損失の期待値は、ベイズ推測を実施する前の自由エネルギー（選択情報量のような量）の期待値がどれだけ増加するかに等しいことになる。

…なるほど、もし仮にその事前分布と確率モデルの下で $X^n$ を観測する事象の選択情報量も $X^{n+1}$ を観測する事象の選択情報量も常に変わらないというなら、 $n+1$ 個目のデータには全く「新たな情報」「意外さ」がない、 $X^n$ を観測したらそれがどんな $X^n$ であっても次に観測される $X_{n+1}$ が確実にわかってしまう、そんな状況だね。そんな状況では汎化損失の期待値もゼロだ。確実に $X_{n+1}$ がわかるんだから誤差は生じない。でも、 $n+1$ 個目のデータ $X_{n+1}$ に僅かでも「新たな情報」があれば、その「新たな情報」はベイズ推測に誤差を生じさせる。 $X_{n+1}$ を観測するまで得られない情報がある状況なんだから、完璧な推測はできない。汎化損失の期待値はゼロにならない。…といったけど、正確には汎化損失は「誤差」って感じじゃないね。KL情報量じゃなくて交差エントロピーだから、完璧に予測分布を $q(x)$ にしても汎化損失 $G_n$ はゼロにならない。 $q(x)$ のエントロピー（連続分布なので微分エントロピー）が理論下限だ。

え、えっと？ともかく、何が何の変数で、何が確率変数なのかややこしいですね。改めて整理します。10ページの上・中・下3箇所の数式は、それぞれ「確率密度」「選択情報量のような量」「交差エントロピー」にみえて、2ステップ更新されたようにみえるんです。でも、右辺は2ステップなのですけど、左辺は3ステップあったのですね（中から下への更新に、先に $x_{n+1}$ で積分して、次に残りの変数で積分するという2ステップの更新が含まれています）。つまり、右辺（上側の表）と左辺（下側の表）でそれぞれ主人公が以下のように交代しています。

$Z_n(1)$	事前分布で平均した確率モデルの下でサンプル $X^{n}$ が観測される事象の確率密度 or サンプル $X^{n}$ の下での事前分布で平均した確率モデルの周辺尤度（ $X_n$ により確率的に変動する確率変数）
$F_n(1)$	事前分布で平均した確率モデルの下でサンプル $X^{n}$ が観測される事象の選択情報量のような量 or サンプル $X^{n}$ に対して事前分布で平均した確率モデルを仮定したときの系の自由エネルギー（ $X_n$ により確率的に変動する確率変数）
$\mathbb{E} \bigl[ F_n(1) \bigr]$	「真の分布」と「事前分布で平均した確率モデル」の交差エントロピー or 真の分布の下でのあらゆるサンプル $X_n$ の出方に対する、事前分布で平均した確率モデルを仮定したときの系の自由エネルギーの期待値（確率変数ではない）

$p^\ast(x)$	予測分布の下で点 $x$ が観測される事象の確率密度（ $X_n$ により確率的に変動する $x$ の確率密度関数）
$\displaystyle \log \frac{1}{p^\ast(x)}$	予測分布の下で点 $x$ が観測される事象の選択情報量のような量（ $X_n$ により確率的に変動する $x$ の確率密度関数）
$G_n$	点 $x$ の真の分布と予測分布の交差エントロピー（ $X_n$ により確率的に変動する確率変数）
$\mathbb{E} \bigl[ G_n \bigr]$	真の分布の下でのあらゆるサンプル $X_n$ の出方に対する、点 $x$ の真の分布と予測分布の交差エントロピーの期待値（確率変数ではない）

右辺の表（上側の表）の方のこの色の文字は右辺の主人公を私たちがいくぶんよく知る統計学の言葉でかいたもので、この色の文字はこの本で出てくる言葉にならったものです。意図していなかったんですが自然とこの色の文字の方が「モデルはぶれない、サンプルがぶれる」という、伝統的な頻度論的な統計学に寄った表現になっていますね。比べてこの色の文字は「サンプル所与の下での、モデルの関数」といった感じです。左辺については $G_n$ までは予測分布 $p^\ast (x)$ というものがあるという前提でのこの分布に対する話なので「サンプルとモデルどっちがぶれる？」という要素はありません。ただ、 $G_n$ から $\mathbb{E} \bigl[ G_n \bigr]$ のステップは「モデルがぶれる」という立場ですね。モデル $p^\ast (x)$ がぶれないならば、汎化誤差の「期待値をとる」などという操作は考えられないはずです。

…汎化損失と自由エネルギーにある関係が成り立つのはわかりました。11ページは、なぜ確率密度の対数のマイナスをとるのかという話をしていますね。…これ、「確率密度は $e^{-E}$ の形であることが多いので $E$ を取りたいから」「 $E$ はエネルギーと実感できるから」って、どちらも突拍子もなく感じるんですが。なぜ確率密度がそんな形をしていることが多いなどといえるんです？だいたいどこからエネルギーが出てきたんですか？自由エネルギーとはそのような名前なのだと割り切っていましたが、ここでは紛れもなく物理のエネルギーの話をしていますよね？

まあ確率密度が $e^{-E}$ の形であるかどうかとエネルギーの登場はおいといても、確率 $P$ を対数のマイナス1倍をとって $- \log P$ として「選択情報量」というものさしでみると「大きいほどめずらしい、それが起きたと知ったときの価値が高い出来事だ」って何となくわかりやすかったよね。こっちのものさし方が推測のよさを測るのに感覚に合ってそうだ。もちろん確率 $P$ のままでも推測の誤差を議論することはできると思うけど…でも、元々よく起きる出来事か、レアな出来事かで確率を1%誤る重大さって違う気がするよね。ある年に名古屋では年間に100日くらい雨が降って、10日くらい雪が降ったらしい。真の値より10日多めに110日雨の予報を出してもたぶんあまり怒られないけど、真の値より10日多めに20日雪の予報を出したらたぶんクレームがくるだろうし天気予報を信用してもらえなくなるよね。

いや、雪の日は実際タイヤチェーンなどさまざまな準備が必要ですから、事象自体が誤りの重大さに関係していてその喩えはあまり適切ではないのでは…まあ雰囲気はわかりますが。ただ、それなら確率の誤差ではなく誤差率をみるということもできると思いますが…。

率って扱いにくいし、誤差が小さければ誤差率は対数差分で近似できるしね。 $(P^\ast - P)/P$ がゼロに近ければ $(P^\ast - P)/P \approx \log \bigl( 1 + (P^\ast - P)/P \bigr) = \log P^\ast - \log P$ だよね。

なんと、確かにそうなりますね。

あと確率密度が $e^{-E}$ の形かってのも、何か物理法則にしたがうデータだったらそうなる見込みがあるからね。その辺は統計力学の話になるけど。

やっぱり物理の話じゃないですか…まあ $-\log P$ にさまざまな解釈が与えられるというのはいいです。次節に進みましょう。…事後分布や予測分布を「解析的に計算できない」ことが多い？「解析的に計算できない」って何ですか？ (1.5) 式や (1.8) 式ってそんなに何か困難な要素があるんですか？？

「解析的に解けない」っていうと解が既知の演算や関数でかきあらわせないって意味だね。「5次方程式は一般的に解けない」というのは一般の係数の加減乗除べき根で方程式の解をかきあらわすことができないって意味だし。ただここでいう既知の演算や関数が何かはまだわからないかな。まあこの節に「計算できる例」があるってことなんだからそこから推し量ることはできるんじゃない？

そうですね、この節の例を計算してみましょう。 $\displaystyle p(x|w) = v(x) e^{f(w) \cdot g(x)}$ で $\displaystyle \varphi(w|\phi) = \frac{1}{z(\phi)} e^{\phi \cdot f(w)}$ ということです。 $f( \cdot)$ と $g( \cdot)$ の制約はわかりませんが確率モデルが積分できることは確かですね。このとき、分配関数及び事後分布を計算すると、
$\displaystyle Z_n(\beta) = \int_W \frac{1}{z(\phi)} e^{\phi \cdot f(w)} \prod_{i=1}^n \bigl( v(X_i) e^{f(w) \cdot g(X_i)} \bigr)^{\beta} dw$
$\displaystyle \qquad \; \; \, = \frac{1}{z(\phi)} \left( \prod_{i=1}^n v(X_i)^\beta \right) \int_W e^{\phi \cdot f(w) + \sum_{i=1}^n \beta f(w) \cdot g(X_i)} dw$
$\displaystyle \qquad \; \; \, = \frac{z \bigl( \phi + \sum_{i=1}^n \beta g(X_i) \bigr)}{z(\phi)} \prod_{i=1}^n v(X_i)^\beta$
$\displaystyle p(w|X^n) = \frac{1}{Z_n} \frac{1}{z(\phi)} e^{\phi \cdot f(w)} \prod_{i=1}^n \bigl( v(X_i) e^{f(w) \cdot g(X_i)} \bigr)^{\beta}$
$\displaystyle \qquad \quad \; \; \; = \frac{1}{Z_n} \frac{1}{z(\phi)} \left( \prod_{i=1}^n v(X_i)^\beta \right) e^{\phi \cdot f(w) + \sum_{i=1}^n \beta f(w) \cdot g(X_i)}$
$\displaystyle \qquad \quad \; \; \; = \frac{1}{z \bigl( \phi + \sum_{i=1}^n \beta \cdot g(X_i) \bigr)} e^{\bigl( \phi + \sum_{i=1}^n \beta g(X_i) \bigr) \cdot f(w)}$
$\displaystyle \qquad \quad \; \; \; = \varphi \left(w \middle| \phi + \sum_{i=1}^n \beta g(X_i) \right) \equiv \varphi(w | \hat{\phi})$
これは…ハイパーパラメータ $\phi$ が更新された形になりましたね。

…確率モデル $\displaystyle \prod_{i=1}^n p(X_i|w)^\beta$ と事前分布 $\varphi(w|\phi)$ の積が $a \varphi(w|\hat{\phi})$ の形にかける（ $a$ は $w$ に依存しない係数）なら $Z_n(\beta) = a$ になって事後分布が $\varphi(w|\hat{\phi})$ になるってことか。逆に事後分布をこの形にしたいなら…やっぱり $\displaystyle \prod_{i=1}^n p(X_i|w)^\beta$ と $\varphi(w|\phi)$ の積が $a \varphi(w|\hat{\phi})$ の形にかけないといけない。 $\displaystyle \prod_{i=1}^n p(X_i|w)^\beta$ と $\varphi(w|\phi)$ の積の $w$ への依存性はこの形でないといけない。「事前分布を事後分布にする」ことを「ハイパーパラメータ $\phi$ を更新する」形で達成することは、以下を満たす確率モデルと事前分布を選ぶことと同じように感じる。

$\varphi(w|\phi) \displaystyle \prod_{i=1}^n p(X_i|w)^\beta = a \varphi(w|\hat{\phi})$ 指数型分布と共役な事前分布はばっちりこれを満たす。指数型分布でない $p(x|w)$ であったら絶対にこれを満たさないのかはすぐわからないな…でも、14ページをみるとどうもそうみたいだね。

「解析的に解ける」というのは結局「パラメータの更新で済む」という意味だったんでしょうか….？もちろん、「パラメータの更新で済む」ような確率モデルと事前分布を選ばなければならないという意味ではないと14ページにありますね。あくまで人間の都合です。…1.3節はさまざまな推測方法ということですが…最尤推測って $\beta=\infty$ のベイズ推測なんですか？？

最尤推測では $w$ の分布でなくただ1つの $w$ を求めることになる。 $p(w|X^n)=\delta(w - w_{ML})$ って感じかな。15ページの $\varphi(m|\hat{\phi_1}, \hat{\phi_2})$ は $\beta \to \infty$ でデルタ関数に近づくと思う。近づく先は $\sum_{i=1}^n X_i / n$ だからパラメータ $m$ の最尤推定値だね。一般の場合にも最尤推定値にたつデルタ関数っていうのを示すのはどうやるのかすぐ思いつかないけど…。

2019-05-07 追記
＞一般の場合にも最尤推定値にたつデルタ関数っていうのを示すのはどうやるのかすぐ思いつかない
以下のような感じな気がします。
https://twitter.com/CookieBox26/status/1125552313572020224

事後確率最大化というのは、ベイズ推測をした上で事後分布の最頻値を採用する方式なのですかね。他方、平均値を採用するのが平均プラグイン？なぜプラグインというのか…。他の方式としては…平均場近似？？

148ページからかいてあるよ。これはどちらかというと、ベイズ推測において事後分布を平均場近似によって得ることを変分ベイズ法というのかな。

17ページの (3) にある場合にベイズ推測と最尤推測が同じといえるのかいえないのかといった気になることがかいてありますね…まあいまはとばしますか。18ページからは解析的に解ける場合で数値実験しているようですね。…節変わって23ページに、「確率モデルが仮のものである場合」とありますね。もちろん確率モデルがわかっている場合とわかっていない場合とあるのはわかるんですが、だから何を言いたかったんでしょうか…。

「1000人の試験結果のときと100万人の試験結果のときで推測されるパラメータがずれていてもいい」って感じに読めるね。仮に100万人の試験結果が正規分布を $J$ 個混合した分布で上手くフィッティングできても「中学生全体が $J$ 個のグループに分けられると結論されたのではない」というのは、 $J$ は「100万人の試験結果の解析に有用なパラメータ」に過ぎないってことだね。

じゃあ、中学生全体の試験結果が何個にクラスタリングされるかを研究することはできないってことですか？

やっちゃいけないのは、特に根拠なく混合正規分布を仮定して、100万人の試験結果から推測した $J$ でもって「中学生というものは $J$ 個のグループに分けられます」と結論付けることだ。「この100万人の中学生は」と付けるなら間違ってはないだろうけどね。本当に「中学生は $J$ 個のグループに分けられるのではないか」という仮説を検証したいなら、きっと人数を増やして収束していくかとか確認すべきで、現実にはそこで頓挫するんじゃないかな。

（ノート3があれば）つづく

2019-03-28

カーネル法入門：ノート2

本読み機械学習

以下の本を読みます。キャラクターは架空のものです。解釈の誤りは筆者に帰属します。お気付きの点がありましたらコメント等でご指摘いただけますと幸いです。

カーネル法入門―正定値カーネルによるデータ解析 (シリーズ多変量データの統計科学)

作者: 福水健次
出版社/メーカー: 朝倉書店
発売日: 2010/11/01
メディア: 単行本
クリック: 19回
この商品を含むブログ (11件) を見る

以前カーネル法を勉強したときの記事：パターン認識と機械学習下（第6章：その6〈終〉） - クッキーの日記

前回：ノート1 ／次回：まだ

前回のあらすじをまとめると、

データを分析するには（つまり、ある変数とある変数が依存し合っているかを調べたり、それに基づいて回帰したり分類したりするには）、線形な依存性のみを調べるのでは不十分なこともあるので、非線形な依存性も調べたい。
言い換えると、データを $\Phi : \Omega \to H$ なる特徴写像で特徴空間 $H$ に写像することで非線形な特徴も引き出して、この特徴空間で分析をしたい。
しかし、あらゆる高次項を特徴に追加したのでは、組み合わせが爆発してしまい、主成分分析のために標本共分散行列の固有値分解をしたり、回帰分析のためにデータ行列の擬似逆行列を求めたりすることができない。
が、少なくとも特徴空間における主成分分析は、行列の各成分が $k(X_i, X_j) = \langle \Phi(X_i), \Phi(X_j) \rangle$ であるような行数と列数がデータサイズの（有限の大きさの）行列の2乗の固有値分解をすることで達成できる。ので、 $k(X_i, X_j)$ の値さえ評価できればよいことになる。ここでポイントとなったのは、

特徴空間における第1主軸 $f$ は $\displaystyle f = \sum_{i=1}^N a_i \Phi(X_i)$ の形でかける（もしかけない成分があったら、かけない成分をゼロにしてその分かける成分を引き伸ばした方が各データの $f$ との内積の分散が大きくなる）。
そうすると解くべき最大化問題の中にデータは $k(X_i, X_j)$ の形でしか出てこない。

…普通に考えれば $\Phi(X_i)$ の標本分散共分散行列を固有値分解しなければならないところを、 $f$ は $H$ の部分空間 $H_0$ から探せば十分ということに着目して $f$ を表現し、最適化問題から $\Phi(X_i)$ を排除してしまうとは、なかなかトリッキーですね…。

カーネルトリックというくらいだからね。

しかし、 $f$ は $H_0$ から探せば十分、というのは通常の主成分分析でも成り立つ話ですよね？通常の主成分分析の言葉でいうと、第1主軸 $u$ を求めるのに、 $u$ をあらゆる向きの $m$ 次元単位ベクトルから探すのではなくて、 $\displaystyle u = \sum_{i=1}^N a_i X_i$ でかけるものを探せば十分、ということになりますよね。

まさしくそうだけど、通常の主成分分析では固有値分解すべき標本分散共分散行列が $m \times m$ の大きさだったのが、カーネルPCAとして解くと固有値分解すべき行列が $N \times N$ の大きさになるから、どちらがありがたいのかという話だね。

あー、元データの次元数 $m$ よりデータサイズ $N$ の方が大きければかえって計算量が大きくなってしまいますね…。8ページからは主成分分析からまた変わって、リッジ回帰という例がありますね。リッジ回帰とは、線形モデルを最小2乗フィッティングするのにパラメータベクトルのL2ノルムの2乗に比例する正則化項を付けて、「誤差2乗和を小さくするパラメータベクトルを求めよ。ただし、パラメータベクトルのL2ノルムの2乗は小さく抑えること」といったもののようですね。その解 $\hat{a}$ は…「容易に求められるように」？何気ない「容易に」が読者を傷つけるのでは？

本だと尺の都合があるからね。最小化したいのは以下だよね。

$\displaystyle \sum_{i=1}^N || Y_i - a^\top X_i ||^2 + \lambda || a ||^2$ これを $a$ で偏微分しようか。ベクトルで微分するときの公式をつかうよ。 $\begin{split} \displaystyle \frac{\partial}{\partial a} \left( \sum_{i=1}^N || Y_i - a^\top X_i ||^2 + \lambda || a ||^2 \right) &= 2 \sum_{i=1}^N ( Y_i - a^\top X_i ) X_i + 2 \lambda a \\ &= 2 \sum_{i=1}^N Y_i X_i - 2 \sum_{i=1}^N (a^\top X_i) X_i + 2 \lambda a \\ &= 2 X^\top Y - 2 X^\top Xa + 2 \lambda a \end{split}$ あとはこれが $0$ になるように $a$ を解けばいい。8ページ中ほどの式になるのはわかるんじゃないかな。

まあそこまで説明してもらえばわかります。その式の下の、リッジ回帰はどんなときに利用されるかというので、「 $X^\top X$ が退化していたり」というのは何ですか？

退化というのは $X^\top X$ が正則でない、逆行列をもたないという意味だね。もし今の回帰問題で正則化項がなかったら、つまり、 $\lambda = 0$ だったら、回帰係数の解は $\hat{a} = (X^\top X)^{-1} X^\top Y$ になるけど、 $X^\top X$ が逆行列をもたなかったらこれは求まらない。

確かに $X^\top X$ が逆行列をもっていなかったら $\hat{a}$ は求まりませんね。いえでも、解けないというのはどういうことなんでしょう？それに、解けないから正則化項 $\lambda || a||^2$ を導入して解決、というのもどういうことなのか…。

$X^\top X$ が正則でないというのは、 $X^\top X$ の階数が $m$ に満たない、さらに言い換えると $X$ の階数が $m$ に満たないって状況だね。そうなっちゃうのは以下のようなときかな。

データ数 $N$ を上回るほど、データの次元 $m$ が巨大であるとき。
データセットが線型従属なとき（例えば、同じ点が重複して含まれているとか）。

こんな状況だと、誤差2乗和を最小にする $\hat{a}$ は無数に出てきちゃうんだよね。そうだな…以下の PRML の 146 ページに Figure 3.4 があるよね。

Christopher Bishop at Microsoft Research

Figure 3.4 の左側の図を例にすると、これはちゃんと $X^\top X$ が退化していないケースで、誤差2乗和（青い等高線）が最小になる点があるんだよね（青い点）。誤差2乗和はこの青い点を底にしたお椀になってる。でも、 $X^\top X$ が退化しているとこれみたいなお椀にならない。スノボのハーフパイプみたいにずーっと谷底が続いてたり、最悪まっさらな平面で「全体的に底」って感じになっちゃう。そうなると最小点を選べなくて困る。その解決策として、「じゃあ別のお椀をもってくればいいじゃん」というのが正則化項だ。Figure 3.4 の左側の図の赤い等高線が原点を底にした「L2ノルムのお椀」だね。このお椀を重ね合わせれば最小点が求まるって寸法だね（※）。もちろん、元々の $X^\top X$ が退化していない場合でも正則化のお椀を重ねていい。Figure 3.4 は誤差2乗和の青いお椀と正則化の赤いお椀を重ねて、重ねてできた新しいお椀の底になるのが $w^\ast$ ってことだね。

※	厳密には、もし $X^\top X$ が固有値 $-\lambda$ をもっていたら、正則化項を導入しても $X^\top X + \lambda I_N$ は正則になりません。が、ピンポイントでそうなることはあまりないのではないかと思っています。

なるほど…誤差2乗和のみでは回帰係数が選びきれないときでも選ぶ基準を新たに与えたんですね。原点を底にしたお椀ということは「選べなかったら原点になるべく近い解にしといて」って感じですよね…？

そうだね。ちなみに、L2ノルムのお椀は Figure 3.4 の左側の赤い等高線のように丸いけど、L1ノルムのお椀はその右側のようにひし形をしてるね。お椀っていうかひっくり返した四角錐だね。だから、L1ノルムで正則化するラッソ回帰では一部の成分がゼロになるような解が求まりやすい。この図でも $w_1 = 0$ の点が解になってるね。

じゃあラッソ回帰は「なるべく平面 $w_1 = 0, \; w_2 = 0, \; \cdots$ に貼り付く解にしといて」って感じですかね。…まあラッソ回帰は置いておいて、今回はリッジ回帰を特徴空間でやりたいわけです。特徴空間でのリッジ回帰の目的関数は式 (1.9) ですね。また特徴空間 $H$ での内積とノルムが出てきていますね…。

データ空間でやっていた「 $X_i$ の各成分を回帰係数 $a$ で足し合わせてそれを予測値としようとする」というのは、特徴空間 $H$ でやろうとすると「 $\Phi(X_i)$ とベクトル $f$ との内積をとってそれを予測値をしようとする」ということだからね。

しかし $\Phi(X_i)$ というベクトルはいま得体の知れないもので、カーネルPCAのときは最終的に $\Phi(X_i)$ は排除されたのですよね。リッジ回帰ではどうなるのでしょう。リッジ回帰でも解は $\displaystyle f = \sum_{j=1}^N c_j \Phi(X_j)$ の形でかけるとありますね。

もし $f$ に直交補空間の成分 $f_\bot$ があったらそれを丸々削った方が正則化項の $|| f||^2$ が抑えられるからね。

全くですね。ということは…今回もまた最小化したい式の中にデータは $k(X_i, X_j)$ の形でしか出てきませんね。 $\Phi(X_i)$ は出てきません。

最小化したい式の中に出てこないだけじゃなく、予測値を計算するときにも $\Phi(X_i)$ が出てこないことが肝心だね。 $\Phi(X_i)$ を直接取り扱いたくないから。でも未知の点 $x$ に対して「 $\Phi(x)$ とベクトル $f$ との内積をとるということをするとやっぱりデータは内積 $k(X_j, x) = \langle \Phi(X_j), \Phi(x) \rangle$ の形でしか出てこない。

えっと、まとめると、特徴空間でそのままリッジ回帰する（予測したい未知データの特徴 $\Phi(x)$ と内積を取るべき特徴空間内のベクトル $f$ そのものをもとめる）のではなく、 $f$ が $\displaystyle f = \sum_{j=1}^N c_j \Phi(X_j)$ の形でかけることに注目して $c_j$ を求める問題にしたのですよね。うーん、カーネルPCAでもカーネルリッジ回帰でも「解はこの形でかけるじゃん」というところがどうも降って湧いたようなんですが…。

最適化問題中の $\Phi(X_j)$ を意地でも $\langle \Phi(X_i), \Phi(X_j) \rangle$ の形にするために、 $f$ から $\Phi(X_i)$ を絞り出したって考えればそんなに降って湧いてもないんじゃないかな。

絶対に $\Phi(X_i)$ を直接取り扱いたくないのだという強い意思といったところですか…。ところで、カーネルPCAとカーネルリッジ回帰ではデータを特徴写像 $\Phi$ で特徴空間 $H$ にとばしましたが、 $\Phi$ や $H$ の詳細には立ち入っていませんでした。 $H$ についてわかっているのは、ベクトル空間であることと、内積 $\langle \cdot, \cdot \rangle$ をもつということだけですね。 $\Phi$ もどんな写像なのかまるでわかりません。直接取り扱っていないので。常にカーネル関数 $k$ を通して取り扱っています。…この $H$ と $\Phi$ は適当に決めていいんでしょうか。

ではないね。まずカーネル関数 $k(X_i, X_j) = \langle \Phi(X_i), \Phi(X_j) \rangle$ の値が評価できないといけないからこれが発散するような $\Phi$ だと困ると思うし…そもそも私たちは直接 $\Phi$ を決めないし…まあ [tex\Phi] とは何かは次章だね。

そもそも非線形な依存性を調べたくて $\Phi$ などという特徴写像を持ち出したのですから、どのような $\Phi$ が可能なのかが肝心な気がしますが…。ところで今回カーネルリッジ回帰を学びましたが、python の scikit-learn でカーネルリッジ回帰ができるようなので実行してみます。iris データの全品種ごちゃ混ぜで、がく片の長さからがく片の幅を予測するという何がしたいのかわからない例ですが…カーネル関数 $k$ には色々なものが指定できるし、自分で実装したものを渡すこともできそうです。組み込みのものを利用するときはおそらく以下の pairwise.kernel_metrics にあるものが指定できると思います。以下のスクリプトではいじっていませんが、パラメータも変更できます。

sklearn.kernel_ridge.KernelRidge ‐ scikit-learn 0.20.3 documentation

sklearn.metrics.pairwise.kernel_metrics ‐ scikit-learn 0.20.3 documentation

f:id:cookie-box:20190328234405p:plain:w480

もともとが線形関数でも非線形関数でもフィッティングできなさそうなのでそもそも例として適当ではありませんが、線形フィッティングは駄目そうですね。 $k$ にpolyカーネル $k(x, y) = (x^\top y)^3$ を選んだときよりもrbfカーネル $k(x, y) = \exp(-\gamma || x-y ||)^2$ を選んだときの方が点の密度が大きいところを通れているといえばいるような。

（ノート3があれば）つづく

スクリプト

import numpy as np
from sklearn.kernel_ridge import KernelRidge
from sklearn.datasets import load_iris

%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 8, 6
rcParams['font.size'] = 16

iris = load_iris()
#x = np.array(iris['data'][:,0][iris['target']==0]) # setosa の sepal length (cm)
#y = np.array(iris['data'][:,1][iris['target']==0]) # setosa の sepal width (cm)
x = np.array(iris['data'][:,0]) # sepal length (cm)
y = np.array(iris['data'][:,1]) # sepal width (cm)

print(x.shape)
print(y.shape)

x = x.reshape(-1, 1)
y = y.reshape(-1, 1)

print(x.shape)
print(y.shape)

# ----- オリジナル -----
plt.scatter(x, y, label='original')
# ----- linearカーネル -----
kr = KernelRidge(kernel='linear')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(linear)')
# ----- polyカーネル -----
kr = KernelRidge(kernel='poly')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(poly)')
# ----- rbfカーネル -----
kr = KernelRidge(kernel='rbf')
kr.fit(x, y)
plt.scatter(x, kr.predict(x), label='predict(rbf)')

plt.legend()
plt.show()

(150,)
(150,)
(150, 1)
(150, 1)

2019-03-22

カーネル法入門：ノート1

カーネル法

カーネル法入門―正定値カーネルによるデータ解析 (シリーズ多変量データの統計科学)

作者:健次, 福水
朝倉書店

Amazon

以前カーネル法を勉強したときの記事：パターン認識と機械学習下（第6章：その6〈終〉） - クッキーの日記

次回：まだ

まえがきから読み始めると冒頭からカーネル法とは何かきっぱり書いてありますね。つまり、「正定値カーネルの定める再生核ヒルベルト空間を用いたデータ解析の方法論」とのことです。…データ解析に「空間を用いる」とは何でしょう？「ほにゃららモデルを用いる」、ならわかるんですが…まえがきの続きによるとこの再生核ヒルベルト空間とは関数空間であるらしいですが、何が何やら…まあ置いておきましょう。冒頭にはカーネル法とは「データの非線形性や高次モーメントを扱うための方法論」ともありますね。データの非線形性や高次モーメントを扱うとは？

おそらく1章の2～3ページで挙げられているのが最も簡単な例だね。「データの線形性しか扱わない」というのは、「データの空間のどこかにデータを仕切るような超平面を置いて、その超平面から表側に垂直に離れるほどプラスのスコア、裏側に垂直に離れるほどマイナスのスコア、のようなスコア付けでもってデータを分類したり回帰したりしようとする」ということだろう。でもまっすぐな仕切りで仕切れないのが図1.1の例だね。「高次モーメントを扱う」というのは「2つの変数 $X, \, Y$ に対して $X$ の $n$ 次かつ $Y$ の $m$ 次のモーメントからわかる依存関係を捉える」って感じなんだろうか。図1.2は $X$ の方を2乗しないと相関係数が小さくなる例だね。相関係数がゼロであることは独立であることを意味しないというのは有名だね（ウィキペディアの「相関係数」の頁の絵）。

相関係数 - Wikipedia

そこまでごり押されなくても「線形の解析では不十分」であることに異存はないですよ…。確かに、図1.1で分類しようとしたり、図1.2で $X$ がどうなると $Y$ がどうなるかという関係を見出したりするには線形関数では足りませんね…もっと柔軟なモデルでなければ…解決策を思い付きました。決定木を使いましょう。

解決できるのはできると思うけど、本のタイトルが「決定木入門」になっちゃうかな。

図1.1の方は図の2次元空間の上に $f(x_1, x_2)$ という「その場所では＋か ○ のどちらか」という関数を見つけたいのですよね。それを $f(x_1, x_2) = \hat{f}(\phi(x_1, x_2)) = \hat{f}(z_1, z_2, z_3)$ とすることで $\hat{f}$ を線形関数で足りるようにしたんですね。図1.2の方は、 $f_1(x)$ と $f_2(y)$ の相関を最大化するのに、 $f_1(x) = \hat{f_1}(\phi(x)) = \hat{f_1}(x^2)$ とすることで $\hat{f_1}$ を線形関数で足りるようにしたという感じでしょうか。…なんというか、これらをみるとカーネル法とはむしろ「何が何でも線形関数で解く。線形関数で解けるようにデータの方が動くべき」というものにみえますね。ニューラルネットワークの方が人としての器が大きいのではないですか？あちらは「データはそのままでどんな関数でも表現する」でしょう？

人じゃないからね。カーネル法はカーネル法で（あるクラスの最適化問題に対する正定値カーネルの定める再生核ヒルベルト空間内での）最適解が線形関数の形でかけるという理論的裏付けと、さまざまな拡張があるんだよ。「線形関数で解けるようにデータの方が動くべき」というようにデータの方が一旦動いてくれたら、その上に線形な依存性に基づいたさまざまな分析手法を応用することができるんだから。

はあ。まあ何にせよカーネル法ではデータを都合よく動かさなければならないということですよね。しかし4ページをみると、あらゆる高次のモーメントを考慮するには組合せが爆発してしまうようですが。前途多難では？

にもかかわらず効率的に高次モーメントを抽出できますよというのがカーネル法の売りだからね。

どうやってそんなことを…まあそれがこれからの話題ですか。1.1.2節に移りましょう。データの空間が $\Omega$ で、それを移す実ベクトル空間が $H$ ですか。図1.1や図1.2の左側が $\Omega$ 、右側が $H$ ということですよね？ …実ベクトル空間というのは？

集合であって、和やスカラー倍が入ってて、和やスカラー倍が結合律やら分配律やらを満たすやつ。「実」ベクトル空間といったときはスカラーは $\mathbb{R}$ の元だね。

それは結局、ベクトルを集めた集合ではないんですか？ 3次元実数ベクトルが入っている $H$ は結局 $\mathbb{R}^3$ ということではないんですか？

どちらかというと実ベクトル空間の条件を満たす集合があったとき、その元をベクトルというかな。あと、3次元実数ベクトルを含む $H$ でも、 $\mathbb{R}^3$ 内のある平面だけとかある直線だけとかの実ベクトル空間はありうるね。

あー平面や直線の中で和やスカラー倍はできますね…まあ、データを和やスカラー倍がちゃんとできる空間にもってくるのだと理解しました。そうですね、元のデータが単語や文章などであったら、和やスカラー倍をしようとしても和やスカラー倍とは何か定かではありませんね。あれでも、 $H$ は内積ももたなければならないと続きにありますね。

データ分析では点がどのように密集（分布）しているかに興味があるからね。「ここに仕切りを入れられそうだ」「これが主成分だ」とかいうには点どうしの近さが決まっていないと始まらない。

それって決めないと決まらないものですか？内積って $x_1 y_1 + x_2 y_2 + x_3 y_3$ のように授業で習った気がするのですけど。

それは実数ベクトルの標準内積だね。それに決めてもいいよ。でもそれとは違う内積に決めてもいい。例えば、 $100 x_1 y_1 + x_2 y_2 + x_3 y_3$ みたいに、1つ目の軸方向だけやたら重視する内積でもいいよ。そうしたかったらね。内積が満たすべき性質は、付録の193～194ページだね。

なんと、そういう内積もありなんですか？一旦193ページにとびますか。…って、内積より前にまずノルムというものが出てきますね。その後に内積はこういうものというのも出てきて…内積が定義できればその $1/2$ 乗をノルムにできるってことですね。なるほど、副部長がいう定義も内積の性質を満たしますね。

対称行列 $A$ について、 $x^{\top}Ay$ は内積になるからね。

まあ何にせよ何かしら内積を決めたとして、でも $H$ が高次元だと内積計算にはコストがかかると…それはそうかもしれませんが。でもカーネル法ではその内積計算がカーネル関数 $k$ の評価で済む？えっと、話を整理しましょう。いま、手元のデータを分析するのに、線形の依存性を調べるのでは不十分なので、高次元の依存性を調べたくて、つまり高次元の項も含めた内積を計算したくて、それにはカーネル値 $k(X_i, X_j)$ を計算すればよい。…ちょっと最後が飛躍しすぎですね。

それはおいおいかな。そこがカーネル法の基本原理だからね。5ページからは実用例みたいな話が始まってるね。

カーネル主成分分析、ですか。通常の主成分分析もわからないのですが…python の scikit-learn でも使えるようなのでまずは実行してみますか。

sklearn.decomposition.KernelPCA — scikit-learn 1.0 documentation

データはそうですね、非線形性があるものがいいのですよね。以下にスイスロールのデータがあったのでこれをお借りすることにします。3次元で、y軸に垂直に切ると渦巻きになっているようなデータですね。

a Swiss Roll Dataset

上のデータに sklearn.decomposition の PCA と KernelPCA を適用したのが下図です（スクリプトは記事の最下部）。左上が元のデータをxy平面に落とした図、右上が元のデータをxz平面に落とした図ですね。左下が通常の主成分分析の第1主成分-第2主成分、右下がカーネル主成分分析の第1主成分-第2主成分です。上段の実線は通常の主成分分析の第1主軸の方向、点線が第2主軸の方向ですね。データ点は上のリンク先の gif と同じ4色の塗分けをしています。よくみると、スイスロールはこの4つのクラスタに切り分けられるのかもしれませんね。色が同じ部分がやや密集しているようです。
f:id:cookie-box:20190321180703p:plain:w540

f:id:cookie-box:20190321180703p:plain:w540

f:id:cookie-box:20190321180731p:plain:w540

…通常の主成分分析はほとんどxy平面に落としただけだね。データがよく広がっている方向をとらえられたという感じはしない。元々がスイスロールなので無理はあるんだけど。それに比べてカーネル主成分分析は、最初の2つの主成分でデータが広がっている方向をとらえることができたのかもしれない。非線形性を扱えるから、丸めたカーペットのように丸まったスイスロールを少し開くことができたって感じかな。

えっと、5ページの通常の主成分分析の復習をみると、通常の主成分分析というのは、データの空間の中で、すべてのデータをその単位ベクトルに正射影したときに正射影した後のデータたちの分散が最大になるような単位ベクトルの向きを求めて、それを第1主軸とするのですね。続けて第2主軸以降も求めるときは、第1主軸と直交する空間の中でそのようなベクトルを探していくと…それが「標本分散共分散行列 $V$ の大きい方から $d$ 個の固有値に対応する固有ベクトル」？ $u^{\top}Vu$ を最大化する $u$ を求めたいというのはいいですが、それが固有ベクトルというのはどういうことですか？？

落ち着いて追いかけようか。 $u^{\top}Vu$ を $u^\top u = 1$ の下で最大化したいわけだ。じゃあラグランジュ関数を $\displaystyle L(u, \lambda) = u^{\top}Vu - \lambda \left(u^\top u - 1\right)$ とおこう。これを $u$ で偏微分する。微分公式 $\displaystyle \frac{\partial}{\partial u} u^\top A u = (A + A^\top) u$ をつかえば $\displaystyle \frac{\partial}{\partial u} L(u, \lambda) = 2Vu - 2\lambda u$ だ。つまり、 $\displaystyle \frac{\partial}{\partial u} L(u, \lambda) = 0 \; \Leftrightarrow \; Vu = \lambda u$ になる。 $u$ は $V$ の固有ベクトルでなければならないね。それも、 $u^{\top}Vu = \lambda u^{\top} u =\lambda$ を最大化したいんだから、最大の固有値 $\lambda$ に対応する $u$ を第1主軸に選ぶべきだ。第2主軸以降はそれまでの軸と直交するように選びたいわけだけど、いま $V$ は対称行列だから固有ベクトルどうしが直交している（証明は以下などを参照してね）。

対称行列の固有値と固有ベクトルの性質の証明 | 高校数学の美しい物語

なるほど、ありがとうございます。6ページからは主成分分析を今度は特徴空間 $H$ で実行しようとしているんですね。先ほどは各データ $X_i$ を $u$ 方向に正射影した $u^\top X_i$ たちの分散を最大化しようとしていましたが、今度は $\langle f, \, \Phi(X_i)\rangle$ なるものの分散を最大化しようとしていますね。 $||f||=1$ の下で。

内積 $\langle \cdot, \cdot \rangle$ とノルム $||\cdot||$ さえ決まっている空間であれば、主成分分析はできるということだね。

それで、解は $\displaystyle f = \sum_{i=1}^N a_i \tilde{\Phi} (X_i)$ の形で探せば十分？？え、えっと、なぜこの形で十分かもわからないですが、 $f$ とは何なんですか？先ほど通常の主成分分析でデータの分散が最大となる向きの単位ベクトル $u$ を探したように、 $f$ も $H$ 内の単位ベクトルというわけではないんですか？これ、ベクトルなんですか？？

無論ベクトルだよ。でも部長はまだベクトルを「成分を縦に並べたもの」とイメージしているかな。ベクトルというのはあくまで「和やスカラー倍ができる集合の元」だ。それでしかない。そして「データの分散が最大となる向き」というのは、「各データとその単位ベクトルとの内積をとったものの分散が最大になるような単位ベクトル（ノルムが1のベクトル）」だ。この「向き」は、成分を縦に並べたものである必要はないよね。

そういうことですか…ではイメージを捨てて、内積で考えるようにしましょう。…すみません、部分空間と直交補空間とは何でしょう。

うーん、よく知っているベクトルのイメージに戻っちゃうけど、 $\mathbb{R}^3$ というベクトル空間を考えて、この空間の中のある平面だけみると、この平面は和とスカラー倍で閉じているから $\mathbb{R}^3$ の部分空間だ。直交補空間というのは、その平面に垂直な直線の中で閉じている空間だね。 $\displaystyle \left( \begin{array}{c} x \\ y \\ 0 \end{array}\right)$ とかけるベクトルたちの空間と $\displaystyle \left( \begin{array}{c} 0 \\ 0 \\ z \end{array}\right)$ とかけるベクトルたちの空間って感じかな。それでいま考えている話に戻ると、 $H$ という特徴空間の中には $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ の和やスカラー倍ではつくれない「ベクトル」がある。 $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ の和やスカラー倍でつくれない、というのは、 $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ たちがつくる部分空間に射影したベクトルを差し引いてもゼロでないベクトルが残っているってことだね。ただ、もし $f$ というベクトルがそのようなゼロでないベクトルを含むようなものであったとしても、その差し引き残ったベクトルってあってもなくても $\tilde{\Phi}(X_i)$ との内積の値に影響しない。いま「データの分散が最大となる向き」 $f$ を探したいわけだけど、もし $f = f_0 \oplus f_{\bot}$ の $f_{\bot}$ がゼロでなかったとしよう。このとき $||f_0||^2 = ||f|| - || f_\bot ||^2 =1-|| f_\bot ||^2 < 1$ だけど、 $f' = 1/(1-|| f_\bot ||^2)^{1/2} f_0 \oplus 0$ とすれば、 $f'$ への射影の方が $f$ への射影よりデータの分散が単純に $1/(1-|| f_\bot ||^2)^{1/2}$ 倍大きくなるんだよね。いま $f$ のノルムは1という制約があるんだから、 $\tilde{\Phi}(X_i)$ との内積の値に影響しない成分に割り振ってる余裕はないってとこかな。

だから $\{\tilde{\Phi}(X_i)\}_{i=1}^N$ たちがつくる部分空間内で探せば十分、つまり、(1.6) 式の形で探せば十分ということですか。次元がよくわからない、というか次元という概念があるのかどうかもよくわからないベクトル空間でも、有限個の係数を探せばよいということになるんですね。そうなると、さっきとそっくりな解き方ができるんですね。しかし、先ほどと違って $a^\top \tilde{K}^2 a$ の $a$ というのはもはや「ベクトルの向き」を意味しませんね。そして、7ページの下部の式のように第 $p$ 主成分まで求まるんですね。これがカーネルPCA…。

（ノート2があれば）つづく

スクリプト

import numpy as np
import pandas as pd

df = pd.read_csv('swissroll.dat', sep='\s+', header=None)
df = df.apply(lambda x: x / x.std(), axis=0) # 標準化
df[1] = df[1] - df[1].mean() # y軸だけ中央がゼロでないのでシフト
print(df.head())
c = np.array(['tomato', 'dodgerblue', 'gold', 'mediumseagreen']).repeat(400) # 4色に色分け

# PCA
from sklearn.decomposition import PCA
pca = PCA()
pca.fit(df)
df_pca = pca.transform(df)

# KernelPCA
from sklearn.decomposition import KernelPCA
kpca = KernelPCA(kernel='rbf', gamma=1.5)
kpca.fit(df)
df_kpca = kpca.transform(df)

# プロット
%matplotlib inline
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['font.size'] = 16

x = np.array([-10, 10])

fig, axes = plt.subplots(ncols=2, figsize=(12, 5), sharex=True, sharey=True)
axes[0].scatter(df.iloc[:,0], df.iloc[:,1], c=c)
axes[0].plot(x, x * pca.components_[0][1] / pca.components_[0][0], c='black', linestyle='solid', linewidth=2) # 第1主成分
axes[0].plot(x, x * pca.components_[1][1] / pca.components_[1][0], c='black', linestyle='dotted', linewidth=2) # 第2主成分
axes[0].set_xlim([-3, 3])
axes[0].set_ylim([-3, 3])

axes[1].scatter(df.iloc[:,0], df.iloc[:,2], c=c)
axes[1].plot(x, x * pca.components_[0][2] / pca.components_[0][0], c='black', linestyle='solid', linewidth=2) # 第1主成分
axes[1].plot(x, x * pca.components_[1][2] / pca.components_[1][0], c='black', linestyle='dotted', linewidth=2) # 第2主成分

plt.subplots_adjust(wspace=0.15)
plt.show()

fig, (axL, axR) = plt.subplots(ncols=2, figsize=(12, 5))
axL.scatter(df_pca[:,0], df_pca[:,1], c=c)
axR.scatter(df_kpca[:,0], df_kpca[:,1], c=c)
plt.show()