Extended Logistic Regression

ExtendedLogisticRegression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logistic Regression (Multi-class and Regularization)\n",
    "- 「データ解析のための統計モデリング入門」(久保,2012) 6章\n",
    "- [Simplified Cost Function and Gradient Descent (Cousera)](https://www.coursera.org/lecture/machine-learning/simplified-cost-function-and-gradient-descent-MtEaZ)\n",
    "- [Multiclass Classification: One-vs-all (Cousera)](https://www.coursera.org/lecture/machine-learning/multiclass-classification-one-vs-all-68Pol)\n",
    "- [Regularized Logistic Regression (Cousera)](https://www.coursera.org/lecture/machine-learning/regularized-logistic-regression-4BHEy)\n",
    "\n",
    "[ここ](https://gist.github.com/shotahorii/63b10aa9337a4e3893ed3dac28af84e5)でロジスティック回帰について書いた。ただ、かなりさらっと書いただけだったので、もう少し整理して書いてみたいと思う。\n",
    "  \n",
    "ロジスティック回帰はGLMの一つのインスタンスであると言える。すなわち、目的変数の**誤差構造**を想定し、説明変数を線形結合して**線形予測子**を作り、**リンク関数**で目的変数(あるいは目的変数の元になる分布のパラメータ)と線形予測子をつなぐ。\n",
    "\n",
    "最小二乗法やポアソン回帰もGLMの一つのインスタンスであると言える。例えばポアソン回帰であれば、誤差構造にポアソン分布を想定し、リンク関数は一般に対数関数を用いるのだった。ここで、ポアソン分布のキモは誤差構造にポアソン分布を仮定する点である。例えば目的変数が0より小さい値を取らないカウントデータなどである場合、誤差がポアソン分布に従う可能性が考えられ、ポアソン回帰が手法の候補となる。一方で、リンク関数は誤差構造ほど厳密ではなく、対数関数でなくidentity関数でも(式の範囲においてλ<0になり得るという点で最適とは言い難いが)問題はない。\n",
    "\n",
    "では、ロジスティック回帰はどうか。ロジスティック回帰では、誤差関数に二項分布を想定し、リンク関数にロジット関数を用いる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**二項分布**  \n",
    "ではここで、二項分布とはどんな分布だったか。\n",
    "二項分布は、成功確率qで結果が成功か失敗のいずれかである試行（ベルヌーイ試行）を独立にn回行った時の成功回数を確率変数とする分布である。例えば、n個中y個xxxであった、というようなデータが考えられる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "二項分布の確率質量関数は以下である。\n",
       "<br><br>\n",
       "$$p(y|n,q)= \\left(\\begin{array}{cc} n \\\\ y \\end{array} \\right) q^y (1-q)^{n-y} \n",
       "= \\frac{n!}{y!(n-y)!} q^y (1-q)^{n-y}$$\n",
       "<br><br>\n",
       "ここで、平均と分散は以下である。\n",
       "<br><br>\n",
       "$$E(y)=nq$$<br>\n",
       "$$Var(y)=nq(1-q)$$\n",
       "<br><br>\n",
       "平均が$$nq$$になるというのは直感的にとてもわかりやすい。例えば表が出る確率$$q=0.5$$のコインを\n",
       "$$n=20$$回投げたら、表が出る回数の期待値は$$n*q=20*0.5=10$$になるだろう、と直感的にわかる。\n",
       "<br><br>\n",
       "分散については、まず$$n$$が大きい時分散が大きくなる、というのは直感的にほぼ自明だ。\n",
       "(表が出る確率0.5のコインを2回投げた時よりも、200回投げた時ばらつきが大きいのは当然。)\n",
       "<br>\n",
       "$$q$$についてはどうだろうか。$$q$$が0.5に近いほど分散は大きくなり、0と1に近いほど分散は小さくなる。\n",
       "これも直感的にわかりやすい。例えば絶対に表がでるコイン$$q=1$$や絶対に裏しかでないコイン$$q=0$$の分散は0であろう。\n",
       "<br>\n",
       "そこまで極端でなくても、例えば99.99%の確率で表が出るコインを複数回投げた時の分散はかなり小さそうだ、とわかる。"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "二項分布の確率質量関数は以下である。\n",
    "<br><br>\n",
    "$$p(y|n,q)= \\left(\\begin{array}{cc} n \\\\ y \\end{array} \\right) q^y (1-q)^{n-y} \n",
    "= \\frac{n!}{y!(n-y)!} q^y (1-q)^{n-y}$$\n",
    "<br><br>\n",
    "ここで、平均と分散は以下である。\n",
    "<br><br>\n",
    "$$E(y)=nq$$<br>\n",
    "$$Var(y)=nq(1-q)$$\n",
    "<br><br>\n",
    "平均が$$nq$$になるというのは直感的にとてもわかりやすい。例えば表が出る確率$$q=0.5$$のコインを\n",
    "$$n=20$$回投げたら、表が出る回数の期待値は$$n*q=20*0.5=10$$になるだろう、と直感的にわかる。\n",
    "<br><br>\n",
    "分散については、まず$$n$$が大きい時分散が大きくなる、というのは直感的にほぼ自明だ。\n",
    "(表が出る確率0.5のコインを2回投げた時よりも、200回投げた時ばらつきが大きいのは当然。)\n",
    "<br>\n",
    "$$q$$についてはどうだろうか。$$q$$が0.5に近いほど分散は大きくなり、0と1に近いほど分散は小さくなる。\n",
    "これも直感的にわかりやすい。例えば絶対に表がでるコイン$$q=1$$や絶対に裏しかでないコイン$$q=0$$の分散は0であろう。\n",
    "<br>\n",
    "そこまで極端でなくても、例えば99.99%の確率で表が出るコインを複数回投げた時の分散はかなり小さそうだ、とわかる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q=0.05の時 var= 0.0475\n",
      "q=0.2の時 var= 0.16000000000000003\n",
      "q=0.5の時 var= 0.25\n",
      "q=0.999の時 var= 0.000999000000000001\n"
     ]
    }
   ],
   "source": [
    "#具体的に計算してみる。(ここではn=1とする)\n",
    "v = lambda q:q*(1-q)\n",
    "print('q=0.05の時 var=',v(0.05))\n",
    "print('q=0.2の時 var=',v(0.2))\n",
    "print('q=0.5の時 var=',v(0.5))\n",
    "print('q=0.999の時 var=',v(0.999))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ロジスティック回帰**  \n",
    "さて、話をロジスティック回帰に戻そう。ロジスティック回帰では、誤差構造に二項分布を想定するのであった。これはつまり(話を簡単にするためここでは説明変数を**x**ただ一つとする)、**x**がある値をとるとき、その値に基づいて二項分布のパラメータ**q**が決まり、その**q**をパラメータとする二項分布にしたがってその**x**の値における**y**が生成される、という想定である。\n",
    "\n",
    "**n**はどこにいったんじゃい、という話であるが、これは扱う問題(=用いるデータ)に依る値である。例えば、久保本のデータのように、あるデータ(1レコード)がある植物の1個体を表し、「植物の体サイズ**x**」と「植物の種子数**n(常に8)**」と「その種子中生き残った種子数**y**」であるとすると、これを回帰する際の二項分布の**n**は8となる。   \n",
    "~~例えばこのデータで、植物の種子数**n**が個体ごとに異なる場合(n=常に8という前提がない時)には、**y'=y/n**、 **n'=1**とし、回帰分析の二項分布の**n**を1として考える必要がある。~~(←二項分布を想定しているのでyが整数以外になってはまずい。nの数はデータレコードごとに異なっていてもロジスティック回帰は可能。)  \n",
    "\n",
    "さて、であるので、回帰で求めたいパラメータは二項分布の**q**である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "で、である。いま、線形予測子を以下のように考える。\n",
       "<br><br>\n",
       "$$\\beta_0+\\beta_1x$$\n",
       "<br><br>\n",
       "この回帰では、xの値によって決まるパラメータqを持つ二項分布にyが従う、と仮定していると上で書いた。\n",
       "すなわち、この線形予測子の値によって$$q$$が決まる、ということ。$$x$$はすでにデータとして与えられているので、\n",
       "ここで行うのはこの$$\\beta_0,\\beta_1$$の値をいじくることで各$$x$$の値における$$q$$を決めることである。\n",
       "<br><br>\n",
       "さて、シンプルなのは以下のように考えることである。\n",
       "<br><br>\n",
       "$$q_i = \\beta_0+\\beta_1x_i$$\n",
       "<br><br>\n",
       "でも本当にこれで良いのか？全然良くない。なぜなら、そもそも$$q$$がなんであったかを考えると、$$q$$\n",
       "は二項分布における成功確率のパラメータであったので、$$0 \\le q \\le 1$$である。<br>\n",
       "上のような式だと、$$\\beta_0,\\beta_1$$の値によって$$q$$は容易に0~1の範囲を超えてしまう。<br>\n",
       "そこで登場するのがロジット関数である。\n",
       "<br><br>\n",
       "$$logit(q)=log(\\frac{q}{1-q})$$\n",
       "<br><br>\n",
       "で、これをリンク関数として用いて、\n",
       "<br><br>\n",
       "$$logit(q_i)=log(\\frac{q_i}{1-q_i}) = \\beta_0+\\beta_1x_i$$\n",
       "<br><br>\n",
       "とすると、線形予測子が$$-\\infty, \\infty$$の範囲でどのような値をとっても必ず$$q_i$$は0~1の間に収まる。\n",
       "<br><br>\n",
       "で、これを$$q=...$$の形に書き換える。ロジット関数の逆関数であるロジスティック関数を使って以下のように書き換えられる。\n",
       "<br><br>\n",
       "$$q_i = logistic(\\beta_0+\\beta_1x_i) = \\frac{1}{1+exp(-(\\beta_0+\\beta_1x_i))}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "で、である。いま、線形予測子を以下のように考える。\n",
    "<br><br>\n",
    "$$\\beta_0+\\beta_1x$$\n",
    "<br><br>\n",
    "この回帰では、xの値によって決まるパラメータqを持つ二項分布にyが従う、と仮定していると上で書いた。\n",
    "すなわち、この線形予測子の値によって$$q$$が決まる、ということ。$$x$$はすでにデータとして与えられているので、\n",
    "ここで行うのはこの$$\\beta_0,\\beta_1$$の値をいじくることで各$$x$$の値における$$q$$を決めることである。\n",
    "<br><br>\n",
    "さて、シンプルなのは以下のように考えることである。\n",
    "<br><br>\n",
    "$$q_i = \\beta_0+\\beta_1x_i$$\n",
    "<br><br>\n",
    "でも本当にこれで良いのか？全然良くない。なぜなら、そもそも$$q$$がなんであったかを考えると、$$q$$\n",
    "は二項分布における成功確率のパラメータであったので、$$0 \\le q \\le 1$$である。<br>\n",
    "上のような式だと、$$\\beta_0,\\beta_1$$の値によって$$q$$は容易に0~1の範囲を超えてしまう。<br>\n",
    "そこで登場するのがロジット関数である。\n",
    "<br><br>\n",
    "$$logit(q)=log(\\frac{q}{1-q})$$\n",
    "<br><br>\n",
    "で、これをリンク関数として用いて、\n",
    "<br><br>\n",
    "$$logit(q_i)=log(\\frac{q_i}{1-q_i}) = \\beta_0+\\beta_1x_i$$\n",
    "<br><br>\n",
    "とすると、線形予測子が$$-\\infty, \\infty$$の範囲でどのような値をとっても必ず$$q_i$$は0~1の間に収まる。\n",
    "<br><br>\n",
    "で、これを$$q=...$$の形に書き換える。ロジット関数の逆関数であるロジスティック関数を使って以下のように書き換えられる。\n",
    "<br><br>\n",
    "$$q_i = logistic(\\beta_0+\\beta_1x_i) = \\frac{1}{1+exp(-(\\beta_0+\\beta_1x_i))}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "と、いうことで、「何回中何回成功」だとか、「この人はxxxする確率がどれくらい、この人はどれくらい、...」とか、そんな感じのデータであればロジスティック回帰を用いて予測できる可能性がある。後者の例については、0.5を閾値として0/1に変換すれば「この人はxxxした、この人はしなかった、この人は、...」というようなデータについての分類問題に対応できる。実際にはこの分類としての使い方が圧倒的にメジャーだと思う。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "さて、ここまでで誤差構造、線形予測子、リンク関数が揃った。ここからパラメータ推定したい。<br>\n",
       "上で書いたように、$$x_i$$の値から、その値における$$q_i$$が計算されるのであった。で、$$q_i$$\n",
       "が決まった時、その$$q_i$$をパラメータとする二項分布の元である$$y_i$$\n",
       "が出現する確率は、二項分布の確率質量関数より以下のように計算できるのであった。\n",
       "<br><br>\n",
       "$$p(y_i|q_i)=\\left(\\begin{array}{cc} n_i \\\\ y_i \\end{array} \\right) q_i^{y_i} (1-q_i)^{n_i-y_i}$$\n",
       "<br><br>\n",
       "パラメータを最尤推定したいので、あるパラメータ$$\\beta_0,\\beta_1$$をとった時の各$$y_i$$\n",
       "の出現する確率の積(=尤度)を最大化するようなパラメータが推定パラメータとなる。\n",
       "<br><br>\n",
       "$$L(\\beta) = \\prod_i\\left(\\begin{array}{cc} n_i \\\\ y_i \\end{array} \\right) q_i^{y_i} (1-q_i)^{n_i-y_i}$$\n",
       "<br><br>\n",
       "対数をとって対数尤度とする。\n",
       "<br><br>\n",
       "$$logL(\\beta) = \\sum_i\\{ log\\left(\\begin{array}{cc} n_i \\\\ y_i \\end{array} \\right) + y_ilog(q_i)\n",
       "                        + (n_i-y_i)log(1-q_i)\\}$$\n",
       "<br><br>\n",
       "で、これは最大化問題であったので、パラメータ$$\\beta$$について定数である項は式から排除できる。\n",
       "<br><br>\n",
       "$$argmax_\\beta logL(\\beta) = argmax_\\beta \\sum_i\\{ y_ilog(q_i) + (n_i-y_i)log(1-q_i)\\}$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "さて、ここまでで誤差構造、線形予測子、リンク関数が揃った。ここからパラメータ推定したい。<br>\n",
    "上で書いたように、$$x_i$$の値から、その値における$$q_i$$が計算されるのであった。で、$$q_i$$\n",
    "が決まった時、その$$q_i$$をパラメータとする二項分布の元である$$y_i$$\n",
    "が出現する確率は、二項分布の確率質量関数より以下のように計算できるのであった。\n",
    "<br><br>\n",
    "$$p(y_i|q_i)=\\left(\\begin{array}{cc} n_i \\\\ y_i \\end{array} \\right) q_i^{y_i} (1-q_i)^{n_i-y_i}$$\n",
    "<br><br>\n",
    "パラメータを最尤推定したいので、あるパラメータ$$\\beta_0,\\beta_1$$をとった時の各$$y_i$$\n",
    "の出現する確率の積(=尤度)を最大化するようなパラメータが推定パラメータとなる。\n",
    "<br><br>\n",
    "$$L(\\beta) = \\prod_i\\left(\\begin{array}{cc} n_i \\\\ y_i \\end{array} \\right) q_i^{y_i} (1-q_i)^{n_i-y_i}$$\n",
    "<br><br>\n",
    "対数をとって対数尤度とする。\n",
    "<br><br>\n",
    "$$logL(\\beta) = \\sum_i\\{ log\\left(\\begin{array}{cc} n_i \\\\ y_i \\end{array} \\right) + y_ilog(q_i)\n",
    "                        + (n_i-y_i)log(1-q_i)\\}$$\n",
    "<br><br>\n",
    "で、これは最大化問題であったので、パラメータ$$\\beta$$について定数である項は式から排除できる。\n",
    "<br><br>\n",
    "$$argmax_\\beta logL(\\beta) = argmax_\\beta \\sum_i\\{ y_ilog(q_i) + (n_i-y_i)log(1-q_i)\\}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "つまりこの上の式の右辺の最大化問題を解くことでパラメータを求めるのである。  \n",
    "この時点でもうわかると思うが、よくある0/1のロジスティック回帰は上の式においてn=1、すなわちyが0か1かしか取らないケースと考えれば良い。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さて、久保本ではここから先(実際に最大化問題を解いて最尤推定を行う部分)についてはライブラリで計算しましょう、として詳しい手法については触れられていなかった。\n",
    "\n",
    "この最大化問題については最小二乗法の二乗誤差最小化のように解析的に解けない(はじパタp94)ので、最適化手法を用いて数値的に解く必要がある。勾配降下法(Gradient Descent)でこれを実装する方法がAndrew Ng先生の[講座](https://www.coursera.org/lecture/machine-learning/simplified-cost-function-and-gradient-descent-MtEaZ)で説明されている。\n",
    "\n",
    "**ロジスティック回帰の拡張**  \n",
    "勾配降下法とそれを使ったロジスティック回帰の実装については別の記事で書きたいと思うが、この記事ではもう少し、ロジスティック回帰の拡張について書いていきたいと思う。\n",
    "\n",
    "ちょっとした注意点として、ここまでの記事で書いてきた数式の記述(久保本に準じている)とAndrew講座でいくつか数式の書き方に違いがあるので以下に書いておく。\n",
    "1. ここまでは目的変数についてn回中y回、というnが任意の二項分布を想定して書いてきたが、Andrew講座ではn=1、y=0or1の分類問題を想定して書いている。\n",
    "2. ここまではパラメータをβで表したが、Andrew講座ではθで表す。\n",
    "3. ここまでは、二項分布のパラメータqを単にq<sub>i</sub>と書いてきたが、Andrew講座では、xを引数にとるパラメータθが与えられた関数によって定義される値であるということで、h<sub>θ</sub>(x<sub>i</sub>)と書かれる。\n",
    "4. 最大化する対数尤度関数について、ここまでは単純に全データ数のSUMとしているが、Andrew講座ではデータ数mで割られている。(mは推定パラメータに依らない値なので、別にあってもなくても解は変わらない)\n",
    "5. 最後に、Andrew講座では最大化問題の式の符号を正負転換して最小化問題としている。つまり尤度の負をとったものをコスト関数ととらえ最小化問題に変えている。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "上にリストした違いを上の尤度の数式に反映させると以下になる。\n",
       "<br><br>\n",
       "<b>久保本に準じた書き方</b>\n",
       "<br><br>\n",
       "$$argmax_\\beta \\sum_i\\{ y_ilog(q_i) + (n_i-y_i)log(1-q_i)\\}$$\n",
       "<br><br>\n",
       "<b>Andrew講座に準じた書き方</b>\n",
       "<br><br>\n",
       "$$argmin_\\theta -[\\frac{1}{m}\\sum_i^m\\{ y_i log h_\\theta(x_i) + (1-y_i)log(1-h_\\theta(x_i))\\}]$$\n",
       "<br><br>\n",
       "注意点として、久保本では$$n_i$$は任意の正の整数であるが、Andrew講座では$$1$$と固定されているので\n",
       "完全に同一の式ではないが、$$n_i=1$$となっている以外は全く同じ式である。\n",
       "<br><br>\n",
       "また、yの確率についても以下のような違いがある。\n",
       "<br><br>\n",
       "<b>久保本に準じた書き方</b>\n",
       "<br><br>\n",
       "$$p(y|q)$$\n",
       "<br><br>\n",
       "<b>Andrew講座に準じた書き方</b>\n",
       "<br><br>\n",
       "$$p(y|x;\\theta)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "上にリストした違いを上の尤度の数式に反映させると以下になる。\n",
    "<br><br>\n",
    "<b>久保本に準じた書き方</b>\n",
    "<br><br>\n",
    "$$argmax_\\beta \\sum_i\\{ y_ilog(q_i) + (n_i-y_i)log(1-q_i)\\}$$\n",
    "<br><br>\n",
    "<b>Andrew講座に準じた書き方</b>\n",
    "<br><br>\n",
    "$$argmin_\\theta -[\\frac{1}{m}\\sum_i^m\\{ y_i log h_\\theta(x_i) + (1-y_i)log(1-h_\\theta(x_i))\\}]$$\n",
    "<br><br>\n",
    "注意点として、久保本では$$n_i$$は任意の正の整数であるが、Andrew講座では$$1$$と固定されているので\n",
    "完全に同一の式ではないが、$$n_i=1$$となっている以外は全く同じ式である。\n",
    "<br><br>\n",
    "また、yの確率についても以下のような違いがある。\n",
    "<br><br>\n",
    "<b>久保本に準じた書き方</b>\n",
    "<br><br>\n",
    "$$p(y|q)$$\n",
    "<br><br>\n",
    "<b>Andrew講座に準じた書き方</b>\n",
    "<br><br>\n",
    "$$p(y|x;\\theta)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "で、だ。少し横道にそれたが、ロジスティック回帰の拡張の話であった。\n",
    "ここでは二つの拡張について触れる。\n",
    "1. 多クラス分類\n",
    "2. 正則化\n",
    "\n",
    "**多クラス分類**  \n",
    "One-vs-all あるいは One-vs-rest と呼ばれる手法でロジスティック回帰を多クラス分類に拡張する。\n",
    "\n",
    "1,2,3の三つのクラスを持つ目的変数yを例に多クラス分類問題を考えていく。  \n",
    "まず最初のステップとして、1,2,3それぞれのクラスについて通常のロジスティック回帰モデルを作る。つまり以下のように、クラス数分の回帰モデルができる。  \n",
    "\n",
    "(モデル1) y=1の時y'=1, それ以外の時をy'=0として、y'を目的変数とした回帰モデル  \n",
    "(モデル2) y=2の時y'=1, それ以外の時をy'=0として、y'を目的変数とした回帰モデル  \n",
    "(モデル3) y=3の時y'=1, それ以外の時をy'=0として、y'を目的変数とした回帰モデル"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "ここで、それぞれのモデルが何をしているのか、という意味を考えてみよう。\n",
       "<br><br>\n",
       "ロジスティック回帰では、yが二項分布に従って確率的に出現すると仮定するのであった。つまり二項分布のパラメータ$$n,q$$\n",
       "が与えられたとき、$$y=k$$である確率は以下のように表される。\n",
       "<br><br>\n",
       "$$p(y=k|n,q)= \\left(\\begin{array}{cc} n \\\\ k \\end{array} \\right) q^k (1-q)^{n-k}$$\n",
       "<br><br>\n",
       "ここでモデル1~モデル3の各モデルでは、$$y=1 \\, or \\, 0$$の二値分類を考えている。すなわち上の式では\n",
       "$$n=1, y=1 \\, or \\, 0$$である。\n",
       "<br><br>\n",
       "ここで、$$y=1$$の場合を考えてみよう。\n",
       "<br><br>\n",
       "$$p(y=1|n,q)=\\left(\\begin{array}{cc} n \\\\ y \\end{array} \\right) q^{y} (1-q)^{n-y}\n",
       "=\\left(\\begin{array}{cc} 1 \\\\ 1 \\end{array} \\right) q^1 (1-q)^{(1-1)}=q$$\n",
       "<br><br>\n",
       "つまり、だ。モデル1~モデル3の各モデルが何をしているのかというと、(ざっくり言うと)\n",
       "あるxが与えられた時に各クラスに分類される確率、である。\n",
       "<br>\n",
       "モデル1が出すのは$$p(y=1|x;\\theta)$$で、モデル2が出すのは$$p(y=2|x;\\theta)$$、モデル3は\n",
       "$$p(y=3|x;\\theta)$$である。\n",
       "<br><br>\n",
       "つまりクラスごとにロジスティック回帰モデルを作ることで、そのクラスに分類される確率が算出された。\n",
       "単純に各xについて、最も確率の高いクラスに分類すれば良いだけである。"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "ここで、それぞれのモデルが何をしているのか、という意味を考えてみよう。\n",
    "<br><br>\n",
    "ロジスティック回帰では、yが二項分布に従って確率的に出現すると仮定するのであった。つまり二項分布のパラメータ$$n,q$$\n",
    "が与えられたとき、$$y=k$$である確率は以下のように表される。\n",
    "<br><br>\n",
    "$$p(y=k|n,q)= \\left(\\begin{array}{cc} n \\\\ k \\end{array} \\right) q^k (1-q)^{n-k}$$\n",
    "<br><br>\n",
    "ここでモデル1~モデル3の各モデルでは、$$y=1 \\, or \\, 0$$の二値分類を考えている。すなわち上の式では\n",
    "$$n=1, y=1 \\, or \\, 0$$である。\n",
    "<br><br>\n",
    "ここで、$$y=1$$の場合を考えてみよう。\n",
    "<br><br>\n",
    "$$p(y=1|n,q)=\\left(\\begin{array}{cc} n \\\\ y \\end{array} \\right) q^{y} (1-q)^{n-y}\n",
    "=\\left(\\begin{array}{cc} 1 \\\\ 1 \\end{array} \\right) q^1 (1-q)^{(1-1)}=q$$\n",
    "<br><br>\n",
    "つまり、だ。モデル1~モデル3の各モデルが何をしているのかというと、(ざっくり言うと)\n",
    "あるxが与えられた時に各クラスに分類される確率、である。\n",
    "<br>\n",
    "モデル1が出すのは$$p(y=1|x;\\theta)$$で、モデル2が出すのは$$p(y=2|x;\\theta)$$、モデル3は\n",
    "$$p(y=3|x;\\theta)$$である。\n",
    "<br><br>\n",
    "つまりクラスごとにロジスティック回帰モデルを作ることで、そのクラスに分類される確率が算出された。\n",
    "単純に各xについて、最も確率の高いクラスに分類すれば良いだけである。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "これをAndrew講座に準じた書き方で書くと以下のようになる。(クラス数を$$k$$とする。)\n",
       "<br><br>\n",
       "$$h_\\theta^{(i)}(x) = p(y=i|x;\\theta) $$<br>\n",
       "$$(i=1,2,3,...,k)$$\n",
       "<br><br>\n",
       "で、予測に際しては以下のように最大確率のクラスを選ぶ。\n",
       "<br><br>\n",
       "$$max_i h_\\theta^{(i)}(x)$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "これをAndrew講座に準じた書き方で書くと以下のようになる。(クラス数を$$k$$とする。)\n",
    "<br><br>\n",
    "$$h_\\theta^{(i)}(x) = p(y=i|x;\\theta) $$<br>\n",
    "$$(i=1,2,3,...,k)$$\n",
    "<br><br>\n",
    "で、予測に際しては以下のように最大確率のクラスを選ぶ。\n",
    "<br><br>\n",
    "$$max_i h_\\theta^{(i)}(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**正則化ロジスティック回帰**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "上で書いたように、ロジスティック回帰の尤度最大化は符号の正負を逆転して以下のようなコスト関数の最小化問題と捉えられる。\n",
       "(データ数$$m$$とする)\n",
       "<br><br>\n",
       "$$argmin_\\theta -[\\frac{1}{m}\\sum_i^m\\{ y_i log h_\\theta(x_i) + (1-y_i)log(1-h_\\theta(x_i))\\}]$$\n",
       "<br><br>\n",
       "このコスト関数に正則化項を追加してパラメータが過大になるのを防ぐ。<br>\n",
       "ここで、仮に線形予測子が以下のようであるとしよう。(説明変数の数を$$k$$とする)\n",
       "<br><br>\n",
       "$$logit(q_i) = \\theta_0 + \\theta_1x_1 + \\theta_2x_2 + ... + \\theta_kx_k$$\n",
       "<br><br>\n",
       "ちなみに、重要なのは切片項を除いた推定すべきパラメータ数がk個であるという点。以下のような線形予測子でも同じこと。\n",
       "<br><br>\n",
       "$$logit(q_i) = \\theta_0 + \\theta_1x_1x_2x_3^2 + \\theta_2x_2^3 + ... + \\theta_kx_1^2$$\n",
       "<br><br>\n",
       "この時、正則化項は以下となる。\n",
       "<br><br>\n",
       "$$\\frac{\\lambda}{2m}\\sum_{j=1}^k\\theta_j^2$$\n",
       "<br><br>\n",
       "つまり、最小化したいコスト関数は以下となる。\n",
       "<br><br>\n",
       "$$J(\\theta)=-[\\frac{1}{m}\\sum_i^m\\{ y_i log h_\\theta(x_i) + (1-y_i)log(1-h_\\theta(x_i))\\}]\n",
       "+ \\frac{\\lambda}{2m}\\sum_{j=1}^k\\theta_j^2$$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%latex\n",
    "上で書いたように、ロジスティック回帰の尤度最大化は符号の正負を逆転して以下のようなコスト関数の最小化問題と捉えられる。\n",
    "(データ数$$m$$とする)\n",
    "<br><br>\n",
    "$$argmin_\\theta -[\\frac{1}{m}\\sum_i^m\\{ y_i log h_\\theta(x_i) + (1-y_i)log(1-h_\\theta(x_i))\\}]$$\n",
    "<br><br>\n",
    "このコスト関数に正則化項を追加してパラメータが過大になるのを防ぐ。<br>\n",
    "ここで、仮に線形予測子が以下のようであるとしよう。(説明変数の数を$$k$$とする)\n",
    "<br><br>\n",
    "$$logit(q_i) = \\theta_0 + \\theta_1x_1 + \\theta_2x_2 + ... + \\theta_kx_k$$\n",
    "<br><br>\n",
    "ちなみに、重要なのは切片項を除いた推定すべきパラメータ数がk個であるという点。以下のような線形予測子でも同じこと。\n",
    "<br><br>\n",
    "$$logit(q_i) = \\theta_0 + \\theta_1x_1x_2x_3^2 + \\theta_2x_2^3 + ... + \\theta_kx_1^2$$\n",
    "<br><br>\n",
    "この時、正則化項は以下となる。\n",
    "<br><br>\n",
    "$$\\frac{\\lambda}{2m}\\sum_{j=1}^k\\theta_j^2$$\n",
    "<br><br>\n",
    "つまり、最小化したいコスト関数は以下となる。\n",
    "<br><br>\n",
    "$$J(\\theta)=-[\\frac{1}{m}\\sum_i^m\\{ y_i log h_\\theta(x_i) + (1-y_i)log(1-h_\\theta(x_i))\\}]\n",
    "+ \\frac{\\lambda}{2m}\\sum_{j=1}^k\\theta_j^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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
}