chapter9_bayes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9.1 ベイズの定理\n",
    "2つの壺 $H_1$, $H_2$ がある. \n",
    "\n",
    "$H_1$ には赤玉が $30$ 個, 青玉が $70$ 個ある. \n",
    "\n",
    "$H_2$ には赤玉が $70$ 個, 青玉が $30$ 個ある. \n",
    "\n",
    "2つの壺のから<u>ランダム</u>に壺を選んで, 12回玉を復元抽出した結果, 赤玉8個, 青玉4個という結果だった. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **問題：選ばれた壺は$H_1$, $H_2$ のどちら？**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 回答\n",
    "\n",
    "求めるもの：壺 $H_1$ が選ばれたと言える比率。\n",
    "$$\n",
    "\\dfrac{\\textrm{壺 $H_1$ を選び, 抽出結果が赤玉8個, 青玉4個となる確率}}{\\textrm{抽出結果が赤玉8個, 青玉4個となる確率}}\n",
    "$$\n",
    "ところで, \n",
    "\\begin{align*}\n",
    "&\\textrm{抽出結果が赤玉8個, 青玉4個となる確率}\n",
    "\\\\\n",
    " &= \\textrm{壺 $H_1$ を選び, 抽出結果が赤玉8個, 青玉4個となる確率} \n",
    "\\\\ \n",
    " &\\quad + \\textrm{壺 $H_2$ を選び, 抽出結果が赤玉8個, 青玉4個となる確率}   \n",
    "\\end{align*}\n",
    "である. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$P(H_1)$：$H_1$の壺を選ぶ確率. \n",
    "\n",
    "$A$ : 抽出結果が赤玉8個, 青玉4個となる事象\n",
    "\n",
    "とすると, 求めるものは\n",
    "\\begin{equation}\n",
    "\\dfrac{P(H_1\\cap A)}{P(A)} \n",
    "= \\dfrac{P(H_1\\cap A)}{P(H_1\\cap A)+P(H_2\\cap A)}\n",
    "\\end{equation}\n",
    "である. \n",
    "さらに条件付き確率\n",
    "$$\n",
    "P(H_1 | A) = \\dfrac{P(H_1\\cap A)}{P(A)} \n",
    "$$\n",
    "という記法を採用すると, 求めるものは"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "P(A|H_1)\n",
    "&= \\dfrac{P(H_1\\cap A)}{P(A)} \\\\\n",
    "&= \\dfrac{P(H_1\\cap A)}{P(H_1\\cap A)+P(H_2\\cap A)}\\\\\n",
    "&= \\dfrac{P(H_1)P(H_1| A)}{P(H_1)P(H_1| A)+P(H_2)P(H_2| A)}.\n",
    "\\end{align}\n",
    "ここで, ランダムに壺を選んだという仮定から, \n",
    "$$\\begin{split}\n",
    "P(H_1\\cap A)=\\dfrac12 {}_{12} C_8 (0.7)^8(0.3)^4,\\\\\n",
    "P(H_2\\cap A)=\\dfrac12 {}_{12} C_4 (0.7)^4(0.3)^8\n",
    "\\end{split}$$\n",
    "なので, "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "&\\dfrac{P(H_1)P(H_1| A)}{P(H_1)P(H_1| A)+P(H_2)P(H_2| A)}\\\\\n",
    "&= \\dfrac{\\dfrac12 {}_{12} C_8 (0.7)^8(0.3)^4}\n",
    "{\\dfrac12 {}_{12} C_8 (0.7)^8(0.3)^4 +\\dfrac12 {}_{12} C_4 (0.7)^4(0.3)^8}\\\\\n",
    "&= \\dfrac{7^4}{7^4+3^4}\\\\\n",
    "&= \\dfrac{2401}{2401+81}\n",
    "\\end{align}\n",
    "思ったより1 に近い.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一方で, $H_2$ が選ばれたと言える確率は, \n",
    "$$\n",
    "\\dfrac{81}{2401+81}\n",
    "$$\n",
    "である. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## この計算から感じること\n",
    "1. 結果から原因の確率分布がわかる\n",
    "1. ランダムに壺を選んだという事前の設定がある\n",
    "\n",
    "$P(H_i)$ を事前確率, $P(H_i|A)$ を事後確率という. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  問題9.1 \n",
    "## 「3つのつぼ $H_1$, $H_2$, $H_3$ があり, それぞれ $3:1$, $1:1$, $1:2$ の比率で, S, F と記された玉が入っている. ランダムに指定されたつぼから, 玉を復元抽出したとき, 玉はSであった. 」\n",
    "つまり, \n",
    "$$P(H_i)= 1/3,$$\n",
    "$$\n",
    "P(S|H_i) = \\dfrac{9}{12}, \\dfrac{6}{12}, \\dfrac{4}{12}.\n",
    "$$\n",
    "## (i) その玉がつぼ $H_i$, $(i=1,2,3)$ である事後確率 $P(H_i|S)$ を求めよ. \n",
    "$$\n",
    "P(H_i|S)\n",
    "= \\dfrac{P(H_i)P(S|H_i)}{P(S)}\n",
    "$$\n",
    "ここで, \n",
    "\\begin{align}\n",
    "P(S) \n",
    "&= \\sum_{i=1}^3 P(H_i)P(S|H_i)\\\\\n",
    "&= \\dfrac13 \\times \\left(\\dfrac{9}{12} +\\dfrac{6}{12}+ \\dfrac{4}{12}\\right)\\\\\n",
    "&= \\dfrac{19}{36}.\n",
    "\\end{align}\n",
    "一方で, \n",
    "$$\n",
    "P(H_i)P(S|H_i) = \\dfrac{9}{36}, \\dfrac{6}{36}, \\dfrac{4}{36}\n",
    "$$\n",
    "により, \n",
    "$$\n",
    "P(H_i|S) = \\dfrac{9}{19}, \\dfrac{6}{19}, \\dfrac{4}{19}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (ii) 指定された同じつぼから, もう一度玉を抜いたとき, 再び S だった. 回を添字で指定して S$_1$, S$_2$ と表そう. この時点でのつぼ $H_i$ の事後確率 $P(H_i|S_1,S_2)$ を求めよ. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "まず次を示す: \n",
    "$$\n",
    "P(H_i|S_1,S_2) = \\dfrac{P(H_i|S_1)P(S_2|H_i)}{P(S_2|S_1)}\n",
    "$$\n",
    "（証明）\n",
    "\\begin{align}\n",
    "&P(H_i|S_1,S_2) P(S_2|S_1)\n",
    "\\\\\n",
    "&= \\dfrac{P(H_i,S_1,S_2)}{P(S_1,S_2)}\\times \\dfrac{P(S_2,S_1)}{P(S_1)}\\\\\n",
    "&= \\dfrac{P(H_i,S_1,S_2)}{P(S_1)}.\n",
    "\\end{align}\n",
    "一方で, \n",
    "\\begin{align}\n",
    "&P(H_i|S_1)P(S_2|H_i)\\\\\n",
    "&= \\dfrac{P(H_i,S_1)}{P(S_1)}\\times \\dfrac{P(S_2,H_i)}{P(H_i)}\\\\\n",
    "&= \\dfrac{P(S_1|H_i)P(H_i) P(S_2|H_i)P(H_i)}{P(S_1)P(H_i)}\\\\\n",
    "&= \\dfrac{P(S_1,S_2|H_i)}{P(S_1)}\n",
    "\\end{align}\n",
    "である. 途中で $S_1$, $S_2$ が独立であることを用いた. \n",
    "\n",
    "（終）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "これより, \n",
    "\\begin{align}\n",
    "P(H_i|S_1,S_2) \n",
    "&= \\dfrac{P(H_i|S_1)P(S_2|H_i)}{P(S_2|S_1)}\\\\\n",
    "&= \\dfrac{P(H_i|S_1)P(S_2|H_i)}{\\sum_{j=1}^3P(H_i|S_1)P(S_2|H_i)}.\n",
    "\\end{align}\n",
    "となり, (i)での事後確率を事前確率として計算すれば良い. よって\n",
    "$$\n",
    "P(H_i|S_1,S_2) = \\dfrac{81}{133}, \\dfrac{36}{133}, \\dfrac{16}{133}.\n",
    "$$\n",
    "### この計算からわかること\n",
    "1. ある時点での事後確率を事前確率として用いて新しいデータに対して事後確率を求めることで, 全体の事後確率を求められる. \n",
    "1. データを増やしていくごとに「どちらのつぼが使われたか」を表す確率分布が変動する. \n",
    "##### データから確率分布を予測することをベイズ推論と呼ぶ. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (iii) 各 $H_i$ に対し, $S_1S_2$ を考える方法で事後確率を求めよ. \n",
    "事前確率 $P(H_i) =1/3$ および, \n",
    "$$\n",
    "\\begin{split}\n",
    "P(S_1,S_2|H_i) \n",
    "&= \\left(\\dfrac{9}{12}\\right)^2,\\left(\\dfrac{6}{12}\\right)^2,\\left(\\dfrac{4}{12}\\right)^2\n",
    "\\\\\n",
    "&= \\dfrac{81}{12^2},\\dfrac{36}{12^2},\\dfrac{16}{12^2}\n",
    "\\end{split}$$\n",
    "である. $81+36+15 =133$だから, \n",
    "$$\n",
    "P(H_i|S_1,S_2) = \\dfrac{81}{133}, \\dfrac{36}{133}, \\dfrac{16}{133}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 問9.2「問9.1(i)の事後確率のもとで, 次にSとなる確率 $P(S_2|S_1)$ を求めよ. 」\n",
    "\\begin{align}\n",
    "P(S_2|S_1)\n",
    "&= \\sum_{j=1}^3 P(S_2|H_i)P(H_i|S_1)\\\\\n",
    "&= \\dfrac{9}{12} \\times \\dfrac{9}{19} + \\dfrac{6}{12} \\times \\dfrac{6}{19} + \\dfrac{4}{12} \\times \\dfrac{4}{19}\\\\\n",
    "&= \\dfrac{133}{238}. \n",
    "\\end{align}\n",
    "### この計算からわかること\n",
    "1. 選んだつぼがどれか確定させず(確率分布のまま)に, 次に何が起こるかを確率的に記述できる. \n",
    "1. つぼがどれか選ぶことは, 母数の推定にそのまま相当する. \n",
    "#### 連続分布でも同じことができそう. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9.2 事前確率分布と事後確率分布\n",
    "新しく開発した薬がヒトに効く確率 $\\theta_0 \\in[0,1]$ を知りたい ⇦ 母数が連続になったよ！\n",
    "\n",
    "実験データZ:$n$人のうち$z$人には効いた. \n",
    "\n",
    "このようなことが起こる確率（**尤度**）は, \n",
    "$$\n",
    "f(z|\\theta) = {}_nC_z \\theta^z(1-\\theta)^{n-z}.\n",
    "$$\n",
    "である. \n",
    "$$\\begin{split}\n",
    "&\\dfrac{\\textrm{母数が$\\theta$であり, 実験データZが得られる確率}}{\\textrm{実験データZが得られる確率}}\\\\\n",
    "&=\\dfrac{\\textrm{母数が$\\theta$になる確率}\\times \\textrm{母数が$\\theta$のときに, 実験データZが得られる確率}}\n",
    "{\\textrm{実験データZが得られる確率}}\n",
    "\\end{split}$$\n",
    "という形で事後分布（$\\theta$の確率分布）を求めたい. これを$w'(\\theta|z)$と書く. \n",
    "\n",
    "とはいえ母数が$\\theta$になる確率は不明なので, 事前分布 $w(\\theta)$ は自分で決めてしまおう. \n",
    "\n",
    "実験データを得るまではまだ何もわからないので, \n",
    "$$\n",
    "w(\\theta) = 1 ,\\quad \\theta \\in[0,1]\n",
    "$$\n",
    "とする. このとき\n",
    "$$\n",
    "w'(\\theta|z) = \\dfrac{w(\\theta)f(z|\\theta)}{\\int_0^1 w(\\theta)f(z|\\theta) \\,d\\theta}\n",
    "$$\n",
    "となる. 分母は規格化定数. \n",
    "#### 問9.1(ii)のように逐次的に学習をしたい ⇨ 事前分布と事後分布は同じようになってほしい  ⇨ 事前分布はなんでもいいと言ったが, こうした方が便利\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "では, 実際に実験データ「$5$人のうち$4$人には効いた」が得られた場合を考えよう. \n",
    "\n",
    "$\\theta \\in [0,1],$ に対して, \n",
    "$$\n",
    "w(\\theta)f(4|\\theta) =5\\theta^4(1-\\theta)\n",
    "$$\n",
    "により, \n",
    "$$\n",
    "\\int_0^1 w(\\theta)f(4|\\theta) \\,d\\theta = 5B(4,2).\n",
    "$$\n",
    "よって事後分布は\n",
    "$$\n",
    "w'(\\theta|4)\n",
    "= \\dfrac{1}{B(4,2)} \\theta^4(1-\\theta)\n",
    "$$\n",
    "なるベータ分布となっている. \n",
    "\n",
    "また, 事後分布に対する $\\theta$ の期待値は\n",
    "$$\n",
    "\\int_0^1 \\theta w'(\\theta|4) \n",
    "=\\dfrac{1}{B(4,2)}\\int_0^1 \\theta\\theta^4(1-\\theta) \\,d\\theta = \\dfrac{B(5,2)}{B(4,2)}=\\dfrac54\n",
    "$$\n",
    "となり, 標本平均に一致する. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ここからおまけ. 記法がやや異なる. \n",
    "### [指数型分布族]\n",
    "実数値確率変数 $X$ が未知パラメータ $\\theta$ をもつ(条件付き)確率分布 $p(x|\\theta)$に従うとする\n",
    "(つまり, 観測モデルを $p(x|\\theta)$ と定める). \n",
    "このとき, $X$ が指数型分布族に属するとは, ある値域の次元が等しい関数 $\\phi$, $g$ により, \n",
    "\\begin{align*}\n",
    "p(x| \\theta) \n",
    "&= f(x)\\exp\\langle \\phi(\\theta) , g(x)\\rangle ,\\\\\n",
    "f(x) \n",
    "&= \\left(\\int \\exp\\langle \\phi(\\theta) , g(x)\\rangle\\, d\\theta \\right)^{-1}\n",
    "\\end{align*}\n",
    "が成立することをいう. \n",
    "\n",
    "### [指数型分布族に対する事前分布]\n",
    "ハイパーパラメータ $v \\in \\mathbb{R}^k$ に対して\n",
    "\\begin{align*}\n",
    "p(\\theta| v) \n",
    "&= z(v)\\exp\\langle \\phi(\\theta) , v\\rangle,\\\\\n",
    "z(v)\n",
    "&= \\left(\\int \\exp\\langle \\phi(\\theta) ,v\\rangle\\, d\\theta \\ \\right)^{-1}\n",
    "\\end{align*}\n",
    "を事前分布として定める. \n",
    "このとき, \n",
    "$\\hat{v}= \\sum_{n=1}^N g(X_n)$ により, 事後分布 $p(\\theta|\\bar{X})$ は\n",
    "$$\n",
    "p(\\theta|\\bar{X}) = p(\\theta|v+\\hat{v})\n",
    "$$\n",
    "と表せる. すなわち, 事前分布と事後分布における $\\phi$ は等しい. \n",
    "このように, 分布のうち母数に依存する部分である $\\phi$ がベイズ推論によって不変な事前分布を\n",
    "共役事前分布という. \n",
    "\n",
    "(証明)\n",
    "モデルの設計から, $p(\\bar{X},\\theta|v)=p(\\bar{X}|\\theta)p(\\theta|v)$ である.\n",
    "条件付き確率の定義より, \n",
    "\\begin{align*}\n",
    "p(\\theta|\\bar{X},v)\n",
    "&= \\dfrac{p(\\bar{X}|\\theta,v)}{p(X_\\ast|v)}\\\\\n",
    "&= \\dfrac{p(\\bar{X}|\\theta)p(\\theta|v)}{p(\\bar{X}|v)}\\\\\n",
    "&= \\prod_{n=1}^N p(X_n|\\theta)p(\\theta|v)/p(\\bar{X}|v)\\\\\\\n",
    "&= \\prod_{n=1}^N f(X_n) z(v)p(\\bar{X}|v)^{-1} \\prod_{n=1}^N\\exp\\langle \\phi(\\theta) , g(X_n)\\rangle\n",
    "\\exp\\langle \\phi(\\theta) , v\\rangle\\\\\n",
    "&=z(\\hat{v})\\exp\\langle \\phi(\\theta) , v+\\hat{v}\\rangle.\n",
    "\\end{align*}\n",
    "（終）\n",
    "## 共役事前分布を用いた予測\n",
    "事前分布または事後分布が与えられたとき, 次に観測される $X_\\ast$ の従う確率分布を予測分布と呼ぶ. \n",
    "\n",
    "指数型分布族に対する共役分布については, 次のような公式が与えられる; \n",
    "\n",
    "確率変数 $X$ が指数型分布族に属するとし, 共役事前分布 $p(\\theta|v)$ をと定められたものとする. \n",
    "このとき, 予測分布に対して次が成立する. \n",
    "\\begin{align*}\n",
    "p(X_\\ast|v)\n",
    " &= \\dfrac{f(X_\\ast)z(v)}{z(v+g(X_\\ast))},\\\\\n",
    "p(X_\\ast|\\bar{X},v)\n",
    "&= \\dfrac{f(X_\\ast) z(v+\\hat{v})}\n",
    "{z(v+\\hat{v}+g(X_\\ast))}.\n",
    "\\end{align*}\n",
    "ここで, 1つ目の式を事前分布に対する予測分布, 2つ目の式を事後分布に対する予測分布と呼ぼう. \n",
    "\n",
    "(証明)\n",
    "\n",
    "事前分布 $p(\\theta| v)$ を用いて予測分布 $p(X_\\ast|v)$ を求めよう. モデルの設計により, \n",
    "$$\n",
    "p(\\theta|X_\\ast,v) =\\dfrac{p(X_\\ast|\\theta)p(\\theta|v)}{p(X_\\ast|v)}\n",
    "$$\n",
    "なる関係が成り立つ. これを変形すると, \n",
    "\\begin{align*}\n",
    "p(X_\\ast|v)\n",
    " &=\\dfrac{p(X_\\ast|\\theta)p(\\theta|v)}{p(\\theta|X_\\ast)}\\\\\n",
    " & = \\dfrac{f(X_\\ast)\\exp\\langle \\phi(\\theta),g(X_\\ast)\\rangle z(v)\\exp\\langle \\phi(\\theta) , v\\rangle\n",
    "  }{z(v+g(X_\\ast))\\exp\\langle \\phi(\\theta),g(X_\\ast)+v\\rangle}\\\\\n",
    " &= \\dfrac{f(X_\\ast)z(v)}{z(v+g(X_\\ast))}.\n",
    "\\end{align*}\n",
    "（注：ここでも指数法則が本質的！！）事後分布 $p(\\theta|\\bar{X})$ を用いて予測分布 $p(X_\\ast|\\bar{X})$ を求めよう. \n",
    "再びモデルの設計により, \n",
    "$$\n",
    "p(\\theta|\\bar{X},X_\\ast,v)=\\dfrac{p(X_\\ast|\\theta)p(\\theta|\\bar{X},v)}{p(X_\\ast|\\bar{X},v)}\n",
    "$$\n",
    "が成り立つ. 同様の計算により, \n",
    "\\begin{align*}\n",
    "p(X_\\ast|\\bar{X},v)\n",
    "&= \\dfrac{p(X_\\ast|\\theta)p(\\theta|\\bar{X},v)}{p(\\theta|\\bar{X},X_\\ast,v)}\\\\\n",
    "&= \\dfrac{f(X_\\ast)\\exp\\langle\\phi(\\theta),g(X_\\ast)\\rangle p(\\theta|v+\\hat{v})}\n",
    "{z(v+\\hat{v}+g(X_\\ast))\\exp\\{\\langle\\phi(\\theta),g(X_\\ast)+v+\\hat{v}\\rangle\\}}\\\\\n",
    "&= \\dfrac{f(X_\\ast)\\exp\\langle\\phi(\\theta),g(X_\\ast)\\rangle z(v+\\hat{v})\\exp\\langle\\phi(\\theta),v+\\hat{v}\\rangle}\n",
    "{z(v+\\hat{v}+g(X_\\ast))\\exp\\langle\\phi(\\theta),g(X_\\ast)+v+\\hat{v}\\rangle}\\\\\n",
    "&= \\dfrac{f(X_\\ast) z(v+\\hat{v})}\n",
    "{z(v+\\hat{v}+g(X_\\ast))}.\n",
    "\\end{align*}\n",
    "(終わり)\n",
    "# ここまでおまけ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 問9.4 ある量 $X$ の母集団分布が $N(\\theta,1/4)$ のとき, $\\theta$ の事前確率分布として平均$5$, 分散 $1/2$の正規分布 $N(5,1/2)$ を仮定した. $X$ の観測値が $X=6.5$ のとき, $\\theta$ の事後確率分布を求めよ. \n",
    "### まずは一般に「平均が未知のときのガウス分布」を先ほどの指数型分布族についての議論に持ち込む. \n",
    "分散の逆数である精度 $\\lambda>0$ は既知として, 平均 $\\mu$ が未知の一次元ガウス分布を考える. \n",
    "$$\n",
    "N(x|\\mu,\\lambda^{-1})\n",
    "= \\sqrt{\\dfrac{\\lambda}{2\\pi}}e^{-\\frac{\\lambda}{2}x^2}\\exp\\{-\\dfrac{\\lambda}{2}\\mu^2 +\\lambda x\\mu\\}\n",
    "$$\n",
    "により, \n",
    "\\begin{align*}\n",
    "f(x)&=\\sqrt{\\dfrac{\\lambda}{2\\pi}}e^{-\\frac{\\lambda}{2}x^2}, \\\\\n",
    "\\phi_1(\\mu)&= -\\dfrac{1}{2}\\mu^2,\\quad \\phi_2(\\mu) =\\mu,\\\\\n",
    "g_1(x) &= \\lambda, \\quad g_2(x)=\\lambda x\n",
    "\\end{align*}\n",
    "なる指数型分布族である. したがって, 事前分布はパラメータ $v=(v_1,v_2)$ より, \n",
    "$$\n",
    "p(\\mu|v ) \n",
    "=\\sqrt{\\dfrac{v_1}{2\\pi}}e^{-\\frac{v_2^2}{2v_1}}\\exp\\{-\\dfrac{v_1}{2}\\mu^2 + v_2\\mu\\}\n",
    "$$\n",
    "と置けば事後分布は \n",
    "$$\n",
    "p(\\mu|\\bar{X})=p(\\mu|v+\\hat{v})\n",
    "$$\n",
    "である. 事前分布は平均が$\\frac{v_2}{v_1}$, 精度が$v_1$である正規分布となっているので, \n",
    "$a=\\frac{v_2}{v_1}$, $b = v_1$ と表すと, \n",
    "事前分布は平均$a$, 分散$b$ の正規分布であり, \n",
    "\\begin{align*}\n",
    "\\hat{a} &= \\dfrac{(v+\\hat{v})_2}{(v+\\hat{v})_1}\\\\\n",
    "&= \\dfrac{ab+\\lambda x}{\\hat{b}},\\\\\n",
    "\\hat{b} &= (v+\\hat{v})_1\\\\\n",
    "&= b+\\lambda\n",
    "\\end{align*}\n",
    "と置くと, 事後分布は平均$\\hat{a}$, 分散$\\hat{b}$ の正規分布である. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\lambda = 4$, $a=5$, $b=2$を代入すれば,\n",
    "$$\n",
    "\\hat{b} = 2+4=6,\n",
    "$$\n",
    "$$\n",
    "\\hat{a}= \\dfrac{5\\times 2 + 4\\times 6.5}{8}=6\n",
    "$$\n",
    "となり, 事後分布は $N(6,1/6)$ である. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9.3 ベイズ推定 & 9.4 統計的決定理論\n",
    "さっきまでの議論を一旦忘れる。\n",
    "\n",
    "データから真の母数 $\\theta_0$ の推定を行う. 推定とは実現値から母数のよい近似値をみつけること. \n",
    "点推定とは違うアプローチを考えよう. \n",
    "\n",
    "尤度関数 $f(z|\\theta)$ は既知であるものとする. \n",
    "\n",
    "すなわち, 母数が $\\theta$ であったならば実現値 $z$ が得られる確率が $f(z|\\theta)$ という形で得られることを仮定する. \n",
    "\n",
    "ここでの目標は, 実現値 $z$ による推定値 $a(z)$ を損失関数の最小化によって構成することだ. \n",
    "\n",
    "損失関数の具体的な形は後で与えるとして, $L(\\theta-a(z))$ の期待値を下げることを考えていこう. \n",
    "\n",
    "まず, $\\theta$ を固定する. このとき$z$ の実現確率は $f(z|\\theta)$ なので, その期待値は, \n",
    "$$\n",
    "\\int f(z|\\theta)L(\\theta-a(z))\\,dz.\n",
    "$$\n",
    "母数 $\\theta_0$ は未知であるから, 何らかの意味での $\\theta$ についての加重平均を取る必要がある. \n",
    "これを事前確率$w(\\theta)$と呼ぶことにすると, この密度関数による期待値\n",
    "$$\n",
    "\\int w(\\theta)\\int f(z|\\theta)L(\\theta-a(z))\\,dzd\\theta.\n",
    "$$\n",
    "を最小化すればよいことが分かる. \n",
    "積分の順序を入れかえることで, 各$z$ について\n",
    "\n",
    "$$\n",
    "\\int w'(\\theta|z)L(\\theta-a(z))\\,d\\theta\n",
    "$$\n",
    "を最小化するものが求めるものだとわかる. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9.5 ベイズ判別\n",
    "有限個のクラスへの分類を行う. つまり母数は離散値(つぼの例に帰ってきた). \n",
    "\n",
    "統計的決定理論において, 損失関数を正解なら0, 不正解なら1を与えるように設定する."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 問9.7 アイリスデータのベイズ判別をせよ\n",
    "全データ使います. \n",
    "\n",
    "平均の推定値はそれぞれの標本平均, 分散の推定値は全体の標本分散に等しい.  \n",
    "（仮定なのか結論なのか分からなかった）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal length (cm)</th>\n",
       "      <th>sepal width (cm)</th>\n",
       "      <th>petal length (cm)</th>\n",
       "      <th>petal width (cm)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   sepal length (cm)  sepal width (cm)  petal length (cm)  petal width (cm)\n",
       "0                5.1               3.5                1.4               0.2\n",
       "1                4.9               3.0                1.4               0.2\n",
       "2                4.7               3.2                1.3               0.2\n",
       "3                4.6               3.1                1.5               0.2\n",
       "4                5.0               3.6                1.4               0.2"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn import datasets\n",
    "iris = datasets.load_iris()\n",
    "df = pd.DataFrame(iris.data, columns = iris.feature_names)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sl</th>\n",
       "      <th>sw</th>\n",
       "      <th>pl</th>\n",
       "      <th>pw</th>\n",
       "      <th>answer</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    sl   sw   pl   pw  answer\n",
       "0  5.1  3.5  1.4  0.2  setosa\n",
       "1  4.9  3.0  1.4  0.2  setosa\n",
       "2  4.7  3.2  1.3  0.2  setosa\n",
       "3  4.6  3.1  1.5  0.2  setosa\n",
       "4  5.0  3.6  1.4  0.2  setosa"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.columns = ['sl', 'sw', 'pl', 'pw']\n",
    "df['answer'] = iris.target_names[iris.target]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(sl    6.588\n",
       " sw    2.974\n",
       " pl    5.552\n",
       " pw    2.026\n",
       " dtype: float64, sl    5.006\n",
       " sw    3.418\n",
       " pl    1.464\n",
       " pw    0.244\n",
       " dtype: float64, sl    5.936\n",
       " sw    2.770\n",
       " pl    4.260\n",
       " pw    1.326\n",
       " dtype: float64)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_vi = df[df['answer']=='virginica']\n",
    "df_se = df[df['answer']=='setosa']\n",
    "df_ve = df[df['answer']=='versicolor']\n",
    "\n",
    "mu_vi = df_vi.mean()\n",
    "mu_se = df_se.mean()\n",
    "mu_ve = df_ve.mean()\n",
    "mu_vi, mu_se, mu_ve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.00000000e+00, -1.11022302e-16,  1.24900090e-16,\n",
       "         1.11022302e-16],\n",
       "       [ 0.00000000e+00,  1.00000000e+00, -3.46944695e-17,\n",
       "         1.38777878e-16],\n",
       "       [-2.22044605e-16, -1.11022302e-16,  1.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [-1.38777878e-17, -1.38777878e-17,  0.00000000e+00,\n",
       "         1.00000000e+00]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V = (df_vi.cov() + df_se.cov() + df_ve.cov() ).values\n",
    "V_inv = np.linalg.inv(V)\n",
    "\n",
    "np.dot(V,V_inv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sl</th>\n",
       "      <th>sw</th>\n",
       "      <th>pl</th>\n",
       "      <th>pw</th>\n",
       "      <th>answer</th>\n",
       "      <th>pred_vi</th>\n",
       "      <th>pred_se</th>\n",
       "      <th>pred_ve</th>\n",
       "      <th>pred</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-3.517584</td>\n",
       "      <td>60.615621</td>\n",
       "      <td>27.609643</td>\n",
       "      <td>pred_se</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-6.330993</td>\n",
       "      <td>49.630676</td>\n",
       "      <td>23.197758</td>\n",
       "      <td>pred_se</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-8.392017</td>\n",
       "      <td>50.724460</td>\n",
       "      <td>21.678906</td>\n",
       "      <td>pred_se</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-7.744714</td>\n",
       "      <td>45.428420</td>\n",
       "      <td>20.872961</td>\n",
       "      <td>pred_se</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-4.120693</td>\n",
       "      <td>60.622414</td>\n",
       "      <td>27.026360</td>\n",
       "      <td>pred_se</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    sl   sw   pl   pw  answer   pred_vi    pred_se    pred_ve     pred\n",
       "0  5.1  3.5  1.4  0.2  setosa -3.517584  60.615621  27.609643  pred_se\n",
       "1  4.9  3.0  1.4  0.2  setosa -6.330993  49.630676  23.197758  pred_se\n",
       "2  4.7  3.2  1.3  0.2  setosa -8.392017  50.724460  21.678906  pred_se\n",
       "3  4.6  3.1  1.5  0.2  setosa -7.744714  45.428420  20.872961  pred_se\n",
       "4  5.0  3.6  1.4  0.2  setosa -4.120693  60.622414  27.026360  pred_se"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['pred_vi'] = 2*(df['sl']*np.dot(V_inv,mu_vi)[0] + df['sw']*np.dot(V_inv,mu_vi)[1] + df['pl']*np.dot(V_inv,mu_vi)[2] + df['pw']*np.dot(V_inv,mu_vi)[3] ) - np.dot(mu_vi,np.dot(V_inv,mu_vi))\n",
    "df['pred_se'] = 2*(df['sl']*np.dot(V_inv,mu_se)[0] + df['sw']*np.dot(V_inv,mu_se)[1] + df['pl']*np.dot(V_inv,mu_se)[2] + df['pw']*np.dot(V_inv,mu_se)[3] ) - np.dot(mu_se,np.dot(V_inv,mu_se))\n",
    "df['pred_ve'] = 2*(df['sl']*np.dot(V_inv,mu_ve)[0] + df['sw']*np.dot(V_inv,mu_ve)[1] + df['pl']*np.dot(V_inv,mu_ve)[2] + df['pw']*np.dot(V_inv,mu_ve)[3] ) - np.dot(mu_ve,np.dot(V_inv,mu_ve))\n",
    "df['pred'] = df.loc[:,['pred_vi','pred_se','pred_ve']].idxmax(axis=1)\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sl</th>\n",
       "      <th>sw</th>\n",
       "      <th>pl</th>\n",
       "      <th>pw</th>\n",
       "      <th>answer</th>\n",
       "      <th>pred_vi</th>\n",
       "      <th>pred_se</th>\n",
       "      <th>pred_ve</th>\n",
       "      <th>pred</th>\n",
       "      <th>is_correct</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-3.517584</td>\n",
       "      <td>60.615621</td>\n",
       "      <td>27.609643</td>\n",
       "      <td>pred_se</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-6.330993</td>\n",
       "      <td>49.630676</td>\n",
       "      <td>23.197758</td>\n",
       "      <td>pred_se</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-8.392017</td>\n",
       "      <td>50.724460</td>\n",
       "      <td>21.678906</td>\n",
       "      <td>pred_se</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-7.744714</td>\n",
       "      <td>45.428420</td>\n",
       "      <td>20.872961</td>\n",
       "      <td>pred_se</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>setosa</td>\n",
       "      <td>-4.120693</td>\n",
       "      <td>60.622414</td>\n",
       "      <td>27.026360</td>\n",
       "      <td>pred_se</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    sl   sw   pl   pw  answer   pred_vi    pred_se    pred_ve     pred  \\\n",
       "0  5.1  3.5  1.4  0.2  setosa -3.517584  60.615621  27.609643  pred_se   \n",
       "1  4.9  3.0  1.4  0.2  setosa -6.330993  49.630676  23.197758  pred_se   \n",
       "2  4.7  3.2  1.3  0.2  setosa -8.392017  50.724460  21.678906  pred_se   \n",
       "3  4.6  3.1  1.5  0.2  setosa -7.744714  45.428420  20.872961  pred_se   \n",
       "4  5.0  3.6  1.4  0.2  setosa -4.120693  60.622414  27.026360  pred_se   \n",
       "\n",
       "   is_correct  \n",
       "0        True  \n",
       "1        True  \n",
       "2        True  \n",
       "3        True  \n",
       "4        True  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['is_correct'] = (df['answer'].str[:2]==df['pred'].str[5:])\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True     147\n",
       "False      3\n",
       "Name: is_correct, dtype: int64"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['is_correct'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}