MLE_LeastSquares

MLE_LeastSquares.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MLE and least squares \n",
    "\n",
    "Maximum likelihood estimation(最尤推定)を用いて曲線フィティングを行う。これはつまり、ある観測データが与えられた時に、そのデータがあるパラメータを持つ分布に従って生成されたと仮定し、その観測データを鑑みて、分布のパラメータは何であったと推定するのが尤もらしいか？というのを考えるということ。そしてその尤もらしさ（尤度）を最大化するようなパラメータをもった関数によって曲線フィッティングする。\n",
    "- https://qiita.com/yryrgogo/items/e2774a5a891700c3916a\n",
    "- https://mathtrain.jp/mle\n",
    "- https://mathwords.net/saisyonijoho\n",
    "- https://www.slideshare.net/tanutarou/ss-80824894"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "ある観測データが与えられた時のパラメータの尤度は、以下の様に記述できる。<br><br>\n",
       "$$L(\\theta|y_1,y_2,...,y_n)$$ <br>\n",
       "$$=\\alpha P(y_1,y_2,...,y_n|\\theta)$$<br>\n",
       "$$=\\alpha \\prod_{i=0}^n P(y_i|\\theta)$$<br>\n",
       "ここで、最尤推定における曲線フィッティングの目的は、尤度を最大化するパラメータを求めること、つまり上の式の最大化問題と考えられる。αは正の定数なので、\n",
       "大小関係を変えないため最大化パラメータ算出の目的においては無視できる。つまり以下の最適化問題を解くことになる。<br><br>\n",
       "$$argmax(L(\\theta|y_1,y_2,...,y_n)) = argmax(\\prod_{i=0}^n P(y_i|\\theta))$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "ある観測データが与えられた時のパラメータの尤度は、以下の様に記述できる。<br><br>\n",
    "$$L(\\theta|y_1,y_2,...,y_n)$$ <br>\n",
    "$$=\\alpha P(y_1,y_2,...,y_n|\\theta)$$<br>\n",
    "$$=\\alpha \\prod_{i=0}^n P(y_i|\\theta)$$<br>\n",
    "ここで、最尤推定における曲線フィッティングの目的は、尤度を最大化するパラメータを求めること、つまり上の式の最大化問題と考えられる。αは正の定数なので、\n",
    "大小関係を変えないため最大化パラメータ算出の目的においては無視できる。つまり以下の最適化問題を解くことになる。<br><br>\n",
    "$$argmax(L(\\theta|y_1,y_2,...,y_n)) = argmax(\\prod_{i=0}^n P(y_i|\\theta))$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "ここで、推定したいパラメータθについて考える。パラメータθと簡単に記述したが、実際にはどんな分布のどんなパラメータを推定しようとしているのか？  \n",
       "まず、フィッティングさせたい関数を以下の様に設定する。（次数は任意）<br><br>\n",
       "$$f(x,w)=w_0+w_1 x+w_2 x^2+...+w_n x^n$$<br>\n",
       "\n",
       "ここでは、観測データには真の分布に誤差が加えられていると仮定し、平均0,分散βの正規分布に従う誤差項を加える。  \n",
       "すなわち、各観測データy_iは以下によって生成されていると考える。(ノイズに正規分布を仮定した尤度最大化問題となる。)<br><br>\n",
       "$$f(x_i,w) + \\epsilon _i(\\sim N(0,\\beta))$$<br>\n",
       "\n",
       "これは、以下のように変形できる。<br><br>\n",
       "$$N(f(x_i,w),\\beta)$$<br>\n",
       "イメージとしては、誤差項無しのf(x)のグラフの各x時点ごとに、縦(y軸方向)に正規分布が乗っているようなイメージ。<br>  \n",
       "(ここがわかりやすい: https://www.hellocybernetics.tech/entry/2016/04/24/055832)<br><br>\n",
       "\n",
       "つまり、最初の質問「どんな分布のどんなパラメータを推定しようとしているのか」の答えは、上記の分布のパラメータ<b>(x,w,β)</b>を推定しようとしている、ということになる。\n",
       "ただし、x_1,x_2,...x_nは、y_1,y_2,...y_nと同時に与えられるので、実際にはxとyを用いてwとβを推定する。\n",
       "また、後に出るがβはどのように設定しても尤度を最大化するwの値に影響を与えないため、特に推定する必要はない。"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "ここで、推定したいパラメータθについて考える。パラメータθと簡単に記述したが、実際にはどんな分布のどんなパラメータを推定しようとしているのか？  \n",
    "まず、フィッティングさせたい関数を以下の様に設定する。（次数は任意）<br><br>\n",
    "$$f(x,w)=w_0+w_1 x+w_2 x^2+...+w_n x^n$$<br>\n",
    "\n",
    "ここでは、観測データには真の分布に誤差が加えられていると仮定し、平均0,分散βの正規分布に従う誤差項を加える。  \n",
    "すなわち、各観測データy_iは以下によって生成されていると考える。(ノイズに正規分布を仮定した尤度最大化問題となる。)<br><br>\n",
    "$$f(x_i,w) + \\epsilon _i(\\sim N(0,\\beta))$$<br>\n",
    "\n",
    "これは、以下のように変形できる。<br><br>\n",
    "$$N(f(x_i,w),\\beta)$$<br>\n",
    "イメージとしては、誤差項無しのf(x)のグラフの各x時点ごとに、縦(y軸方向)に正規分布が乗っているようなイメージ。<br>  \n",
    "(ここがわかりやすい: https://www.hellocybernetics.tech/entry/2016/04/24/055832)<br><br>\n",
    "\n",
    "つまり、最初の質問「どんな分布のどんなパラメータを推定しようとしているのか」の答えは、上記の分布のパラメータ<b>(x,w,β)</b>を推定しようとしている、ということになる。\n",
    "ただし、x_1,x_2,...x_nは、y_1,y_2,...y_nと同時に与えられるので、実際にはxとyを用いてwとβを推定する。\n",
    "また、後に出るがβはどのように設定しても尤度を最大化するwの値に影響を与えないため、特に推定する必要はない。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "最大化問題の式に戻る。<br><br>\n",
       "$$argmax(\\prod_{i=0}^n P(y_i|\\theta))$$<br>\n",
       "\n",
       "ここで、推定したいパラメータはx,w,βなので、この式を直すと以下となる。<br><br>\n",
       "$$argmax(\\prod_{i=0}^n P(y_i|x_i,w,\\beta))$$<br>\n",
       "\n",
       "ここで、パラメータx_i,w,βが与えられた際に<b>y_i</b>は<b>N(f(x_i,w),β)</b>によって生成されるので、<b>P(y_i|x_i,w,β)</b>は、\n",
       "分布<b>N(f(x_i,w),β)</b>の確率密度関数によって表すことができる。<br>\n",
       "ここで、正規分布$$N(\\mu,\\sigma^2)$$の確率密度関数は以下である。<br><br>\n",
       "$$f(x)=\\frac{1}{\\sqrt{2\\pi\\sigma^2}} exp(-\\frac{(x-\\mu)^2}{2\\sigma^2}) $$<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",
    "最大化問題の式に戻る。<br><br>\n",
    "$$argmax(\\prod_{i=0}^n P(y_i|\\theta))$$<br>\n",
    "\n",
    "ここで、推定したいパラメータはx,w,βなので、この式を直すと以下となる。<br><br>\n",
    "$$argmax(\\prod_{i=0}^n P(y_i|x_i,w,\\beta))$$<br>\n",
    "\n",
    "ここで、パラメータx_i,w,βが与えられた際に<b>y_i</b>は<b>N(f(x_i,w),β)</b>によって生成されるので、<b>P(y_i|x_i,w,β)</b>は、\n",
    "分布<b>N(f(x_i,w),β)</b>の確率密度関数によって表すことができる。<br>\n",
    "ここで、正規分布$$N(\\mu,\\sigma^2)$$の確率密度関数は以下である。<br><br>\n",
    "$$f(x)=\\frac{1}{\\sqrt{2\\pi\\sigma^2}} exp(-\\frac{(x-\\mu)^2}{2\\sigma^2}) $$<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": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "ここで、このまま最適化問題を解くのは困難なので、対数をとる。(大小関係は変わらないので問題ない。)<br><br>\n",
       "$$L(x,w,\\beta|y_1,y_2,...,y_n)=\\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(x_i,w))^2}{2\\beta})$$<br>\n",
       "$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}=log{(\\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(x_i,w))^2}{2\\beta}))}\n",
       "=-\\frac{n}{2}\\log{(2\\pi\\beta)}-\\frac{1}{2\\beta}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$<br><br>\n",
       "\n",
       "-----------------<br>\n",
       "ここで一点注意: PRMLのp29の式(1.62)では、分散を$$\\beta$$ではなく$$\\beta^{-1}$$と置いているため、\n",
       "上の式と若干異なっている点に注意。(βを1/βとおきかえると同じ式になる。)<br>\n",
       "-----------------<br><br>\n",
       "\n",
       "ここで、β(誤差項の分散)はコントロール不能なパラメータであるので、定数として扱う。（β=ρ^2であるので、βは常に正の値である）<br>\n",
       "定数部分をCとおくと以下のような式になる。<br><br>\n",
       "$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}=C_1-C_2\\sum_{i=0}^n (y_i-f(x_i,w))^2$$<br>\n",
       "ここで、<br><br>\n",
       "$$C_1 = -\\frac{n}{2}\\log{(2\\pi\\beta)}$$<br>\n",
       "$$C_2 = \\frac{1}{2\\beta}$$<br>\n",
       "である。c_2は必ず正の値となるため、この最大化問題はすなわち、以下の最小化問題となる。<br><br>\n",
       "$$argmin(\\sum_{i=0}^n (y_i-f(x_i,w))^2)$$<br><br>\n",
       "つまり、解法は最小二乗法と同様。別の言い方をすれば、最小二乗法は、ノイズに正規分布を仮定した最尤推定である。"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "ここで、このまま最適化問題を解くのは困難なので、対数をとる。(大小関係は変わらないので問題ない。)<br><br>\n",
    "$$L(x,w,\\beta|y_1,y_2,...,y_n)=\\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(x_i,w))^2}{2\\beta})$$<br>\n",
    "$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}=log{(\\prod_{i=0}^n \\frac{1}{\\sqrt{2\\pi\\beta}} exp(-\\frac{(y_i-f(x_i,w))^2}{2\\beta}))}\n",
    "=-\\frac{n}{2}\\log{(2\\pi\\beta)}-\\frac{1}{2\\beta}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$<br><br>\n",
    "\n",
    "-----------------<br>\n",
    "ここで一点注意: PRMLのp29の式(1.62)では、分散を$$\\beta$$ではなく$$\\beta^{-1}$$と置いているため、\n",
    "上の式と若干異なっている点に注意。(βを1/βとおきかえると同じ式になる。)<br>\n",
    "-----------------<br><br>\n",
    "\n",
    "ここで、β(誤差項の分散)はコントロール不能なパラメータであるので、定数として扱う。（β=ρ^2であるので、βは常に正の値である）<br>\n",
    "定数部分をCとおくと以下のような式になる。<br><br>\n",
    "$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}=C_1-C_2\\sum_{i=0}^n (y_i-f(x_i,w))^2$$<br>\n",
    "ここで、<br><br>\n",
    "$$C_1 = -\\frac{n}{2}\\log{(2\\pi\\beta)}$$<br>\n",
    "$$C_2 = \\frac{1}{2\\beta}$$<br>\n",
    "である。c_2は必ず正の値となるため、この最大化問題はすなわち、以下の最小化問題となる。<br><br>\n",
    "$$argmin(\\sum_{i=0}^n (y_i-f(x_i,w))^2)$$<br><br>\n",
    "つまり、解法は最小二乗法と同様。別の言い方をすれば、最小二乗法は、ノイズに正規分布を仮定した最尤推定である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "では、最小二乗法としてではなく最尤法として考えることにメリットはあるのか？一つは、wの推定だけでなく、ノイズ項の分散βも推定可能となること。\n",
       "また、ノイズをガウス分布以外の分布で仮定しても上記の枠組みで推定できるようになること。<br>\n",
       "ノイズ項の分散βを推定するには、対数をとった尤度関数をwでなくβで偏微分して0とおけばよい。\n",
       "<br><br>\n",
       "$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}\n",
       "=-\\frac{n}{2}\\log{(2\\pi\\beta)}-\\frac{1}{2\\beta}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$<br><br>\n",
       "\n",
       "これをβで偏微分し0と置く。<br>\n",
       "\n",
       "$$\\frac{\\partial L}{\\partial \\beta} = \n",
       "-\\frac{n}{2\\beta}+\\frac{1}{2\\beta^2}\\sum_{i=0}^n (y_i-f(x_i,w))^2 = 0$$<br>\n",
       "$$\\beta = \\frac{1}{n}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "では、最小二乗法としてではなく最尤法として考えることにメリットはあるのか？一つは、wの推定だけでなく、ノイズ項の分散βも推定可能となること。\n",
    "また、ノイズをガウス分布以外の分布で仮定しても上記の枠組みで推定できるようになること。<br>\n",
    "ノイズ項の分散βを推定するには、対数をとった尤度関数をwでなくβで偏微分して0とおけばよい。\n",
    "<br><br>\n",
    "$$\\log{L(x,w,\\beta|y_1,y_2,...,y_n)}\n",
    "=-\\frac{n}{2}\\log{(2\\pi\\beta)}-\\frac{1}{2\\beta}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$<br><br>\n",
    "\n",
    "これをβで偏微分し0と置く。<br>\n",
    "\n",
    "$$\\frac{\\partial L}{\\partial \\beta} = \n",
    "-\\frac{n}{2\\beta}+\\frac{1}{2\\beta^2}\\sum_{i=0}^n (y_i-f(x_i,w))^2 = 0$$<br>\n",
    "$$\\beta = \\frac{1}{n}\\sum_{i=0}^n (y_i-f(x_i,w))^2$$"
   ]
  },
  {
   "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
}