Created
February 25, 2019 00:15
-
-
Save shotahorii/e34f065c43a2a5da81a736fab0af7e3a to your computer and use it in GitHub Desktop.
MAP_curve_fitting
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Curve Fitting with MAP estimation\n", | |
| "Related gists:\n", | |
| "- [MAP Estimation](https://gist.github.com/shotahorii/386f75acb7efac0f6518c096293e08b9)\n", | |
| "- [MLE and least squares](https://gist.github.com/shotahorii/be3bda5e42eb11fe7b721254768408ab)\n", | |
| " \n", | |
| "References:\n", | |
| "- https://www.iwanttobeacat.com/entry/2018/02/14/231841\n", | |
| "- http://aidiary.hatenablog.com/entry/20100404/1270359720\n", | |
| "- http://www.geocities.co.jp/kashi_pong/FSM2017_6.pdf (pdfダウンロードリンク)\n", | |
| " \n", | |
| "MAP推定で曲線フィッティングを行ってみる。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "#### その前に...最尤推定おさらい" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 68, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "最尤推定とは、最も尤度の大きいパラメータを見つけ出し、それを真のパラメータであるとする推定法である。\n", | |
| "尤度とは、ある観察データが得られた時、その背後のパラメータθの取るある値がどれほど尤もらしいか、\n", | |
| "という度合いのことである。ある事象yが観測された時、その背後にあるパラメータθの尤度は$$P(y|\\theta))$$と表される。\n", | |
| "複数(n個)の観測データが与えられた時、θの尤度は以下のように表される。\n", | |
| "これを最大化するようなθを求めるのが最尤推定である。<br><br>\n", | |
| "$$\\prod_{i=0}^n P(y_i|\\theta)$$<br>\n", | |
| "\n", | |
| "Related gistsの\"MLE and least squares\"で書いたように、最尤推定による曲線フィッティングでは、\n", | |
| "観測データは真の分布+等分散正規分布に従う誤差で表されると仮定した。つまり以下の式である。\n", | |
| "(ここで、以下の式では説明変数をxのみのn次式としているが、x1, x2, ... と説明変数の数は任意。次数も任意である。)\n", | |
| "<br><br>\n", | |
| "$$y_i = w_0+w_1 x+w_2 x^2+...+w_n x^n + \\epsilon _i(\\sim N(0,\\beta)) \n", | |
| "= f(X_i,W) + \\epsilon _i(\\sim N(0,\\beta))$$<br>\n", | |
| "つまり、$$y_i = N(f(X_i,W),\\beta)$$<br><br>\n", | |
| "\n", | |
| "ここで、上で書いたようにパラメータθについて尤度は$$P(y_i|\\theta)$$と表され、ここで推定するパラメータはX,W,βであるので、\n", | |
| "尤度最大化問題は以下の式で表される。<br><br>\n", | |
| "$$argmax(\\prod_{i=0}^n P(y_i|X_i,W,\\beta))$$<br><br>\n", | |
| "\n", | |
| "これを、正規分布の確率密度関数で置き換えて以下の最大化問題とした。<br><br>\n", | |
| "$$argmax(\\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(X_i,W))^2}{2\\beta}))$$<br>" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "最尤推定とは、最も尤度の大きいパラメータを見つけ出し、それを真のパラメータであるとする推定法である。\n", | |
| "尤度とは、ある観察データが得られた時、その背後のパラメータθの取るある値がどれほど尤もらしいか、\n", | |
| "という度合いのことである。ある事象yが観測された時、その背後にあるパラメータθの尤度は$$P(y|\\theta))$$と表される。\n", | |
| "複数(n個)の観測データが与えられた時、θの尤度は以下のように表される。\n", | |
| "これを最大化するようなθを求めるのが最尤推定である。<br><br>\n", | |
| "$$\\prod_{i=0}^n P(y_i|\\theta)$$<br>\n", | |
| "\n", | |
| "Related gistsの\"MLE and least squares\"で書いたように、最尤推定による曲線フィッティングでは、\n", | |
| "観測データは真の分布+等分散正規分布に従う誤差で表されると仮定した。つまり以下の式である。\n", | |
| "(ここで、以下の式では説明変数をxのみのn次式としているが、x1, x2, ... と説明変数の数は任意。次数も任意である。)\n", | |
| "<br><br>\n", | |
| "$$y_i = w_0+w_1 x+w_2 x^2+...+w_n x^n + \\epsilon _i(\\sim N(0,\\beta)) \n", | |
| "= f(X_i,W) + \\epsilon _i(\\sim N(0,\\beta))$$<br>\n", | |
| "つまり、$$y_i = N(f(X_i,W),\\beta)$$<br><br>\n", | |
| "\n", | |
| "ここで、上で書いたようにパラメータθについて尤度は$$P(y_i|\\theta)$$と表され、ここで推定するパラメータはX,W,βであるので、\n", | |
| "尤度最大化問題は以下の式で表される。<br><br>\n", | |
| "$$argmax(\\prod_{i=0}^n P(y_i|X_i,W,\\beta))$$<br><br>\n", | |
| "\n", | |
| "これを、正規分布の確率密度関数で置き換えて以下の最大化問題とした。<br><br>\n", | |
| "$$argmax(\\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(X_i,W))^2}{2\\beta}))$$<br>" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "#### では、MAP推定による曲線フィッティングとは" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 31, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "観測データから考えて最も尤もらしいパラメータを選んだ最尤推定に対して、MAP推定は、もともとパラメータそのものに\n", | |
| "確率分布が存在している(事前確率)と想定し、観測データによってそれが更新されて事後確率が得られたと考え、この\n", | |
| "事後確率を最大化するパラメータを真のパラメータとする推定法である。<br>\n", | |
| "データDが与えられた時のパラメータθの事後確率は、以下のように表される。<br><br>\n", | |
| "$$P(\\theta|D)$$<br>\n", | |
| "ここで、ベイズの定理より、以下のように書ける。<br>\n", | |
| "$$P(\\theta|D) = \\frac{P(D|\\theta)P(\\theta)}{P(D)}$$<br><br>\n", | |
| "ここで、P(D)はパラメータθによらない値であるので、事後分布の最大化問題は以下のように書き換えられる。<br><br>\n", | |
| "$$argmax_\\theta P(\\theta|D) = argmax_\\theta P(D|\\theta)P(\\theta)$$" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "観測データから考えて最も尤もらしいパラメータを選んだ最尤推定に対して、MAP推定は、もともとパラメータそのものに\n", | |
| "確率分布が存在している(事前確率)と想定し、観測データによってそれが更新されて事後確率が得られたと考え、この\n", | |
| "事後確率を最大化するパラメータを真のパラメータとする推定法である。<br>\n", | |
| "データDが与えられた時のパラメータθの事後確率は、以下のように表される。<br><br>\n", | |
| "$$P(\\theta|D)$$<br>\n", | |
| "ここで、ベイズの定理より、以下のように書ける。<br>\n", | |
| "$$P(\\theta|D) = \\frac{P(D|\\theta)P(\\theta)}{P(D)}$$<br><br>\n", | |
| "ここで、P(D)はパラメータθによらない値であるので、事後分布の最大化問題は以下のように書き換えられる。<br><br>\n", | |
| "$$argmax_\\theta P(\\theta|D) = argmax_\\theta P(D|\\theta)P(\\theta)$$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 72, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "いよいよこれを曲線フィッティングにあてはめていく。最尤推定同様、観測データは真の分布+等分散正規分布に従うと仮定する。\n", | |
| "つまり、尤度P(D|θ)は以下のように表される。<br><br>\n", | |
| "$$P(y|X,W,\\beta) = \\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(X_i,W))^2}{2\\beta})$$\n", | |
| "<br><br>\n", | |
| "ここで、PRMLの表記に合わせて少し書き方をかえる。この上の式も含め、ハイパーパラメータβについて、分散=βとおいて\n", | |
| "式を書いていたが、ここではPRMLの表記に合わせて分散=1/βとして式を書き換える。\n", | |
| "<br><br>\n", | |
| "$$P(y|X,W,\\beta) = \\prod_{i=0}^n \\sqrt{\\frac{\\beta}{2\\pi}} exp(-\\frac{\\beta}{2}(y_i-f(X_i,W))^2)$$\n", | |
| "<br><br>\n", | |
| "で、これは最大化問題であり、この尤度は必ず正の値であることから、対数をとる。以下はPRMLのp29の式(1.62)と同じ。<br>\n", | |
| "$$logP(y|X,W,\\beta) = -\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 \n", | |
| "+ \\frac{n}{2}log\\beta - \\frac{n}{2}log(2\\pi)$$\n", | |
| "<br><br>\n", | |
| "で、次に事前分布P(θ)について。ここではパラメータwに事前分布があると考える。観測データが正規分布に従うので、\n", | |
| "この事前分布は正規分布の共役事前分布、つまり正規分布をとる。<br>\n", | |
| "本来であれば平均、分散共にパラメータであり、かつ全体を通して等分散でなければならないわけではないが、簡便のため、\n", | |
| "平均0、分散α(これがハイパーパラメータ)の等分散正規分布であると仮定する。<br>\n", | |
| "つまり、ハイパーパラメータαが与えられた時のwの事前分布は以下のようになる。\n", | |
| "<b>多項式(複数の説明変数$$x_1,x_2,...$$)を前提としているため、\n", | |
| "wは一つの数値ではなくベクトル($$(w_0, w_1,w_2,...)$$)であることに注意。</b>$$w_0$$は切片項。\n", | |
| "説明変数の数Mに対してベクトルwのエレメント数はM+1となる点に注意。\n", | |
| "<br><br>\n", | |
| "$$P(W|\\alpha) = N(W|0,\\alpha I) \n", | |
| "= \\prod_{m=0}^{M+1} \\frac{1}{\\sqrt{2\\pi\\alpha}} exp(-\\frac{w_m^2}{2\\alpha})$$<br>\n", | |
| "$$= (2\\pi\\alpha)^{-(M+1)/2} exp(-\\frac{1}{2\\alpha}W^T W)$$\n", | |
| "<br><br>\n", | |
| "と、ここまで書いたが、分散をαではなく1/αとおく方が一般的っぽいので、そのように書き直したのが以下の式である。\n", | |
| "これはPRMLのp30の式(1.65)と全く同じもの。<br><br>\n", | |
| "$$P(W|\\alpha) = N(W|0,\\alpha^{-1} I) = (\\frac{\\alpha}{2\\pi})^{(M+1)/2} exp(-\\frac{\\alpha}{2}W^T W) $$" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "いよいよこれを曲線フィッティングにあてはめていく。最尤推定同様、観測データは真の分布+等分散正規分布に従うと仮定する。\n", | |
| "つまり、尤度P(D|θ)は以下のように表される。<br><br>\n", | |
| "$$P(y|X,W,\\beta) = \\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(X_i,W))^2}{2\\beta})$$\n", | |
| "<br><br>\n", | |
| "ここで、PRMLの表記に合わせて少し書き方をかえる。この上の式も含め、ハイパーパラメータβについて、分散=βとおいて\n", | |
| "式を書いていたが、ここではPRMLの表記に合わせて分散=1/βとして式を書き換える。\n", | |
| "<br><br>\n", | |
| "$$P(y|X,W,\\beta) = \\prod_{i=0}^n \\sqrt{\\frac{\\beta}{2\\pi}} exp(-\\frac{\\beta}{2}(y_i-f(X_i,W))^2)$$\n", | |
| "<br><br>\n", | |
| "で、これは最大化問題であり、この尤度は必ず正の値であることから、対数をとる。以下はPRMLのp29の式(1.62)と同じ。<br>\n", | |
| "$$logP(y|X,W,\\beta) = -\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 \n", | |
| "+ \\frac{n}{2}log\\beta - \\frac{n}{2}log(2\\pi)$$\n", | |
| "<br><br>\n", | |
| "で、次に事前分布P(θ)について。ここではパラメータwに事前分布があると考える。観測データが正規分布に従うので、\n", | |
| "この事前分布は正規分布の共役事前分布、つまり正規分布をとる。<br>\n", | |
| "本来であれば平均、分散共にパラメータであり、かつ全体を通して等分散でなければならないわけではないが、簡便のため、\n", | |
| "平均0、分散α(これがハイパーパラメータ)の等分散正規分布であると仮定する。<br>\n", | |
| "つまり、ハイパーパラメータαが与えられた時のwの事前分布は以下のようになる。\n", | |
| "<b>多項式(複数の説明変数$$x_1,x_2,...$$)を前提としているため、\n", | |
| "wは一つの数値ではなくベクトル($$(w_0, w_1,w_2,...)$$)であることに注意。</b>$$w_0$$は切片項。\n", | |
| "説明変数の数Mに対してベクトルwのエレメント数はM+1となる点に注意。\n", | |
| "<br><br>\n", | |
| "$$P(W|\\alpha) = N(W|0,\\alpha I) \n", | |
| "= \\prod_{m=0}^{M+1} \\frac{1}{\\sqrt{2\\pi\\alpha}} exp(-\\frac{w_m^2}{2\\alpha})$$<br>\n", | |
| "$$= (2\\pi\\alpha)^{-(M+1)/2} exp(-\\frac{1}{2\\alpha}W^T W)$$\n", | |
| "<br><br>\n", | |
| "と、ここまで書いたが、分散をαではなく1/αとおく方が一般的っぽいので、そのように書き直したのが以下の式である。\n", | |
| "これはPRMLのp30の式(1.65)と全く同じもの。<br><br>\n", | |
| "$$P(W|\\alpha) = N(W|0,\\alpha^{-1} I) = (\\frac{\\alpha}{2\\pi})^{(M+1)/2} exp(-\\frac{\\alpha}{2}W^T W) $$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 91, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "ここまでで、尤度$$P(y|X,W,\\beta)$$と、事前分布$$P(W|\\alpha)$$を式で表すことができた。<br>\n", | |
| "また、ベイズの定理より、$$P(W|X,W,\\alpha,\\beta) \\propto P(y|X,W,\\beta)P(W|\\alpha)$$である。<br>\n", | |
| "さらに、これらは正の値であることから、両辺に対数をとり、最大化問題は以下のように書ける。<br><br>\n", | |
| "\n", | |
| "$$argmax_W logP(W|X,W,\\alpha,\\beta) = argmax_W logP(y|X,W,\\beta)P(W|\\alpha)$$<br><br>\n", | |
| "ここで左辺は以下のように変形される。<br><br>\n", | |
| "$$argmax_W -\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 + \\frac{n}{2}log\\beta - \\frac{n}{2}log(2\\pi)\n", | |
| "+ \\frac{M+1}{2}log(\\frac{\\alpha}{2\\pi}) - \\frac{\\alpha}{2}W^T W$$<br><br>\n", | |
| "Wに対して定数である項を除くと以下となる。<br><br>\n", | |
| "$$argmax_W (-\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 - \\frac{\\alpha}{2}W^T W)$$<br><br>\n", | |
| "符号を変えて最小化問題に書き換える。<br><br>\n", | |
| "$$argmin_W (\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 + \\frac{\\alpha}{2}W^T W)$$<br><br>\n", | |
| "これは、PRMLのp30の式(1.67)である。" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "ここまでで、尤度$$P(y|X,W,\\beta)$$と、事前分布$$P(W|\\alpha)$$を式で表すことができた。<br>\n", | |
| "また、ベイズの定理より、$$P(W|X,W,\\alpha,\\beta) \\propto P(y|X,W,\\beta)P(W|\\alpha)$$である。<br>\n", | |
| "さらに、これらは正の値であることから、両辺に対数をとり、最大化問題は以下のように書ける。<br><br>\n", | |
| "\n", | |
| "$$argmax_W logP(W|X,W,\\alpha,\\beta) = argmax_W logP(y|X,W,\\beta)P(W|\\alpha)$$<br><br>\n", | |
| "ここで左辺は以下のように変形される。<br><br>\n", | |
| "$$argmax_W -\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 + \\frac{n}{2}log\\beta - \\frac{n}{2}log(2\\pi)\n", | |
| "+ \\frac{M+1}{2}log(\\frac{\\alpha}{2\\pi}) - \\frac{\\alpha}{2}W^T W$$<br><br>\n", | |
| "Wに対して定数である項を除くと以下となる。<br><br>\n", | |
| "$$argmax_W (-\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 - \\frac{\\alpha}{2}W^T W)$$<br><br>\n", | |
| "符号を変えて最小化問題に書き換える。<br><br>\n", | |
| "$$argmin_W (\\frac{\\beta}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 + \\frac{\\alpha}{2}W^T W)$$<br><br>\n", | |
| "これは、PRMLのp30の式(1.67)である。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 93, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "ここで、βは定義上必ず正の値であるので、式をβで割っても同等の最小化問題である。\n", | |
| "つまり以下のように書き換えられる。<br><br>\n", | |
| "$$argmin_W (\\frac{1}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 + \\frac{\\alpha}{2\\beta}W^T W)$$<br><br>\n", | |
| "これは、q=2の正則化項を加えた正則化最小二乗法（Ridge回帰）と全く同じ式である。<br>\n", | |
| "ちなみに、事前分布を正規分布でなくラプラス分布とするとLasso回帰となる。" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "ここで、βは定義上必ず正の値であるので、式をβで割っても同等の最小化問題である。\n", | |
| "つまり以下のように書き換えられる。<br><br>\n", | |
| "$$argmin_W (\\frac{1}{2}\\sum_{i=0}^n (y_i-f(X_i,W))^2 + \\frac{\\alpha}{2\\beta}W^T W)$$<br><br>\n", | |
| "これは、q=2の正則化項を加えた正則化最小二乗法（Ridge回帰）と全く同じ式である。<br>\n", | |
| "ちなみに、事前分布を正規分布でなくラプラス分布とするとLasso回帰となる。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "kernelspec": { | |
| "display_name": "Python 3", | |
| "language": "python", | |
| "name": "python3" | |
| }, | |
| "language_info": { | |
| "codemirror_mode": { | |
| "name": "ipython", | |
| "version": 3 | |
| }, | |
| "file_extension": ".py", | |
| "mimetype": "text/x-python", | |
| "name": "python", | |
| "nbconvert_exporter": "python", | |
| "pygments_lexer": "ipython3", | |
| "version": "3.5.2" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 0 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment