3_Estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 估计"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 骰子问题"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4 0.0\n",
      "6 0.392156862745\n",
      "8 0.294117647059\n",
      "12 0.196078431373\n",
      "20 0.117647058824\n",
      "==================\n",
      "4 0.0\n",
      "6 0.0\n",
      "8 0.943248453672\n",
      "12 0.0552061280613\n",
      "20 0.0015454182665\n"
     ]
    }
   ],
   "source": [
    "from thinkbayes import Suite\n",
    "class Dice(Suite):\n",
    "    def Likelihood(self, data, hypo):\n",
    "        if hypo < data:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1.0 / hypo\n",
    "\n",
    "suite = Dice([4, 6, 8, 12, 20])\n",
    "# 如果转动得到6\n",
    "suite.Update(6)\n",
    "suite.Print()\n",
    "print(\"==================\")\n",
    "# 多摇几次\n",
    "for roll in [6, 8, 7, 7, 5, 4]:\n",
    "    suite.Update(roll)\n",
    "suite.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 火车头问题\n",
    "铁路上以 $1$ 到 $N$ 命名火车头。有一天你看到一个标号 $60$ 的火车头，请估计铁路上有多少火车头？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "151.849587959\n",
      "164.305586423\n",
      "171.338181092\n"
     ]
    }
   ],
   "source": [
    "class Train(Suite):\n",
    "    def Likelihood(self, data, hypo):\n",
    "        if hypo < data:\n",
    "            return 0\n",
    "        else:\n",
    "            return 1.0 / hypo\n",
    "        \n",
    "def Mean(suite):\n",
    "    total = 0\n",
    "    for hypo, prob in suite.Items():\n",
    "        total += hypo * prob\n",
    "    return total\n",
    "\n",
    "for num in [500, 1000, 2000]:\n",
    "    hypos = xrange(1, num+1)\n",
    "    suite = Train(hypos)\n",
    "    for data in [60, 30, 90]:\n",
    "        suite.Update(data)\n",
    "    # 使用后验概率的平均值来作为估计值会减少从长远来看的均方差\n",
    "    print(Mean(suite))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在数学上，幂律表示公司的数量与公司规模成反比\n",
    "$$PMF(x) \\propto \\left(\\frac{1}{x}\\right)^\\alpha $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(130.70846986256004, (91, 235))\n",
      "(133.27523137503127, (91, 242))\n",
      "(133.99746308073063, (91, 243))\n"
     ]
    }
   ],
   "source": [
    "from thinkbayes import Pmf\n",
    "from thinkbayes import Percentile\n",
    "class InvTrain(Train):\n",
    "    # 构造一个服从幂律分布的先验\n",
    "    def __init__(self, hypos, alpha=1.0):\n",
    "        Pmf.__init__(self)\n",
    "        for hypo in hypos:\n",
    "            self.Set(hypo, hypo**(-alpha))\n",
    "            \n",
    "for num in [500, 1000, 2000]:\n",
    "    hypos = xrange(1, num+1)\n",
    "    suite = InvTrain(hypos)\n",
    "    for data in [60, 30, 90]:\n",
    "        suite.Update(data)\n",
    "    interval = Percentile(suite, 5), Percentile(suite, 95)\n",
    "    print(Mean(suite), interval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "贝叶斯当中，有两种途径选择先验分布。\n",
    "- 选择最能代表问题相关背景资料的先验概率，在这种情况下，先验被认为是“**信息**”。问题是，人们可能会使用不同的背景信息（或者进行不同的诠释），所以基于信息的先验往往显得主观。\n",
    "- 另一种方法是所谓的“**无信息参考的先验**”，其目的是为了让数据说话，越没有约束越好。\n",
    "\n",
    "有大量的数据时，先验的选择不是特别关键；如果你没有太多的参考数据，那么采用相关的背景信息就有很大区别了。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}