はじめに
前回の記事
の続きってほど滑らかに繋がってはいませんが、少し突っ込んだ話に行きます。
ここでは前回、データ$x$から$y$を予測する場合のモデルが
$$ y = f(x) $$
と表せるようなケースを想定して話を進めました。その際のポイントとして、
$$ y = f(x) = w\cdot \phi(x) $$
と、（非線形）変換$phi(x)$を噛ませた後に$w$を重みとした線型結合を考えるという手順を踏みました。
$\phi(x)$の選び方や$w$の求め方に特に具体的な制限を設けずに前回は話しましたが、今回は少し具体的な話に入っていきます。
これはまだ途中書きです。
一般線形モデル
まず、$y=f(x)$の回帰の問題を考えましょう。これまでに話したとおり、
$$ y = f(x) = w\cdot \phi(x) $$
と表してしまうことにします。 $x$に対して上手い加工$\phi$を施すことができれば、$w$を使った線型結合で$y$を上手く表せるであろうという仮定をします。
これは所謂線形モデルなのですが、$\phi(x)$は非線形関数でも構いません。 ここで重要なのは、（$\phi$は固定のもとで）求めなければいけないのは$w$の方であり、$w$に対してはモデルは線形になっています。
最小二乗法
いま、このように仮定した線形回帰モデルを学習させるために下記のような損失関数を与えましょう。
$$ L(w) = \frac{1}{2}\sum_i (y_i - w\cdot \phi(x_i))^2 $$
今手元に、$D=\{ (x_1, y_1), ... , (x_N, y_N) \}$というN個のデータがあるとしましょう。 私達は今勝手に、$y = w\cdot \phi(x)$と表せるだろうと目論んでいるわけですから、本当にピッタリとそうならば上記の損失関数はゼロになってくれるはずです（理想的には$y-w\cdot \phi(x)=0$のはずだから）。
もしも少しでもズレがあるならば、損失関数の右辺は正の値を持ちます。この値が小さければ小さいほど、上手く$w$を選ぶことができたと言えるわけです。 このような指針で$w$を決めることを最小二乗法と言います。
ここで非常においしいのが、モデルが$w$の線形関数になっている点です。二乗したところで損失関数は凸関数になっており、必ず最小値を求めることが可能です。 具体的な方法は極値を見つけさえすればそこが最小値になっていることが保証できるので、微分して$0$となる点が最小値になります。
ここで線形代数のテクニックを使うと、極めて簡潔に問題を解くことができます。 まず$x_i$一つ一つはベクトルです。そしてそれを変換した$z_i=\phi(x_I)$なる新たなベクトル考えることにします（ここは固定なので、$z_i$で置いたほうが楽です）。 また、もしかしたら$z=\phi(x)$によらない定常的なスカラー値$w_0$が影響するかもしれないので、モデルを改めて
$$ y = w\cdot z + w_0 $$
のように表しておくことにしましょう。これは、形式的に$z = (1, z_1, z_2, ..., z_D)$とベクトルの最初の成分に$1$をパディングし、$w=(w_0, \cdots, w_D)$としてしまうことで （ここでの下付き添字はベクトルの成分のindexです）、相変わらず
$$ y = w\cdot z $$
と表すことができますのでそれほど気にしなくていいでしょう。さて、線形代数を使った最小二乗法の解法に入ります。
ここで$z_i$を格納（ここでの下付き添字はデータ数のindex）した行列
$$ Z = (z_1, \cdots, z_N)^T $$
を考えます。これはひとつのデータが横ベクトルの形で格納されている状態です。
同じように$Y=(y_1,\cdots, y_N)^T$というベクトル（$y_i$ひとつひとつはスカラー）を作っておけば、全てのデータの線形関係を下記のように表現しておくことができます。
$$ Y = Zw $$
損失関数は簡単に下記で書き換えられます。
$$ L(w) = \frac{1}{2}|Y-Zw|2 $$
これの微分が$0$となればいいわけですから、
$$ \frac{\partial L}{\partial w} = Z^T(Y-Zw) = 0 $$
を$w$について解けばよく、
$$ \begin{eqnarray} Z^TY-Z^TZw & = & 0 \\ Z^TZw &=& Z^TY \\ w &=& (Z^TZ)^{-1}Z^TY \end{eqnarray} $$
と直ちに求まります（$Y$も$Z$も手持ちのデータからすぐ決まる！）。
線形モデルは数学的に容易に扱うことができ、しかも線形代数と組み合わせると至ってシンプルに解けるようです。 しかし、実際にはこのように線形モデルを最小二乗法によってフィッティングすることは、（機械学習を必要とするような）多くの場合では上手く行きません。
最小二乗法の背後に潜む仮定
いま、私達は勝手に
$$ L(w) = \frac{1}{2}\sum_i (y_i - w\cdot \phi(x_i))^2 $$
という損失関数を考えてきました。これは、モデルの予測$w\cdot \phi(x_i)$が$y_i$から外れるほど損失を大きくするという、一見素晴らしい評価指針に見えます。 しかし、果たしてそうでしょうか。モデルの予測が実際の値より小さくなってしまった場合にこそ損失を与えたい（つまり上下のブレに対して非対称に損失を与えたい）だとか、 モデルの予測がちょっとずれるくらいならどうでも良いとか、状況に応じて本当は様々なケースが考えられるはずです。
また、そもそもデータの性質上、モデルの出力が実際に得られている値から上下に均等にブレると言えるのかは怪しいはずです。
実はここで紹介した最小二乗法というのは、データが発生する背後に「値が上下に均等にブレるようなノイズが混入している」という仮定を設けた場合の 統計的な最尤推定の結果に一致しているのです。
より具体的には「データ$y$はデータ$x$がとある変換$f(x)$を受けた後に、ガウス分布から生起するノイズ$\epsilon$が加算され観測される」という仮定を設けた場合の最尤推定に一致します。
$$ y = f(x) + \epsilon $$
最尤推定とは、「手元のデータが手に入った」という事実が既に存在するわけですから、「手元にそのようなデータが得られる確率は高いはず」という考えのもと、 「手元のデータが得られる確率が最大となるような$w$を求める」方法です。
上記のモデルでは確率的に変動するのはノイズ$\epsilon$のみです。コイツは、平均的には$0$ですが、とある分散でランダムに値を取ります。 ノイズが無ければ、$y$は直ちに$y=f(x)$と得られるはずですから、$y$は平均的には$f(x)$という値を取るものの、$\epsilon$の分散によって散らばる確率変数になります。
$$ y〜 N(y \mid f(x), \sigma^2) = N(y \mid w\cdot \phi(x) , \sigma^2) $$
とサンプリングされるという仮定を設けるのです。 ガウス分布は平均からずれた値は生じにくいです。すなわち、平均$f(x)$からずれるような値を$y$は取りにくいはずです。 しかし事実$y$と$f(x)$は全てのデータに対して厳密には当てはまりません。ですのでノイズが生じていると考えるのは妥当でしょう。
そして最尤推定「手元にデータが得られるような確率を最大化する」を実行した場合には、$y$が$f(x)$からずれるようなことは極力起こりにくいはずなので、 ノイズが小さく見えるような$w$が求まることになります。実際に計算をするとすぐに分かります。
手元に$D=\{(x_1, y_1), ..., (x_N, y_N) \}$が得られるような確率（尤度）は
$$ \prod_i N(y_i\mid w\cdot \phi(x_i) , \sigma^2) ∝\prod_i \exp \left\{ - \frac{(y_i-f(x_i))^2}{2\sigma^2} \right\} $$
で表され、これを$w$に関して最大化することは、最小二乗法を実行することに他なりません（$\exp$の中身の分子に注目！）。 これをしっかり解きたい場合は、尤度の最大化を行えばよく、実際には$\exp$の中身で、パラメータ$w$に関連のある部分が最大化されれば良いので
$$ {\rm maximize}_w \left[ - \sum_i \frac{(y_i-f(x_i))^2 }{2\sigma^2}\right] $$
を解けば良いことになります。これは尤度の対数を取って最大化することと全く同じです。 そして最大化問題を最小化問題（符号を反転）すると、最小二乗法とやはり全く同じ問題になります。
$$ {\rm minimize}_w \left[ {\sum_i (y_i - w \cdot \phi (x_i) )^2} \right] $$
ここで関係のない分母だとか、係数部分は無視しました。 いま私達は最小二乗法を使うことが、データの生成過程でガウスノイズが混入すると考えていると等価であることを見ました。
確率的なモデルであると考え、手元にこのようなデータが手に入る確率（尤度）を最大化すると 損失関数なるものを考えなくとも、自然と解くべき問題が導かれます。
それは裏を返せば、確率モデルの設計に対して解くべき問題が直接付随するわけで、 どのようなモデルを仮定するかが問題の結果に直接響いてくるというわけです。
最小二乗法が、今考えている問題に対して適当な解法になっているかの保証は当然ありませんし、 今回のようなガウスノイズを仮定した最尤推定が本当に正しい結果を導くかも誰も保証してくれません。
ただ、そう仮定した場合においては、（その仮定が正しいとした時）とある結果が得られるだけです。
しかし確率モデルを考えることひとつのメリットは、データがどのように生まれてきているのかを意識できる点です。
最小二乗法による解法では、とりあえずモデルの予測結果が手元のデータに一致するように帳尻を合わせる損失関数を決めて それを最小化するという手順を踏みました。しかし、もっと違う損失関数を考えてみたくなった時、その設計はいかにも抽象的です。
確率モデルならば、例えば、出力は２値である（男か女である）などとした場合には、出力がベルヌーイ分布であるとか、 ノイズは値を増やす方向にしか働かないと分かっていれば、ポアソン分布などのノイズ項を考えることができます。 ある種、モデルの設計手段に対して柔軟性持たせることができるのです。
まとめ
かなり駆け足で来てしまいましたが、以下のまとめの事実を知っておくだけでも有益でしょう。
- 最小二乗法とガウスノイズを仮定した最尤推定は同じ結果を出す（しかし、使っている人の考えていることや面持ちは違うかもしれない）。
- いずれの方法も、平均から上下に均等値がずれるような場面を想定している（ガウスノイズや、最小二乗法の損失関数を見れば何となく分かる）。
- これらの手法で導かれる結果は、上記の仮定が正しいとしたもとで導かれる（仮定が正しくなくとも、もしそうだとしたら、という観点で何かしらの結果は出る）。
- ノイズに仮定する分布などをいろいろ変えると、最尤推定の式がそれに応じて変わる（その式は、単に損失関数をいじろうと思った場合は非自明に見える）。
実際にどのような手法を使うかは場面に応じてマチマチだとして、これらの考え方があることを知っておくのは重要だと思います。
ちょっとした例題
- 私の年齢が分かっている
- 私の誕生日は4月
- 数年分の月単位の収入データがある
- ただし、ブログなどの副収入が含まれている。
- 勤務先は歴史ある企業で年功序列
（※フィクションです）
という時に、私の年齢に応じた会社からの収入を回帰して予測したい。
このような場合にどのように分析をするでしょうか。 月の手取りのデータを合算して年収を算出してもいいでしょう。つまり、ある年齢の時の年収という形式にしてモデルを作ってもいいですし、 もっと細分化して、4月が誕生日ですから入社から月単位で回帰しても良いはずです。
ここで重要なのは、私の会社からの収入が知りたいということです。 ブログの収入は含めません（そしてそれがどれほどかは分からない）。 しかし、ブログの収入を含んだ収入のデータが手元にあります。
このとき、私の年齢（あるいは月齢？）を元に収入を回帰する問題を最小二乗法で解くのは得策でしょうか。 解いたとしたら、それはどのようなデータになるだろうか。何を何で回帰し、何をノイズと見なすのか。 場合によってはある種の信号源推定にも見えます。
ひとつ言えるのは、会社からの収入を年齢（月齢）に応じて平均的に獲得し、 ブログの収入がノイズとして得られるようなモデルを考えた場合、最小二乗法（あるいはノイズがガウス分布）という方法は得策ではありません。 なぜでしょうか。
一般化線形モデル
そこで、$y = f(x) = w\cdot \phi(x)$なる関係を基本に、これを大幅に拡張したものが一般化線形モデルと呼ばれるモデル群です。 一般化線形モデルでは、ノイズの性質やあるいは$y$という出力自体にも何かしらの加工を施します。
出力を加工する
まず出力を加工するという概念を見ましょう。
例えば出力$y$が収入や体重など絶対に負の値を取らないようなものであったとしましょう。 もしも何も考えず線形モデル
$$ y = f(x) = w\cdot x $$
を使ってフィッティングをした場合おかしなことが起こりそうです。 仮にフィッティングの結果が下の図のようになったとしましょう。
手持ちのデータに対してはそれらしいフィッティングになっていそうですが、引いた直線は左のほうで負の値を取っています。 私達はいま、体重であるとか収入であるとか、正の値を予測しようとしているのにこれは明らかに変です。
$x$という値が増えれば$y$も増えそうな傾向はあるのですが、$y$が負になるようなことは絶対に無い！と分かっているのであれば、 仮定するモデルにこのことを反映させたいと考えます。その最も単純な方法が
$$ \log y = w\cdot \phi(x) $$
というモデルを考えてしまうことです。 なんと$y$の方に適当な関数を被せてしまいました。 $\log$の定義域（引数の取れる範囲）は$0$以上です。 このようなモデルを考えることは、$y$が$0$より大きくあれという要請です。
しかしそんなことして良いのか？と思いますが、そもそも線形モデルで表して良いのか？ ってことだって分かりませんし、だから、このようにして良いのかも分かりません。 ただ、できる限り知りうる知識を反映した形でモデルを表したいというだけの話なのです。
ここで、上記のやり方がわかりにくいと思うのであれば、上記と等価の以下の式を見てください。
$$ y = \exp(w\cdot \phi(x)) $$
なるほど確かに$y$は$w\cdot x$がどんな値を取ろうとも正になってくれそうです。 特に$\log$である必然性など無いので、状況に応じて適したものを用いると良いでしょう。 それは例えば
$$ \frac{1}{y} = \exp(w\cdot \phi(x)) $$
でも良いですし、
$$ y^2 = \exp(w\cdot \phi(x)) $$
でも良いですし、
$$ \frac{y}{1-y^2+2y} = \exp(w\cdot \phi(x)) $$
でもいいかもしれません（これがデータを上手く表せるかは知りませんが）。
中でも最も有名だと言っても過言ではないのが、
$$ \log \frac{y}{1-y} = \exp(w\cdot \phi(x)) $$
のような形式です。本当に何でもありという感じですね。 しかし、これは後に見るようにちゃんと意味のあるモデル化が行われています。
ノイズについて考えなおす
更にノイズについて考えなおしてみましょう。最小二乗法では
$$ \begin{align} y &= \exp(w\cdot \phi(x)) + \epsilon \\ \epsilon &\sim {\cal N}(0, \sigma^2) \end{align} $$
というような状況を想定していました。 ノイズはガウス分布に従っていると思っているわけですから、 プロットされた点たちの真ん中くらいを通るような線を引く最小二乗法は確かに妥当そうですし、 実際、最尤推定の解は最小二乗法と全く同じになります。
ところで、ノイズなど諸々を含めて、$x$を受け取った際の出力$y$は以下の確率分布に従います。
$$ y \sim {\cal N}(y \mid w\cdot \phi(x)), \sigma^2) $$
平均的には$w\cdot \phi(x))$でありながら、ノイズが乗って観測されるというような雰囲気です。 ノイズという確率的な変動を考えたがために、その影響を受けている$y$の方も確率的に観測されるということです。 当然ノイズにガウス分布以外を想定した場合、$y$もそれに応じた確率分布に従うと想定するわけであり、 どのようなモデルを考えるのかの重要な要素となります。
それまでの線形モデルと言えば、最小二乗法を使うのが主流であり、それはガウス分布を想定していることにほかなりませんでした。 一般化線形モデルでは柔軟に、データに合わせていろいろな分布を想定していくことになります。
ロジスティック回帰モデル
さて、一般化線形モデルというのは、左辺側をゴニョゴニョしつつ、確率分布についても真剣に考えなおしてみた。 というような、線形モデルの発展版です。
線形モデルをいじったロジスティック回帰モデルの考え方
今やニューラルネットワークの文脈で、中間層を持たない２クラス分類の単純パーセプトロンとして紹介されることの多いロジスティック回帰は 一般化線形モデルの例にもなっています。それはどのようなものかというと、
$$ y = w\cdot \phi(x) $$
で、$x$に対する加工$\phi(x)$は行わずに（いや、実際は行っても構わない）、
$$ y = w\cdot x $$
からスタートして、左辺側に細工をして以下のようにしてしまいましょう。
$$ \log \frac{y}{1-y} = w\cdot x $$
先ほど現れた有名な加工です。何じゃこれは！！？と思ってしまいますね。
ここで、"$\log y = w\cdot \phi(x)$"を"$y = \exp(w\cdot \phi(x))$"と変えてみると観やすかったように、 今回も$y$について整理してみることで、これの意味を把握してみましょう。
$$ \begin{align} &&\log \frac{y}{1-y} &= w\cdot x \\ \Leftrightarrow &&\frac{y}{1-y} &= \exp(w\cdot x) \\ \Leftrightarrow &&\frac{1-y}{y} &= \frac{1}{\exp(w\cdot x)} \\ \Leftrightarrow &&\frac{1}{y} - 1 &= \frac{1}{\exp(w\cdot x)} \\ \Leftrightarrow &&\frac{1}{y} &= \frac{1}{\exp(w\cdot x)} + 1 \\ \Leftrightarrow &&\frac{1}{y} &= \frac{1 + \exp(w\cdot x)}{\exp(w\cdot x)} \\ \Leftrightarrow && y &= \frac{\exp(w\cdot x)}{1 + \exp(w\cdot x)} \\ \Leftrightarrow && y &= \frac{1}{1 + \exp(-w\cdot x)} \\ \Leftrightarrow && y &= {\rm sigmoid}(w\cdot x) \end{align} $$
上記のような変形の末、わけの分からん加工を左辺に施したがために 実は、$y$を$w \cdot x$をシグモイド関数で変換した出力を得ていることになったのです。 「やっぱり単純パーセプトロンの2クラス分ついと同じだ！！」と思うわけですが、実際には話が逆で、 単純パーセプトロンがロジスティック回帰と同じだということです。
現実の話の流れでは、そもそも$y$に謎の加工を施した時点でしっかりと意図があります。 まず予測したいものが$1$か$0$（あるいは男・女でもいい）ならば、それを表すモデルはベルヌーイ分布が適していそうです。
要するに確率$p$で$1$が生起するようなベルヌーイ分布を${\rm Bern}(p)$と表すとして、 私達が作りたい予測モデルの出力は
$$ {\rm output} \sim {\rm Bern}(p) $$
という形式になっていて欲しいと考えます。さて、問題は、$p$がどんな値になっているかということです。 この$p$という値はどうも、手持ちのデータ$x$から見積もられるはずです。 それは例えば
$$ p = w \cdot x $$
という形になっているかもしれません。 しかし、これは記事の最初の方で見たとおり、確率なのにも関わらず負の数を取る可能性が出てきます。 それではまずいです。そこで、最初に見たテクニックで、
$$ \log p = w \cdot x $$
と被せてやることで、$p$自体が負の値を取ることは一切無くなりました。しかし、まだ問題があります。 なぜかというと、"$p = \exp(w\cdot x)$"のグラフを想像してください。$\exp$は増加が世界一早くて有名です。 $w\cdot x$が少し大きくなるだけで、$p$というのはまたたく間に無限の彼方へ飛んでいくでしょう（確率$p$を考えているのに１を超えてしまう）。
$\log$を被せれば負に行くことは避けられます。 あとは無限の彼方へ飛んでいく状況を考慮しなければなりません。そこで、左辺に新たなる加工を施しましょう。 今、$\log$の中身が無限の彼方へ飛んでいってもいいような$p$の形が必要です。 それが
$$ \frac{p}{1-p} $$
です。$p$はベルヌーイ分布から$1$が生起する確率ですから、$1-p$は$0$が生起する確率です。 つまりこれは、「$0$が生起する確率と$1$が生起する確率の比」です。オッズ比と呼ばれています。
これはいくらでも大きな値を取ることがありえます。しかし負にはなりません。 条件は整いました。こいつを$\log$で包んでやれば、今まで可笑しいと感じていたことが全て避けられます。
$$ \log \frac{p}{1-p} = w \cdot x $$
とすることに決めたのです。
私達は確率($0 \sim 1$)をオッズ比($0\sim \infty$)にし、更に$\log$を被せることで全区間で値を取るように加工し、 その値を$w\cdot x$で回帰することに決めたのです（応用上2クラス分類に使われるが、本質的に回帰です）。
これが冒頭で見た、
$$ \log \frac{y}{1-y} = w\cdot x $$
の意味になります。そうです。ここで$y$というのはベルヌーイ分布のパラメータであり、 これがロジスティックシグモイド関数によって$0 \sim 1$に抑えられることになったのは必然だったのです （ロジスティックシグモイド関数によって抑えられたというより、$0 \sim 1$に抑えるような上手い加工をした結果得られた関数をロジスティックシグモイド関数と呼んでいるのでしょう）。
最後に
さて、ロジスティック回帰モデルは対数オッズ比を線形回帰するという形で書けるようです。
$$ \log \frac{y}{1-y} = w\cdot x $$
そして、ここでの$y$をベルヌーイのパラメータとした
$$ {\rm output} \sim {\rm Bern}({\rm output} \mid y) = {\rm Bern}(\rm{output} \mid {\rm sigmoid}(w\cdot x)) $$
によって予測&{\rm output}&を得ようということです。 こんな芸当はガウス分布にこだわっているうちは絶対に出てこないでしょう（と言っても、最小二乗法とタメを張れるくらいド定番ですが）。
要するに確率分布について、あるいは左辺の加工（リンク関数とか言われるらしいです）についてもっとよくしれば、いろいろなモデルを考えることができるはずです。
あとは、これを手持ちのデータが$(x, {\rm output})$の形で沢山得られているのであれば、
$$ \prod_i {\rm Bern}({\rm output} \mid {\rm sigmoid}(w\cdot x)) $$
で最尤推定ができるという話になります。 ここで是非、この最尤推定の式については、ちゃんとベルヌーイ分布の式を入れて対数尤度を取って計算してみてください。
そこでよくご存知の交差エントロピーが現れてくるはずです。
確率モデルを考えることができるようになれば、あとは事前分布を置いて事後分布を計算し、 ベイズ予測分布を利用するなどすれば、更に発展的な手法を得ることができると思われます。
最後の最後に、この記事で紹介した流れは私の妄想なので、歴史的にそうであるかは知りません。 本当はパーセプトロンの方が先かもしれないしね。
