bfarzin/Bayesian CR.ipynb

## Bayesian CR.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.stats import norm,invgamma\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data1 = {0: 46.0, 1: -556.0, 2: 637.0, 3: -11992.0, 4: -163.0,\n",
    " 5: -281.0, 6: 348.0, 7: -1508.0, 8: 1009.0, 9: -1056.0,\n",
    " 10: 1580.0, 11: 2941.0, 12: -1180.0, 13: -21527.0, 14: -16924.0,\n",
    " 15: -2151.0, 16: 2606.0, 17: -11803.0, 18: -714.0, 19: 2342.0,\n",
    " 20: 2517.0, 21: 1447.0, 22: 1627.0, 23: -1306.0, 24: -2715.0,\n",
    " 25: -2279.0}\n",
    "data2 = {0: -3601.0, 1: -318.0, 2: -4577.0, 3: 862.0, 4: 490.0,\n",
    " 5: -302.0, 6: 878.0, 7: 719.0, 8: 935.0, 9: 645.0,\n",
    " 10: -1554.0, 11: 627.0, 12: 8083.0, 13: 2360.0, 14: 843.0,\n",
    " 15: 2398.0, 16: 2929.0, 17: 4212.0, 18: 2755.0, 19: 4223.0,\n",
    " 20: -897.0, 21: 1412.0, 22: 164.0, 23: -5466.0, 24: 223.0,\n",
    " 25: -1201.0, 26: 1191.0}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ex1_data = pd.Series(data1)\n",
    "ex2_data = pd.Series(data2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def draw_mus_and_sigmas(data,m0,k0,s_sq0,v0,n_samples=10000):\n",
    "    # number of samples\n",
    "    N = np.size(data)\n",
    "    # find the mean of the data\n",
    "    the_mean = np.mean(data) \n",
    "    # sum of squared differences between data and mean\n",
    "    SSD = np.sum( (data - the_mean)**2 ) \n",
    "\n",
    "    # combining the prior with the data - page 79 of Gelman et al.\n",
    "    # to make sense of this note that \n",
    "    # inv-chi-sq(v,s^2) = inv-gamma(v/2,(v*s^2)/2)\n",
    "    kN = float(k0 + N)\n",
    "    mN = (k0/kN)*m0 + (N/kN)*the_mean\n",
    "    vN = v0 + N\n",
    "    vN_times_s_sqN = v0*s_sq0 + SSD + (N*k0*(m0-the_mean)**2)/kN\n",
    "\n",
    "    # 1) draw the variances from an inverse gamma \n",
    "    # (params: alpha, beta)\n",
    "    alpha = vN/2\n",
    "    beta = vN_times_s_sqN/2\n",
    "    # thanks to wikipedia, we know that:\n",
    "    # if X ~ inv-gamma(a,1) then b*X ~ inv-gamma(a,b)\n",
    "    sig_sq_samples = beta*invgamma.rvs(alpha,size=n_samples)\n",
    "\n",
    "    # 2) draw means from a normal conditioned on the drawn sigmas\n",
    "    # (params: mean_norm, var_norm)\n",
    "    mean_norm = mN\n",
    "    var_norm = np.sqrt(sig_sq_samples)/kN\n",
    "    mu_samples = norm.rvs(mean_norm,scale=var_norm,size=n_samples)\n",
    "\n",
    "    # 3) return the mu_samples and sig_sq_samples\n",
    "    return mu_samples, sig_sq_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count       26.000000\n",
       "mean     -2271.346154\n",
       "std       6195.052781\n",
       "min     -21527.000000\n",
       "25%      -1990.250000\n",
       "50%       -418.500000\n",
       "75%       1337.500000\n",
       "max       2941.000000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex1_data.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count      27.000000\n",
       "mean      667.888889\n",
       "std      2725.633297\n",
       "min     -5466.000000\n",
       "25%      -310.000000\n",
       "50%       719.000000\n",
       "75%      1886.000000\n",
       "max      8083.000000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex2_data.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu_samp_1, sig_samp_1 = draw_mus_and_sigmas(ex1_data.values,0,10,6195,100)\n",
    "mu_samp_2, sig_samp_2 = draw_mus_and_sigmas(ex2_data.values,0,10,2725,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def my_hist(data,label=None):\n",
    "    weights = np.ones_like(data)/float(len(data))\n",
    "    plt.hist(data, weights=weights,label=label,alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1640.8321913423276\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFXCAYAAABZQMyNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X90VPWd//HX/CAhPwbJkMloEdGmbFizYRvbYt0oWE8U\nU6WHrxUaVFDPqfbbGKUmdaGBDXRXGqtt+eLm9MvhWNz1fD2mKs13pd9SSrdflm+bYGJRxLhdS9rG\nVCRh8oMkM4Exmfv9wzoSSTI3MpPJhzwf53hO7tz53Pue91x8zb0z916HZVmWAACAkZzJLgAAAHx8\nBDkAAAYjyAEAMBhBDgCAwQhyAAAMRpADAGAwW0FeU1Oj0tJSrV69WkePHh31Od///ve1Zs2aCY0B\nAADnJ2aQNzc3q62tTXV1dXr00Ue1devWc57T2tqqV155RQ6Hw/YYAABw/mIGeWNjo4qLiyVJubm5\n6uvrUzAYHPGc7373u6qsrJzQGAAAcP5iBnkgEJDX641OZ2VlKRAIRKfr6+t1zTXX6JJLLrE9BgAA\nxMeEf+x29hVdT506pX/7t3/T3XffrfGu9MpVYAEASAx3rCfk5OSM2Jvu7OyUz+eTJB06dEhdXV26\n4447dObMGbW3t+uxxx5TTk6OTp48OeqYsViWFf2OHQAA2BMzyIuKilRbW6tVq1appaVFfr9f6enp\nkqRly5Zp2bJlkqR33nlH3/rWt7Rhwwa9+uqrqq2t1Ve+8pVzxozF4XDo5Mn+OLykC5vP56FPNtEr\ne+iTPfTJPnplj8/nictyYgZ5YWGh8vPzVVpaKpfLperqatXX18vj8UR/0GZnDAAAiD/HVLqNKZ/g\nYuOTrn30yh76ZA99so9e2ROvPXKu7AYAgMEIcgAADEaQAwBgMIIcAACDEeQAABgs5ulnAAAkWiQS\nUW9vT1yXOXt2lpzOsfdXa2v/h/7rv/5T3d1dGhwc1KWXztOsWbP06KOPj7vcUCiolpaj+tznPq9d\nu3Zq9uws3XbbyrjWPhEEOQAg6Xp7e/TiL48qPXNWXJYXGujT7cUF8nrnjPmc8vJvSJL27v2p/vjH\nVpWVrbO17N/97j/V1PSyPve5z8el1vNFkAPABSoRe7l2OJ1hdXePfh75eHvJ6ZmzlOmZncjSYnr1\n1d+qru5/aXBwUGVlD+mb33xIP/3pLyVJmzat15e/vErbtj2uUCikyy6bL0l6663f6ZFH1un48Xe0\nbt03tXjx5AY8QQ4AF6h47+XalZ6eolAofM7jdvaSp4I//KFVzz33E7ndbkkj7wHicDh0xx1r9cc/\ntmr58hXatWungsGgnnhiu5qaDqm+/kWCHAAQP8nYy83ISJXTdWZS1xlPn/rUgr+EuD2LFn1akpSd\n7VMoFExUWWPiV+sAAJzF7Z5x1tSHVzEfHh4a9fkul+vDZyfhqufskQMAMAaHw6kzZ87Isiy99dZ/\n/eUxh4aGRg/1ZCDIAQBTQmigb8ota8WKL+v+++/R5ZdfoYUL/1qS9Fd/tVA7dtQqJ8cvh8MRYwmJ\nx93PDMNdheyjV/bQJ3tM7FN3d5d+dqgtKd+RB4Pnfkc+0N+rL35+/qg/dkvGeeTJNmn3IwcAINGc\nTueU/zX7VDV1P6oAAICYCHIAAAxGkAMAYDCCHAAAgxHkAAAYjF+tAwCSLlmnn5048a7Wri2NniNu\nWZYcDoe2bn1CHs/I08M6Ok6ou7tLf/3X+XGt83wR5ACApOvt7dH/PvpTZczKjMvygn0DWlFwq61T\n2ubPn68nn9wR83mHD7+iUChEkAMAMJqMWZnKvCg+F0k5X83Nh7Rz5//UzJkz5fV69fDD67Vr1065\n3W5dfPElKiq6LtklRhHkAIBpbbTrm+7e/bwefPBhLVr0aR08eECWFVFJya2aPTtrSoW4RJADAKa5\nt99u00MP/ffo9+OXXTZfN9xwo5544ju66aYvqrj4JmVleZNd5pgIcgDAtDbWd+SLF1+jgwf/rzZs\nqNA//dNjSajMHk4/AwBMa6MdWv+Xf3lKLpdLX/rSf9MNN9yoP/3pT3I6nVPq9qUfYI8cADAlBPsG\nkrKs9vb3D61LH55+dtVVn9U3vlEmj8ejWbMuUmnpXUpLS9N3vvNtZWVl6cYbb45breeL25gaxsRb\nKSYLvbKHPtljYp+4jSm3MQUAYFJwG9OPz1aQ19TU6MiRI3I4HKqqqlJBQUF03vPPP6/du3fL5XJp\n4cKFqq6uVlNTk9atW6cFCxbIsizl5eVp06ZNCXsRAABMVzGDvLm5WW1tbaqrq1Nra6s2btyouro6\nSdLp06e1d+9ePffcc3I6nbr77rv12muvSZIWL16s7du3J7Z6AACmuZhfHjQ2Nqq4uFiSlJubq76+\nPgWDQUnSzJkz9fTTT8vpdGpwcFADAwPKzs6W9P4PBgAAQGLFDPJAICCv98MT4bOyshQIBEY8Z+fO\nnbrppptUUlKiSy+9VJLU2tqqsrIy3XnnnWpoaIhz2QAAQPoYP3YbbU/7/vvv1z333KOvfvWr+sxn\nPqPLL79c5eXlKikpUXt7u9auXav9+/fL7R5/dfH6Bd+Fjj7ZR6/soU/2mNYnpzOs9PQUZWSkTvq6\nR1tnZDhF2dkezZljVh+nuphBnpOTM2IPvLOzUz6fT5LU29urt956S4sXL1ZKSoqWLFmiw4cPq7Cw\nUCUlJZKkefPmKTs7Wx0dHZo7d+646zLt1I5kMPEUmGShV/bQJ3tM7FN3d79CobCcrnNPBUuksU4/\nC4XCCgT6FYmkTGo9U9WknX5WVFSk2tparVq1Si0tLfL7/UpPT5ckDQ8Pq6qqSnv27FFaWppef/11\nrVixQnv27FFbW5vKy8vV1dWl7u5u+f3+uBQMADBTJBJRT098zxWPh6l+vnksMYO8sLBQ+fn5Ki0t\nlcvlUnV1terr6+XxeFRcXKzy8nKtWbNGbrdbCxcu1A033KBgMKjKykqtXr1almVpy5YtMQ+rAwAu\nbIOhfv2f3wTkzc5JdilRoYE+3V5cYPQ57FzZzTAmHt5LFnplD32yx8Q+TbUru3W82yaXK1XZORdP\naj3jGe9qc4kWr0Pr5h5LAAAABDkAACYjyAEAMBhBDgCAwQhyAAAMRpADAGAwghwAAIMR5AAAGIwg\nBwDAYAQ5AAAGI8gBADAYQQ4AgMEIcgAADEaQAwBgMIIcAACDEeQAABiMIAcAwGAEOQAABiPIAQAw\nGEEOAIDBCHIAAAxGkAMAYDCCHAAAgxHkAAAYjCAHAMBgBDkAAAYjyAFgmohEIhro71UkEkl2KYgj\nghwApolQsE+vvbtfoWBfsktBHBHkADCNpKbNTHYJiDO3nSfV1NToyJEjcjgcqqqqUkFBQXTe888/\nr927d8vlcmnhwoWqrq6OOQYAAMRHzCBvbm5WW1ub6urq1Nraqo0bN6qurk6SdPr0ae3du1fPPfec\nnE6n7r77br322mt67733xhwDAADiJ+ah9cbGRhUXF0uScnNz1dfXp2AwKEmaOXOmnn76aTmdTg0O\nDmpgYEDZ2dnjjgEAAPETM8gDgYC8Xm90OisrS4FAYMRzdu7cqZtuukklJSW69NJLbY0BAADnb8I/\ndrMs65zH7r//fv37v/+7Dh48qMOHD9saAwAAzl/M78hzcnJG7E13dnbK5/NJknp7e/XWW29p8eLF\nSklJ0ZIlS3T48OFxx4zH5/N8nNcw7dAn++iVPfTJHtP65HSGlZ6eooyMVElSZDhF7gH3iMcSZbTl\np6WlyOWekfB1T0RkOEXZ2R7NmWPWe3u2mEFeVFSk2tparVq1Si0tLfL7/UpPT5ckDQ8Pq6qqSnv2\n7FFaWppef/11rVixQllZWWOOGc/Jk/3n/4oucD6fhz7ZRK/soU/2mNin7u5+hUJhOV1nJEmhUFhD\nQ0MjHkuEjIxUBYPnLn9wMCyXyzHqvGQJhcIKBPoViaRM+rrj9cEwZpAXFhYqPz9fpaWlcrlcqq6u\nVn19vTwej4qLi1VeXq41a9bI7XZr4cKFuuGGGyTpnDEAACD+bJ1HXlFRMWI6Ly8v+veKFSu0YsWK\nmGMAAED8cWU3AAAMRpADAGAwghwAAIMR5AAAGIwgBwDAYAQ5AAAGI8gBADAYQQ4AgMEIcgAADEaQ\nAwBgMIIcAACDEeQAABiMIAcAwGAEOQAABiPIAQAwGEEOAIDBCHIAAAxGkAMAYDCCHAAAgxHkAAAY\njCAHAMBgBDkATCOngyEFB04luwzEEUEOAIDBCHIAAAxGkAMAYDCCHAAAgxHkAAAYjCAHAMBgBDkA\nAAZz23lSTU2Njhw5IofDoaqqKhUUFETnHTp0SNu2bZPL5dIVV1yhrVu3qqmpSevWrdOCBQtkWZby\n8vK0adOmhL0IAACmq5hB3tzcrLa2NtXV1am1tVUbN25UXV1ddP7mzZv1zDPPyO/3a926dTp48KBm\nzpypxYsXa/v27QktHgCA6S7mofXGxkYVFxdLknJzc9XX16dgMBidv3v3bvn9fkmS1+tVb2+vJMmy\nrETUCwAAzhIzyAOBgLxeb3Q6KytLgUAgOp2ZmSlJ6uzsVENDg5YuXSpJam1tVVlZme688041NDTE\nu24AACCb35GfbbQ97a6uLn3961/Xli1bdNFFF2n+/PkqLy9XSUmJ2tvbtXbtWu3fv19u94RXBwAA\nxhEzWXNyckbsgXd2dsrn80WnBwYGdN9996myslLXXHONJMnv96ukpESSNG/ePGVnZ6ujo0Nz584d\nd10+n+djvYjphj7ZR6/soU/2mNYnpzOs9PQUZWSkSpIiwylyu1xKS/vwsUQZbflpaSlyuWckfN0T\nERlOUXa2R3PmmPXeni1mkBcVFam2tlarVq1SS0uL/H6/0tPTo/Mfe+wx3XvvvSoqKoo+tmfPHrW1\ntam8vFxdXV3q7u6Ofo8+npMn+z/my5g+fD4PfbKJXtlDn+wxsU/d3f0KhcJyus5IkkKhsIaGhzU4\nGFYweCZh683ISB11+YODYblcjoSue6JCobACgX5FIimTvu54fTCMGeSFhYXKz89XaWmpXC6Xqqur\nVV9fL4/Ho2uvvVYvvfSS3n77bT3//PNyOBxavny5brnlFlVUVGj16tWyLEtbtmzhsDoAAAlgK10r\nKipGTOfl5UX/fv3110cds2PHjvMoCwAA2MGV3QAAMBhBDgCAwQhyAAAMRpADAGAwghwAAIMR5AAA\nGIwgBwDAYAQ5AAAGI8gBADAYQQ4AgMEIcgAADEaQAwBgMIIcAACDEeQAABiMIAcAwGAEOQAABiPI\nAQAwGEEOAIDBCHIAAAxGkAMAYDCCHAAAgxHkAAAYjCAHAMBgBDkAAAYjyAEAMBhBDgCAwQhyAAAM\nRpADAGAwghwAAIMR5AAAGMxt50k1NTU6cuSIHA6HqqqqVFBQEJ136NAhbdu2TS6XS1dccYW2bt0a\ncwwAAIiPmEHe3NystrY21dXVqbW1VRs3blRdXV10/ubNm/XMM8/I7/dr3bp1OnjwoNLS0sYdAwAA\n4iPmofXGxkYVFxdLknJzc9XX16dgMBidv3v3bvn9fkmS1+tVb29vzDEAACA+YgZ5IBCQ1+uNTmdl\nZSkQCESnMzMzJUmdnZ1qaGjQ0qVLY44BAADxYes78rNZlnXOY11dXfr617+uLVu26KKLLrI1ZjQ+\nn2ei5UxL9Mk+emUPfbLHtD45nWGlp6coIyNVkhQZTpHb5VJa2oePJcpoy09LS5HLPSPh656IyHCK\nsrM9mjPHrPf2bDGDPCcnZ8TedGdnp3w+X3R6YGBA9913nyorK3XNNdfYGjOWkyf7J1T8dOTzeeiT\nTfTKHvpkj4l96u7uVygUltN1RpIUCoU1NDyswcGwgsEzCVtvRkbqqMsfHAzL5XIkdN0TFQqFFQj0\nKxJJmfR1x+uDYcxD60VFRdq3b58kqaWlRX6/X+np6dH5jz32mO69914VFRXZHgMAAOIj5h55YWGh\n8vPzVVpaKpfLperqatXX18vj8ejaa6/VSy+9pLffflvPP/+8HA6Hli9frpUrV+rKK68cMQYAAMSf\nre/IKyoqRkzn5eVF/3799ddHHVNZWXkeZQEAADu4shsAAAYjyAEAMBhBDgCAwQhyAAAMRpADAGAw\nghwAAIMR5AAAGIwgBwDAYAQ5AAAGI8gBADAYQQ4A00AkElFw4JRs3lUaBpnw/cgBAOYJBft09MQB\nOd3sv11oeEcBYJpITZuZ7BKQAAQ5AAAGI8gBADAYQQ4AgMEIcgAADEaQAwBgMIIcAACDEeQAABiM\nIAcAwGAEOQAABiPIAQAwGEEOAIDBCHIAAAxGkAMAYDCCHAAAgxHkAAAYjCAHAMBgbjtPqqmp0ZEj\nR+RwOFRVVaWCgoLovHA4rH/4h3/QsWPHtHv3bklSU1OT1q1bpwULFsiyLOXl5WnTpk2JeQUAAExj\nMYO8ublZbW1tqqurU2trqzZu3Ki6urro/Mcff1yLFi1Sa2vriHGLFy/W9u3b418xAACIinlovbGx\nUcXFxZKk3Nxc9fX1KRgMRudXVlbq+uuvP2ecZVnxqxIAAIwqZpAHAgF5vd7odFZWlgKBQHQ6LS1t\n1HGtra0qKyvTnXfeqYaGhjiUCgAAPsrWd+Rns7OnPX/+fJWXl6ukpETt7e1au3at9u/fL7d7wqsD\nAADjiJmsOTk5I/bAOzs75fP5xh3j9/tVUlIiSZo3b56ys7PV0dGhuXPnjjvO5/PYqXnao0/20St7\n6JM9pvXJ6QwrPT1FGRmpigynyO1yyYo4lZb2/mOJNNry09JS5HLPSPi6JyIynKLsbI/mzDHrvT1b\nzCAvKipSbW2tVq1apZaWFvn9fqWnp494jmVZI/bU9+zZo7a2NpWXl6urq0vd3d3y+/0xizl5sv9j\nvITpxefz0Ceb6JU99MkeE/vU3d2vUCgsp+uMQqGwhoaHNTQc0eBgWMHgmYStNyMjddTlDw6G5XI5\nErruiQqFwgoE+hWJpEz6uuP1wTBmkBcWFio/P1+lpaVyuVyqrq5WfX29PB6PiouLde+99+rEiRN6\n9913tXz5ct1zzz0qKSlRRUWFVq9eLcuytGXLFg6rAwCQALbStaKiYsR0Xl5e9O+nn3561DE7duw4\nj7IAAIAdXNkNAACDEeQAABiMIAcAwGAEOQAABiPIAQAwGEEOAIDBCHIAAAxGkAMAYDCCHACmESti\nKTTQp0gkkuxSECcEOQBMI+HTp3Us+IpCwb5kl4I4IcgBYJpJSZuZ7BIQRwQ5AAAGI8gBADAYQQ4A\ngMEIcgAADEaQAwBgMIIcAACDEeQAABiMIAcAwGAEOQAABiPIAQAwmDvZBQDAhSASiai3tyfZZYzQ\n09Mjy7KSXQYSjCAHgDjo7e3Ri788qvTMWckuJSpw4s/KvChbnqlTEhKAIAeAOEnPnKVMz+xklxEV\nHDiV7BIwCfiOHAAAgxHkAAAYjCAHAMBgBDkAAAYjyAEAMJitIK+pqVFpaalWr16to0ePjpgXDoe1\nfv163X777bbHAACA+IgZ5M3NzWpra1NdXZ0effRRbd26dcT8xx9/XIsWLZrQGAAAEB8xg7yxsVHF\nxcWSpNzcXPX19SkYDEbnV1ZW6vrrr5/QGAAAEB8xgzwQCMjr9Uans7KyFAgEotNpaWkTHgMAAOJj\nwj92+zjX7eVavwAAJEbMS7Tm5OSM2Jvu7OyUz+eL+xhJ8vk8MZ8D+jQR9Moe+mTPeH1yOsNKT09R\nRkbqJFY0vrS0FLncM5SRkarIcIrcLpfcLqfcLmfCax1t2WfXM1VEhlOUne3RnDnm/huIGeRFRUWq\nra3VqlWr1NLSIr/fr/T09BHPsSxrxF63nTGjOXmy/2O8hOnF5/PQJ5volT30yZ5Yferu7lcoFJbT\ndWYSqxrf4GBYLpdDweAZhUJhDQ0Pa2g4oqHhSEJrzchIVTB47rLPrmeqCIXCCgT6FYmkTPq64/UB\nOmaQFxYWKj8/X6WlpXK5XKqurlZ9fb08Ho+Ki4t177336sSJE3r33Xe1fPly3XPPPfryl7+sK6+8\ncsQYAAAQf7buflZRUTFiOi8vL/r3008/PeqYysrK8ygLAADYwZXdAAAwGEEOAIDBCHIAAAxGkAMA\nYDCCHAAAgxHkAAAYjCAHAMBgBDkAAAYjyAEAMBhBDgCAwQhyAAAMRpADAGAwghwAAIPZuvsZcKGL\nRCLq7e1JdhkjzJ6dJaeTz9oAxkeQA5J6e3v04i+PKj1zVrJLkSSFBvp0e3GBvN45yS4FwBRHkAN/\nkZ45S5me2ckuAwAmhON2AAAYjCAHAMBgBDkAAAYjyAEAMBhBDgCAwQhyAAAMRpADAGAwghwAAIMR\n5AAAGIwgBwDAYAQ5AAAGI8gB4AIXiUQUHDgly0p2JUgEbpoCABe4ULBPR08ckNPNvtuFiHcVAKaB\n1LSZyS4BCWJrj7ympkZHjhyRw+FQVVWVCgoKovMaGhq0bds2uVwuLVmyRGVlZWpqatK6deu0YMEC\nWZalvLw8bdq0KWEvAgCA6SpmkDc3N6utrU11dXVqbW3Vxo0bVVdXF52/detW7dq1Szk5Obrrrru0\nbNkySdLixYu1ffv2xFUOAABiH1pvbGxUcXGxJCk3N1d9fX0KBoOSpPb2ds2ePVt+v18Oh0NLly7V\noUOHJEkWv6oAACDhYgZ5IBCQ1+uNTmdlZSkQCIw6z+v1qrOzU5LU2tqqsrIy3XnnnWpoaIh33QAA\nQB/jV+vj7Wl/MO/yyy9XeXm5SkpK1N7errVr12r//v1yu8dfnc/nmWg50xJ9ss9ur5zOsNLTU5SR\nkZrgiuyJDKcoO9ujOXMm571mm7JnvD5NtW1IktLSUuRyz1BaWorcLpesiFPDLqfcLmfCax1t2R/U\nM5V6NNn/1hIhZpDn5ORE98AlqbOzUz6fLzrv5MmT0XkdHR3KyclRTk6OSkpKJEnz5s1Tdna2Ojo6\nNHfu3HHXdfJk/8d6EdOJz+ehTzZNpFfd3f0KhcJyus4kuCp7QqGwAoF+RSIpCV8X25Q9sfo01bYh\nSRocDMvlcsiyXBoaHtbQcCT6XyJrzchIVTB47rI/qGe0eckymf/WPipeH6BjHlovKirSvn37JEkt\nLS3y+/1KT0+XJM2dO1fBYFDHjx/X0NCQDhw4oGuvvVZ79uxRbW2tJKmrq0vd3d3y+/1xKRgAAHwo\n5h55YWGh8vPzVVpaKpfLperqatXX18vj8ai4uFibN29WRUWFJOnWW2/V/PnzlZ2drcrKSq1evVqW\nZWnLli0xD6sDAICJs5WuHwT1B/Ly8qJ/f/aznx1xOpokZWRkaMeOHXEoDwAAjIcruwEAYDCCHAAA\ngxHkAAAYjF+gTQORSES9vT3JLmOE2bOz5HTyORIAzhdBngC/+n9N6hoYSsiyMzNSNTDBczC7O48r\nkuqVZ9bshNQ0UaGBPt1eXCCvd06ySwEA4xHkCTAUcWpG5sUJWbY7I1UzHBMLcqu3X2npmcr0TI0g\nB5BcVsRScOCU0jNmcWTsAsA7CADTTPj0ab3Z/WuFgn3JLgVxQJADwDSUmjYz2SUgTghyAAAMRpAD\nAGAwghwAAIMR5AAAGIwgBwDAYAQ5AAAGI8gBADAYQQ4AgMEIcgAADEaQAwBgMIIcAACDEeQAABiM\nIAcAwGAEOQAABiPIAQAwGEEOAIDBCHIAAAxGkAPABSwSiSg4cEqWlexKkCjuZBcAAEic06EBvXWq\nUU43+20XKt5ZALjApabNTHYJSCCCHACmIStiKThwSpFIJNml4DzZCvKamhqVlpZq9erVOnr06Ih5\nDQ0NWrlypUpLS/XDH/7Q1hgAQHKdGRzUm92/VijYl+xScJ5ifkfe3NystrY21dXVqbW1VRs3blRd\nXV10/tatW7Vr1y7l5OTorrvu0rJly9Td3T3uGABA8nHI/cIQM8gbGxtVXFwsScrNzVVfX5+CwaAy\nMjLU3t6u2bNny+/3S5KWLl2qxsZGdXd3jzkGAADET8xD64FAQF6vNzqdlZWlQCAw6jyv16uTJ0+O\nOwYAMDkikWGdOTOY7DKQYBM+/cwa52TEseaNN+ZCZEXCCp06npiFh1MUCoYnNOS9UI8sxwy53VPj\nbMPQQJ96enoSvh6nM6zu7n5bz+3p6VFoYOp8VzhZPZIm1qfpLFafpto2JEnH/3xMf+p4UxmXZGro\nvff/vxEePC3nDJfc7hkaHh7SyeH2uK83MpyiUOjc/08NDvTL6T6jgf6pc0h/qr1nH0fM/7Pn5OSM\n2Jvu7OyUz+eLzjt58mR0XkdHh3JycjRjxowxx4zH5/NMqPip6q5Vy5JdAv5izpw5Np95uT7/+cKE\n1jKV2e/T9DZ+n6biNrQ02QVgEsQ8tF5UVKR9+/ZJklpaWuT3+5Weni5Jmjt3roLBoI4fP66hoSEd\nOHBA11577bhjAABA/DgsG8e9f/CDH6ipqUkul0vV1dV688035fF4VFxcrFdeeUXf+973JEk333yz\n7rnnnlHH5OXlJfSFAAAwHdkKcgAAMDVxZTcAAAxGkAMAYDCCHAAAgyXsxOLh4WFt3LhRb7/9tiKR\niP7+7/9eV111ldasWaPTp09r5syZcjgc2rBhg6688ko99dRT2rdvn5xOp8rKyrR06VINDAyosrJS\n/f39ysjI0Pe//33NmjVLDQ0N2rZtm1wul5YsWaKysrJEvYxJMVavfve732nLli1yOp3Ky8vT5s2b\nJWla9+rll1/Www8/rJqaGi1d+v6pNWxToxutV2xTY6uvr9f27dt12WWXSXr/jJ2vfe1rcenZdFFT\nU6MjR47I4XCoqqpKBQUFyS5p0jU1NWndunVasGCBLMtSXl6evvrVr+qRRx6RZVny+Xx6/PHHNWPG\nDL300kt65pln5HK5tHLlSt1+++0aGhrShg0bdPz4cblcLtXU1OjSSy8df6VWguzevdvavHmzZVmW\n9fvf/966/fbbLcuyrLvuuss6duzYiOe2t7dbt912mzU0NGR1dXVZN998sxWJRKx//ud/tn70ox9Z\nlmVZP/7xj63vfe97lmVZ1he/+EXrxIkTViQSse64445zlmeasXq1Zs0a64033rAsy7IqKiqsgwcP\nTutetbUZe9auAAAFmklEQVS1WQ888ID14IMPWgcOHIg+zjZ1rrF6xTY1tp/85CfWd7/73XMeP5+e\nPfHEE5P6GpKpqanJ+trXvmZZlmUdO3bM+spXvpLkipLj5Zdfth566KERj23YsMHat2+fZVmW9YMf\n/MB67rnnrFAoZC1btswaGBiwTp8+bd16663WqVOnrPr6eusf//EfLcuyrF//+tfWN77xjZjrTNih\n9S996Uv61re+Jen9S7eeOnXq7A8PI5778ssva8mSJXK5XPJ6vZo7d65+//vf69ChQ7rxxhslSV/4\nwhf0m9/8ZsT13R0Oh5YuXapDhw4l6mVMitF69d577+nPf/6z8vPzJUk33HCDGhoapnWvLr74YtXW\n1o56zX62qZFG69V7772nd955h21qAs63Zw0NDcksf1KNdV+O6eij/z9qamrSF77wBUkfbhdHjhzR\nokWLlJGRodTUVF111VX67W9/O6KPf/d3f6fDhw/HXF/CDq273e7oJUH/9V//VcuXL4/Oe/LJJ9Xd\n3a3c3FxVVVWdc232OXPmRK/ZnpWVNeKxrq6uc67v3t4e/0sMTqbRetXT06PZs2dHn+P1etXZ2ams\nrKxp26uUlJQx57FNjTRar3p6enTRRRdFp9mmztXU1KT77rtPQ0NDWr9+vbxe73n1bDrdYyIQCOhv\n/uZvotMf3GNjOt4sq7W1VWVlZTp16pQeeOABnT59WjNmzJD0/nbR2dk56r+lj96rxOFwyOl0amho\naNxLbMclyF944QW9+OKLcjgcsixLDodDDz74oIqKivTss8/qzTff1I4dOyRJd999t/Ly8jRv3jx9\n+9vf1rPPPnvO8ka70f0Hyx3tcZPY7VVXV5et5V2ovRqvTx/FNmW/V3ZcyL36wGg9u+WWW/Tggw9q\n6dKleu211/TII4/oRz/6ka3XOFbPprPp+vrnz5+v8vJylZSUqL29XWvXrtXQ0FB0/lh9Gevx0bat\nj4pLkK9cuVIrV6485/EXXnhBBw4c0A9/+EO5XC5Jih4ykKTrr79eP//5z3X11VfrD3/4Q/Txjo4O\n+f3+6HXeMzMzo9dxH+v67qaw2yuv1zviphln92Q69GqsPo2Gbcper6b7NnW2WD379Kc/rZ6eHmVl\nZam3tzf6+Mfp2XQx3n05phO/36+SkhJJ0rx585Sdna033nhD4XBYKSkpI7aVj/5bKiwsjPYxLy8v\n+gEg1g2vEvYdeXt7u3784x+rtrY2ekhBev8Xxh+82a+88ooWLFigq6++Wv/xH/+hoaEhdXR0qLOz\nU5/61KdUVFSkvXv3SpJ+8Ytf6LrrrtMnPvGJUa/vbrLReuV2u/XJT34y+v3IB69/uvfqA2d/emWb\nGt8HvWKbGt9TTz2lF154QZJ07Ngxeb1ezZgx47x7Nl1wj4337dmzR7W1tZKkrq4udXV16bbbbtPP\nf/5zSdK+fft03XXXadGiRXrjjTc0MDCgYDCoV199VZ/5zGdUVFQUfe6vfvUrXX311THXmbBLtG7b\ntk0/+9nPdMkll0QPXe3atUv79+/Xzp07lZmZqZycHH3nO99Ramqqnn32Wb300ktyOBx6+OGHdfXV\nVysUCumRRx5Rb2+vZs2apSeeeEKZmZljXt/dVGP1qq2tTdXV1bIsS3/7t3+r9evXS9K07dX+/fv1\n5JNPqrOzUxkZGcrKytLu3bu1d+9etqmPGKtXra2tbFNj6Ojo0De/+U1J7x/O3LBhgwoKCuLSs+mC\ne2xIwWBQlZWVOnXqlCzL0gMPPKCFCxdq/fr1CofD+sQnPqGamhq5XC794he/0FNPPSWn06k1a9bo\nlltuUSQS0caNG9XW1qbU1FQ99thj8vv9466Ta60DAGAwruwGAIDBCHIAAAxGkAMAYDCCHAAAgxHk\nAAAYjCAHAMBgBDkAAAYjyAEAMNj/B8PjROCy6KFuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f93fa00cfd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_hist(ex1_data.values,label='Truth')\n",
    "my_hist(mu_samp_1,label='Est')\n",
    "plt.legend()\n",
    "print(np.mean(mu_samp_1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "487.1483416209812\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFXCAYAAAC7nNf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X98VPWd7/H3zIRAJplAhkxGixQt0tBi6EYLSiNi2RRM\nK31wFWkQpfh4VPeKsdikFUxo4HbLBrHKpaYty0Nw17YPI8ojLXRrUXdrqU0osWi0sdY2tjFCSTL5\nQZKZQCBz7h9eRiIhMyEzTPLl9fwr58f35PPhhHnPOWfmHJtlWZYAAMCoZ493AQAAIDoIdQAADEGo\nAwBgCEIdAABDEOoAABiCUAcAwBAJkaxUVlam2tpa2Ww2FRcXKysr66x1Hn30Ub3++uv68Y9/rIMH\nD2r16tWaNm2aLMtSZmam1q1bF/XiAQDAh8KGek1NjRoaGlRRUaH6+nqVlJSooqKi3zr19fV69dVX\nNWbMmNC82bNna+vWrdGvGAAADCjs6ffq6mrl5uZKkqZOnarOzk75/f5+6zz88MMqKirqN4972gAA\ncGGFDXWfzye32x2aTktLk8/nC01XVlZqzpw5uvTSS/uNq6+v16pVq7R8+XJVVVVFsWQAADCQiK6p\nn+nMI/Bjx47p5z//uXbu3KkjR46Elk2ZMkUFBQXKy8tTY2OjVqxYoRdffFEJCUP+dQAAIEJhj9Qz\nMjL6HZk3NzfL4/FIkg4cOKDW1lbdfvvtuv/++/WnP/1JmzZtktfrVV5eniRp8uTJSk9PV1NT06C/\nh9P1AAAMT9hD55ycHJWXl2vp0qWqq6uT1+uV0+mUJC1cuFALFy6UJB0+fFgPPfSQ1q5dq71796qh\noUEFBQVqbW1VW1ubvF7voL/HZrOppaUrCi2NTB6Py9j+TO5Nor/Rjv5GL5N7kz7oL9rChnp2drZm\nzJih/Px8ORwOlZaWqrKyUi6XK/QBuo+aP3++ioqKtGzZMlmWpQ0bNnDqHQCAGLONpEevmv6OzNT+\nTO5Nor/Rjv5GL5N7k2JzpM4d5QAAMAShDgCAIQh1AAAMQagDAGAIQh0AAEPwPTMAQNwFg0F1dLT3\nm2e396qt7fw//T5hQprs9nMfu5aX/1/9+c9/Ultbq3p6enTZZZOVmpqq735386DbDQT8qqt7U7Nm\nXaedO7drwoQ03XLLbeddZzQR6gCAuOvoaNdzL70pZ0pqaJ7TmahAoPe8thfo7tSS3Cy53RPPuU5B\nwQOSpOef/4X+9rd6rVq1OqJtv/32n3Tw4O81a9Z151VbLBHqAIARwZmSqhTXhNB0cvJY2R0nLmgN\nr732B1VU/EQ9PT1aterr+uY3v65f/OIlSdK6dWt0661LtWXLZgUCAX3841MkSe+887a+9a3VOnLk\nsFav/qZmz45f2HNNHQCAM7z7br0ee6xc06d/WpKt3zKbzabbb1+hf/7nL2jRosWSJL/fr0ce2arV\nq7+pysrn4lDxhwh1AADOcOWV04Z0a/OZM/9JkpSe7lEg4I9VWREh1AEAOENCwpgzpj68k3pf36kB\n13c4HB+uHec7r3NNHQCAc7DZ7Dpx4oQsy9I77/z5/8+z6dSpgQM+3gh1AMCIEOju7Dcd7Bvep9+j\nYfHiW3XPPSt1+eVXaPr0T0mSPvnJ6dq2rVwZGV7ZbLYwW7iweErbBWLy04ZM7k2iv9GO/kaHgb6n\nnp7uks8Xu++px1tcnqcOAECs2e32s75TPnGiS8FgYpwqGp1G7lsYAAAwJIQ6AACGINQBADAEoQ4A\ngCEIdQAADMGn3wEAcRePR69K0tGj/9CKFfmh76BbliWbzaaNGx+Ry9X/K2dNTUfV1taqT31qxnnX\nFGuEOgAg7jo62vWzN3+h5NSU0LyklkT1nOfNZ/yd3VqcdfOgj149bcqUKfr+97eFXe/QoVcVCAQI\ndQAAwklOTVHK+A+Pjp3JY+UYc2EfvXpaTc0Bbd/+I40bN05ut1vf+MYa7dy5XQkJCbrkkkuVkzM3\nLnWFQ6gDAC5qA91XdffuXbr//m9o5sx/0v79L8uygsrLu1kTJqSN2ECXCHXgnE5f4xvudb1YGOm3\nvwRGk/fea9DXv/6/Q9fTP/7xKZo//wt65JF/04IFX1Ru7gKlpbnjXWZEIgr1srIy1dbWymazqbi4\nWFlZWWet8+ijj+r111/Xj3/844jHACNZR0e7nnvpTaVnpJ/3QyViIdDdqSW5WRFdKwQQ3rmuqc+e\nPUf79/9aa9cW6l//dVMcKhu6sKFeU1OjhoYGVVRUqL6+XiUlJaqoqOi3Tn19vV599VWNGTMm4jHA\naOBMSZUrNU12R3yu6wGIvYFOv//HfzyhW2/9ir785f+l9vY2/f3vf5fdbh+xj1w9LWyoV1dXKzc3\nV5I0depUdXZ2yu/3Kzk5ObTOww8/rKKiIn3/+9+PeAwAAGfyd3b3m+47eWJYn36PVGPjB6ffpQ+/\n0nb11Z/VAw+sksvlUmrqeOXn36GkpCT927/9H6WlpekLX7jpvOqKtbCh7vP5dNVVV4Wm09LS5PP5\nQgFdWVmpOXPm6NJLL414DAAAZ5owIU2Ls27uNy8aj14N55JLLtW+fb8ZcNnKlV/rNz1r1rWqrPzl\neddzIQz5g3JnPn792LFj+vnPf66dO3fqyJEjEY0BAOCjePRqdIQN9YyMDPl8vtB0c3OzPB6PJOnA\ngQNqbW3V7bffrhMnTqixsVGbNm1SRkaGWlpaBhwzmFg8MH4kMbk/E3uz23vldH7wgpKcPDbO1Xwo\n2Jeo9HSXJk6M3r+5ifvvTPQ3epncWyyEDfWcnByVl5dr6dKlqqurk9frldPplCQtXLhQCxculCQd\nPnxYDz30kNauXavXXntN5eXl+spXvnLWmMG0tIysrw1Fk8fjMrY/U3tra+tSINArV6rk94+cD8oF\nAr3y+bqidgRj6v47jf5GL5N7k2LzhiVsqGdnZ2vGjBnKz8+Xw+FQaWmpKisr5XK5Qh+Gi2QMAACI\nrYiuqRcWFvabzszMPGudSZMm6amnnjrnGAAAEFvckgoAAEMQ6gAAGIJQBwDAEIQ6AACGINQBADAE\noQ4AgCEIdQAADEGoAwBgCEIdAABDEOoAABiCUAcAwBCEOgAAhiDUAQAwBKEOAIAhCHUAAAxBqAMA\nYAhCHQAAQxDqAAAYglAHAMAQhDoAAIYg1AEAMAShDgCAIQh1AAAMQagDAGAIQh0AAEMkRLJSWVmZ\namtrZbPZVFxcrKysrNCyXbt2affu3XI4HJo+fbpKS0t18OBBrV69WtOmTZNlWcrMzNS6deti1gQA\nAIgg1GtqatTQ0KCKigrV19erpKREFRUVkqTjx4/r+eef19NPPy273a6vfvWrev311yVJs2fP1tat\nW2NbPQAACAl7+r26ulq5ubmSpKlTp6qzs1N+v1+SNG7cOD355JOy2+3q6elRd3e30tPTJUmWZcWw\nbAAA8FFhQ93n88ntdoem09LS5PP5+q2zfft2LViwQHl5ebrsssskSfX19Vq1apWWL1+uqqqqKJcN\nAAA+KqJr6mca6Aj8nnvu0cqVK/W1r31N11xzjS6//HIVFBQoLy9PjY2NWrFihV588UUlJAz51wEA\ngAiFTdmMjIx+R+bNzc3yeDySpI6ODr3zzjuaPXu2EhMTdcMNN+jQoUPKzs5WXl6eJGny5MlKT09X\nU1OTJk2aNOjv8nhcw+llxDO5PxN7s9t75XQmSpKSk8fGuZoPBfsSlZ7u0sSJ0fs3N3H/nYn+Ri+T\ne4uFsKGek5Oj8vJyLV26VHV1dfJ6vXI6nZKkvr4+FRcXa+/evUpKStIbb7yhxYsXa+/evWpoaFBB\nQYFaW1vV1tYmr9cbtpiWlq7hdzRCeTwuY/sztbe2ti4FAr1ypUp+/4l4lxMSCPTK5+tSMJgYle2Z\nuv9Oo7/Ry+TepNi8YQkb6tnZ2ZoxY4by8/PlcDhUWlqqyspKuVwu5ebmqqCgQHfeeacSEhI0ffp0\nzZ8/X36/X0VFRVq2bJksy9KGDRs49Q4AQIxFlLSFhYX9pjMzM0M/L168WIsXL+63PDk5Wdu2bYtC\neQAAIFLcUQ4AAEMQ6gAAGIJQBwDAEIQ6AACGINQBADAEoQ4AgCEIdQAADEGoAwBgCEIdAABDEOoA\nABiCUAcAwBCEOgAAhiDUAQAwBKEOAIAhCHUAAAxBqAMAYAhCHQAAQxDqAAAYglAHAMAQhDoAAIYg\n1AEAMAShDgCAIQh1AAAMQagDAGAIQh0AAEMkRLJSWVmZamtrZbPZVFxcrKysrNCyXbt2affu3XI4\nHJo+fbpKS0vDjgEAANEXNtRramrU0NCgiooK1dfXq6SkRBUVFZKk48eP6/nnn9fTTz8tu92ur371\nq3r99dd18uTJc44BAACxEfb0e3V1tXJzcyVJU6dOVWdnp/x+vyRp3LhxevLJJ2W329XT06Pu7m6l\np6cPOgYAAMRG2FD3+Xxyu92h6bS0NPl8vn7rbN++XQsWLFBeXp4uu+yyiMYAAIDoGvIH5SzLOmve\nPffco//+7//W/v37dejQoYjGAACA6Ap7TT0jI6PfUXZzc7M8Ho8kqaOjQ++8845mz56txMRE3XDD\nDTp06NCgYwbj8bjOp4dRw+T+TOzNbu+V05koSUpOHjuksV2d7ZIkV2pa1OsK9iUqPd2liROj929u\n4v47E/2NXib3FgthQz0nJ0fl5eVaunSp6urq5PV65XQ6JUl9fX0qLi7W3r17lZSUpDfeeEOLFy9W\nWlraOccMpqWla/gdjVAej8vY/kztra2tS4FAr1ypkt9/YkhjA4FeSZLdMbRxkW7b5+tSMJgYle2Z\nuv9Oo7/Ry+TepNi8YQkb6tnZ2ZoxY4by8/PlcDhUWlqqyspKuVwu5ebmqqCgQHfeeacSEhI0ffp0\nzZ8/X5LOGgMAAGIrou+pFxYW9pvOzMwM/bx48WItXrw47BgAABBb3FEOAABDEOoAABiCUAcAwBCE\nOgAAhiDUAQAwBKEOAIAhCHUAAAxBqAMAYAhCHQAAQxDqAAAYglAHAMAQhDoAAIYg1AEAMAShDgCA\nIQh1AAAMQagDAGAIQh0AAEMQ6gAAGIJQB6IsGAzK331MwWAw3qUAuMgQ6kCUBfydevPoy+oJdMW7\nFAAXGUIdiIGxSePiXQKAixChDgCAIQh1AAAMQagDAGAIQh0AAEMQ6gAAGCIhkpXKyspUW1srm82m\n4uJiZWVlhZYdOHBAW7ZskcPh0BVXXKGNGzfq4MGDWr16taZNmybLspSZmal169bFrAkAABBBqNfU\n1KihoUEVFRWqr69XSUmJKioqQsvXr1+vp556Sl6vV6tXr9b+/fs1btw4zZ49W1u3bo1p8QAA4ENh\nT79XV1crNzdXkjR16lR1dnbK7/eHlu/evVter1eS5Ha71dHRIUmyLCsW9QIAgHMIG+o+n09utzs0\nnZaWJp/PF5pOSUmRJDU3N6uqqkrz5s2TJNXX12vVqlVavny5qqqqol03AAD4iIiuqZ9poCPw1tZW\n3XvvvdqwYYPGjx+vKVOmqKCgQHl5eWpsbNSKFSv04osvKiFh8F/n8biGWs6oYnJ/JvZmt/fK6UyU\nJCUnj414XLAvUQkOh5KSEoc0bijbT093aeLE6P2bm7j/zkR/o5fJvcVC2FDPyMjod2Te3Nwsj8cT\nmu7u7tbdd9+toqIizZkzR5Lk9XqVl5cnSZo8ebLS09PV1NSkSZMmDfq7WlrMvVe2x+Mytj9Te2tr\n61Ig0CtXquT3n4h4XCDQq1N9ferp6R3SuKFs3+frUjCYGJXtmbr/TqO/0cvk3qTYvGEJe/o9JydH\n+/btkyTV1dXJ6/XK6XSGlm/atEl33XWXcnJyQvP27t2r8vJySR8cxbe1tYWuuwMAgNgIe6SenZ2t\nGTNmKD8/Xw6HQ6WlpaqsrJTL5dL111+vPXv26L333tOuXbtks9m0aNEifelLX1JhYaGWLVsmy7K0\nYcOGsKfeAQDA8ESUtIWFhf2mMzMzQz+/8cYbA47Ztm3bMMoCAABDxR3lAAAwBKEOAIAhCHUAAAxB\nqAMAYAhCHQAAQxDqAAAYglAHAMAQhDoAAIYg1AEAMAShDgCAIQh1AAAMQagDAGAIQh0AAEMQ6gAA\nGIJQBwDAEIQ6AACGINQBADAEoQ4AgCEIdQAADEGoAwBgCEIdAABDEOoAABiCUAcAwBCEOgAAhiDU\nAQAwREIkK5WVlam2tlY2m03FxcXKysoKLTtw4IC2bNkih8OhK664Qhs3bgw7BgAARF/YUK+pqVFD\nQ4MqKipUX1+vkpISVVRUhJavX79eTz31lLxer1avXq39+/crKSlp0DEAACD6wp5+r66uVm5uriRp\n6tSp6uzslN/vDy3fvXu3vF6vJMntdqujoyPsGAAAEH1hQ93n88ntdoem09LS5PP5QtMpKSmSpObm\nZlVVVWnevHlhxwAAgOiL6Jr6mSzLOmtea2ur7r33Xm3YsEHjx4+PaMxAPB7XUMsZVUzuz8Te7PZe\nOZ2JkqTk5LERjwv2JSrB4VBSUuKQxg1l++npLk2cGL1/cxP335nob/QyubdYCBvqGRkZ/Y6ym5ub\n5fF4QtPd3d26++67VVRUpDlz5kQ05lxaWrqGVPxo4vG4jO3P1N7a2roUCPTKlSr5/SciHhcI9Ork\nyVNqbfEpJdUruz26XzIJBHrl83UpGEyMyvZM3X+n0d/oZXJvUmzesIR9tcnJydG+ffskSXV1dfJ6\nvXI6naHlmzZt0l133aWcnJyIxwCm6z3eo7/6X1XA3xnvUgBcRMIeqWdnZ2vGjBnKz8+Xw+FQaWmp\nKisr5XK5dP3112vPnj167733tGvXLtlsNi1atEi33XabPv3pT/cbA1xsEpPGxbsEABeZiK6pFxYW\n9pvOzMwM/fzGG28MOKaoqGgYZQEAgKHijnIAABiCUAcAwBCEOgAAhiDUAQAwBKEOAIAhCHUAAAxB\nqAMAYAhCHQAAQxDqAAAYglAHAMAQhDoAAIYg1AEAMAShDgCAIQh1AAAMQagDAGAIQh0AAEMQ6gAA\nGIJQBwDAEIQ6AACGINQBADAEoQ4AgCEIdQAADEGoAwBgCEIdAABDEOoAABgiIZKVysrKVFtbK5vN\npuLiYmVlZYWW9fb26tvf/rb++te/avfu3ZKkgwcPavXq1Zo2bZosy1JmZqbWrVsXmw4AAICkCEK9\npqZGDQ0NqqioUH19vUpKSlRRURFavnnzZs2cOVP19fX9xs2ePVtbt26NfsUAAGBAYU+/V1dXKzc3\nV5I0depUdXZ2yu/3h5YXFRXpxhtvPGucZVnRqxIAAIQVNtR9Pp/cbndoOi0tTT6fLzSdlJQ04Lj6\n+nqtWrVKy5cvV1VVVRRKBQAAg4nomvqZIjkCnzJligoKCpSXl6fGxkatWLFCL774ohIShvzrAABA\nhMKmbEZGRr8j8+bmZnk8nkHHeL1e5eXlSZImT56s9PR0NTU1adKkSYOO83hckdQ8apnc33B6CwaD\nam9vj2I10WG3n1RSUqIkKTl5bMTjgn2JctgdSnDY5XQmDmlspNtPT3dp4sTo/T2Z/Lcp0d9oZnJv\nsRA21HNyclReXq6lS5eqrq5OXq9XTqez3zqWZfU7gt+7d68aGhpUUFCg1tZWtbW1yev1hi2mpaXr\nPFoYHTwel7H9Dbe3trZWPffSm3KmpEaxquHzHX1fKePTlTo+TX7/iYjHBQK96gv26VRfUIFAr+yO\nyMdGun2fr0vBYGJUtmfy36ZEf6OZyb1JsXnDEjbUs7OzNWPGDOXn58vhcKi0tFSVlZVyuVzKzc3V\nXXfdpaNHj+of//iHFi1apJUrVyovL0+FhYVatmyZLMvShg0bOPWOQTlTUpXimhDvMvrxdx+LdwkA\nMCQRJW1hYWG/6czMzNDPTz755IBjtm3bNoyyAADAUHFHOQAADEGoAwBgCEIdAABDEOoAABiCUAcA\nwBCEOgAAhiDUAQAwBKEOAIAhCHUAAAxBqAMAYAhCHQAAQxDqAAAYglAHAMAQhDoAAIYg1AEAMASh\nDkRRMBiUv/uYLCvelQC4GBHqQBQF/J168+jLOnXqZLxLAXARItSBKBubNC7eJQC4SBHqAAAYglAH\nAMAQhDoAAIZIiHcBAIYmGAyqvb09atuz23vV1tYVlW1NmJAmu51jBSBeCHVglOkJdOm/fueTOz0j\nKttzOhMVCPQOezuB7k4tyc2S2z0xClUBOB+EOjAKOZNTleKaEJVtJSePld1xIirbAhBfnCcDAMAQ\nEYV6WVmZ8vPztWzZMr355pv9lvX29mrNmjVasmRJxGMAAED0hQ31mpoaNTQ0qKKiQt/97ne1cePG\nfss3b96smTNnDmkMAACIvrChXl1drdzcXEnS1KlT1dnZKb/fH1peVFSkG2+8cUhjAABA9IUNdZ/P\nJ7fbHZpOS0uTz+cLTSclJQ15DAAAiL4hf1DOOo/HT53PGAAAMDRhv9KWkZHR7yi7ublZHo8n6mMk\nyeNxhV1nNDO5v+H0Zrf3yulMVHLy2ChWNHxJSYlyJIyRpIhrC/YlKsHhUILDrgSHPSZ9na4rmtuN\nxraCfYlKT3dp4sSR93du8v89yez+TO4tFsKGek5OjsrLy7V06VLV1dXJ6/XK6XT2W8eyrH5H45GM\nGUhLS3TuajUSeTwuY/sbbm9tbV0KBHpH3Hele3p65XDYJEl+f2S1BQK9OtXXp1N9QZ3qC8akr9N1\nRVpTOMnJY6OyrUCgVz5fl4LBxChUFT0m/9+TzO7P5N6k2LxhCRvq2dnZmjFjhvLz8+VwOFRaWqrK\nykq5XC7l5ubqrrvu0tGjR/WPf/xDixYt0sqVK3Xrrbfq05/+dL8xAAAgtiK6o1xhYWG/6czMzNDP\nTz755IBjioqKhlEWAAAYKu4oBwCAIQh1AAAMQagDAGAIQh0AAEMQ6gAAGIJQBwDAEIQ6AACGINQB\nADAEoQ4AgCEIdQAADEGoAwBgCEIdAABDRPRAFwAIJxgMqr29Pd5l9BMMBiUdV1ubP96lnGXChDTZ\n7RxXIbpGfai/9sZbqj98LN5lnGVCsk25N1wX7zKAC6Yn0KX/+p1P7vSMeJcS4jv6vpJdKUpKnhDv\nUvoJdHdqSW6W3O6J8S4Fhhn1oX7i5CklpHjjXcZZTgVb410C4swKWvJ3H5MzOfWiOSJzJqcqxTVy\nAtTffUzOlBQlp4ycmoBYujheaYA46D1+XG+1vaKAvzPepQC4SBDqQAyNTRoX7xIAXEQIdQAADEGo\nAwBgCEIdAABDEOoAABiCUAcAwBCEOgAAhiDUAQAwBKEOAIAhIrpNbFlZmWpra2Wz2VRcXKysrKzQ\nsqqqKm3ZskUOh0M33HCDVq1apYMHD2r16tWaNm2aLMtSZmam1q1bF7MmAABABKFeU1OjhoYGVVRU\nqL6+XiUlJaqoqAgt37hxo3bu3KmMjAzdcccdWrhwoSRp9uzZ2rp1a+wqBwAA/YQ9/V5dXa3c3FxJ\n0tSpU9XZ2Sm//4PHGDY2NmrChAnyer2y2WyaN2+eDhw4IEmyLCuGZQMAgI8KG+o+n09utzs0nZaW\nJp/PN+Ayt9ut5uZmSVJ9fb1WrVql5cuXq6qqKtp1AwCAjxjyo1cHOwI/vezyyy9XQUGB8vLy1NjY\nqBUrVujFF19UQsKof9IrAAAjVtiUzcjICB2ZS1Jzc7M8Hk9oWUtLS2hZU1OTMjIylJGRoby8PEnS\n5MmTlZ6erqamJk2aNGnQ3+XxuIbcwPgJSeqyjR3yuFhLtpLO6ud8+hsthtOb3d4rpzNRyckjaz8m\nJSXKkTBGkiKuLdiXqASHQwkOu+x2mxIc9qj3drquaG4zGtuKRV3DlZSUKCk6/UVTsC9R6ekuTZwY\nndcEXltwWthQz8nJUXl5uZYuXaq6ujp5vV45nU5J0qRJk+T3+3XkyBFlZGTo5Zdf1qOPPqq9e/eq\noaFBBQUFam1tVVtbm7xeb9hiWlq6htzAsY4e+QOOIY+Ltb5gT79+PB7XefU3Ggy3t7a2LgUCvbI7\nTkSxquHr6emVw2GTJPn9kdUWCPTqVF+fTvUF5XDYdaovGPXeTtcVaU3hJCePjcq2ol1XNPT09CrF\nlTiiapI++Dvx+boUDCYOe1u8toxesXjDEjbUs7OzNWPGDOXn58vhcKi0tFSVlZVyuVzKzc3V+vXr\nVVhYKEm6+eabNWXKFKWnp6uoqEjLli2TZVnasGEDp94BAIixiJL2dGiflpmZGfr5s5/9bL+vuElS\ncnKytm3bFoXyAABApLijHAAAhiDUAQAwBKEOAIAhCHUAAAzBR9JjJBgMqq2tNTRtt/eqrS2+X80I\nBoOSbLLbbVHd7nB7a29v57bCwAgQDAbV0dEe7zJCznxtmTAhTXY7x6HhEOox4vd36bmX3pQzJVWS\n5HQmKhDojWtNvqPvy56QKHd6RlS3O9zefEffV8r4dLlSo1gUgCHr6Gjv97oVb6dfWwLdnVqSmyW3\ne2K8SxrxCPUYcqakKsU1QdIHN/iI981V/N3H5HCMDdUULcPtzd99LIrVABiOM1+34m0kvG6ONpzL\nAADAEIQ6AACGINQBADAEoQ5ESTAYlL/7mPggP4B44YNyQJQE/J168+jLsifwXhlAfPDqA0TR2KRx\n8S4BwEWMUAcAwBCEOgAAhiDUAQAwBKEOAIAhCHUAAAxBqAMAYAhCHQAAQxDqAAAYglAHYsgKWvJ3\nH1MwGIx3KQAuAoQ6EEMnenr0VtsrCvg7410KgIsAoQ7EGLeOBXChEOoAABgioqe0lZWVqba2Vjab\nTcXFxcrKygotq6qq0pYtW+RwOHTDDTdo1apVYccAAIDoCxvqNTU1amhoUEVFherr61VSUqKKiorQ\n8o0bN2rnzp3KyMjQHXfcoYULF6qtrW3QMQAAIPrChnp1dbVyc3MlSVOnTlVnZ6f8fr+Sk5PV2Nio\nCRMmyOv1SpLmzZun6upqtbW1nXMMYCp/9zGdPNkb7zIAXMTCXlP3+Xxyu92h6bS0NPl8vgGXud1u\ntbS0DDoGAADERkTX1M9kWdaQlw02ZrgS7DYF2o/EbPvnKyF4UoHeD7/GFOxLVCAQ36O4nu4u2RNO\nqLsrup+JYmxKAAAKLElEQVTGHm5vsapruE7X1dXZHra/9vYm/f1vdQqM6ZIk9fYcl32MQwkJYxTo\n6lJLX2PU64rWv1e0/jZH4n7s6e6S3WbJshzxLqWfQHen2tvbo7Itu71XbW1dUdlWe3u7At0j5+uX\np/82R1JNI13YUM/IyOh3lN3c3CyPxxNa1tLSElrW1NSkjIwMjRkz5pxjBuPxuIZUvCTlLfic8oY8\nCgDMMXHixCht6XJdd112lLaFeAh7+j0nJ0f79u2TJNXV1cnr9crpdEqSJk2aJL/fryNHjujUqVN6\n+eWXdf311w86BgAAxIbNiuDc+GOPPaaDBw/K4XCotLRUb731llwul3Jzc/Xqq6/qe9/7niTppptu\n0sqVKwcck5mZGdNGAAC42EUU6gAAYOTjjnIAABiCUAcAwBCEOgAAhhjy99TP144dO7R3716NGTNG\n69ev11VXXaW3335bGzZskN1uV2ZmptavXy9JeuKJJ7Rv3z7Z7XatWrVK8+bNU3d3t4qKitTV1aXk\n5GQ9+uijSk1NvVDlR8Tn8+mLX/yifvCDH2jWrFnG9NfX16eSkhK99957CgaDevDBB3X11Vcb099g\nRvMzDDZv3qxDhw6pr69P99xzj7KysvStb31LlmXJ4/Fo8+bNGjNmjPbs2aOnnnpKDodDt912m5Ys\nWaJTp05p7dq1OnLkiBwOh8rKynTZZZfFu6WznDhxQjfffLPuu+8+XXfddUb1t2fPHu3YsUMJCQn6\n+te/rszMTGP6CwQCWrNmjY4dO6aTJ0/qvvvu05VXXjnq+3v77bd1//33a+XKlVq+fLmOHj067J7O\n9Tp7TtYF8Je//MW69dZbrWAwaL311lvW448/blmWZd15553WH//4R8uyLKuwsNDav3+/1djYaN1y\nyy3WqVOnrNbWVuumm26ygsGg9fjjj1s7duywLMuynnnmGeuRRx65EKUPyYMPPmjdcsst1sGDBy3L\nMqe/3bt3W+vXr7cs64N9uWTJEsuyzOnvXA4ePGj9y7/8i2VZlvXXv/7V+spXvhLniiJ34MAB6+67\n77Ysy7La29utG2+80Vq7dq31q1/9yrIsy3rsscesp59+2goEAtbChQut7u5u6/jx49bNN99sHTt2\nzKqsrLS+853vWJZlWa+88or1wAMPxK2XwTz22GPWkiVLrMrKSmvt2rXWvn37QvNHc3/t7e3WggUL\nrEAgYLW0tFjf/va3jervJz/5ifXYY49ZlmVZTU1N1k033TTq/z4DgYC1cuVKa/369dZPfvITy7Ks\nqOyzgV5nB3NBTr//+te/Vl5enmw2mz71qU+poKBAJ0+e1OHDhzVjxgxJ0vz581VVVaXf//73uuGG\nG+RwOOR2uzVp0iT95S9/0YEDB/SFL3xBkvT5z39eVVVVF6L0iB04cEAul0uf/OQnJcmo/r785S/r\noYcekvTBrYBPv7t+//33jejvXM713IPRYNasWdq6daskKTU1VYFAQDU1NZo/f76kD/dBbW2tZs6c\nqeTkZI0dO1ZXX321/vCHP/Tr/XOf+5wOHToUt17O5d1339Xf/vY3zZs3T5ZlqaamRp///Ocljf7+\nqqqqlJOTo6SkJKWnp+s73/mODh48aEx/brc7dEe9Y8eOye12j/q/z7Fjx+rf//3flZ6eHpo3nH32\n2muvnfN1djAXJNQPHz6sI0eO6Gtf+5ruuusuvf3222pvb9f48eND67jdbjU3N6u1tbXffeMnTpwY\nup98WlpaaN5Iupf8yZMn9aMf/UgPPPBAaJ5J/SUkJGjs2LGSpP/8z//UokWL1N7ergkTJoTWGc39\nnctofoaB3W5XUlKSJOm5557TjTfeqJ6eHo0ZM0bSB/tgoP010PMbbDab7Ha7Tp06deEbGcTmzZu1\ndu3a0LRJ/R0+fFg9PT269957dccdd6i6ulrHjx83pr+8vDwdPXpUCxYs0IoVK7RmzZpRv//sdrsS\nExP7zRtOTzabTT6f76zX2TPv4jqQqF9Tf/bZZ/Xcc8/JZrNJ+uC+762trZo7d66eeOIJ/eEPf9C6\ndev0wx/+MKJ7wgeDwbPmRTIuVs7sz7Is2Ww2XX/99Vq2bJlSUlL6rWtKf/fff79ycnL005/+VG+9\n9Za2bdum1tbWiLY30vobjtFY90svvaTdu3drx44dWrBgQWj+uXo51/yB9mM8/exnP9OsWbP0sY99\nbMDlo70/y7LU0dGhH/zgBzp8+LBWrFjRr/bR3t+ePXt0ySWXaPv27frzn/+skpKSfstHe38DGWpP\np19/h/q6E/VQv+2223Tbbbf1m1deXq5PfOITkqRrrrlGR44c0cSJE9XR0RFap6mpSV6vVxkZGXr3\n3XcHnO/z+ZSSkhK6x3w8DNTfsmXL9Morr+jJJ5/Ue++9pzfffFPf+973dOzYsdA6o7k/6YOwf/nl\nl/XDH/4wdGr9zAdSjJb+hmKw5x6MBr/97W+1fft27dixQykpKUpOTlZvb68SExP77ZePPr8hOzs7\n1HtmZmboCCgh4YJ9rjas3/zmN3r//ff1wgsvqKmpSWPGjJHT6TSmv/T0dGVnZ8tut2vy5MlKTk5W\nQkKCMf0dOnRIc+fOlSRlZmaqqalJSUlJxvR32nD+z1n//8N1H83JcK+dF+T0+9y5c/Xb3/5WklRf\nX69LLrlEDodDn/jEJ0LXQl544QXNnTtX1157rX7zm9/o1KlTampqUnNzs6688krl5OTo+eef77fu\nSPH000+roqJCzzzzjG688UatX79e06dP1xVXXGFEf42NjXrmmWdUXl4eOpWUkJBgzP47l9H8DIPu\n7m498sgj2rZtm1yuDx6UNGfOnFA/+/bt09y5czVz5kz98Y9/VHd3t/x+v1577TVdc801ysnJ0a9+\n9StJ0v/8z//o2muvjVsvA9myZYueffZZPfPMM1qyZInuu+8+zZkzJ1TzaO8vJydHv//972VZ1gdP\nTgsEjOpvypQpev311yV9cKnB6XTqc5/7nDH9nTbc/3PnysnBXLDbxD7++OP63e9+J0l66KGH9JnP\nfEb19fUqLS2VZVn6zGc+ozVr1kiSfvrTn2rPnj2y2Wz6xje+oWuvvVaBQEDf+ta31NHRodTUVD3y\nyCNnne4eCR566CHdcsstmjVrljH9bdmyRb/85S916aWXhk4J7dy5Uw0NDUb0N5jR+gyDXbt2qby8\nXJdffnlonz388MMqKSlRb2+vPvaxj6msrEwOh0MvvPCCnnjiCdntdt1555360pe+pGAwqJKSEjU0\nNGjs2LHatGmTvF5vvNsaUHl5uS677DJdf/31evDBB43pb9euXXr22Wdls9m0atUqXXXVVcb0FwgE\nVFxcrNbWVvX19emBBx7QFVdcoTVr1oza/mpra7Vu3Tq1tbXJ4XBo/Pjx2rFjh9auXTusns6VI+fC\nvd8BADAEd5QDAMAQhDoAAIYg1AEAMAShDgCAIQh1AAAMQagDAGAIQh0AAEMQ6gAAGOL/AVAxJn/u\nuO5EAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f94453cf290>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_hist(ex2_data.values,label='Truth')\n",
    "my_hist(mu_samp_2,label='Est')\n",
    "plt.legend()\n",
    "print(np.mean(mu_samp_2))\n"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}