ざっくりML

機械学習に興味ある大学院生によるブログです.

交差検定

今回はモデル選択の時に使われる手法の、交差検定について軽くまとめて実装して見たいと思う.

交差検定

訓練とテストに使えるデータには限りがあり、良いモデルを選択するために得られたデータをできるだけたくさん訓練に使いたい. しかし、確認用に使うテストデータが小さいとたまたま良い精度が出てしまったのかもしれないなど、うまく評価ができない. そこでよく用いられる手法に交差検定 (cross-validation) がある. 得られたデータのうち\frac{S-1}{S}の割合のデータを訓練に使い、残りをテストデータに使う. 例えば全データ数が[N = 100]とし、S = 4とする. この時訓練に使うデータは全体の\frac{4-1}{4}、つまり7.5割(75個)を使い残りの2.5割(25個)をテストデータに使う. この場合全データをS分割したことになる. したがって、テストに使える訓練データはS個あるので、学習と評価をS回繰り返す. S回分の精度の平均を学習器の精度として採用することで、比較を行うことができる. 1回目の学習では最初の25個をテストデータ、残りを訓練データに、2回目は次のブロックのテストデータを使う. これを繰り返す.
f:id:linearml:20170923033620p:plain 欠点としては、分割数と訓練数が比例することである. 一回の訓練に時間のかかる学習器に対して、交差検定をしようとすると、S回訓練を行うことになるので、膨大な時間がかかることは想像つく.
以下に私が書いたソースコード載せる. 分類器はサポートベクタマシン、データセットはirisを使う.

交差検定

16,17,25行目で与えられたデータをシャッフルするためのマスクを定義している. 28,29行目でs回目のテストデータ、30,31行目では学習に使う訓練データを作成している. setdiff1d(a,b)はaとbの差集合を求めることができるので、それを利用して訓練データとテストデータを分割することにした. 36,37,38行目で精度(この指標についてはまた今度)を計算して保存しておく. 最後にその平均を出力することにしている

結果は以下のように得られた.
micro precision : 0.96
micro recall : 0.96
micro F1 : 0.96

多項式曲線フィッティング

ビショップ本では多項式曲線フィッティングから始まっている. この単純な回帰から多くの重要な概念を説明したいらしい. まず、実数値の入力変数xから、実数値の目的変数tを予測することを考える. 今回は人工的に生成したデータを使って例で考え、後で実装を試みる. ここで訓練データとしてN個のxを並べた、\boldsymbol{x} = (x_1,\cdots,x_N)^{T}とそれぞれに対する観測値である\boldsymbol{t} = (t_1,\cdots,t_N)^{T}が与えられたとする. 例えば、N=10個のデータからなる訓練集合をプロットしたものを以下に載せる. f:id:linearml:20170919065043p:plain ただし、ここでは関数sin(2\pi x)にランダムなノイズを加えて生成したものをプロットしている. この訓練データを用いて、新たな入力変数\hat{x}に対する観測値\hat{t}を予測することが目標である. そのためのアプローチの一つに曲線フィッティングがある. ビショップ本ではあまり形式ばらない形で話を進めている. ここでは以下のような多項式を使ってフィッティングすることを考える.

\displaystyle y(x,\boldsymbol{w}) = \sum_{i=0}^{M}w_{i}x^{i} 

ただし、M多項式の次数である. \boldsymbol{w}に関して線形であることから、線形モデルと呼ばれている. (3章、4章で詳しく議論しているのでそこらへんはまたの機会に) フィッティングでは訓練データを多項式に当てはめて係数の値、つまり\boldsymbol{w}を求める. これは、\boldsymbol{w}をランダムに固定したときの多項式の値と、訓練データの目的変数の値との差を最小化することで達成することができる. 全ての訓練データの入力変数と目的変数との差が小さくなるように作られた曲線が求めるべき曲線ということになるので直感的にも理解しやすい. 誤差関数にも色々あるが、今回は単純で広く用いられている、二乗和誤差で考える.

\displaystyle E(\boldsymbol{w}) = \frac{1}{2}\sum_{n=1}^{N}\{y(x_{n},\boldsymbol{w}) - t_{n}\}^{2}

これを最小にする\boldsymbol{w}を求めるのが目標なので、\boldsymbol{w}偏微分して0と置き、方程式を解くことになる. その微分動作の時に降りてくる指数部の2と相殺させるために係数を1/2にしている. 誤差関数を\boldsymbol{w}偏微分すると次のようになる(演習問題となっている).

\displaystyle \sum_{j=0}^{M}A_{ij}w_{j} = T_{i}
\displaystyle A_{ij} = \sum_{n=1}^{N}(x_{n})^{i+j}
\displaystyle T_{i} = \sum_{n=1}^{N}(x_{n})^{i}t_{n}

この式を\boldsymbol{w}に関して解けば良い. 求まった係数を\boldsymbol{w}^{*}とすると、結果としてy(x,\boldsymbol{w}^{*})が求まった多項式として表される. あとは、次数Mを幾つに選択すれば良いかという問題が残っているが、この問題をモデル選択と呼ぶ. 例として、M = 0,1,3,9に対してフィッテイングした結果と実装したソースを以下に示す.

f:id:linearml:20170919063200p:plainf:id:linearml:20170919063204p:plainf:id:linearml:20170919063208p:plainf:id:linearml:20170919063218p:plain
M = 0,1,3,9に対してのフィッテイング結果

多項式曲線フィッティング

M=3のときが最もよく再現できているように見える. 次数をM=9にした場合、全ての点を通っているのでE(\boldsymbol{w}^{*})=0になっているが曲線は元の形からは程遠い. このような振る舞いを過学習と呼ばれている. 訓練データに対してはうまく表現できていても、新しいデータに対しては全く使い物にならない、汎用性がない. さて、この汎化性能とMの関係を定量的に評価したいときの指標として平均二乗平方根誤差(Root Mean Square error; RMS)というものがある.

\displaystyle E_{RMS} = \sqrt{2E(\boldsymbol{w}^{*})/N}

平均をとることで、サイズが異なるデータに対して比較ができる. 訓練データとテストデータのRMSをグラフ化したものを以下に載せる. f:id:linearml:20170919063356p:plain Mが小さい時は、誤差がかなり大きく、対応する多項式は柔軟性に欠け、元の関数をうまく表現できていない. Mが3から7の間では誤差が小さくうまく表現できているように見える. Mが8を超えると訓練データとの誤差は0に限りなく近くなるがテストデータとの誤差は極端に増える. よってMが大きければ良いわけではない. 今回の場合はM=3で十分であると言える. テイラー展開などの級数展開を思い浮かべると次数を増やすほど良い結果が得られると期待してしまうから困ってしまう. また、訓練データ数を増やしたときの結果を以下に示す.

f:id:linearml:20170919063419p:plainf:id:linearml:20170919063426p:plain

訓練データ数を増やしたとき
図の通り、訓練データ数が多いほど過学習が起きにくい. また、過学習を避けるためにベイズ的アプローチを採用すれば良いがこの話は3章の内容なのでまた今度. この章では、過学習を抑制するためのテクニックとして正則化を紹介している. 誤差関数に対し罰金項を付加することで、係数が大きな値になることを防ごうとしている. 罰金項のうち最も単純なものは係数の2乗和である.

\displaystyle \hat{E}(\boldsymbol{w}) = \frac{1}{2} \sum_{n=1}^{N}\{y(x_{n},\boldsymbol{w})-t_{n}\}^{2} + \frac{\lambda}{2} \boldsymbol{w}^{T}\boldsymbol{w}

で表される. 係数\lambdaは罰金項と二乗誤差の和の項との相対的な重要度を調節している. この場合も先ほどと同様に\boldsymbol{w}偏微分して0と置いて方程式を解けば\boldsymbol{w}^{*}が閉じた形で求まる. 統計学の分野では縮小推定、ニューラルネットワークでは荷重減衰と呼ばれている. \lambdaは交差検定などを行なって決定すると良い.

テイラー展開

最適化数学でテイラー展開がたくさん出てきたのでとりあえず実装してみることにしました. 今回はサインカーブを使おうと思います. まず、関数のある一点での導関数の無限和で表されたものをテイラー級数と言います. その、テイラー級数を得るためにテイラー展開を行います. この時、ある一点をx=aとすると、x=a周りのテイラー展開と表現します. 要は、x=a近傍における導関数の無限和をとることを意味します. ちなみに、原点回りのテイラー展開のことをマクローリン展開といいます(ややこしい...). ここで、このテイラー展開マクローリン展開について気軽に学べる良書として数学ガールを紹介しておきます. 私は、学部時代はテイラー展開の有難みが全く分かりませんでしたが、この数学ガールを読んでテイラー展開の使いどころを理解できるようになりました.

数学ガール (数学ガールシリーズ 1)

数学ガール (数学ガールシリーズ 1)

さて、関数f(x)x=a周りのテイラー展開は次の式で表されます.

\displaystyle f(x) = \sum_{k=0}^{\infty}\frac{f^{(k)}(a)}{k!}(x-a)^{k}

無限級数で近似していることになりますので、例えばk=5で止めれば関数を5次関数で近似できることになります. kが大きいほど元の関数に近づきます.

また、マクローリン展開は原点回りのテイラー展開のことだったので、上式のa0を代入すればよいので、

\displaystyle f(x) = \sum_{k=0}^{\infty}\frac{f^{(k)}(0)}{k!}x^{k}

と表されます.

ここで、サインカーブに対してマクローリン展開してみたいと思います. \sin x微分していくと\sin x \rightarrow \cos x \rightarrow -\sin x \rightarrow -cosxとなります. また、\sin 0 = 0なのでサインの項が消えます. よってサインカーブマクローリン展開は、

\displaystyle \sin x = \sin 0 + \frac{cos 0}{1!}x + \frac{-sin 0}{2!}x^{2} + \frac{-cos 0}{3!}x^{3} + \frac{sin 0}{4!}x^{4} + \cdots = x - \frac{x^{3}}{3!}+\cdots

と規則性が見えてきます. コサインカーブマクローリン展開も同様に導くことができます.

今回は、サインカーブに限定して実装したので、導関数のところは事前にこちらで用意することにしました(4通りしかないので). 一般の関数に対しての実装はまた別の機会にあげられたらと思います. また、sympyというのを使えば、テイラー展開もn次導関数も求めることができますが、今回はnumpyだけを用いています.

テイラー展開展開(sin)

f:id:linearml:20170915200418p:plain

かなり見づらいですが、元のサインカーブが青色で他は各次元で近似したものを図示しています. 例えば緑色は1次元で近似しているので、直線になっています. その後次元数を上げると元のサインカーブに近づいていき、11次元の時は大分サインカーブに近づいていることが分かります. これなら分かる最適化数学ではn変数のテイラー展開について触れていましたが、それもまた別の機会で...

勾配法

「これなら分かる最適化数学」を読み始めたので、最適化についても紹介していこうと思います.

 

これなら分かる最適化数学―基礎原理から計算手法まで

これなら分かる最適化数学―基礎原理から計算手法まで

関数の最適化の代表の手法の「勾配法」を軽く説明する.

  1. 勾配法

1変数の場合と多変数の場合があるのでまずは1変数の場合を説明する. そもそも最適化とは、ある制約条件のもとで関数の値を最大(最小)にする変数の値を求めることを言う. 勾配法のアルゴリズムは単純であるが昔からよく使われる手法の一つである. まず初期値 x_0を与える. この初期値の与え方についても多くの研究がされているが今回は割愛. f^{\prime}(x_0) = 0ならそこで最大(最小)値をとる. ある点での勾配が正、つまりf^{\prime}(x_t ) > 0のとき(下図の青)はその点から負の方向に、勾配が負のとき、つまりf^{\prime}(x_0)  < 0のとき(下図の赤)はその点から正の方向に少し進むことで、目的地(最大、最小となる点)に近づいていく. 簡単にまとめたアルゴリズムを以下に示す. f:id:linearml:20170903230043p:plain

  1. xの初期値x^{0}を与える
  2. もしf^{\prime}(x^{0}) = 0のとき、そのときのx_tが最大(最小)値をとる地点である
  3.  x^{t} = x^{t-1} - \alpha f^{\prime}(x^{t-1})に更新する (ただし、\alphaは学習率と呼ばれ1ステップでどれだけ更新するかを表すハイパーパラメータである)
  4. 収束するまで2,3を繰り返す

計算機で実装する上でちょうどf^{\prime}(x)=0となるxを求めるのは現実的ではないので、ある小さい値\epsilonに対して、|f^{\prime}(x)|< \epsilonになる点に到達した時に収束したとみなすことにする. また導関数 \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}を用いて数値微分する. 以下にソースと結果を載せる.

勾配法

f:id:linearml:20170903232812p:plain 初期値はx=5で53回の繰り返しで収束している. またオレンジの線が移動した経緯を示している.

多変数の場合も基本的な流れは同じで各点での導関数を求め少しずつ動かしていけば良い.

勾配法は単純で有効な手法であるが問題点もいくつかある. 一つ目は初期値の与え方によって局所解に陥ってしまうことである. 二つ目は勾配を計算しなければならないが、微分できない関数であったり、できたとしても複雑すぎる場合がある. また変数が数百、数千のときも勾配を求めるのが大変なことである. なので色々な最適化の手法を知ることが必要であり、それぞれの関数に対して適切な手法を選ぶことが求められる.

ガウス分布と条件付きガウス分布

今回は、ビショップ本の2章のメインとも言えるガウス分布についてです. ガウス分布については結構な量があるのでいくつかに分けてまとめていきたいと思います.

  1. ガウス分布について ガウス分布 (Gaussian distribution)は、正規分布 (normal distribution)とも呼ばれていて、連続変数の分布をモデル化する際によく使われます. 1変数のときのガウス分布は、
{\displaystyle  N(x|\mu,\sigma^{2}) = \frac{1}{(2\pi\sigma^{2})^{\frac{1}{2}} } \exp(-\frac{1}{2\sigma^{2}}(x-\mu)^{2} )}

と表されます. ここで \muは平均、\sigma^{2}は分散です. D次元のときは、

{\displaystyle  N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{\frac{D}{2}}}\frac{1}{|\boldsymbol{\Sigma}|^{\frac{1}{2}}}  \exp\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\} }

と表され、\boldsymbol{\mu}はD次元の平均ベクトル、\boldsymbol{\Sigma}D \times Dの共分散行列です.

1次と2次の時のガウス分布を実装してみたので以下にソースコードと結果をあげておきます.

ガウス分布 (1次の場合)

f:id:linearml:20170828234310p:plain

ガウス分布 (2次の場合)

f:id:linearml:20170828234316p:plain

ビショップ本では、ガウス分布の幾何的な形状について考え、平均と共分散行列について解釈することから始めています. このことについてはまた別の機会にまとめられたらと思っています…

さて、多変量ガウス分布の重要な特性の一つに、2つの変数の同時分布がガウス分布に従うならば、それらの変数を用いた条件付き分布及び、周辺分布もガウス分布になるというのがあります. 今回は、その内の条件付き分布についてのお話です.

\boldsymbol{x}N(\boldsymbol{x}|\boldsymbol{\mu},\boldsymbol{\Sigma})に従うD次元ベクトルとし、これを互いに素な部分集合\boldsymbol{x_a}\boldsymbol{x_b}に分割します. つまり、

\displaystyle \boldsymbol{x} = \begin{bmatrix}\boldsymbol{x_{a}} \\ \boldsymbol{x_{b}} \end{bmatrix}

とし、これに対する平均ベクトルと共分散行列もそれぞれ、

\displaystyle \boldsymbol{\mu} = \begin{bmatrix}\boldsymbol{\mu_a} \\ \boldsymbol{\mu_b}\end{bmatrix}
\displaystyle \boldsymbol{\Sigma} = \begin{bmatrix}\boldsymbol{\Sigma_{aa}} && \boldsymbol{\Sigma_{ab}} \\ \boldsymbol{\Sigma_{ba}} && \boldsymbol{\Sigma_{bb}} \end{bmatrix}

で与えれます. 共分散行列が与えられたが、ビショップ本をはじめとし、共分散の逆行列を考えた方が便利になることが多いため、\boldsymbol{\Lambda} = \boldsymbol{\Sigma}^{-1}を考えています. これを精度行列と呼んでいます. 以降では精度行列をベースに話が進んでいきます.

\displaystyle \boldsymbol{\Lambda} = \begin{bmatrix} \boldsymbol{\Lambda_{aa}}&&\boldsymbol{\Lambda_{ab}} \\ \boldsymbol{\Lambda_{ba}} && \boldsymbol{\Lambda_{bb}}\end{bmatrix}

となるが、 \boldsymbol{\Lambda_{aa}}\boldsymbol{\Sigma_{aa}}逆行列となっていないことに注意しなくてはいけません. (このことについては最後に述べます.)

まず、一般のガウス分布の指数部は、 \displaystyle -\frac{1}{2}(\boldsymbol{x} - \boldsymbol{\mu})^{T} \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) = - \frac{1}{2} \boldsymbol{x}^{T} \boldsymbol{\Sigma}^{-1} \boldsymbol{x} +  \boldsymbol{x}^{T} \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu} + ( \boldsymbol{x}とは独立な項) と変形できます. 右辺の\boldsymbol{x}の2次の項の係数行列が精度行列(共分散の逆行列)と等しく、\boldsymbol{x}の線形項の係数は\boldsymbol{\Sigma^{-1}\mu}と等しくなります. ここから平均ベクトルを計算することができます. このことを利用して、条件付き分布の平均と共分散行列を求めますが、基本的な考えとしては分割した変数が従うガウス分布を、\boldsymbol{x_b}を固定した上での\boldsymbol{x_a}の関数とみなします. さて、分割した変数が従うガウス分布の指数部は、

\displaystyle -\frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu} )^{T} \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) = \frac{1}{2}(\boldsymbol{x_a} - \boldsymbol{\mu_a})^{T}\boldsymbol{\Lambda_{aa}}(\boldsymbol{x_a} - \boldsymbol{\mu_a}) - \frac{1}{2}(\boldsymbol{x_a} - \boldsymbol{\mu_a})^{T}\boldsymbol{\Lambda_{ab}}(\boldsymbol{x_b} -  \boldsymbol{\mu_b}) -\\  \frac{1}{2}(\boldsymbol{x_b} - \boldsymbol{\mu_b})^{T}\boldsymbol{\Lambda_{ba}}(\boldsymbol{x_a} - \boldsymbol{\mu_a}) -
\frac{1}{2}(\boldsymbol{x_b} - \boldsymbol{\mu_b})^{T}\boldsymbol{\Lambda_{bb}}(\boldsymbol{x_b} - \boldsymbol{\mu_b})

となることはすぐ確認できると思います. この中の\boldsymbol{x_a}の2次の項を全て取り出すと、

\displaystyle -\frac{1}{2}\boldsymbol{x_a}^{T}\boldsymbol{\Lambda_{aa}}\boldsymbol{x_a}

となるので、\boldsymbol{\Sigma_{a|b}} = \boldsymbol{\Lambda_{aa}}^{-1}となります. 次に、\boldsymbol{x_a}について線形の項を全て取り出すと、 \boldsymbol{x_a}^{T}
{\boldsymbol{\Lambda_{aa}}\boldsymbol{\mu_a} - \boldsymbol{\Lambda_{ab}}(\boldsymbol{x_b} - \boldsymbol{\mu_b}) }が得られます. この式の\boldsymbol{x_a}の係数が  \boldsymbol{\Sigma_{a|b}}^{-1}\boldsymbol{\mu_{a|b}}になることは説明しました. これらより、\boldsymbol{\mu_{a|b} = \boldsymbol{\mu_a} - \boldsymbol{\Lambda_{aa}}^{-1} \boldsymbol{\Lambda_{ab}}(\boldsymbol{x_a} - \boldsymbol{\mu_b})}となります. 今までの結果は精度行列を用いての表現方法でしたがもちろん共分散行列での表現も可能です. 分割された行列の逆行列に関する計算はビショップ本では公式のように取り上げられています. (演習問題にも取り上げられており、東京大学大学院試でも出題されていました笑)

\displaystyle \begin{bmatrix}A && B \\ C && D\end{bmatrix}^{-1} = \begin{bmatrix} M && -MBD^{-1} \\ -D^{-1}CM && D^{-1} + D^{-1}CMBD^{-1}\end{bmatrix}

となります. ここでM = (A - BD^{-1}C)^{-1}となります. このMの逆行列には部分行列Dに関するシューア補行列という名前がついているようです. この公式と、

\displaystyle \begin{bmatrix}\Sigma_{aa} && \Sigma_{ab} \\ \Sigma_{ba} &&\Sigma_{bb}\end{bmatrix}^{-1} = \begin{bmatrix} \boldsymbol{\Lambda_{aa}} && \boldsymbol{\Lambda_{ab}} \\ \boldsymbol{\Lambda_{ba}} && \boldsymbol{\Lambda_{bb}}\end{bmatrix}

より、

\displaystyle \boldsymbol{\Lambda_{aa}} = (\boldsymbol{\Sigma_{aa}} - \boldsymbol{\Sigma_{ab}}\boldsymbol{\Sigma_{bb}}^{-1}\boldsymbol{\Sigma_{ba}})^{-1}
\displaystyle \boldsymbol{\Lambda_{ab}} = -(\boldsymbol{\Sigma_aa} - \boldsymbol{\Sigma_{ab}}\boldsymbol{\Sigma_{bb}}^{-1}\boldsymbol{\Sigma_{ba}})^{-1}\boldsymbol{\Sigma_{ab}}\boldsymbol{\Sigma_{bb}}^{-1}

が得られます. これらから、

\displaystyle \boldsymbol{\mu_{a|b}} = \boldsymbol{\mu_a} +  \boldsymbol{\Sigma_{ab} \boldsymbol{\Sigma_{bb}}^{-1}( \boldsymbol{x_b} -\boldsymbol{ \mu_b}})
\displaystyle \boldsymbol{\Sigma_{a|b}} = \boldsymbol{\Sigma_{aa}} +  \boldsymbol{\Sigma_{ab}} \boldsymbol{\Sigma_{bb}}^{-1}\boldsymbol{\Sigma_{ba}}

となります.

いやあ、まとめるの非常に下手くそですね。。。 このブログを通してまとめる技術もついたらいいなと思ってます笑

ベータ分布

ベルヌーイ分布と二項分布のパラメータ\mu最尤推定で求めた結果は\mu_{ML} = \frac{m}{N}で、データ集合中でx = 1となる観測値の割合になる. しかし、例えば3回の試行中3回とも表が出た場合に\mu_{ML} = 1となってしまい今後もずっと表が出るように予測することになる.これは過学習してしまっている. そこで、この問題をベイズ主義的に扱うことを考える. ベイズの定理は、

{\displaystyle p(x|y) = \frac{p(y|x)p(x)}{p(y)}}

であり、p(x)を事前分布、p(x|y)を事後分布と呼んだ. したがって、ベイズ主義的立場でパラメータを推定するときには、パラメータの事前分布p(\mu)を考え、事後分布p(\mu|x)を最大にしたり(MAP推定)利用したりすることをする. 今回はベルヌーイ分布や二項分布のパラメータの事前分布についてまとめる. まず尤度は\mu^{x}(1-\mu)^{1-x}の形をしている. したがって、事前分布に \mu 1-\muのべき乗に比例するように選べば、事後分布は事前分布と同じ形になる. なぜなら事後分布は尤度と事前分布に比例するからである. この性質は共役性と呼ばれている. なので、事前分布には事後分布と共役なものを選ぶと良い. このことから、ベルヌーイ分布や二項分布のパラメータの事前分布には次のベータ分布を選ぶ.

{\displaystyle Beta(\mu|a,b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1}}

ただし、\Gamma(x)はガンマ関数と呼ばれ、\Gamma(x) \equiv \int_{0}^{\infty}u^{x-1}e^{-u}duで定義されている. aとbはパラメータの分布p(\mu)を決めるためのパラメータで超パラメータ(hyperparameter)と呼ばれる. ここで、ベータ分布を私が実装したコードを以下に示す.

gist163dce714aebb6093bc8e65586caf82b

また、超パラメータを色々な値に設定した時のベータ分布の図をいかに載せる. f:id:linearml:20170822212717p:plain

\muに依存する部分だけを取り出すと事後分布は、 p(\mu|m,l,a,b) \propto \mu^{m+a-1}(1-\mu)^{l+b-1}になるのは容易に確かめられる. ただし、l = N -mには注意. この式は事前分布の共役性によって事前分布と同じ形になっていることがわかる.

ベルヌーイ分布

第2章の確率分布の一番最初に紹介されているベルヌーイ分布について

二値確率変数 (binary random variable, 以下二値r.v.) {x \in \{0,1\}}が1つの場合について考える. 例えば、 x = 1で「表」、x = 0で「裏」を表現するとする. これを使えばコイン投げの結果が二値r.v.で表現できる. 一般の話を扱うので「表」が出る確率は0.5ではなくパラメータ\muで表すことにする. つまり、

{\displaystyle  p(x = 1 | \mu) = \mu }

となる. 確率なので  0 \leq \mu \leq 1 を満たし、 p(x = 0 | \mu) = 1 - \muとなる. これより xが従う分布は

{\displaystyle Bern(x | \mu) = \mu ^{x} (1 - \mu)^{1-x}}

となり、これがベルヌーイ分布である. 平均、分散はそれぞれ、E[x] = \mu, V[x] = \mu (1-\mu)である.

ベルヌーイ分布をpythonとmatplotlibを使って可視化してみようと思う. 私が書いたソースは以下の通りである.

パラメータをそれぞれ0.1,0.5,0.8に設定したものを載せておく.

f:id:linearml:20170818103828p:plain f:id:linearml:20170818103833p:plain f:id:linearml:20170818103839p:plain

次に、観測データ集合 \{x_{1},\cdots,x_{N}\}がベルヌーイ分布から独立に与えられたとする. その時の尤度は

{\displaystyle \prod_{n=1}^{N}\mu^{x_{n}}(1-\mu)^{1-x_{n}}}

となり、この尤度を最大にするようにしてパラメータを推定するのが最尤推定である. 積の形では扱いにくいので尤度の対数をとって和の形にしてしまう. この対数尤度をパラメータ \mu微分して0とおき、方程式を解けば最尤推定量が求まり以下のように得られる.

{\displaystyle \mu_{ML} = \frac{1}{N} \sum_{n=1}^{N} x_{n}}

ここで表が出た回数、つまりx = 1となった回数をmとすると

{\displaystyle \mu_{ML} = \frac{m}{N}}

となる. 例えば試行3回中3回とも表が出た場合に、最尤推定量は1となるので、今後もずっと表が出続けると予測してしまっている. これは過学習の極端な例であるが、これではまずい. 今後はパラメータが従う事前分布を導入して、もっと常識的な結果を得る方法について考える.