SML p set 2 solution

sml-ps2-solution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "theta = 0.7\n",
    "coins = np.random.uniform(size=(3, 1000)) <= theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True,  True, False, ..., False, False,  True],\n",
       "       [False, False,  True, ...,  True,  True,  True],\n",
       "       [False,  True,  True, ..., False, False,  True]], dtype=bool)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# MLE estimator of theta\n",
    "theta_hat = coins.sum(axis=0)/3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "theoretical bias 0.0 observed bias 0.00133333333333\n",
      "theoretical variance 0.07 observed variance 0.0741315555556\n",
      "theoretical mse 0.07 observed mse 0.0741333333333\n"
     ]
    }
   ],
   "source": [
    "# E(theta_hat)\n",
    "e_theta_hat = theta_hat.mean()\n",
    "\n",
    "# bias\n",
    "print 'theoretical bias', 0.,\n",
    "print 'observed bias', (theta_hat - theta).mean()\n",
    "\n",
    "# variance\n",
    "print 'theoretical variance', theta*(1-theta)/3,\n",
    "print 'observed variance', ((theta_hat - e_theta_hat)**2).mean()\n",
    "\n",
    "# mean squared error\n",
    "print 'theoretical mse', theta*(1-theta)/3,\n",
    "print 'observed mse', ((theta_hat - theta)**2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# the sample (a, b) from problem 2\n",
    "ab = [(1., 1.), (2., 4.), (20., 20.), (25., 10.)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEaCAYAAABEsMO+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYVOXZx/EvVQJIEwUEdLGgYEcFpK6KIqhYECsaBI3t\ntaDGWJKAeWOJb4y9K2CJQAJREUH6KqCiGKSDgKCAAiJKC5Gy+/5xz7izy8zO2dlTZ36f65qL2Z2z\nc+7Vufc552k3iIiIiIiIiIiIiIiIiIiIiIiIiIiIiIj8okrQAYivCrD/53MCjkMkagpQ7viuctAB\nCI2BMcBaoBA4KM3xecA0YDuwGDi9HOcqij0q6kVgCbAH+HWaY/cBhgCbge+AgS6cXwTgbGAG8CP2\n2XoJqF3G8XkEnztxV2H5PqCMY3I+d9RABa8QGAf0dnj8cOBzoAFwHzAKaOhNaCl9AdwI/Jv0STsY\nOBRreE8F7gK6exmc5Iw6wJ+AJkAroCnwf2UcH4bcAagP3AssoOz8GYxyRyrgbmA5sAVYCJxfgfeq\nSvo7qJbAf4FaCd/7ALjO4TmmAQ8Cs7CrsrexZMnUdOxKsCxrgW4JX9+P/aGQ3OZm7sRdAMxL8VqY\ncud54IbYe/Yv47iczx3dQVXMcqATdiV3P/AG0Cj2Wies6yHVo0MG5zsK+ArrooibG/u+E5WwBuVq\n7KpzN/Bkwus/lRHvXRnEWz92nrkJ35tXjngle3mRO12xu5JkwpI7bYE2WCNVFuWOuG4O0CvDn3Vy\nB3Ul8HGp7/0ZGOrwHPGrwLhWwM9Y8mUi3R1Uc+x3qp7wvTOAlRmeT7JXRXIH7HO1CTgsxethyJ0q\nwGdYIxV/z1R3UModdAdVUVdhiRW/Ujoa2M/D823DrjgT1cO6SZxanfD8G6Aa3vXDb4v9mxhzXWCr\nR+eT6HAzd9oDf8fGcZenOCYMuXMjdhf0acL3UjVwyh3UQFXEwdhstpuwQdf6WPdC/APXGfswpXp0\nzOCcC4FDKDlT6bjY9506qNTzXcDG2Nfbyoj37gzijc+uOr5UvKm6YSQ3uJk7JwDvAP2wO5JUwpA7\np2HjZN/FHh2ARynZVRin3JEKaQ3swAZfq2B907soe9AzlRpY4hTG3q9GGcd+jM1UqgFciH2Q41ee\neZTdTViAXQW2AmoC/8T6/surWuz8M4FrYs9TXQk+FDtvvdh5vwPOzOCckj3cyp2jgfVAH4fHB507\ndYEDYo9GWP7cBuyb4njljlTIn4EfgO+xK6F0s3JSKYw99iT8G/dc7BF3cOw8/8HWcpyW8FpnbCA4\n1QLsacADFM9Eege7gi2vglLxFgJdYq9dQcmrvOrAK7HzrcMSUsSN3BmCTVZIvFuZn/B6GHOn9Hsm\n/s7KnQwMwa5S5pdxzJPAMmzGyQkJ36+HrTVYDCzC+orFO/cB1wYdhKTVHPvjtBD7g3RLkmPysT9M\nc2KP3/sVXI5S7kRUZ6zRSdVA9cQWmgK0Az5JeO1Viq8QqmK3uCK5rjHFYwu1gaVYF06ifGyHERFJ\nI4/UDdTzwCUJXy/B+lfrYrfMIlK2t9l725184F3/QxEJDzdm8TWl5PTLNUAzoAXWvzwU2xLnJWxw\nUUSK5WE9FLNKfb8Im+U1F+uhaO1vWCLBc2uaeekZXEVYl14b4NnYv9tJMlX50EMPjW/CqIceYX+k\nWmOTqdrYGO2tFK97ifs3NlZ1HPAUdpe1F+WPHhF5ZJQ7bjRQa7FEimsW+96a2OOz2PdHYQ1VCStW\nrKCoqCjwx6BBgwKPIWyxKI6SD2zjTrdUA0ZjU5WTNT5bsdlmAONjx+81a0z5oziiEEemueNGAzWG\n4u1u2mN7Uq3HpkWuxtY6gG16WJ5FcSLZqhI2fXgR8HiKYxpR3DPRNvZ8k/ehiYRHVQfHDMc2YWyI\nNTiDsKs5gBew/vGe2C3cdmzRXdzN2BYk1YEVpV4TyVUdgb7YtjfxAnj3UrxI9AXgImzH693YndSl\nPscoEjgnDdRlDo75nxTfnwuc7Dyc4OTn5wcdwi/CEovi8MwM0vdePBN7REJY/h8pjpLCEkemMt3F\n2k1FsT5KkVCrVKkShCNnEil/JPQyzR1tFisiIqGkBkpEREJJDZSIiISSGigREQklNVAiIhJKaqBE\nRCSU1ECJiEgoOWmgKlKwEKxC5RxUOkBERMrBSQM1FDirjNd7AocBhwO/oWSJZbCdmhdhO9qKiIg4\n4mSro+lYzZpUemGVc8Fq2tTDNrpcj+1s3hN4ALg94yhFJHK+/x6WLIFNm+Dnn6FWLWjSBA4/HPbd\nN+joJAqcNFDpJCtY2BRroB4DfgvUceE8ItmiOfAacADWs/Ai1k1e2pNAD2yz2H4UbywbSkVF8OGH\n8Oab8P77sGULtGoFDRvCPvvAtm3w7bewbBkccgiceir07g2dO0OVKkFHL2HkRgMFe++xVAk4B9iA\nJVV+WT88ePDgX57n5+dHfoNDyQ4FBQUUFBR48da7gIHAF1jRws+BScDihGMSu87bYV3n7b0IpqKK\nimDUKBg82J737w+33QZHHgmVkuy+tns3fPEFTJwIAwfCDz/A9dfDddfBfvv5Hr6EmNPN+/KwSQ7H\nJHnteaAAGBH7egnWIN0CXImVC6iB3UWNprh2VJw2u5RI8HCz2LexqrlTEr73PDANGBn7eglW9mZ9\nqZ8NNH+WLrUGaccOeOghOPPM5I1SWf79b3jqKRgzBm64Ae66C+qozyWrBLlZbLKCheuw+jbNgRZY\nLZup7N04ieS6PGzm66xS30/Wdd7Mp5gcefFF6NgRLrsMZs+G7t3L3zgBtGkDQ4fae6xZA0ccYV/r\nulW8LliYSB83kZJqA6Owma7bkrxe+s99KHJo92645RaYNg0++ghatkz/M060aAHDhllDdf318MYb\n8NJLNl4lucnrgoVxH8QeImKqYV3eb2BdfKWtxXog4prFvrcXP8dwd+60O6Zt2+CTT6BuXffPcdJJ\n9t6PPQbt2sFf/wpXXZXZ3ZkEw63x2zD8L9cYlESCi2NQlbClGT9gkyWS6Yld+PXEus4fJ/kkCd/y\nZ9cum3VXuTKMHGkz87w2b541iCedBM89BzVren9OcZ8KFopER0egL3AqNst1Djad/LrYA6zr/Cus\n6/wF4Eb/wyxWVAS/+Q3s2QP//Kc/jRPAscfCp59CYSGccgqsWuXPeSUcdAcl4lAul3z/059g7Fgb\nd6pVy/PT7aWoCJ58Eh5+GP7xD1s7JdGRae6EIdnUQEkk5GoD9d57tkZp9mxo3NjTU6U1cSL07WuN\n1aWXBhuLOKcGSsRjudhArVplExX+9S+bUh4G8+bBOefYYuDbtYFaJKiBEvFYrjVQe/bAaafB2Wfb\n4tkwWb3aFgVfcAE88IBm+IWdJkmIiKsee8zGfu64I+hI9ta8OUyfbl1+t9xikygk+4ThukN3UBIJ\nuXQHtWyZzZr77DNbQBtWmzdDjx5w1FHwwgs2BV7Cx8s7qEwLFjbH9hJbCCzA9uYTkZArKoKbboJ7\n7gl34wS2UHjCBCvrcd11upPKNl4WLIzv2HwUtsDwJqBVxpGKiC9GjoT1663rLAr23RfGjYPFi+HG\nG7WHXzZx0kBNB34s4/VUBQvXYeUEwPYZWwwcmFmYIuKHHTtsQsSzz0K1aumPD4t4IzVnjs3sUyOV\nHdzosXWy63IeyXdsFpEQefxxaNs2PFPKy6NOHSuUWFAAgwYFHY24wauChYnXL+l2bFbBQgklDwsW\nhtKGDfDoo/Dxx0FHkrn69W1MqksXG58K4wxEcc6rgoXxwmrVgLHAeGyzy2Q0i08iIdtn8d1+u20I\n+9RTrrxdoFavtu2QBg+Gfv2CjkYyzR037qDGYLsuj6C4YOH6WDCvAItI3TiJSAisW2e1mBYsCDoS\ndzRvbt19+flWRv7cc4OOSDLhZAxqOPARcAQ21tQfZ7suJ9uxuazZgCK5JN3yjXxgM8W583svg3nk\nEbjySjgwi6YxHXmklZHv39/qS0n0hKG7Ql18Egkud/F1xsZkXyN513k+cDs2S7YsFc6fDRvsj/mC\nBdnVQMWNG2eN1AcfWDl58Z+2OhKJlnTLN8CnC8inn4aLL87OxgmgZ0948EHbcWLduqCjkfJwaxaf\niLirCOiA7c6yFrgTG8911fbtVql25ky33zlc+veHb76xXdALCqB27aAjEifUQImE07+x7cL+g1Xb\nfRtomezAiizTGDYMOnWClknfObsMGmSN1GWXwdtvQ5UqQUeUvdxaoqExKBGHPJhmnkfq5RulrQRO\nBDaV+n7G+VNYaGMyw4ZFc2FuJnbtsi6/I46w6fQq0+EPjUGJZJdGFCd029jz0o1ThUycaFsEdejg\n5ruGW7VqMGqUdfM98UTQ0Ug66uITCcZwbEF7Q2z5xiBsYTvYco2LgBuA3Vg3n+sFzp991jZXzbW7\niLp1YexYa5hbtIDzzgs6IkklDB9NdfFJJGTTThKrVsGJJ9qOCzVruh9UFHz2mXX3TZgAbdoEHU12\nUxefiDj24ou2MDdXGyeAk0+G55+HXr1gzZqgo5FkvCxYCLZzxJLYa7/LMEYRcdHu3fDqq3DttUFH\nErzeva3u1bnnwrakW1lLkLwsWFgFeDr2s62By1DBQpHATZxoe9UddVTQkYTDb39rXXxXXAF79gQd\njSTyqmBhY2zm0XJgFVZddwSg4UiRgA0ZYgtXxVSqZIuVt26F36mfJ1S8KljYFKuem+z7IhKQjRth\n8mS45JKgIwmX6tVt+vmYMTY+J+HgVcHC8v1wpcEJX+XHHiJBK4g9ssfIkXD22TbVWkpq0ADee8/q\nSB16KJx+etARiRsN1FpsS5a4ZtjdUrVS328e+/5eiooGuxCGiNvySbxYqlTp/qACcc2bb8J99wUd\nRXgdfrg14hdfbLufH3lk0BHlNje6+MYAV8WeJxYsnI1NnMgDqgOXxI4VkQCsXAlffglnnBF0JOHW\ntSv85S+2sezGjUFHk9uc3EGlW/E+DpvJtxzYDlwde203Vml3Ajaj7xVgsVuBi0j5DB8OffrYdj9S\ntn79YOlSOP98mDIF9tkn6IhyUxhWxWsnCYmEqO8kcfTRtjC1UyePI8oShYXW1VejBrz+eu5tCeUm\n7SQhIiktWgSbN+fWxrAVVbkyvPaadYv+6U9BR5Ob1ECJBKMiO7SU2+jRtmtCZWV8udSsaVPPhw6F\nv/896Ghyjz6uIsHIdIeWjIwaZQ2UlF/jxrb7+cCBMH160NHkFjVQIsHIZIeWRpmcaNky2LBB3XsV\ncfTRdgfVp49NnhB/qIESCadkO7Q0y+SNRo+GCy9UifOKOuMMeOABK9GxYUPQ0eQGFSwUCa/Ss56S\nTtcbPHjwL8/z8/PJz88v8frbb8Of/+xyZDlqwACrpdWrF0ydmtvlSspSUFBAQUFBhd8nDBMnNc1c\nIsGDaeZ5wLvAMUleex7bZ2lE7Osl2HrE9aWOKzN/1q2DVq1g/Xrbb04qrqgIfv1r2LLF7k51Z5qe\nl9PM09V0qg+8hc00mgUkbuJ/D7AQm6n0JqDlbiLOpNqhpVzGjoXu3dU4ualSJXj5Zdi+HW6+2Ros\n8Ua6BspJTad7gX8Dx2EJ9UTs+3nAtUAb7AqxCnCpG0GLZIHhwEfAEdhYU3/gutgDbIeWr7AdWl4A\nbszkJGPGWHeUuKt6dbt7+ugjePDBoKPJXunGoBJrOkFxTafELYtaAQ/Hni/FGqb9gS1YHaiawJ7Y\nv2tdiFkkG1zm4Jj/qcgJ/vMfKCiw6rnivjp1YNw46NgRmjRRjS0vpLuDSlXrKdFc4MLY87bAwdhs\no03Ao8A3wLdYF8XkCsYrIg5NmQInngj16wcdSfY68ECYMMF2iH/33aCjyT7pGignvasPY2s05mBX\nfHOwO6ZDgduwO6oDgdrAFZkGKiLlM368TYkWb7VsaV2pAwZoIa/b0nXxla71lKym01as/zxuJdZ3\nfjbWx/5D7Pv/AjoAe20Ykm6arEgQ3JoqG4SiImugdFXvj5NPtoW8vXvDxIlw/PFBR5Qd0k37q4qN\nK52OddN9ivWdJ45B1QV2ADuxSREdgX7A8cAbwMnAf4FhsZ9/ptQ5NM1cIiFKu5kvXWoVYVev1i7c\nfho1Cm65BaZNgyOOCDqa8Mg0d9LdQaWq6RSfafQCNrtvGNYduAAYEHvtC+A1rHBhITbT78XyBigi\n5Td+PPToocbJbxddBFu3wplnWkXevLygI4q2MHx8dQclkRClO6izzoJrr9UGsUF5+ml47DH48ENo\nWnpaWQ7KNHfCkGxqoCQSotJA7dgBBxwAa9ZA3boBRSU88gi88ordSTVuHHQ0wfKqi09EImbmTDj2\nWDVOQbvrLti5E047zcakGmW0F31uUwMlkmUmT4Zu3YKOQgB+/3vYvdsaqalT1UiVl8ptiGSZKVNs\nBp+Ew+DBcPHFkJ8P334bdDTREob+dI1BSSREYQxq0yabObZxozaIDZsHH7TS8VOmwEEHBR2Nv7zc\nzVxE3JeuSkA+sBnbmWUO8Hsnbzptmu0Np8YpfO69F268Ebp0gS+/DDqaaNAYlIj/4lUCumG7tXyG\nlddYXOq4D7DS745NmaLxpzAbONA2mc3Pt41mteNE2XQHJeK/xCoBuyiuElBaubtEpk61AXkJrwED\n4KmnbDFvRHfS8o3XBQvrAaOwK8NFWOE1kVznpEpAEbZ35VysNlTrdG+6bp1Vzj32WLfCFK/07g0j\nR9rkiZEjg44mvNJ18TnpiogXLLwAK772TOx4sOKF44CLYueq5VbgIhHmZFbQv7HNmf8D9ADeBlom\nOzC+2fKCBdCqVT5VquS7EqR469RTbUnA2WfD11/Db3+bPVtTubXRcrr/HKcAg7C7KIC7Y/8+nHDM\n2NjXM2JfL4/93E5scPeQNOfQLD6JBBdn8bUHBlOcV/dg+1X+pYyfWQmciNVZS/RL/tx4Ixx2GNx+\nuwsRim/WrIFzz4U2beC557JzgotXs/gqUrCwBfA9MBS7GnwJq6orkutmA4djtdKqA5dgPROJGlGc\n0G1jz0s3TiV88AF07epqnOKDZs2sjtTGjTbBZcOGoCMKj3RdfE4LFj6B3S3Np7hgYXWgDbYb+mfA\n49gd2B9Lv4HqQUkYeVgPykmVgIuAG2LH/ge4tKw33LAB1q7VrLCoql0b3noL/vhHaNvWynacdFLQ\nUQUv3S1Xpl0Rx2AVdD/G7qQAOmEN1DmljlcXn0RCmBfqjhoFw4bB2LFBhyMVNXo0XH+9Ley95prs\nGJfyqovPSVdE3dhrYAULPwC2Aeuw7sH4wG43YGF5AxSR9KZPh86dg45C3NC7t/3/fPJJ6NvX6kvl\nqnQNVGJXxCJgJMVdEfHuiNZY194SoDtwa8LP34yVeJ8LHAs86FbgIlJsxgzo1CnoKMQtRx4Js2ZB\nzZpwwgn2PBeF4eZRXXwSCWHt4tu6tYhGjeCHH6BGjaDDEbeNHm0zNK+/Hu67L5qz/LQXn0iOmjXL\nrrLVOGWn3r1hzhyYPRvatYPPPw86Iv+ogRKJuBkzbINYyV4HHmgTYAYOhJ494c47c2NsSg2USMTN\nnKkGKhdUqgRXXQXz59uaqdat4c03IZtHSMLQn64xKImEsI5B1alTxIoV0LBh0KGIn2bOhFtvhWrV\n4JFHwj2LU2NQIjmqSRM1TrmoY0f49FObQHHllban32efBR2Vu9RAiUTcKacEHYEEpXJla5yWLrUG\n6sIL4ayzrHBlNnRMqYESibj2KmKT8/bZx+6kli+HPn3ghhtsZueQIbB9e9DRZS4M/ekag5JICOsY\n1BdfFHHccUGHIWFSWAgTJ8Izz9hYVZ8+dqfVoYPddfnNyzGoihQsBNsMcw7wbnmD85NHm4JmJCyx\nKA5PpcsrgCdjr88FTkj1RkeVzrgAhOX/keIwlStbV98ddxQwbx60aAHXXQcHH2wTK6ZOhV27Ag3R\nkXQNVLxg4VnYlkaXAa1KHRMvWHgccBW2s3miW7FtkkJ9mxT0BypRWGJRHJ5xklc9gcOwvTB/AzyX\n6s2qpqtJ4IOw/D9SHCUVFBTQrBncfTcsXAgTJtiEmrvvtn/PPRcefdQmW+zcGXS0e0vXQLXFChCu\nAnYBI4DzSh3TCpgWe74U21h2/9jXzbBEe5nwdY2IBMVJXvUCXo09nwXUw2pEiWSsdWv4wx+sQVqx\nwjaj/eor2zW9fn0r9XHttfDEE9aYrVgR7J1WumuvZAUL25U6Jl6wcAYlCxZ+DzwG/Bao40awIlnC\nSV4lO6YZsN7b0CRXNGwIl1xiD4Bt22DuXHssWgRjxljjtXYt7L+/7WbRqJH9XIMGUKeO1bGqVcu2\n2dpnH7ubr1rVuhjjY11elgvpjVXCjesLPFXqmH2BIdg402vAp1h33znAM7Fj8kk9BrUc6/7TQ4+w\nP5bjDid59S6QuD/EZKwAaGnKHz2i8Mgod9LdQa0Fmid83Ry7kku0Feif8PVK4CusdlQvrIuvBnYX\n9Ro2TpXosPKFLBJ5TvKq9DHNYt8rTfkjOasqsILigoVfsPdgbumChcOSvE9XQj6LT8RHTvKqJzAu\n9rw98IlfwYlESQ9s8sNyrOQ7lCxYeErs9SXAKKzBKq0re1fiFcll6fIKbKbfcmycN1n3noiIiIhk\nM9cWJnocxxWx888DZmKl6oOII+5kYDc2UzKoOPKxSTALgAKP4nASS0PgfaxLbAHQz4MYhmAz5eaX\ncYwfn9NEYckdJ7Eof/aWj/f5o9ypgCpYV0UeUI30fe7t8KbP3Ukcp1DcTXlWgHHEj5sKjMVmfgUR\nRz1gITZID/ZB94KTWAYDDyXE8QPpJ/qUV2cscVIlmR+f00RhyR2nsSh/SvIjf7I2d/zalSksCxOd\nxPExsDkhjma4z0kcADdj43rfexCD0zguB0ZTPMtsY4CxfEfxmro6WJLtdjmO6cCPZbzu9wLasOSO\n01iUPyX5kT9Zmzt+NVDJFh02dXCM2x9uJ3EkGkBxi+93HE2xD1l8i5uigOI4HGiA7RYyG7jSgzic\nxvISttfjt1gXwa0exVIWPz6n6c4XRO44jSWR8sef/Mna3PFrFy+nH47Sa47d/lCV5/1OxdZ3eVFM\n20kcjwN3x46thDdbRTmJoxo2g+x0oCZ2hfwJ1o/sdyz3Yt0X+cChwCRsUfhWl2NJx+vPaSbv7UdM\nyp/yx+FH/mRt7vjVQLm5MNHrOMAGdl/C+tDLumX1Mo4TsVt1sD7jHtjtu5vT9Z3EsRrrltgRe3yI\nfbDdbqCcxNIBeCD2fAW2KPwI7MrUL358Tss6X1C54zQWUP4k8iN/lDsVFJaFiU7iOAjrz/WyDJyT\nOBINxZtZSE7iOBLbZqcKdgU4H9uBO4hY/gYMij1vhCVhAw9iycPZQK8fC2jDkjtOY1H+lORH/ih3\nXBCWhYnp4ngZG0CcE3t8GlAcibxKMKdx3InNRJoP3OJRHE5iaYjtSDI3FsvlHsQwHOun34ld/fYn\n+AW0YckdJ7Eof4LJH+WOiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiISHJ5QCH+\n1QATyRZ5KHcCof/g3jsbmIGVHfgOK0NQO+H1fYAhWBXS74CBad7vcuBrYBvwFlDf5XjT+TW2Rf9m\nbEPIv2A7Ncc1iMW1DavweVma9xuI/d6bgVew3ZhFIH3uDAN+xmoabQW2UHbdp6Bz52hgAlbhtzDJ\n68od8d1lwJlADazE8TiKq3wCPAR8ANTFtub/Duie4r2OwpKwE1AL+Du2g7ATebhzFXg9VoSuKnAg\n1lj9LuH14bFHzdhxP5G6vEB3YB1WGqAeVnX0oQrGJ9kjXe4MBf7k8L3CkDstgaux0ufJGijljpTL\n3djW8Fuw7fLPd+E9LwDmJXy9FuiW8PX9pE6cB4E3Er4+BLuCrOXgvHlYUlwbO+e3wB2OIi7bQIqL\nwNWKxXNYwuuvkjpx3gT+nPD1qVgDLdHnR+4MBf7X4c+GKXcOY+8GSrmThLr4yrYcu+KqgzUcb2DF\nvoh9/8cyHh1SvGdXYEHseX2gCVYbJW4edrWXTOtSx36FfahbOv2FsJLPh2FXpr/DSlGDdX+k+l02\nYdUvk0n8fVoCu7H/bnFzcf77zMP++/rd9SLu8zp34m7E6k/NpuyaT2HMnUTKHamwOdjteabOwD6w\n8auk5tiVVPVSx6xM8fOTgd+U+t4aoIuDc+fFzpWYkH/BCsxlqj/wDcWVOTuz91XctVj3QzLLsWSP\nqxaL8aAKxCTh5HbuAJyA/UGujBXs20Lqxi1MuZPsDkq5k4TuoMp2FZZY8auho4H9Mnyv9li/d2+K\nr5K2xf6tk3BcXWzAN5ltsdcTlXV8MqsTnn+DjSNl4nys26QH9ocjHl+dUsel+31K/+6UcbxEh9e5\nQ8L7FwLjY8ekuosKU+4ko9xJQg1UagcDLwI3YXcI9bHuhfgsoc4Uzx5K9uiY8F4nAO8A/Sh5RRSf\nnXR8wveOY+9ujLiFsdfjDsXuvr4sx+91UKnna2PPr0jyOyTOjkrspjgL+29zTiymuC+xyROJV7np\nfp/Sv/t67L+LRJcfuVNeYcmdVJQ7Ui6tgR3YbX0VbPbNLqxbqzyOxj44fVK8/hBQgM3EaYU1WGem\nOLY1NqU0PhPpzdgjbjCpkzgPu9J8HfgV1re9npITNJw4Devz75Ti9eGxmGrGjvkJ+72S6Y79vq2w\nP2IF2F2ZRJtfuXMRNu28MpYzW0jdZReG3AGbkdg69n77xB5xyh0plz9jf4y/Bx7FPsDlTbIh2OBn\n4lXV/ITXq2NrGDZj00ZvK/Xzpa8oL6PkWo56Ca+9QupZTXnAHuAa7MrvO+DOcv4uAFOBnZT8fd5L\neL0+JddyXJrw2kGx4xOvKAdiv3d8LUe1DGKS8PEjdz7E/ohvxrr7Li7182HLnTysYSqMvV8hNlkj\nTrkTgHrAKGAxsAjrTxZvzEGzeLKN8scfyp0c9SrFV05V2XugUkRSU/6IeKQuJW9hRcQ55Y/kNK9n\n8bXA+qCHAv/G9tKq6fE5RbKF8kdymtcNVFWgDfBs7N/t2BYovzj00EOLAD30iMIjcQ2OH5Q/emTL\nw+/ccaRpsum0AAAWrklEQVQxJXdF6ASMLXVMURgMGjQo6BB+EZZYFEdJWKL5SflTToqjpLDEQYa5\n4/Ud1Dps9XV8i5BulFzYKSKpKX8kp1X14Rw3Y1uQVAdWYIv2RMQZ5Y/kLD8aqLnAyT6cp0Ly8/OD\nDuEXYYlFcYSC8qccFEdJYYkjU2VVn/RLrItSJNwqVaoE4ciZRMofCb1Mc0ebxYqISCipgRIRkVBS\nAyUiIqGkBkpEREJJDZSIiISSGigREQklNVAiIiG2Zw889hj07g0XXggrV6b/mWzhx0JdsOqQW7Aq\nkruAtj6dVyTqVqHcyVl79sDVV8OqVXDjjfZv27bw6qvQs2fQ0XnPrwaqCMgHNvl0PpFsodzJYXfe\nCWvWwPvvQ81YoZVOnaBPH1iwAPbbL9j4vOZnF1/YVuCLRIVyJwfNnQtvvgmjRxc3TmAN1MUXw+23\nBxebX/xqoIqAycBs4Fqfzplzvv4a7rgDWrSw/urJk4OOSFyg3MlBRUVw220weDDUr7/36w88ANOm\nwaxZvofmK7+6+DoC3wH7A5OAJcD0+IuDBw/+5cD8/PzIb3AYhNWroXNnu7L6179g9mzo1w/++le4\n9NKgo4umgoICCgoKgg6jzNwB5U82eu892LgRrk1xSVK7NtxyCzz9NLRr529sTriVO0F0HQwCtgGP\nxr7WZpcV9NNP0LEj9O9vd1Bx8+dDt24wfDicdlpw8WWLEGwWWzp3QPmTlbp1s3y+/PLUx2zaBIcc\nAl9+CQcc4F9smQjzZrE1gX1jz2sBZwLzfThvzrj/fmjfvmTjBHDMMfDaa3DNNfDzz8HEJhWi3MlB\nixfbBIjevcs+rkEDO+bll/2JKwh+XA22AN6KPa+KFV97KOF1XQFWwPLl1jgtWpT6KqpXL+jSxWYE\nSeYCuINKlzug/Mk6N98M9erB//5v+mM//9waqZUroVKIp9Jkmjth+JWUYBVw0UVw4olwzz2pj1my\nxManFi+Ghg39iy3bhKCLLxnlTxbZtg0OOgjmzYNmzdIfX1QELVvCiBH2dyCswtzFJx5ZuhRmzLDZ\nPmU58kg47zx47jl/4hKRzLz7rvWIOGmcwO6aLrzQJkZlIzVQEfbsszBgAPzqV+mPvflmeOEF2L3b\n+7hEJDPDh5d/1u0FF8Bbb6U/LorC0F2hLooMbNsGBx8Mc+ZYl4ATnTrZ4r4LL/Q2tmylLj7x0o8/\nQl6eLRmpU8f5zxUWQvPmMGWK9ZaEkbr4csybb9q4ktPGCeCmm+CZZ7yLSUQy99ZbNr28PI0TQOXK\ncP752XkXpQYqooYMgeuuK9/P9O5tg6+rVnkSkohUwIgRcMklmf3s2WfDhAnuxhMGYeiuUBdFOa1c\naavH166FatXK97PXX2+L++66y5vYspm6+MQrP/1k3XTffWe7RJTXtm3QuDGsXw+1arkfX0Wpiy+H\njBxpd0PlbZzABmBHjHA/JhHJ3Pvv21rFTBonsJ874QSb1ZtN/GqgqgBzgHd9Ol9WGzEi8/31OneG\ndetsirpEgnInB4wZY0tBKqJbN5sokU38aqBuBRZhOzNLBSxeDN9/bzPyMlGlim0oO3Kku3GJZ5Q7\nWW7XLruDOuecir3P6adnXwUDPxqoZkBP4GXC138fOf/6l3XvVamS+Xv06QOjRrkXk3hGuZMDpk+H\nww6DAw+s2Pu0awcrVtgu6NnCjwbqMeC3QKEP58p6Y8bY3noV0b69dfN9/bU7MYlnlDs5YOzYit89\ngY1Jd+oEH3xQ8fcKC68bqHOADVgfuq4AKyg+dtSlS8Xep0oV6NnTtlWR0FLu5Ijx4y0f3dCpU3ZN\nlPC6YGEHoBfWTVEDqAO8BlyVeJAKrjnz3nvQvTtUr17x9+rVy7Y++p//qfh7ZauACxY6yh1Q/kTZ\nqlXwww/Qpo0779ep095ld4IQxYKFXYE7gXNLfV/rOBw67zwbP+rbt+LvtW2b9XmvWVP+leu5KsB1\nUKlyB5Q/kfb88zBzJrz+ujvv99//wn77wYYN4VoPFZV1UMqkDO3YAdOmQY8e7rxf7dpWhXfiRHfe\nTzyn3MlC77/vXk4D1KgBxx8Ps2a5955B8rOB+gDrspAMfPghHHusXR25pWdP6/+W0FPuZKGdO+2i\n88wz3X3fbBqH0k4SETFhApx1lrvv2b27va96iET899FHcMQR7hcR7dTJug2zgRqoiHj/fWtQ3HT4\n4TbhYuFCd99XRNKbONH9nAbo0AE++QT27HH/vf2mBioCvvnGdo9wu6RzpUrFd1Ei4q+JE+GMM9x/\n3/32gyZNYNEi99/bb2qgImDCBPsgV/bg/1b37nZ3JiL+2bgRli2zRfNeaNfO7qKiTg1UBHgx/hR3\n2mn2Qd6+3Zv3F5G9TZkCXbu6s6Yxmfbt1UCJD/bsgalTvekKAFsDdcIJth+YiPhj4kT3Z+8lUgMl\nvvj8c2ja1PqUvXLGGdm3C7JIWBUVeTf+FHfMMbbX5k8/eXcOP6iBCrlJk6zOi5e6dbPziIj3liyx\n8eSWLb07R9Wqtn3SZ595dw4/+NFA1QBmAV9gdW0e8uGcWWPSJG+vtABOPtlmCq5f7+15pNyUO1ko\nntOVPN40Kxu6+fxooP4LnAocDxwbe55hub3csn07zJ5d8d3L06laFfLzs68aZxZQ7mQhr8ef4rJh\nJp9fXXz/if1bHSthvcmn80bahx/a2qfatb0/l7r5Qku5k0V27rQJSaef7v252rWDTz+N9k4xfjVQ\nlbFuivXANKy7QtKYPNmfDzJYl8OkSdH+MGcp5U4W+fhjG3tyc0/NVJo2tSKGq1Z5fy6veF0PKq4Q\n66aoC0wA8oGC+IuqZ5Pc5Mm2Hb8fDj/cBm6//NL2B5PA60HFlZk7oPyJEr+698DGuNq1s53NW7Tw\n55xxUawHFfcHYAfw19jXqmeTxPr11lBs3GhjRH4YMMBm/tx0kz/ni5oA60HFlc4dUP5EykknwaOP\n2iJdPzz0kG2T9re/+XO+VMJcD6ohUC/2/FfAGVgZaynDlCk2ccGvxgk0DhVCyp0s8v33tr3RKaf4\nd874OFRU+fHnrwnwKtYYVgZeBzRfLI3Jk71f/1TaaafBDTfA7t3+NoySknIni0yaBKee6t32Rsmc\neCLMmQO7dtl4VNT48WdoPtDGh/NkjaIia6Duusvf8zZqBAcfbFPbvdrEUspFuZNFJkzwprxGWerW\ntZxesMC2NIsa7SQRQsuWQWFhMJMV1M0n4r7CwmAaKCieKBFFaqBCKL5Pl9crzZOJTzcXEffMmwf7\n7guHHOL/udVAiasmTfJvKmppXbpYn/XWrcGcXyQbjR/vXcmcdNRAiWt27YKCAv8nSMTVrAlt21oM\nIuKOcePg7LODOfcxx9hem1Hc2VwNVMjMmgWHHQb77x9cDGeead2MIlJxP/4Ic+faspEgRHlnczVQ\nIePnSvNU1ECJuGfCBFuYW6NGcDFEtZtPDVTIeF3IzInjjrOrvq+/DjYOkWwwbhz07BlsDO3bq4GS\nCvrhB1i0CDp2DDaOypWtkXz//WDjEIm6PXssj4JuoOJ3UFHbFcuPBqo5tgvzQmABcIsP54ykiROt\nn3qffYKOxBJq/Pigo8h5yp2I++QTaNzYFssGqVkz28Fi5cpg4ygvPxqoXcBA4CigPXAT0MqH80bO\n+PHQo0fQUZju3WHaNKtfI4FR7kTcO+/A+ecHHYXp0AE++ijoKMrHjwZqHVbPBmAbsBg40IfzRkph\noXUFhKWBatgQjjwSZswIOpKcptyJsKIieOutcDVQM2cGHUX5+D0GlQecAERwuM5bn39ujUJeXtCR\nFFM3X6jkodyJlMWL4eefw7MHXseO0buD8nPP6trAKOBW7GrwFyq4Fq7uvbgePeDqq+H//i/oSIIR\nkoKFUEbugPInrN5+G847L5gty5I5/nhYsQI2b7ZNZL0UtYKF1YCxwHjg8VKvqeAaVsjskUes5EVY\nFBbCgQfaVVcQe4iFTUAFC8vKHVD+hNbJJ8PDD8PppwcdSbGuXeG++/xfaxnmgoWVgFeARSRPsJy3\nZg189RV07hx0JCVVrgznnGMDvRII5U5ErVhh2wv5VTnXqY4dozUO5UcD1RHoC5yKVQOdAwS0bWI4\njRlj4z1hLCh23nlqoAKk3Imof/wDevcOX+HPqE2UCEPvaM53UXTvDtdcA336BB3J3nbssHUcX30F\n++0XdDTBCqiLL52cz58wOuEEeOyx4PbfS+XHH21N1g8/+HtBHOYuPinDli02xhPUVvzp/OpXNi42\nblzQkYhEw5dfwvr14euyB6hf3zajjsrGsWqgAjZuHHTqZMXMwur88209h4ikN3w4XHQRVKkSdCTJ\n5edHp5yOGqiAjRwJF18cdBRl69ULpkyxuz0RSa2wEIYNg379go4kNTVQ4siWLTB1KlxwQdCRlK1+\nfau0O2ZM0JGIhNuHH1pvSFgW5ybTpQt8/HE0tjFTAxWgd96xaaj16gUdSXqXXAIjRgQdhUi4DR1q\nd09hWZybTL160LJlNMah1EAFaORI+8MfBb16wfTpNgtIRPa2datddPbtG3Qk6eXnW+9N2KmBCsjG\njfYHv1evoCNxpk4dqxE1alTQkYiE02uv2a4RBxwQdCTpnXmmVfoNOzVQAXnjDWucwjx7r7Rf/xqG\nDAk6CpHwKSqCp5+Gm28OOhJnunaFefNg06agIymbHw3UEGA9MN+Hc0VCURG88goMGBB0JOXTo4eV\ngV+0KOhIcoZyJyKmTLFdI8K2tVEqNWpYrBMnBh1J2fxooIai7VlKmD3bdmiIyoc5rmpV3UX5TLkT\nEU89ZXdPYZ4cUVrPnuFfgO9HAzUd0NB6gldegf79o/VhjuvfH15/PRpTVLOAcicCFi600u5RmByR\nqEcPK5JaWBh0JKlpDMpnmzbZRpJXXx10JJk5/HA45hibgSgi8MADMHAg1KwZdCTlk5dnRVI//TTo\nSFILxV67uVRw7YUXbIfwJk2CjiRzd9wB995rV4xRvAt0KkQFC8uUS/kTNkuXwqRJltdRdNFFdsHc\nvr277xu1goV5wLvAMUley5ndmHfuhBYtrHrusccGHU3mCgvh6KNt1lKYCix6LaDdzPNInTuQQ/kT\nRpdfDq1awR/+EHQkmVm0yJaPrF5t9d+8ot3MI+DNN6F162g3TmAf5IEDc7cUvAjYdkEffmi5EFWt\nW1s334wZQUeSnB8N1HDgI6AlsBqI6OhLxezcCfffD3/8Y9CRuOOqq2Dx4mgVP4sg5U5IFRZaw/Tg\ng1C7dtDRVMyll9oO7GEUhhGEnOiiePppm9IZ9mmd5fHqq/Dyy3YVmc1jUXEqWChxL79s406zZnnb\nNeaHlSuhbVvr5qtRw5tzqIsvxLZutSutBx8MOhJ39e1rsxLfey/oSET88803cM89th4w6o0T2Lj4\niSeGczPoLPjPG3733Wdl3Y8/PuhI3FWlCvztb3DLLbB9e9DRiHivsNB2gBk40JZbZItbb4UnnrBd\nbsIkDN0VWd1F8fHH0Ls3LFgADRoEHY03rrzSNsh89NGgI/GWuvjkD3+ADz6wncCrhmKRjjsKC202\n4ksvWb0ot6mLL4S2brUFuY8/nr2NE8Bjj8Hf/25jUSLZatQo27F81KjsapzAuipvuw0efjjoSEpS\nA+WRwkLbt65Ll/CXdK+ohg1twsSll8KaNUFHI+K+cePgxhvhrbeiUU4jEwMGwLJl4SrDoQbKI4MG\nwbp1tolkLuje3fqxzz8fNm8OOhoR97zzjlXJHTMG2rQJOhrvVK9u3fQDB8KuXUFHY9RAuayoyPqp\n33rLHvvsE3RE/rnrLtsy5YwzVHlXoq+w0Gbe3nSTzVR1ezugMDr3XGjePDwzjsMw4Js1g7zbt9uH\n+YsvbH+u/fcPOiL/FRXZXn3jx1tf/VFHBR2RezRJIncsWwbXXAO7d8M//wkHHhh0RP757js4+WR4\n9ln3Kn6HeZLEWcASYBnwOx/Ol5GKbmw4eTKcdJL9gZ45s2KNU1g2KM0kjkqVbOr5734H+fm2HdLP\nP/sfRxbJifxxS0Xj+PprG2tq3x4uuMAm/mTSOEX5v0eTJjB6tI1JTZ7sfkzl4XUDVQV4Gkuy1sBl\nQCuPz5mRTP5H7toFb79tXVrXX2/b7g8bBrVq+R+LFyoSR79+1lDPmAFHHAF//Sts3Oh/HBGX1fnj\nhUzi2LDB8rZHDxtj2ndf26X8tttsrZ9fcXgh0zjatbMekL59LXeDqv/m9WTJtsByYFXs6xHAecBi\nj8/rup07Ye1aWLHCuvA++cTKPLdubY3TxRfn1niTEy1b2gDzJ59Yd8Ehh9hGuV272h+Cli2tJk3t\n2rmxVVIGsiZ/grZzp10grVtnd0nLllmhwU8/hW+/hW7d4Ior7M4hanWdvNK1K3z0EdxwAzzzjP2d\n69HD/ub5Nc3e69M0xTa5jFsDtCt90LnnehxFTOmu+qKi4seXX9oV/549dmf0889Wln3bNvjpJ/u3\nSZPiP7LnnWf76zVu7E/sUda+vT3++19b5PjxxzYtfflyWLXKjmnQwK5ca9a0/cCqVbMkqFwZvvrK\n/pDEG7HExizLG7ZQ5U9Zli6Fzz937/0SczVZ3hYWFuduYaHl7Z499lkZO7Y4f7dvhy1b7OuGDS1f\nDz4YDj3UloAMHGilY7JtXZNbDjnEpp3PnGnVGHr3tj37mjYtztlq1exOs1Kl5DkaZr2BlxK+7guU\nnni9HCjSQ48IPJbjL+WPHtnyyCh3vL5uWAs0T/i6OXYVmOgwj2MQiSrlj4iHqgIrsKqg1YEvCOkg\nr0gIKX9EPNYDWIrd4t0TcCwiUaP8EREREclVThYcPhl7fS5wQkBxXBE7/zxgJnBsQHHEnQzsBi4M\nMI58YA6wACjwKA4nsTQE3se6uhYA/TyIYQiwHphfxjF+fE4ThSV3nMSi/NlbPt7nj3KnAqpgXRR5\nQDWS96X3BOIF0dsBnwQUxylA3djzswKMI37cVGAsNqMriDjqAQuBZrGvG3oQh9NYBgMPJcTxA+5P\n9OmMJU6qJPPjc5ooLLnjNBblT0l+5E/W5o5fm8UmLjjcRfGCw0S9gFdjz2dh/2MbBRDHx0B8P+5Z\nFH+w/I4D4GZgFPC9BzE4jeNyYDTFs8cy3A/ClVi+A+rEntfBkmy3y3FMB8ra6taPz2misOSO01iU\nPyX5kT9Zmzt+NVDJFhw2dXCM2x9uJ3EkGkBxi+93HE2xD9lzsa+LAorjcKABMA2YDVzpQRxOY3kJ\nOAr4FusiuNWjWMrix+c03fmCyB2nsSRS/viTP1mbO36tn3b64Si9/tjtD1V53u9UoD/Q0eUYnMbx\nOHB37NhKeLOLtpM4qgFtgNOBmtgV8idYP7LfsdyLdV/kA4cCk4DjgK0ux5KO15/TTN7bj5iUP+WP\nw4/8ydrc8auBcrLgsPQxzWLf8zsOsIHdl7A+dC8qGzmJ40TsVh2sz7gHdvs+xuc4VmPdEjtijw+x\nD7bbDZSTWDoAD8SerwBWAkdgV6Z+8eNzWtb5gsodp7GA8ieRH/mj3KkgJwsOEwfQ2uPN4KqTOA7C\n+nO9LE9W3gWYQ/FmFpKTOI4EJmMDsTWxAdDWAcXyN2BQ7HkjLAkbeBBLHs4Ger36nCYKS+44jUX5\nU5If+aPccUGyBYfXxR5xT8den4vdFgcRx8vYAOKc2OPTgOJI5FWCOY3jTmwm0nzgFo/icBJLQ+Bd\n7PMxHxuAdttwrJ9+J3b1259gPqeJwpI7TmJR/gSTP8odERERERERERERERERERERERERERERERER\nEZEA/D/Fq9j03d/mJAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108620250>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "thetas = linspace(0, 1, 101)\n",
    "for i, (a, b) in enumerate(ab):\n",
    "    subplot(2,2,i+1)\n",
    "    plot(thetas, stats.beta.pdf(thetas, a, b))\n",
    "    title('a=%.1f, b=%.1f' % (a, b))\n",
    "    tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# for MAP estimator with beta prior (a, b)\n",
    "# theta is the true theta for the bernoulli distribution\n",
    "def variance(theta, a, b):\n",
    "    return 3*(1-theta)*theta/(1+a+b)**2\n",
    "\n",
    "def bias(theta, a, b):\n",
    "    return ((2-a-b)*theta + a - 1)/(a+b+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a=1.0, b=1.0\n",
      "theoretical bias 0.0\n",
      "empirical bias 0.00133333333333\n",
      "theoretical variance 0.07\n",
      "empirical variance 0.0741315555556\n",
      "theoretical MSE 0.07\n",
      "empirical MSE 0.0741333333333\n",
      "a=2.0, b=4.0\n",
      "theoretical bias -0.257142857143\n",
      "empirical bias -0.256571428571\n",
      "theoretical variance 0.0128571428571\n",
      "empirical variance 0.013616\n",
      "theoretical MSE 0.0789795918367\n",
      "empirical MSE 0.0794448979592\n",
      "a=20.0, b=20.0\n",
      "theoretical bias -0.185365853659\n",
      "empirical bias -0.185268292683\n",
      "theoretical variance 0.000374776918501\n",
      "empirical variance 0.000396897085068\n",
      "theoretical MSE 0.0347352766211\n",
      "empirical MSE 0.0347212373587\n",
      "a=25.0, b=10.0\n",
      "theoretical bias 0.025\n",
      "empirical bias 0.0251111111111\n",
      "theoretical variance 0.000486111111111\n",
      "empirical variance 0.000514802469136\n",
      "theoretical MSE 0.00111111111111\n",
      "empirical MSE 0.00114537037037\n"
     ]
    }
   ],
   "source": [
    "for a, b in ab:\n",
    "    # MAP estimator of theta with beta prior\n",
    "    theta_hat_map = (coins.sum(axis=0) + a - 1) / (a+b+1)\n",
    "    e_theta_hat_map = theta_hat_map.mean()\n",
    "    print 'a=%.1f, b=%.1f' % (a, b)\n",
    "    print 'theoretical bias', bias(theta, a, b)\n",
    "    print 'empirical bias', (theta_hat_map - theta).mean()\n",
    "    print 'theoretical variance', variance(theta, a, b)\n",
    "    print 'empirical variance', ((theta_hat_map - e_theta_hat_map)**2).mean()\n",
    "    print 'theoretical MSE', bias(theta, a, b)**2 + variance(theta, a, b)\n",
    "    print 'empirical MSE', ((theta_hat_map - theta)**2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}