Skip to content

Instantly share code, notes, and snippets.

@joezuntz
Created November 25, 2014 13:08
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save joezuntz/e6d216753d5421f2b2d3 to your computer and use it in GitHub Desktop.
Save joezuntz/e6d216753d5421f2b2d3 to your computer and use it in GitHub Desktop.
Lecture 2 Hard problem
Display the source blob
Display the rendered blob
Raw
{
"metadata": {
"kernelspec": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"display_name": "IPython (Python 2)",
"language": "python",
"name": "python2"
},
"name": "",
"signature": "sha256:88e44dc88edd32f3859751042cf13b0b95258e096edd0d3047b20e44ac570259"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"pylab inline"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First let's generate some simulated data for our problem. This is how I made the data for you to download."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#Simulation input values\n",
"a = 1.0\n",
"b = 3.0\n",
"sigma = 1.5\n",
"s = 0.5\n",
"#Choose some observed points on the x axis\n",
"M_true = np.random.uniform(2.0, 8.0, size=50)\n",
"\n",
"#Generate the intrinsic scatter\n",
"L_true = np.random.normal(a*M_true+b, sigma)\n",
"\n",
"#Add noise from observations\n",
"L_obs = np.random.normal(L_true, s)\n",
"\n",
"#Save the data so the students can use it\n",
"savetxt(\"simulated_data.txt\", transpose([M_true, L_true]), header=\"M_i L_i_obs s\")\n",
"\n",
"#Plot\n",
"errorbar(M_true, L_true, s, fmt='.')\n",
"M_line = arange(0,10,1.0)\n",
"plot(M_line, a*M_line+b, 'r-')\n",
"xlim(0,10)\n",
"ylim(0,15)\n",
"xlabel(\"M\")\n",
"ylabel(\"L\")"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 2,
"text": [
"<matplotlib.text.Text at 0x1095b1890>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHNRJREFUeJzt3XuQpVV57/HvjwHkMjkCSmEETgb0TIVQeBkIEBOkiWKB\nEqHKEwQSbkaPZUAHSIxgEmnqnCQYAyomWhoEZbjqoBwhkkQNexAwBJFbuBxPVAiD5QzioEEkgPPk\nj/ft6bd37717X97rfn+fKqr33r0vq7uZ9bzrWWs9SxGBmZm1z1ZVN8DMzKrhAGBm1lIOAGZmLeUA\nYGbWUg4AZmYt5QBgZtZSW1fdgF4keW2qmdkYIkLDPreWAQBG+yGmmaTZiJituh114N/FPP8u5vl3\nMW/Ui2engMzMWsoBwMyspRwA6q9TdQNqpFN1A2qkU3UDaqRTdQOaSnWsBSQpPAdgZjaaUftOjwDM\nzFrKAcDMrKUcAMzMWsoBwMyspQoLAJIukbRB0n09vvcHkjZL2qWozzczs8GKHAFcChzR/aCkPYHD\ngUcK/GwzM1tCYQEgIr4ObOrxrQuBPyrqc83MbDilzgFIOhpYHxH3lvm5Zma2WGnF4CTtALyfJP2z\n5eEBz5/N3O1ERKeYlpmZNZOkGWBm7NcXuRNY0grg+ojYT9J+wFeBp9Nv7wE8BhwYERu7XuedwGZm\nIxq17yxtBBAR9wG7zd2X9D1g/4j4UVltMDOzeUUuA70KuA1YKelRSad2PaV+RYjMzFrExeDMzKaE\ni8GZmdlQHADMzFrKAcDMrKUcAMzMWsoBwMyspRwAzMxaygHAzKylStsJbGY2iMQM83VtZoBOersT\nseW25cgbwcysdiQiIikW6cAwvFH7TgcAM6udbAAY5nFLeCewmZkNxQHAzKylHADMzFrKAcDMakXi\nU+nXL0vsVHV7ppkDgJnVzcr065HAJ6tsyLRzADCzupk7NvYO4J1zD3pkkD8HADOrmxPSr2+I4MnM\n4x4Z5Mw7gc1soLI3YkXwpJR87fpWz5GBjc8bwcxsaGVtxOr1OWnaZxOwc4/gYIzed3oEYGa10DXS\nQGI2vdmJoDNgZGBj8gjAzIZW5QigzM9vqlqNACRdArwJ2BgR+6WPfQg4CngW+A5wakT8uMh2mFn9\nLTUCKL9F06/QEYCkQ4CngMsyAeBw4GsRsVnS+QARcXbX6zwCMKshjwDqrVYjgIj4uqQVXY99JXP3\nduAtRbbBzJqta2SwziOD/FQ9Cfw24KqK22BmQ8hsxPpnkqWgz1DOstBO5jPqR1oG/BoRt1TdlFFV\nFgAk/THwbERc2ef7s5m7nYjolNEuM+trbiPWQcAjEbw1TcnM5Pkh2R2/wAm1XfUjLQdOAVYDP0I6\njIinB78o7yZoBsb//Re+CihNAV0/NweQPnYK8A7gdRHxTI/XeA7ArGbSDvlIko1Yb0iXZeaek5fo\nAIemdz8XwVvzfP+JSXsApwNvB9YBFwK3UYMllbU/EEbSEcB7gaN7df5mVlv9SjTkrZ47fqX9ka4A\n7gW2Aw4k4i1E3FqHzn8cRa8Cuookkr8Y2ACcC5wDbAv8KH3aNyLi97te5xGAWUlGKfXQfcVf0Aig\nPjt+k/z+UcBZwF7ARcDFRNQyLeUzgc1sbEt16GUEgCLfd4QGZPP7m4ALgC8Q8VxlbRpC7VNAZjYd\nprI8s7QHyf6kh4HDSILAQURcU/fOfxwOAGY2rukpz7wwv789Saff6Pz+MKreB2BmNTdgI9Z26dd6\nTdYOq3d+/7S65veL4DkAM9tilNx7kZO1hc4BLM7vXwhcOw0pnlqVgjCz5hh1A1be5ZkLL/mweP3+\nKdRk/X5VPAIwMwAk7gZemd59HPh4ertvB1z5ap1hSPsDZwJvBNYAFxHxnWobVQwvA7XWK/sIw2nR\ntdP3V4fp2IcNAKX/TRq2fj8vDgBmGY24Qq2JbE4f2JRnAJj0NSO8+Y4kqZ0zmLL8/jAcAMwyHABG\nM/f7GvR7m/RqvpC/SY3r85TJAcAswwFgNMMEgLw+I6c3W0WS5pn6/P4wvBPYzKabtAzpaKR1wHXA\nPcDeRKxuc+c/Di8DNbMF+i0HrXxyveX5/SI4BWRTrS0poLw6Z4kgyaEPrMc/7u81DS7vAG5k2MNe\nnN8fmucAzFJjdTZTYNTOuUfw2B14OfAQ8Gu9fm8TBIAOwx720qL1+3lxADBL1f5kqYJMOuoZpsTD\nBAFg0aliXU9o5fr9vDgAmKWW7GymVB5pr1HPBRjhfXsHF+f3c+EAYJaq1clSJapzAFj0Wuf3c+Vl\noGapuU6/TZ1/GXI5CCapv385Lau/XzdeBmo2xQpautl9EMxwcytpfr/DawG+SJLfP935/eo4AJhN\nkR5r+DuknX6aepnJ4WOeTr8OdxBMV/39T/AuDuXmlzm/Xz3PAdhUa8s+gDmDVj4NkdefYYjRwtBz\nK33y+yI2t+lvUqbaTAJLugR4E7AxIvZLH9sFuAb4JZJDl4+NHsM/BwCbROU7VgsyzM81aOVTnsFw\n4Hsl6/fPStuxBrhIxJ5Ltd0mV6cAcAjwFHBZJgD8JfDDiPhLSe8Ddo6Is3u81gHAbIB+HfCgq/NC\nA4DX79dCbQJA2pgVwPWZAPAQcGhEbJD0EqATEb/c43UOAGYDLFGuuV9wyD8ATPH5uk1U92Wgu0XE\nhvT2BmC3kj/fzHKwO+tBOp8klXsYSRA4iIir3fk3R2WrgCIiJPUdfkiazdztRESn8EaZTalRD3wf\n8Eb7A2fdy84wv37f9XkqImkGxl/ZVUUKaCYifiDpF4GbnAIyG92oKaCJ6iKl+f0f8qLztub5va7j\nmNvfy4d2+CG7fjV9hidya2LUvrPsEcCXgJOBD6Zfryv5881qpc/KnhXp/YfpsWKm19V81/usk5jN\nvoZR1+4nH7Agv/9injgfuPaU+Mxzpwz1BlZ3Ra4CuorkiuPFJPn+DwD/F/gc8N/xMlCrUB2Xiva5\ncs/lan6kukiuz9NYtVoFNC4HACtTXTaLjRAAxqpyuuTP2WP9vvP7zVL3FJCZTe4Ekqv5vp1/vxGO\nxMyCEU7v9funef1+O3gEYK3XtBHAoMcHvW9XGebu9fsXAF/wEs5m8wjAKlfH/HpRGvezLs7vn4Lz\n+63lAGC5K6gCZS3l9bP2Wdkz8dr9NEDx2/r8J6/j8k0/4Re+dx3HbPgLztn4EPvcBxwObBPUMFhZ\n4ZwCskLVJb0ySHcbx72qn/CUrA5dK3smqeyZPmkZcNQ6Xnvdodz8KJn6PE34u9jonAIyG0Gvq+yK\nRjC91umPvnYfFuX3P87vu/6+9eQRgBWq7leaS62pH2eydcx2LFqn3/3YkiOTHuv3t+LnJwdbvQO4\nkUwaqe5/FxuP9wFYrUzYKc5Q8ATrUmvqywoA/V4/ZKonu37/MpL1+99NX9+hR4Abt61VTHo3bqK9\nQg4AVhtpemXR1eeY71XIFetSO2RrGwCGrL/fL8Dl8fusYhThkctgdS8Hbe3SfXh47cx1iJMEJ1g4\nl5AGlWJIy5FOBx4C/hj4OPAyIv6qz+atE9KvQ+8YtvZwALAijTeJWRMjdurFBjtpj676+6eSlGK+\nZtDkbq8AV1qwstpzALAiNf3qc5ROvZBgt4o7QboCuBfYDjiQiLcQccsEm7dqPzKzcjgAWGHySq9U\naJROfaxgJzEjMZuWb14nMbuNnvvsafqbr/2bXv7w3/Gm//wMJ+/6SzzyKRHXzU3uTqjRIzPLjyeB\nrVA5TTbmNpmcec8ZllhZMlIJZXL4WQuqz9Njo9tIP1ef98z9bzLk53oSeACvArJaySkAdBj3NKsJ\nlbIKaMT6+0MGr4HPyWHFUoeS/yZVBZ0m8U5gm0bTmbJI1u+fCbyRpP7+UOfrDrNTOfucglTxN+me\nuyjtQmBaeQRghcppBDBxymKCz853BDDk+v0i2pfH6zKvL/1vMu5BOG3iFJBVroidm1Xlfpf63MzP\nuiL9+nB6u5PenivTsCNJfv8Mko7zQuDavPP7Rb8u7/cY8fMquxBoCgcAm0pldjaTBrCug1cmOl83\nz1x+3oHZO4HrxwHAplKT/uFLRKADWJjfn/h83bFrBRXEAaB+GlEKQtI5ku6XdJ+kKyW9oIp2mOVK\nWoZ0dCdZHPNF4B5gbyJWF3G4unf02qRKDwCSVpAs5VoVEfsBy4Djym6HWW6kHZFOI6nP8yef4F2Q\n1Of5UMGHq3tHr02k9BSQpF2AbwAHA/9BcqX00Yj4auY5TgFZ/csAS7szn9+/GfgwcKuIzQVVLg2S\nOkAz6UPvBnYBvg+8M4Ib8v7MHm2YweWga6sRcwCS/hfJLsefAf8QESd2fd8BwOpLWkWyjLNnfj+H\nJZYz9O7wzu3a0Tv3j9erYgxoQACQ9DLgeuAQ4MfA54G1EXFF5jkOAJaL3K4epa2YX7+/N/Ax4G97\n1N/PdbfqXDDp9b5zAcCTojanCTuBDwBui4gnACR9AXgNcEX2SZJmM3c7EdEpq4E2PSY+33fh+v0n\nSZZxrh2wfr+o3areBWuLSJqBEf+fzqgiADwE/Kmk7YFngNcD/9L9pIiYLbldNqGpytUuzu+fCtw6\nxPr9okokTGc5DJtIemHcmbsv6dxRXl/VHMAfAScDm4FvAW+PzBWVU0DNV/R67XGCzZClGgbm94do\nV667VTMpoF6HxpeeApqqID+Faj8HMAwHgOYreefuUJ814HzdofL7ebdn1Pfquj03JwAVTQJ7U1b9\nNGEOwKweFuf3LyCH+jxFS6/C35B56BaJtfgq3EbkE8GsFbK7Zq/S8fsg/QVJsbbfJMnvH0jE1XXq\n/Pvt9E07+QfSu3cAvxHBrDt/G9VYAUDSGXk3xKxgK1/Nt1jD7x55FDfcDewAHJzD+bpFGrTTt+nn\nLVsNjDUHIOnRiNizgPbMvb/nABquNnMAaX7/Ll516Yt4YpfLOOlR4JA/if/zSM5tmCH/EtgL6t8D\nm7o3glWZg6/6822xUiaBHQBskLKP7uvZEXXl99ez+yf24nuffp5t+k6Y1m2FS2blzwdJlkwvaBNw\nkwOAZTkAWOXKPC92UbBBy0nW7/8emfo8REQp5/vmbFA7Ki4F7fN5ayi3ctCSnpL0H73+A16aS2tt\nWpW5aWklwCruPPIfOfxu4F5ge+qf3286VyKdAn2XgUbE8jIbYlPlBJLURbETlNKyt3Hx8pO4jJV8\n+9kX8J+fBj5WcAlmS3hn8hTwRjArRM6boWbI5OZ3ZeOt53Huq07islduy7MbT2TNq7/PS3e9OV77\nw7zaVNcUUF3mKXw+bz15J7AVYtSOp5AOVNrjfN736Nl88AmS83VHqr/flABQl05+KXUJkjbPO4Gt\nEBNX1ZyEtD/p+brb8zOAgxbU5xnwv3tXZ7pOYja9Pahm0JYNWFQwwZn9XZsVySOAFpr0CnPItMS5\nwGFjX7FKy5ivz7MXaX0eEZu6rzrHuRId9DsAZilpFVOTeQRQP04B2UjG7Dx7vibtVG8CzmPc1MUS\n9Xl6fXYOJ3B1B7QFG7Cc4+7NAaB+HABsJHkGgHHfL33hHszX39+S3+9ewpkpjzxDTnnyHgHAE5x9\nNGV+oq0cAGwklQeATH6fJervF7X5qIhRhVkVctsIZtZLvwqVI77JMqSjkdYBXwTuAfYW8UURJ0rM\nSnTSr7PpVSd485FZrjwCaLlRr3SXKvOwRGG2bH5/E8n5uj3r76enXZ2X3p0hSTX8DvBycs7NewRg\n08IpIBvJGAFg4ARpn8Js3fn9C4HbBpVo6HUSVhG5+X5pJQcAayKngGxoY6Zz+tahX/R+0v5IlzNf\nn+egtD7PMIerLzL3eTlPzDqtZK3ljWDt1t35LbnePYInpb6d8Mqt+DlHccOR/5s//TZJCeOLgNP7\n1ecZsKqkLNmaNmsyG8WG2jQ2Cq+gsbpxCqjFxl3v3q/+/nl84PYTuHLfZ9jupzuz6T178NiaUY5Y\nHHAAes/beahqyadTTFaERswBSNoJuBjYFwjgbRHxz5nvOwCUYNzOb0HnlcnvP8MLbnsdX/ut+9l3\n5ydjp5E700yuf0FensxJWEV0nFV0xg4AVoSmzAF8FPhyROwDvAJ4sKJ2tNpEOfUe+f3t4pk338av\nM07n3yWbmroz+bgJlp2aWU+lzwFIeiFwSEScDBARzwM/LrsdNoa0Pk+H10Kyfn9gfn8C2bz8s8De\nJMHgFgrIzZu1VRWTwHsBj0u6FHglyRXe6oh4evDLrDLSjv/EYX/2Cl504k/Z8ZmrOe7BI7lx7c/Y\nYTnwqsh/4nbLgTLAleljhdblyQSVGTw5ay1R+hyApAOAbwCviYg7JH0E+ElEfCDzHM8BlGSJjVsT\nrd+fpC1Frv1P33+G3ityzi0jN+85ACtCE84DWA+sj4g70vtrgbO7nyRpNnO3ExGd4pvWDkvWyE/q\n85xFknZZQ3f9/cXvcQxs6ZzvGSdF012Df+7xJZadjq1fzX2Jc/P8nF6qPm/ApoekGeb/HY7++opW\nAd0MvD0ivp129NtHxPsy3/cIoGyL6+9fBFw8TH4/j6vZ7hITwLFFrvwZ0I7CP2upchpm42rCCADg\n3cAVkrYFvgOcWlE7TFpOUp9nNUvU5ylY9yHjx5b8+WXygepWC94IVnOF7R4dI78/oI15jAAW5PqL\n3Py1RDvKGAH4vAErRCM2gi3FAaC3XDqnxfn9vvX3S20XfSuAQkkTs3Nt8CSwNZUDwBSb4LStsfP7\nhbZrwPtUUTenqANn+nyWA4DlzgFgivToBA8luUIerhMsKb9fRACoQpmTs1X/rDadmjIJbEPILlVM\n0yNEbFli2fcq+X+y9v7P89v7M5/fP4Ux8/tLmbIljZ6ctVbxCKAhMgGg/1m86AByzu8P0a4OE1w1\n16lEctGTs3X6WW06eQTQNr3r85xWQH2efia6au63IasKRW06y7x/h5r8rGbgEUBjzI0AmJugRM+T\nye+/lat/9RqO27bs9fvTtqTRuXlrMk8CT6m5ALA76/lz3v/QSazZlcz6fRGbq+q4pqnTnKafxdrH\nKaApJPGpVdzJWVzIUdzw82X8vAN8iIjvzj+psuaZWUM5ANRZmt//JquO3ZXHuYj3cBYX3rghdnvX\nwqdN1Uqc0i1ZHM9sSjkFVEdd6/dP52PbfJJ3vup5toEeufaqi4s5bWJWD54DaLI+9XlEvJBkorXn\nMtBxD3efrKle0mhWNw4ATbSwPs9lJOv3v7vwKf33AUzbShwzG48DQFOMWJ9nqI1gTsOYtZpXAdXd\niPV5PEFpZkXxCKAsE9bfX+oK3yMAMxu179yqyMYYSX5fugK4F9ie5HzdtxBxaxHF2czMhuURQBFy\nqr8/ykobjwDMzJPAVZJ2JMnvn0EJ5+t6KaaZZTkAVCHH83XNzMblOYAyJfn9y3F+38waqLIAIGmZ\npLskXV9VG8YiLUM6GmkdSf39u4G9iVhd9OErZmZ5qnIfwGrgAeAXKmzD8JzfN7MpU8kcgJKc+WeA\nPwPOiojf6vp+feYAapDf9wofMxtGU+YAPgy8F9hc0ecvTVrl/L6ZTbPSU0CSjgI2RsRdkmYGPG82\nc7cTEZ2Cm9a9fn9vkvX7p5d4vq6Z2dDSPnRm7NeXfTEr6c+BE4Hnge2A/wZcGxEnZZ5TbgpoYX7/\nSeACCszvj8opoNF5DsXaqFH7ACQdCvxhZXMAi/P7HwZql+JxAJiMf3/WFk2ZA8iqYha63/r9W+rW\n+ZuZFaXSABAR6yLizaV82OL1+/fQgPX72fN+04NfzMxyMf3nAZS8fr8AK9OvRwKfhHLP+zWz6TW9\nAWBxfv8Umlmf5+n06x3AO6tsiJlNl+krBpecr3sm8EZgDcn5urVN8SzF5/1OzpPA1hZNnASe3BTX\n55nr9N35j8dzKGb9NTsF1Pz8vhXPcyhmfTQzAExPft+K5zkUsz6aNQcwZfn9YTmHPT7PoVibNGon\ncD8LfoicztdtGpcyyI8DqLXF9AQAWI7z+5YDBwBri1EDQJ3nAB7B+X0zs8LUOQAc1Ib8vplZVeqb\nAqrLiWDWSJ5DsTaanjkABwAzs5G0cyewmZmNzAHAzKylHADMzFrKAcDMrKUcAMzMWsoBwMyspRwA\nzMxaygHAzKylKgkAkvaUdJOk+yX9q6T3VNEOM7M2q2QnsKSXAC+JiLslLQfuBI6JiAfT73snsJnZ\niBqxEzgifhARd6e3nwIeBF5aRVvMzNqq8jkASSuAVwO3V9sSM7N2qbQcdJr+WQusTkcC2e/NZu52\nIqJTYtPMzGpP0gzzVW9Hf31V1UAlbQPcANwYER/p+p7nAMzMRtSIctCSBHwWeCIizuzxfQcAM7MR\nNSUA/AZwM3AvMNeAcyLi79PvOwCYmY2oEQFgKQ4AZmaja8QyUDMzq54DgJlZSzkAmJm1lAOAmVlL\nOQCYmbWUA4CZWUs5AJiZtZQDgJlZSzkAmJm1lAOAmVlLOQCYmbWUA4CZWUs5AJiZtZQDgJlZSzkA\nmJm1lAOAmVlLOQCYmbWUA4CZWUs5AJiZtZQDgJlZS1USACQdIekhSf9f0vuqaIOZWduVHgAkLQP+\nGjgC+BXgeEn7lN2OppA0U3Ub6sK/i3n+Xczz72J8VYwADgT+LSIejojngKuBoytoR1PMVN2AGpmp\nugE1MlN1A2pkpuoGNFUVAWB34NHM/fXpY2ZmVqIqAkBU8JlmZtZFEeX2x5IOBmYj4oj0/jnA5oj4\nYOY5DhJmZmOICA373CoCwNbA/wNeB3wf+Bfg+Ih4sNSGmJm13NZlf2BEPC/pdOAfgGXAp935m5mV\nr/QRgJmZ1UPtdgJ7k1hC0p6SbpJ0v6R/lfSeqttUJUnLJN0l6fqq21I1STtJWivpQUkPpPNqrSTp\nnPTfyH2SrpT0gqrbVBZJl0jaIOm+zGO7SPqKpG9L+kdJOw16j1oFAG8SW+A54MyI2Bc4GDitxb8L\ngNXAA3gVGcBHgS9HxD7AK4BWplAlrQDeAayKiP1IUsrHVdmmkl1K0ldmnQ18JSJWAl9L7/dVqwCA\nN4ltERE/iIi709tPkfwjf2m1raqGpD2ANwIXA0OvcJhGkl4IHBIRl0AypxYRP664WVX5CcmF0g7p\n4pIdgMeqbVJ5IuLrwKauh98MfDa9/VngmEHvUbcA4E1iPaRXOq8Gbq+2JZX5MPBeYHPVDamBvYDH\nJV0q6VuS/lbSDlU3qgoR8SPgAuDfSVYUPhkRX622VZXbLSI2pLc3ALsNenLdAoCH910kLQfWAqvT\nkUCrSDoK2BgRd9Hyq//U1sAq4OMRsQr4KUsM86eVpJcBZwArSEbHyyX9TqWNqpFIVvgM7FPrFgAe\nA/bM3N+TZBTQSpK2Aa4FLo+I66puT0VeA7xZ0veAq4DflHRZxW2q0npgfUTckd5fSxIQ2ugA4LaI\neCIinge+QPL/S5ttkPQSAEm/CGwc9OS6BYBvAv9D0gpJ2wJvBb5UcZsqIUnAp4EHIuIjVbenKhHx\n/ojYMyL2Ipng+6eIOKnqdlUlIn4APCppZfrQ64H7K2xSlR4CDpa0ffrv5fUkCwXa7EvAyentk4GB\nF46lbwQbxJvEFvh14HeBeyXdlT52TkT8fYVtqgOnCeHdwBXpRdJ3gFMrbk8lIuKedDT4TZL5oW8B\nn6q2VeWRdBVwKPBiSY8CHwDOBz4n6feAh4FjB76HN4KZmbVT3VJAZmZWEgcAM7OWcgAwM2spBwAz\ns5ZyADAzaykHADOzlnIAMFuCpM2S1mTuby3pcZemtqZzADBb2k+BfSVtl94/nKQkgzfRWKM5AJgN\n58vAm9Lbx5PUJXJxOms0BwCz4VwDHJeeOLUf7S3NbVPEAcBsCBFxH0nZ4eOBv6u2NWb5qFUxOLOa\n+xLwVyQFuHatuC1mE3MAMBveJcCmiLhf0kzVjTGblAOA2dICICIeA/4685hXAVmjuRy0mVlLeRLY\nzKylHADMzFrKAcDMrKUcAMzMWsoBwMyspRwAzMxaygHAzKylHADMzFrqvwBlzuaw8PXEDwAAAABJ\nRU5ErkJggg==\n",
"text": [
"<matplotlib.figure.Figure at 0x1095bb690>"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's try to recover the underlying parameters, $a$, $b$, and $\\sigma$ - we will assume we know the measurement errors; this is often true. Let's collect our parameters together as a vector to avoid having to keep writing them out:\n",
"\\begin{equation}\n",
"p \\equiv \\left( a, b, \\sigma\\right)\n",
"\\end{equation}\n",
"\n",
"\n",
"Our model is first that our data points are independent:\n",
"\\begin{equation}\n",
"P(L^\\mathrm{obs}|p) = \\Pi_i P(L_i^\\mathrm{obs}|p)\n",
"\\end{equation}\n",
"and that the intrinsic scatter and noise are both Gaussian:\n",
"\\begin{equation}\n",
"L^\\mathrm{obs}_i \\sim N(L_i, s_i^2)\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"L_i \\sim N(\\mu(M_i), \\sigma^2)\n",
"\\end{equation}\n",
"\n",
"So that the likelihood for a single point is:\n",
"\\begin{equation}\n",
"P(L_i^\\mathrm{obs}|p) = \\int \\mathrm{d}L_i P(L_i^\\mathrm{obs}|p L_i) P(L_i|p)\n",
"\\end{equation}\n",
"\n",
"and the total likelihood for all of them is the product of the $L_i$, since the data points are independent. It's nearly always easier and more accurate to work with log probabilities, since then you can add them together instead of multiplying small numbers. Since logarithms are monotonic increasing you can also always maximize a log-likelihood to get the same answer as maximimizing the likelihood.\n",
"\n",
"If you write in the definition of the Gaussians on the integral and complete the square you can integrate out the dependence on $L_i$ analytically. Because this is a very simple example everything is Gaussian and this is possible; often your likelihoods are more complex and you do your integrals numerically.\n",
"\\begin{equation}\n",
"P(L_i^\\mathrm{obs}|p) \\propto \\frac{1}{(s^2+\\sigma^2) \\sigma s} \\exp{-0.5\\left(\\frac{(L_i^\\mathrm{obs}-(aM_i+b))^2}{s^2+\\sigma^2}\\right)}\n",
"\\end{equation}\n",
"\n",
"\n",
"Let's write a function to compute this likelihood"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#One nice feature of python is that we can write this function so that\n",
"#it can take an array of values for a, b, or sigma.\n",
"def logP_i(L_i_obs, M_i, a, b, sigma, s):\n",
" L_mu = a*M_i+b\n",
" alpha = (sigma**-2 + s**-2)**0.5\n",
" logP = -0.5*(L_i_obs-L_mu)**2/(sigma**2+s**2) - np.log(alpha) -np.log(s) -np.log(sigma)\n",
" return logP\n",
"\n",
"#The total likelihood for all our data points\n",
"def logP(a, b, sigma, s):\n",
" L = np.array([logP_i(L_i_obs, M_i, a, b, sigma, s) for (L_i_obs, M_i) in zip(L_obs, M_true)])\n",
" return np.sum(L,0)\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Remember: \"real\" answers to inference and stats questions are distributions, not single numbers, or even single numbers with error bars. In fact the really real answer is usually a multi-dimensional distribution! But as it's pretty hard to visualize distributions in more than about 2 dimensions we usually have to condense that information down.\n",
"\n",
"The simplest thing you can do is take slices through a distribution. That's what we'll do here - slices in the $a$ and $\\sigma$ parameters:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"a_trial=arange(0.5, 2.0, 0.001)\n",
"plot(a_trial, exp(logP(a_trial, b, sigma, s)))\n",
"xlabel(\"$a$\")\n",
"ylabel(\"Likelihood\")"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 4,
"text": [
"<matplotlib.text.Text at 0x1095e2c10>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEVCAYAAADkckIIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAH7RJREFUeJzt3Xu0nVV57/HvLxfuQgKBAEmQW4jAEYjUiIiDTdUKqaYd\nShVb5MjxFMqQStVhPeVwmnT0KPWcU+VQK00VKMgZMuqlFG28lx1QSyyHJAQCknA5JAESyI1LEBPy\nnD/ed8HaK3vt/e613tva6/cZY4+sd6255/tk7/2uZ8053zmnIgIzM7OGCVUHYGZm9eLEYGZmQzgx\nmJnZEE4MZmY2hBODmZkN4cRgZmZD9ERikHSDpI2SVuVQ12mSfi7pfkkrJX2g6bVjJC2TtEbSrZIm\nd3s+M7Ne0xOJAbgRODenul4EPhwR/yGt8xpJB6avfR7464iYDWwFPprTOc3MekZPJIaIuIvkjfpV\nko6T9D1J90i6U9KcjHWtiYhH0sdPAZuAQyUJOAf4Zlr0JuB3c/tPmJn1iElVB9CFvwcujYi1kt4C\nfBl4x1gqkDQP2CsiHpE0DdgWEbvTlzcAM3KN2MysB/RkYpB0APBW4BvJB30A9kpfex/wF8N82/qI\nOK+pjiOAm4GLio3WzKy39GRiIOkC2xYRc1tfiIhvA98e6ZvTMYXvAldGxC/SpzcDUyRNSFsNM0la\nDWZmfaXwMQZJEyUtl/SdYV4bkLQ9fX25pKuy1BkRzwGPSTo/rUeSTskYz17APwE3p0mkUWcAdwC/\nlz71H4HbstRpZjaelDH4fAWwGmi3jOvSiJibfv334QpI+jrwc2COpHWSLgb+APiopBXA/cCCjPF8\nAHg78JGmhNRIKp8BPilpDTAVuD5jnWZm44aKXHZb0kzgH4DPAp+MiPe2vD4AfKr1eTMzq07RLYYv\nAp8Gdrd5PYAz04lmSySdVHA8ZmY2isISg6T3AJsiYjmgNsXuBWZFxKnA3+A+fTOzyhXWlSTpc8CH\ngV3APsCBwLciou3toZIeA06PiC0tz3ubOTOzDkREuw/mI35T4V/A2cB3hnl+Oq8lp3nA422+P8qI\nM4f/56KqY3CcjtFxOs6mOKOT7ytzHkMASLo0jXYxcD5wmaRdwA7gghLjMTOzYZSSGCJiKbA0fby4\n6fm/Bf62jBjMzCybnlhEr4cMVh1ARoNVB5DRYNUBZDBYdQAZDVYdQEaDVQeQ0WDVARSp0HkMeZEU\n0ckAiplZH+v0vdMtBjMzG8KJwczMhnBiMDOzIZwYzMxsCCcG6ykScyT+SeLUqmMxG698V5L1FIkf\nAhOB/SM4o+p4zOqs0/fOXt3BzfqQxGEkS6ccAayVOCGChysOy2zccVeS9ZLzgB9H8BKwJD02s5w5\nMVgvOZPXZpz+jKT1YGY5c2KwXnIasDx9/O84MZgVwoPP1hMkJgHbgcMjeF5iIrANmBnB9mqjM6sn\nL4lh491s4KkIngeI4BXgYWBOpVGZjUNODNYr3gjc1/Lcw8AJFcRiNq45MVivOB5Y0/LcL3GLwSx3\nhScGSRMlLZf0nTavXytpjaSVkuYWHY/1rGOBR1uec4vBrABltBiuAFaTbu3ZTNJ84PiImA1cAlxX\nQjzWm44FHml5zi0GswIUmhgkzQTmA18FhhsZXwDcBBARy4ApkqYXGZP1rOPYs8XwCHCcNOzflpl1\nqOgWwxeBTwO727w+A1jXdLwemFlwTNZjJPYCDmfo3wokt68GcFDpQZmNY4WtlSTpPcCmiFguaWCk\noi3Hw06skLSo6XAwIga7CtB6yeuBDRHsbH4ygpBYB8wimdNg1tfS99qBbuspchG9M4EF6TjCPsCB\nkm6OiIuaymwguagbZqbP7SEiFhUVqNXeTPZsLTQ8ARwFrCovHLN6Sj8wDzaOJS3spJ7CupIi4sqI\nmBURxwAXAP/akhQAbgcuApB0BrAtIjYWFZP1rBm0+cAAr7YYzCwnZS67HQCSLgWIiMURsUTSfElr\ngReBi0uMx3rHSImh0WIws5yUkhgiYimwNH28uOW1y8uIwXrakcDjbV5bB7yzvFDMxj/PfLZe4K4k\nsxI5MVgvGCkxPAV47otZjpwYrBccCTzZ5rWnSbb6NLOcODFYrUlMIJnc1i4xbAP2ldi3vKjMxjcn\nBqu7Q4HnInh5uBcjCJJWg7uTzHLixGB1N1I3UsPTJK0KM8uBE4PV3XRgtEmPTgxmOXJisLo7FNg0\nShknBrMcOTFY3R0GPDNKGScGsxw5MVjducVgVjInBqs7txjMSubEYHXnFoNZyZwYrO6ytBieIUkg\nZpYDJwaruywthmeBaSXEYtYXnBis7rK0GLYD+0nsXUI8ZuOeE4PVlsR+JHuGPD9SuXRZjGeBQ8qI\ny2y8KzQxSNpH0jJJKyStlnT1MGUGJG2XtDz9uqrImKynHAo8k77xj+ZZPM5glotCd3CLiF9JOici\ndkiaBPxU0lkR8dOWoksjYkGRsVhPyjK+0PAMHmcwy0XhXUkRsSN9uBcwEdgyTDEVHYf1pCzjCw0e\ngDbLSeGJQdIESStIFkK7IyJWtxQJ4ExJKyUtkXRS0TFZzxhri8FdSWY5KLQrCSAidgOnSToI+IGk\ngYgYbCpyLzAr7W46D7gNOKG1HkmLmg4HW+qw8cktBrMxkDQADHRdT0SWcb18SPpvwEsR8b9GKPMY\ncHpEbGl6LiLC3U19RuJ/AJsj+HyGspcDJ0bwseIjM+sNnb53Fn1X0jRJU9LH+wLvApa3lJkuSenj\neSTJarhxCOs/h+IWg1npiu5KOgK4SdIEkiT0tYj4iaRLASJiMXA+cJmkXcAO4IKCY7LecSjJG34W\nHmMwy0mpXUmdcldSf5L4OfDpCH6WoeypwNciOKX4yMx6Qy27ksy6NJXhb28ejlsMZjlxYrA6O5js\niWEzcIjkOTFm3XJisFpK3+APBrZmKR/By8CvgAOLjMusHzgxWF0dALwcwa/H8D1bSLqfzKwLTgxW\nV2PpRmrYkn6fmXXBicHqyonBrCJODFZXY7kjqWEL3pPBrGtODFZXbjGYVcSJwerKicGsIk4MVleZ\nb1Vt4sRglgMnBqsrtxjMKuLEYHXlxGBWEScGq6tO70pyYjDrkhOD1ZVbDGYVcWKwunJiMKuIE4PV\nVSd3JW0FDvYKq2bdKSwxSNpH0jJJKyStlnR1m3LXSlojaaWkuUXFYz1nzC2GCF4CXgH2KyQisz5R\nWGKIiF8B50TEacApwDmSzmouI2k+cHxEzAYuAa4rKh7rHRJ7A5OBFzv4dncnmXWp0K6kiNiRPtwL\nmMienwAXADelZZcBUyRNLzIm6wlTgS0RdLLvrBODWZcKTQySJkhaAWwE7oiI1S1FZgDrmo7XAzOL\njMl6QicDzw2bcWIw68qkIiuPiN3AaZIOAn4gaSAiBluKtQ4UDvspUdKipsPBYeqx8aObxOAWg/Ut\nSQPAQLf1FJoYGiJiu6R/AX4DGGx6aQMwq+l4ZvrccHUsKio+q51O7khqcGKwvpV+YB5sHEta2Ek9\nRd6VNE3SlPTxvsC7gOUtxW4HLkrLnAFsi4iNRcVkPcMtBrMKFdliOAK4SdIEkgT0tYj4iaRLASJi\ncUQskTRf0lqSO1AuLjAe6x1ODGYVKiwxRMQq4E3DPL+45fjyomKwntXJOkkNW4Djc4zFrO945rPV\nkVsMZhVyYrA6cmIwq5ATg9WR70oyq5ATg9VRty2GQ3KMxazvODFYHXWbGKbmGItZ33FisDrq5q6k\nHcDkdCE+M+uAE4PVisRE4EBgeyffny6851aDWRecGKxuDgKei+CVLurwALRZF5wYrG66uSOpYStO\nDGYdazvzWdKnmg6D11ZBDYCI+EKBcVn/6mbgucEtBrMujLQkxutIksAc4M0kC94JeA/wi+JDsz7V\nzcBzg8cYzLrQNjE0lrmWdBfwpoh4Pj1eCCwpJTrrR24xmFUsyxjDYcDOpuOd6XNmRcgjMXiMwawL\nWVZXvRn4haRvk3Ql/S7pPs1mBcirxXBiDrGY9aVRE0NEfFbS94Gz0qc+EhGtG+6Y5eVghu4D3gmP\nMZh1Ievtqq8Au5u+MpE0S9Idkh6QdL+kjw9TZkDSdknL06+rstZv45LHGMwqNmqLQdIVwB8Cja6k\nWyR9JSKuzVD/TuATEbFC0gHA/5X0o4h4sKXc0ohYMNbgbVzK464kjzGYdSHLGMN/Bt4SES8CSPor\n4G5g1MQQEU8DT6ePX5D0IHAk0JoY1Pq91rfcYjCrWNaupN1tHmcm6WhgLrCs5aUAzpS0UtISSSd1\nUr+NG3klBo8xmHUoS4vhRmBZy11JN4zlJGk30jeBKyLihZaX7wVmRcQOSecBtwEnjKV+G1fySAzb\ngIMkJkR09kHGrJ8pIkYvJJ0OvC09vGssdyVJmgx8F/heRFyTofxjwOkRsaXpuQD+oqnYYEQMZo3B\neoOEgJeB10Xwcpd1bQOOieh63SWzniFpABhoemphRIy5qz5LiwGSu5IaGWQsdyUJuB5Y3S4pSJoO\nbIqIkDSPJFnt8YmxMRPbxrX9gZ3dJoVUozvJicH6RvqBebBxnK5UMWajjjGkdyXdAhxKMuP5luFu\nO23jbcCFwDlNt6OeJ+lSSZemZc4HVklaAVwDXDDm/4WNF3nckdTgAWizDo3alSRpFXBG011J+wN3\nR8QbS4ivEUN00hyy3iJxKnBzBKfmUNePgP8ZwQ+7j8ysN3X63lnaXUlmGeSxF0ODWwxmHSrlriSz\njA4GNudUl29ZNetQlrWSviBpKclaSYHXSrLiuMVgVgNZ70paQTKDeRIQko6KiCeKC8v6VB5zGBq2\nAkfkVJdZX8myVtIfAwuBTTBkg/bSBp+tb+SZGLYAJ+dUl1lfydJi+BNgTkTk1fdr1s7BwCM51eUx\nBrMOZbkr6QnguaIDMSP/FoPHGMw60LbFIOlT6cNHgUFJ3wV+nT4XEfGFooOzvpPnBDcvvW3WoZG6\nkl5HchfSEyQ7au2VfonXlscwy5NbDGY1kGkRvap55nN/kHgCOCuCru94k9gX2BLBvt1HZtabOn3v\nHKkr6X9HxBWSvjPMy+Ed16wAuc1jiOAlCUnsG8FLedRp1i9G6kr6WvrvX5cRiPU3ib2BvYHW/Tq6\n0ehO2pBjnWbjXtvEEBH3pP8OlhaN9bOpJF0/efZtOjGYdWCkrqRVI3xfRMQpBcRj/SvPgecGz2Uw\n68BIXUnvLS0Ks2ISg29ZNevASF1JjzceSzoaOD4ifixpP2Bi4ZFZvymqxeDEYDZGWXZwuwT4BrA4\nfWomcFuWyiXNknSHpAck3d9u5zdJ10paI2mlpLlZg7dxxV1JZjWRZUmMj5Esuf0cQEQ8TLLFZxY7\ngU9ExMnAGcDHJJ3YXEDSfJLWyGzgEuC6jHXb+OIWg1lNZEkML0fEq5uzS5pExpnPEfF0RKxIH78A\nPAgc2VJsAXBTWmYZMEXS9Cz127gylfz2YmjwGINZB7IkhqWS/iuwn6R3kXQrDTfpbUTpOMVcYFnL\nSzNIltxoWE/SXWX9xS0Gs5rIkhg+AzwDrAIuBZYAV43lJJIOAL4JXJG2HPYo0nJc/3U6LG8eYzCr\niSz7MSyKiD8H/h5A0kTg/wC/n+UEkiYD3wJuiYjhBq03ALOajmcyzIQkSYuaDgc98W7ccYvBrEuS\nBoCBbuvJkhiOkvRnEXG1pL2BfwQy7fksScD1wOqIuKZNsduBy4FbJZ0BbIuIja2FImJRlnNaz/I8\nBrMupR+YBxvHkhZ2Us+oq6tKmkDSQlgFnAMsiYgvZqpcOgu4E7iP17qHrgSOAoiIxWm5LwHnAi8C\nF0fEvS31eHXVcU7iEeDdEazNsc6pwGMRTMmrTrNe0ul7Z9vEIOl0Xnszn0wyj+HnwFcBWt+8i+TE\nMP5JbAWOi8iv1SAxgWRzqb0jhuxXbtYXikgMgwwdBB6yQU9EnDPWk3XKiWF8k5gIvEwBb+ASW4DZ\nEXjPcus7ue/HEBEDXUVklt0U4PmCPtU3BqCdGMwyGml11Qsj4pZ07+c9Wg7e89lyVMTAc4PvTDIb\no5HuSto//bex97NZUaZSbGLwXAazMRipK2lx+u+i1tckfaLAmKz/FNli8C2rZmOUZebzcD6ZaxTW\n7w7BXUlmtdFpYjDL0zTg2YLqdleS2Rg5MVgdFJ0Y3GIwG4OR7kp6gfaDzvsVE471qWnA/QXVvRU4\nraC6zcalkQafDygzEOtr00hW8C3CZpIxDDPLyF1JVgdFdiU9m9ZvZhk5MVgdFJ0YDi2obrNxyYnB\n6qDIxPAMbjGYjYkTg1VKQiRjAEWtZbQd2Fdi74LqNxt3nBisagcCv4rg5SIqjyDwALTZmDgxWNWK\n7EZqeAaPM5hlVmhikHSDpI2SVrV5fUDSdknL06+riozHaqmMxOABaLMxyLLnczduBP4GuHmEMksj\nYkHBcVh9ldVi8AC0WUaFthgi4i6Smacj8c5s/c0tBrOaqXqMIYAzJa2UtETSSRXHY+Vzi8GsZoru\nShrNvcCsiNgh6TzgNuCE4QpKWtR0OBgRg8WHZyUoKzGcXPA5zConaQAY6LaeShNDRDzf9Ph7kr4s\n6eCI2GNt/uE2DLJxYRrwWMHn8LIY1hfSD8yDjWNJCzupp9KuJEnTJSl9PA/QcEnBxjXfrmpWM4W2\nGCR9HTgbmCZpHbAQmAyvbh16PnCZpF3ADuCCIuOxWvLgs1nNKKLdlgv1ISkiwncvjUMSDwLvj2B1\ngec4AlgeweFFncOsjjp976z6riSzMloMm4FDJP+9m2XhC8UqIzERmEKy/WZhIvg18CJwUJHnMRsv\nnBisStOAbRHsKuFcHoA2y8iJwao0HdhY0rk8AG2WkRODVanMxODZz2YZOTFYlQ4Hni7pXJtIEpGZ\njcKJwapUZovhaZwYzDJxYrAqlZ0YPI/BLAMnBqvSdMrrSnJiMMvIicGqdDhuMZjVjhODVcldSWY1\n5MRgVSo9MUjeMdBsNE4MVol0OYxDSOYXFC6CF0h2DDygjPOZ9TInBqtKYzmMnSWe091JZhk4MVhV\nyuxGanBiMMvAicGqUuas5wYnBrMMCk0Mkm6QtFHSqhHKXCtpjaSVkuYWGY/VShUthqdwYjAbVdEt\nhhuBc9u9KGk+cHxEzAYuAa4rOB6rD3clmdVUoYkhIu4Cto5QZAFwU1p2GTBFktez6Q8zgA0ln9OJ\nwSyDqscYZgDrmo7XAzMrisXK5cRgVlOTqg4A9phwFMMWkhY1HQ5GxGBRAVkpZpJ8ECjT08ARJZ/T\nrDSSBoCBbuupOjFsAGY1Hc+kzafIiFhURkBWmipaDBtwi9TGsfQD82DjWNLCTuqpuivpduAiAEln\nANsiouwBSSuZxASST+5PlnzqTcBBEvuUfF6znlJoi0HS14GzgWmS1gELgckAEbE4IpZImi9pLfAi\ncHGR8VhtHEYy6/nlMk8awW6JJ0laK4+UeW6zXlJoYoiID2Uoc3mRMVgtVTG+0LCOpPvSicGsjaq7\nkqw/zaD6xGBmbTgxWBXa3mRQgvU4MZiNyInBqlCHriQza8OJwapQxa2qDevwLatmI3JisCrMoroW\ng7uSzEbhxGBVOAZ4rKJzuyvJbBSKGHYFilqRFBHhvXrHAYnJwAvAASXv3tY4v4CXgEMieLHs85uV\nqdP3TrcYrGyzgKerSAoAEQTwBPD6Ks5v1gucGKxsVXYjNTwCHFdxDGa15cRgZTsGeLTiGJwYzEbg\nxGBlq0uL4diKYzCrLScGK1sdEsOjuMVg1pYTg5WtDonBLQazETgxWNmOpfrE8ChwdLovhJm18IVh\npZF4HXAAyRablYlgB7CVZGkOM2vhxGBlegPwywh2Vx0ISavB3Ulmwyg8MUg6V9JDktZI+swwrw9I\n2i5pefp1VdExWWXeADxUdRCph4E5VQdhVkdFb+05EfgS8E6S1TT/XdLtEfFgS9GlEbGgyFisFuqU\nGB4ATq46CLM6KrrFMA9YGxGPR8RO4Fbgd4Yp53WQ+oMTg1kPKDoxzCBZzbJhPXsO+AVwpqSVkpZI\nOqngmKw6TgxmPaDQriSSN/3R3AvMiogdks4DbgNOaC0kaVHT4WBEDOYSoZVCYhLJYO+aqmNJrQf2\nkzgkgs1VB2OWB0kDwEC39RSdGDYwdO37PTZoiYjnmx5/T9KXJR0cEVtayi0qMlAr3BuAJyJ4qepA\nIFllVXq11XBn1fGY5SH9wDzYOJa0sJN6iu5KugeYLeloSXsBHwRuby4gabokpY/nkewRsWXPqqzH\nzQWWVx1EC3cnmQ2j0BZDROySdDnwA2AicH1EPCjp0vT1xcD5wGWSdgE7gAuKjMkqU8fEcD/wxqqD\nMKsb7+BmpZC4A7g6gh9WHUuDxFnAFyKYV3UsZkXo9L3TicEKl26nuQWYE8GmquNpkNgfeAaYGsHL\nVcdjljdv7Wl1dhzwQp2SAkC65/PDwKlVx2JWJ04MVoazgJ9WHUQbvwB3JZk1c2KwMryd+iaGZcBb\nqw7CrE6cGKwMZwF3VR1EG4PAOek4iJnhxGAFkzgcmE4yZ6COHgV+TTIBz8xwYrDinQv8OIJXqg5k\nOBEE8BPgHVXHYlYXTgxWtN8Gvlt1EKP4MfDuqoMwqwvPY7DCSOwFbCKZv7Cx6njakZgCPAHMimB7\n1fGY5cXzGKyOfgtYXeekABDBNmAp8N6qYzGrAycGK9LFwD9UHURG/wh8qOogzOrAXUlWCIlpwFrg\n9b3QPZMuj/H/gHkRPFp1PGZ5cFeS1c3HgW/0QlKAV5fHuAG4vOpYzKrmFoPlTmIqyU5tPfXpW2IW\nsAI4JYINVcdj1i2vrmq1IfF3ABH8UdWxjJXE1cCMCC6qOhazbtWyK0nSuZIekrRG0mfalLk2fX2l\npLlFxmPFkziPZO7CsL/vHvA54EyJD1YdiFlVCksMkiYCXyKZ+XoS8CFJJ7aUmQ8cHxGzgUuA64qK\npwzpRty1V1ScEmcANwEX5DG2UMXPM4LngQ8AX5JGnw3d77/zvDnOeiiyxTAPWBsRj0fETuBW4Hda\nyiwgeSMhIpYBUyRNLzCmog1UHUBGA3lWJjFZ4uMk+3l/JIKf5VT1QE71jEkE9wLvB26V+FOJvUco\nPlBOVF0bqDqAjAaqDiCjgaoDKFKRiWEGsK7peH363GhlZhYYk+VEYorEOyT+kmSgeQFwdgRLKg4t\nFxHcCbwZOBtYK/FXEr8pcWjFoZkVblKBdWcd1W4dGCl8NFziD3mt9dJ8/tEej/L6x45N9xHOWl8O\n5+zk8Z/MlJjf4TmnAIeR/J6WA3cD7wOWpwvSjRsRPA78tsRpwO8BfwmcLDER2AxshSumSfwmsAt4\nJf13dx6nz7eOPzpB4vQc6iyY46yDwu5KknQGsCgizk2P/wzYHRGfbyrzd8BgRNyaHj8EnB0RG1vq\nGldvOGZmZenkrqQiWwz3ALMlHQ08CXyQPZccuJ1kQtGtaSLZ1poUoLP/mJmZdaawxBARuyRdDvwA\nmAhcHxEPSro0fX1xRCyRNF/SWuBFkrV1zMysQj0xwc3MzMpTq7WSskyIS8u9WdIuSe8rM7703Fkm\n7Q1IWi7pfkmDJYfYiGHEOCVNk/R9SSvSOD9SQZhIukHSRkmrRihT6STI0WKU9AdpbPdJ+pmkU8qO\nMY1j1J9lWq6y6yc9f5bfeR2uodF+73W5hmZJukPSA2kcH29TLvt1FBG1+CLpbloLHA1MJlmz5sQ2\n5f6VZFew99ctRpK7dh4AZqbH0+r4swQWAVc3YiS5y2ZSBbG+HZgLrGrz+nxgSfr4LcDdNYzxrcBB\n6eNzq4gxS5xNfxuVXD9j+HlWfg1ljLMu19DhwGnp4wOAXw5zvY/pOqpTiyHLhDiAPwa+CTxTZnCp\nLDH+PvCtiFgPEBHPlhwjZIvzKeDA9PGBwOaI2FVijABExF3A1hGKVD4JcrQYI+LfIqIx03sZFc3F\nyfCzhGqvHyBTnHW4hrLEWZdr6OmIWJE+fgF4EDiypdiYrqM6JYZRJ8RJmkHyBtdYOqPsAZIsk/Zm\nAwenTbt7JH24tOhekyXOrwAnS3oSWAlcUVJsY9VrkyA/CvWc5FeD6yerOlxDWdTuGkrvAp1L8gGl\n2ZiuoyJvVx2rLH+k1wD/JSJCkthzclzRssQ4GXgT8A5gP+DfJN0dEWsKjWyoLHFeCayIiAFJxwE/\nknRqRDxfcGydKH0SZCcknQP8J+BtVcfSRtXXT1Z1uIayqNU1JOkAktbgFWnLYY8iLcdtr6M6JYYN\nwKym41kkWa3Z6SRzHiDp0ztP0s6IuL2cEDPFuA54NiJeAl6SdCdwKsmyEWXJEueZwGcBIuIRSY8B\nc0jmn9RJ6/9lZvpcraQDzl8Bzo2I0bpzqlL19ZNVHa6hLGpzDUmaDHwLuCUibhumyJiuozp1Jb06\nIU7SXiQT4ob8wUbEsRFxTEQcQ5IZLyv5j3rUGIF/Bs6SNFHSfiQDPatLjDFrnA8B7wRI+xrnQC03\n1bkdkr0RRpoEWSVJRwHfBi6MiLVVx9NODa6frOpwDWVRi2sobf1dD6yOiGvaFBvTdVSbFkNkmBBX\naYBknrT3kKTvA/eRrJnzlYgo9Y8648/yc8CNklaSfED404jYUmacAJK+TrJQ3TRJ64CFJF0JjZ9n\n5ZMgR4sR+HNgKnBd+ml8Z0TMq2GctZDhd175NZQlTmpyDZF0XV4I3CdpefrclcBRjVjHeh15gpuZ\nmQ1Rp64kMzOrAScGMzMbwonBzMyGcGIwM7MhnBjMzGwIJwYzMxvCicHMzIZwYjAzsyGcGMzMbIja\nLIlh1mskfYhkiYSZwKaI+GrFIZnlwi0Gsw5ImgO8OyJuBl4B7q84JLPcODGYdeZCXlux9lRg+Qhl\nzXqKE4NZZ6YAv0yXNX8d8BsVx2OWG48xmHXmZuC3gJOAR4Ajqg3HLD9edtvMzIZwV5KZmQ3hxGBm\nZkM4MZiZ2RBODGZmNoQTg5mZDeHEYGZmQzgxmJnZEE4MZmY2xP8HbEd7DQFKGeAAAAAASUVORK5C\nYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x1095a7d90>"
]
}
],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"sigma_trial=arange(0.5, 2.0, 0.01)\n",
"plot(sigma_trial, exp(logP(a, b, sigma_trial, s)))\n",
"xlabel(\"$\\sigma$\")\n",
"ylabel(\"Likelihood\")"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 5,
"text": [
"<matplotlib.text.Text at 0x109606850>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEVCAYAAADkckIIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYXGWZ/vHvDQSRRRBRkCRMlEVFiCCLGMA0IGMSBVFw\nYRFFHBAHZVwZGZeoqD91HFFHIcoiogOKsgQMIiINKBJAkrBmSECEEBNkNRBgCHl+f5y3Y3enO13d\nXafec6ruz3XVlaqu0+fc6eTU0+c976KIwMzMrMdauQOYmVm1uDCYmVkfLgxmZtaHC4OZmfXhwmBm\nZn24MJiZWR+1KAySzpS0VNKtTdjXTpKuk3SbpHmS3tnrvZdJmi1pgaTzJI0Z7fHMzOqmFoUBOAuY\n0qR9PQm8JyJ2SPs8RdIL0ntfA74ZEdsCjwJHN+mYZma1UYvCEBHXUnxQryJpa0mXSbpJ0jWSXtHg\nvhZExN3p+V+BB4EXSxKwD/CLtOnZwEFN+0uYmdXEOrkDjMIPgGMjYqGk1wHfB/Ybzg4k7Q6sGxF3\nS9oMeCwiVqa3HwDGNjWxmVkN1LIwSNoQeD1wfvGLPgDrpvfeDnxhgG9bFBFTe+3jpcCPgSPLTWtm\nVi+1LAwUTWCPRcTO/d+IiAuAC9b0zemewqXASRFxQ/ryw8AmktZKVw3jKK4azMw6Sun3GCStLWmO\npEsGef87qRfQPEmrfdAPJCL+DvxZ0iFpH5I0scE86wIXAj9ORaRnnwFcBbwjfem9wEWN7NPMrJ20\n4ubzCcAdwGrTuEqaBmyTegEdA5w60A4knQtcB7xC0v2SjgIOB46WNBe4DTiwwTzvBPYG3pcK1pxe\nReVE4GOSFgAvBM5o9C9pZtYuVOa025LGAT8Cvgx8LCIO6Pf+acBVEfGz9Ho+MDkilpYWyszM1qjs\nK4ZvAZ8EVg7y/ljg/l6vF1G07ZuZWSalFQZJbwEejIg5gNa0ab/XXjnIzCyjMnslTQIOTPcR1gNe\nIOnHEdG7e+gDwPherwfsCSTJxcLMbAQiYk2/mA/6TaU/gMnAJQN8fRowKz3fA7h+kO+PVuRswt9z\neu4MzumMzumcvXLGSL6vleMYAkDSsSntjIiYJWmapIUUcxgd1cI8ZmY2gJYUhoi4Grg6PZ/R773j\nW5HBzMwaU4tJ9GqkO3eABnXnDtCg7twBGtCdO0CDunMHaFB37gAN6s4doEyljmNoFkkRI7mBYmbW\nwUb62ekrBjMz68OFwczM+nBhMDOzPlwYzMysDxcGMzPrw4XBzMz6cGEwa3MSa0m1Xa3RMnBhMGtD\nEltLnC6xBHgGeFTiAomDpTXOdmzmwmDWTiQk8XlgNrAY2B3YAJgAzAQ+B1wo8eJsIa3yPPLZrE2k\n5qLvA68F3hzBaishSjwP+CJwGLBvBAtam9JaaaSfnW53NGsfp1JcGewTwbKBNojgGeBEiQXAFRJ7\nR/RZRdHMhcGsHUi8F9gL2C2CJ4baPoLTJTamKA6vi+Dx0kNabbgpyazmJF5NMdvnPhHcNszvPQ3Y\nIIL3lJHN8vIkemYdSGIt4CzgpOEWheTjwG4ShzY3mdWZC4NZvb0HeA44YyTfHMGTFDeivy0xtpnB\nrL7clGRWUxIbAfOBt0Vwwyj39RVgXARHNiWcVcJIPztdGMxqSuLLFB/m723CvjakKDKHRHD9qMNZ\nJVSyMEhaj2Kt5+cB6wIXR8Sn+23TBVwM3JO+9MuIOLnfNi4MZr1IbAosAHZqVndTifcAHwb2iGBl\nM/ZpeVXy5nNEPA3sExE7AROBfSTtNcCmV0fEzulx8gDvm1lfHwJmNnkMwk+BAN7ZxH1aDZV+8zki\nlqen6wJrA48MsJmvBswaJLE+xW/2X2/mftNVwueAz6beTtahSv/Hl7SWpLnAUuCqiLij3yYBTJI0\nT9IsSduXncms5o4CrovgzhL2/RtgGXBICfu2mmjZzWdJGwOXA/8eEd29vr4R8FxELJc0Ffh2RGzX\n73t9j8GMVeMW7gKOjOC6ko4xjeJqZKLvNdRb5edKiojHJf0K2JVilGbP15f1en6ZpO9L2jQi+jQ5\nSZre62V37+Ji1kH2BZ4A/ljiMS4DpgMHAReUeBxrstSZp2vU+ym5V9JmwIqIeEzS8ymuGL4QEVf2\n2mZz4MGICEm7Az+PiAn99uMrBjNA4mfANRF8r+TjvBP4cAR7l3kcK1dVu6vuCJxNcS9jLeCciPiG\npGMBImKGpH8FjgNWAMuBj0XE9f3248JgHU/iJRTNSBMieKzkY60D3A0cHMFNZR7LylPJwtAsLgxm\nIPFJYPsIjmrR8T5BMU7iiFYcz5rPhcGsjaXlOOcDR5V103mAY25CMfB0hwgWt+KY1lyVHOBmZk3z\nWopxQGXedO4jNVf9D/DBVh3TqsGFwaweDgXOi6DVl/gzgPenew7WIVwYzCoujV14F3Buq48dwa3A\nfcDUVh/b8nFhMKu+PYFHI7g90/F/AByT6diWgQuDWfUdCpyX8fg/ByZJjM+YwVrIhcGswlLb/iFk\nLAwRLE/Hf3+uDNZaLgxm1bYnsChi1XoluZwFHJm6zVqbc2Ewq7a3UixkldufgKcpCpW1ORcGs4pK\nv50fBFyUO0vqJns2jH4ZUas+Fwaz6tqB4hy9JXeQ5KfAwRLPzx3EyuXCYFZdbwUuyjCobUARPADc\nSJHL2pgLg1l1HUQ17i/0dg54Ur1250n0zCpIYixwK/CSCFbkztNDYiNgEcXU34/mzmNr5kn0zNrL\nFOA3VSoKABEsA66kuJqxNuXCYFZNU4Bf5w4xiJ9RzN1kbcpNSWYVk0Y7/41iUZ6/5s7Tn8QGwGJg\n6wgeyp3HBuemJLP2sQdwbxWLAkAET1Jczbw9dxYrhwuDWfVMBS7LHWIIbk5qY6UVBknrSZotaa6k\nOyR9dZDtviNpgaR5knYuK49ZjVT5/kKPy4BdJDbPHcSar7TCEBFPA/tExE7ARGAfSXv13kbSNGCb\niNiWYr73U8vKY1YH6YP25bRwCc+RiOAp4FfAwbmzWPOV2pQUEcvT03Up1qt9pN8mB1LMv0JEzAY2\nkeTfQKyT7Qd0R/Bs7iANcHNSmyq1MEhaS9JcYClwVUTc0W+TscD9vV4vAsaVmcms4t5IMU6gDi4H\ndpTYMncQa65SF/iOiJXATpI2Bi6X1BUR3f0269+VasD+s5Km93rZPcB+zGotzaa6H/D13FkaEcEz\nEjOBdwDfzp3HQFIX0DXq/bRqHIOkzwJPRcR/9vraaRQf8uel1/OByRGxtN/3ehyDtT2JbYCrgXFV\nmThvKBJTgc9EeJ2GKqrcOAZJm0naJD1/PrA/MKffZjOBI9M2ewCP9S8KZh3kjcBv61IUkiuB7SW2\nyB3EmqfMewwvBX6X7jHMBi6JiCslHSvpWICImAXcI2khMAP4UIl5zKpuP+pzfwGACP6Pouuqp+Ju\nI54Sw6wCJNaimAZjYlr3oDYk3gEcHcGU3Fmsr8o1JZnZsOwE/K1uRSH5NTBJYuPcQaw5XBjMqqEL\n+F3uECORpuK+BpiWO4s1hwuDWTV0Ad2ZM4zGhcDbcoew5vA9BrPMJNYGHgJeFcGS3HlGQuLFwAJg\niwiezp3HCr7HYFZfE4EldS0KABH8DZhH0eXWas6FwSy/yRQD2+ruQrzkZ1twYTDLr4t631/ocRFw\nYGoasxpzYTDLKI1feANtcMUQwb3AA+DpMerOhcEsrx0pxi9UchnPEXDvpDbgwmCW1xsoxgC0iwuB\nt6WZYq2mXBjM8poE/CF3iCa6DXiOoqeV1ZQLg1lek4DrcodoljQz7EyK1RmtplwYzDKRGAesTzEw\nrJ1cjGdbrTUXBrN8Xg9cV7P1Fxrxe2BCKnxWQy4MZvnsSRs1I/WIYAUwCzcn1ZYLg1k+bXV/oR/f\nZ6gxT6JnloHE+hQL82wWwVO58zSbxEYUg93GRfD33Hk6lSfRM6uXXYHb2rEowKo1Gn4PXtWtjlwY\nzPJot/ELA5mJeyfVUqmFQdJ4SVdJul3SbZI+MsA2XZIelzQnPT5TZiazimjLG8/9zASmSozJHcSG\nZ52S9/8s8NGImCtpQ+BPkq6IiDv7bXd1RPhGlXWENF3EJODY3FnKFMFiiYXA3tR02dJOVeoVQ0Qs\niYi56fkTwJ3AlgNs6hvL1km2A5ZFsDh3kBbwYLcaatk9BkkTgJ2B2f3eCmCSpHmSZknavlWZzDJp\n526q/c0E3upJ9eql7KYkAFIz0i+AE9KVQ283A+MjYrmkqRSLfWw3wD6m93rZHRHdJcU1K1sn3Hju\ncRuwkmJ68VsyZ2l7krooFn4a3X7KHscgaQxwKXBZRJzSwPZ/BnaJiEd6fc3jGKxtSNwBHB7BnNxZ\nWkHiW8AjEXwpd5ZOU8lxDJIEnAHcMVhRkLR52g5Ju1MUq0cG2tas7iQ2BcYDt+bO0kK+z1AzZTcl\n7QkcAdwiqee3o5OArQAiYgZwCHCcpBXAcuDdJWcyy2kP4IY0n1Cn+D3wMolxESzKHcaG5ikxzFpI\n4mQgIvhs7iytJHEOxUyyp+bO0klG+tk56BWDpI/3ehn8o0tpAETEfw33YGbGJODruUNkcDHwAXBh\nqIM1NSVtRFEEXgHsRtHtTMBbgBvKj2bWXtII4N2A63NnyeBy4EyJF3hSveobsilJ0rXAtIhYll5v\nBMyKiL1bkK8ng5uSrPYkdgHOjmCH3FlykLgMOCuCn+fO0inK7JX0EoqpLXo8m75mZsPTSQPbBnIx\nXqOhFhrplfRj4AZJF1A0JR0EnF1qKrP2NImiSaVTXQJ8RWJMRJ9fNq1iGuqVJGkXYK/08pqIaOnA\nHDclWTuQ+AuwfwR35c6Si8SNwIkRnlSvFcoe4PYcxbD2noeZDYPEOGB9YEHuLJl5sFsNDFkYJJ0A\n/AR4McW9hZ8MtK6Cma3RJIp+/NUfOFSui4EDPaletTXSK+lWYI+IeDK93gC4PiJ2bEG+ngxuSrJa\nkzgF+GsEX8udJadUEO4GDorwpHplK7spaeUgz82sMZ3eIwmAdMXkJT8rrpHCcBYwW9J0SV+gGJxz\nZrmxzNqHxPrAq4GbcmepCN9nqLjh9EraM7281r2SzBonMRn4WgR75M5SBRLrAEuB13hSvXK1oldS\npIebksyGx81IvaSZZWcBB+TOYgNzrySz8rkwrM73GSrMvZLMSpR64TwE7BjB4tx5qkJiI+ABYJwn\n1SuPeyWZVdN2wDIXhb4iWEax7vWbcmex1TUyV1JPr6TecyW5V5JZYyZRfADa6np6J52fO4j1Ndy5\nkgL3SjJrmMTpwJwIvpc7S9VIjKVY+3pzT6pXjqav4NbPXGBJ2j4kbRUR9w33YGYdaBK4KAwkggck\n7gb2Bk+qVyWN9Er6MEWf4yuAS4FfpceQJI2XdJWk2yXdNlhvJknfkbRA0jxJOw8jv1llSWwKjKf4\nrdgG5jUaKqiRK4Z/A14REQ+PYP/PAh+NiLmSNgT+JOmKiLizZwNJ04BtImJbSa+jWBPWA4GsHewB\n3JD67dvAZgIXS3zUEwxWRyO9ku6DkXUni4glETE3PX8CuBPYst9mB5IW/omI2cAmkjYfyfHMKsbj\nF4bWczXVsu7vNrRBrxgkfTw9vQfolnQp8H/paxER/zWcA0maAOwMzO731ljg/l6vFwHjKJqvzOps\nT+js2VSHEkFIq3onebbVilhTU9JGFL2Q7qP44F43PZS+3rDUjPQL4IR05bDaJv1er7Z/SdN7veyO\niO7hZDBrJYkxwK4Uk07aml0MfB34Uu4gdSepC+ga9X4a6a46qgNIYyhuWl8WEacM8P5pFB/056XX\n84HJEbG01zburmq1IrELcHYEO+TOUnWpiC4BJkbwQO487aTpI58lfTv9eckAj5kNhhJwBnDHQEUh\nmQkcmbbfA3isd1EwqynfX2hQGsNwGe6dVBlrako6J/35zVHsf0/gCOAWST2D4k4CtgKIiBkRMUvS\nNEkLgSeBo0ZxPLOqmARcnjtEjVwMHE3RK9EyK70pqRnclGR1I3EfsF8EC3JnqYM0qd5iYKwn1Wue\npo98TrOqDiYiYuJwD2bWCSTGA88HFubOUhcRLJNWTarnuZMyW1NTkhfRMBuZ1wPXecDWsF0EvA0X\nhuwGvfkcEff2PNKXtknPHwRGMgrarFP4xvPIXARMk1gvd5BO18hcScdQVPAZ6UvjKP4BzWxgLgwj\nEMESigk7vUZDZo1MifGvFFNu/x0gIu6iWOLTzPqRWB94NXBT7iw1dT5wSO4Qna6RwvBMRDzT80LS\nOgxz5LNZB9kNuDWCp3IHqakLgLdIPC93kE7WSGG4WtJ/AOtL2p+iol9Sbiyz2nIz0ihE8FeKifX+\nOXeWTtZIYTgR+BvFP9axwCzgM2WGMqsxF4bROx94R+4QnWzIAW6SvhgRn+v1em3gnIg4rOxwvY7p\nAW5WeRJrUfwStWMEi3PnqSuJLYHbgS0ieGao7W1wTZ8rqZetJH06HeR5FG2Adw33QGYd4JXA4y4K\no5N+frcB++fO0qkaKQzvByZKOoliltTuiJheaiqzetoT+H3uEG3CzUkZDdqUJGkX/tH7aAzFOIbr\ngNMBIuLmVgRMWdyUZJUncTbwhwh+kDtL3UmMpbivuUXEqgXCbJhG+tm5psLQTd9uqX0W6ImIfYZ7\nsJFyYbA6kLgbOCCCO3JnaQcSvwe+EsGs3FnqqumFoUpcGKzq0g3TW4EXR7Ayd552IPFvwGsiPBX/\nSJUxu+oREfGTtPbzalcOw13z2azN7UnRjOSi0DznA5+VWC+Cp3OH6SRruvm8Qfpzo36PDdOfZvYP\newF/yB2inaRlPucB03Jn6TQjakqS9NGI+FYJeQY7npuSrNIk/gR8JMLFoZkkPgBMjeDg3FnqqKX3\nGCTdHxHjh/2NI+TCYFWWVh9bAmzqAVnNJbEJ8BfgnyJ4LHeeuilzgJuZrdnrgJtdFJovFYMrwVcM\nrVRqYZB0pqSlgy0TKqlL0uOS5qSH52CyOtoLD2wr00+Bw3OH6CSDFgZJT0haNtAD2LLB/Z8FTBli\nm6sjYuf0OLnR4GYV4sJQrl8BO6VBb9YCa1rac8OI2GiQx9qN7DwirgUeHWIz3zuw2pIYQ9GU5BlV\nS5K6ql4IHJo7S6fIfY8hgEmS5kmaJWn7zHnMhus1wL0RQ/4CZKPzP7g5qWUGHeDWIjcD4yNiuaSp\nFGtJbzfQhpKm93rZHRHd5cczG5KbkVqjG3iJxPaecmRwkrqArlHvp+wpMSRNAC6JiB0b2PbPwC4R\n8Ui/r7u7qlWSxC+ACyP4ae4s7U7im8BTEV4orFG17K4qaXNJSs93pyhUjwzxbWaVICE84rmVfgoc\nnhZEshKV2pQk6VxgMrCZpPuBz1NM4U1EzAAOAY6TtAJYDry7zDxmTbYNsIJiAJaVbw7wJPAGiqYl\nK4lnVzUbIYljgL0jeE/uLJ1C4qPAa/0zb0wtm5LMam4f4KrcITrMOcABaaoMK4kLg9kIpPsLLgwt\nFsFDwBV4TEOpXBjMRuaVwNMR/Dl3kA50JsVa9FYSFwazkfHVQj6/AbaQmJg7SLtyYTAbGReGTCJ4\nDvgRcHTmKG3LvZLMhin1o18K7BLBfbnzdCKJlwOzgXGe7nxw7pVk1jqvBv7uopBPBPcAtwAH5c7S\njlwYzIbPzUjV4JvQJXFhMBu+LlwYquACYFeJCbmDtBsXBrNhSPcXJuPCkF0ETwE/Bj6YO0u7cWEw\nG56JwEMRLM4dxAD4HnC0xPNzB2knLgxmw+P7CxUSwULgRjwSuqlcGMyGx4Wher4LfDhNU2JN4MJg\n1iCJtYG98ZTPVXM5sAGwZ+4g7cKFwaxxOwOLI1iaO4j9QwQrKe41HJ87S7twYTBr3P7Ab3OHsAH9\nCPhniS1zB2kHLgxmjXsTxQRuVjERPA6ci7uuNoXnSjJrgMRGwF+BzSN4MnceW53E9sDvgH/y/EkF\nz5VkVq4u4AYXheqK4A7gVuBdubPUXamFQdKZkpZKunUN23xH0gJJ8yTtXGYes1F4E0XvF6u2bwAn\nphHqNkJl//DOAqYM9qakacA2EbEtcAxwasl5zEbKhaEergCeAg7MHaTOSi0MEXEt8OgaNjkQODtt\nOxvYRNLmZWYyG6409/8LKKZ5tgqLIICvAp/2gLeRy325NRa4v9frRcC4TFnMBvMm4Depv7xV34UU\nhXzf3EHqap3cAWC1qj5gNylJ03u97I6I7rICmfXzZuCc3CGsMRGslPga8Gngytx5WklSF0VHidHt\np+zuqpImAJdExI4DvHcaxYf8een1fGByRCztt527q1oWEusDS4CtIngsdx5rjMQYYCHwjghuyJ0n\nl7p2V50JHAkgaQ/gsf5FwSyzfYGbXRTqJYJngf+kuGqwYSq1KUnSuRSLmmwm6X7g88AYgIiYERGz\nJE2TtBB4EjiqzDxmI/Bm4NLcIWxEzgD+Q2L7NMbBGuSRz2aDSL1a7gP2j2B+7jw2fBKfBnaM4LDc\nWXKoa1OSWZVNBJ4B/jd3EBux/wb2ldghd5A6cWEwG9wBwKWpb7zVUATLgK8DX8qdpU5cGMwG9zaK\nPvFWb6cCu0nsljtIXbgwmA1AYgIwHvh95ig2ShE8BZwMfNWjoRvjwmA2sLcBF0fwXO4g1hRnUMyq\nMOjcbfYPLgxmA3s7cEHuENYcaVzDicA3pErM+FBpLgxm/UhsAexAseiLtY+ZwMPA+zLnqDwXBrPV\nHQRc5lXA2kvqXfZx4IsSG+fOU2UuDGareyfwi9whrPkiuAm4jGIWBhuERz6b9SIxlmJ5yC0jeDp3\nHms+iZcAtwOT232qDI98NmuOdwEXuSi0rwgepOi++l13Xx2YC4NZX4cB5+YOYaX7HvBC4IjcQarI\nTUlmicR2wDXAuAhW5M5j5ZLYlWLm3B0ieCh3njK4Kcls9A4Ffuai0BnSjehzgW/mzlI1vmIwAyTW\noljx610R3Jg7j7WGxIYUnQ2Oj+BXufM0m68YzEanC3gCuClzDmuhCJ6gWCDsBxKb5s5TFS4MZoWj\ngTM8xXbniaCbYtzKf2eOUhluSrKOJ/FC4M/A1hE8nDuPtZ7E+sDNwMkR/CR3nmapbFOSpCmS5kta\nIOnEAd7vkvS4pDnp8ZmyM5n1cxjwaxeFzhXBcuDdwLckts6dJ7dSZxmUtDbF5dkbgQeAGyXNjIg7\n+216dUQcWGYWs4GkAU7HAJ/IncXyimCuxJeA8yT26uS5ssq+YtgdWBgR90bEs8B5wFsH2M7NRJbL\nG4B1gd/mDmKV8F1gEfBfuYPkVHZhGAvc3+v1ovS13gKYJGmepFmSti85k1lvJwDf9U1ng1UzsL4P\n2F/iyMxxsil7wYpGTrabgfERsVzSVOAiYLtyY5mtWr5zMnTuB4CtLoLHJd4OXCVxRxoI11HKLgwP\nUKyb22M8xVXDKhGxrNfzyyR9X9KmEfFI7+0kTe/1sjsiupsf1zrMh4Afpb7sZqtEcJvEMcBFEntE\n9P3cqipJXRRjcka3nzK7q0paB/hfYD9gMXADcGjvm8+SNgcejIiQtDvw84iY0G8/7q5qTZUWarkb\n2D2Ce3LnsWqS+BTFVCl71/EXiJF+dpZ6xRARKyQdD1wOrA2cERF3Sjo2vT8DOAQ4TtIKWNVlzKxs\nxwOzXBRsCN8AtgV+KXFABP+XO1AreICbdZw0P849wBsimJ87j1WbxDrA+cAzwGERrMwcqWGVHeBm\nVkEfBK5yUbBGpNl2DwW2AGakCRfbmq8YrKOkq4UFwJQI5uXOY/UhsRHFetG3A8fV4crBVwxmjfkE\ncKWLgg1XBMuAqcCrgR9KrJ05Uml8xWAdQ2JLirn3XxvBX3LnsXpKV50XAn+nuOdQ2akzfMVgNrQv\nAqe7KNhopG6rbwFWApe34zoOLgzWESR2oziZv5o7i9Vfukp4N8XCTn+U2DZzpKZyYbC2JzEGOB34\neASP5c5j7SGC5yL4BMWa0X+QODh3pmbxPQZrexL/TjFNwFRPlmdlSFekPwMuBT5ZlfsOI/3sdGGw\ntiaxPXANsGsE92aOY21MYhPgTGAr4J1VGFXvm89m/aTlGn8GfMpFwcqWmikPBn4M3CDxkbp2afUV\ng7UtiR8AGwBHuAnJWknilcBpFP//jolgTp4cvmIwW0XiAxT3FT7oomCtlqZb2Qf4PvBriW9KvDBz\nrIa5MFjbkZgCnAy8JY1WNWu5CCKCs4AdgBcAd0mclAbIVZoLg7UViV0p2ngPjuCu3HnMIvhbBP8C\n7AnsCCyQOCHdA6skFwZrGxKTgFnAByL4Q+48Zr1FcFcEhwJTKJqZ/iLx/6Q+q1xWgguDtQWJNwEX\nA0dGMDN3HrPBRDAvgoOA1wPrAfMkfi6xf1V6MblXktWahIBPAh+l6Dt+beZIZsMi8QLgvemxJXAe\n8FPg5tF2nPAAN+s4abbU04CXAm+P4P7MkcxGJXVzPQw4PH3pV+lxdQRPD39/LgzWIdJSi0cDX6Io\nDF+uyhQEZs2QroQnAtOAN6fn11KM4r8WuKmR9acrWRgkTQFOAdYGTo+Irw2wzXcoFr9YDrwvIlYb\nCOLCYLBqMrx3AJ8HFlFMijc3byqz8km8CNgX2AvYG9gOuBn4U/pzDjA/LUPa6/sqNsBN0trAf1Pc\ngd8eOFTSq/ptMw3YJiK2BY4BTi0rTytI6sqdoRF1yynxSokvAPcCHwA+DLyxCkWhbj/LqnPOgUXw\ncATnR3BCBK8FxlKM1VlCcUVxAfC4xM0SP5HYejTHK7NX0u7Awoi4NyKepbih8tZ+2xwInA0QEbOB\nTSRtXmKmsnXlDtCgrtwB1kRic4kD4PCvStwGXAlsTLFO874R/KZCo5m7cgdoUFfuAA3qyh2gQV05\nDx7B4+k8+FoE745gO2AL4FjgtxSry43YOs0IOYix0Odm4CLgdQ1sMw5YWmIuy0jiecCLgE3Tn5sD\nWwMvT39uRzG/zE3w1HLg34Ab67DwullOaZT/jekxKmUWhkZ/o+vf/lX6b4IS/0Lfq5f+GdTge/1e\nf2jrNMiqwe1Hc6zR7PsjW0nsX86+V70eAzyfop/2er2eAzwCPJz+fBC4h6Kt9OfAQuDeCEK6YHoE\nszGzlirMFTsEAAAE2klEQVTt5rOkPYDpETElvf40sLL3DWhJpwHdEXFeej0fmBwRS/vtqyrNBmZm\ntTKSm89lXjHcBGwraQKwGHgXcGi/bWYCxwPnpULyWP+iACP7i5mZ2ciUVhgiYoWk44HLKbqrnhER\nd0o6Nr0/IyJmSZomaSHwJHBUWXnMzKwxtRjgZmZmrVOpSfQkTZE0X9ICSSeuYbvdJK2Q9PZW5kvH\nHjKjpC5JcyTdJqm7xRF7Mqwxp6TNJP1a0tyU830ZYiLpTElLJd26hm2+k/4e8yTt3Mp86fhrzCjp\n8JTtFkl/kDSx1RlTjiF/lmm7bOdPOn4j/+ZVOIeG+nevyjk0XtJVkm5POT4yyHaNn0eRVpPI/aBo\nbloITKDo0TIXeNUg2/0OuBQ4uGoZgU2A24Fx6fVmVfxZAtOBr/ZkpOgltE6GrHsDOwO3DvL+NGBW\nev464PoKZnw9sHF6PiVHxkZy9vq/keX8GcbPM/s51GDOqpxDWwA7pecbAv87wPk+rPOoSlcMjQyI\ng2LU6y+Av7UyXNJIxsOAX0bEIoCIeKjFGaGxnH+lWFWK9OfDEbGCFouIa4FH17BJ9kGQQ2WMiD9G\nxOPp5WyKsTgt18DPEvKeP0BDOatwDjWSsyrn0JKImJuePwHcSTFLa2/DOo+qVBgGGuw2tvcGksZS\nfMD1TJ3R6hskQ2YEtgU2TZd2N0l6T8vS/UMjOX8IvFrSYmAecEKLsg3XYIMgq+poisWCKqcC50+j\nqnAONaJy51DqBbozrDb+Z1jnUZndVYerkf+kpwD/HhEhSaw+sKpsjWQcA7wW2A9YH/ijpOsjYkGp\nyfpqJOdJwNyI6JK0NXCFpNdERBXXSG75IMiRkLQP8H6KJRyrKPf506gqnEONqNQ5JGlDiqvBE9KV\nw2qb9Hs96HlUpcLwAPRZ4m48RVXrbReKMQ9QtOlNlfRsRLRqxa5GMt4PPBQRTwFPSboGeA3Qyv/U\njeScBHwZICLulvRn4BUU40+qpP/fZVz6WqWkG84/BKZExFDNObnkPn8aVYVzqBGVOYckjQF+Cfwk\nIi4aYJNhnUdVakpaNSBO0roUA+L6/IeNiJdHxMsi4mUUlfG4Fv+nHjIjxfKSe0laW9L6FDd67mhh\nxkZzzgfeCJDaGl9BMTVF1cwEjoRVo+kHHASZk6StKGa3PCIiFubOM5gKnD+NqsI51IhKnEPp6u8M\n4I6IOGWQzYZ1HlXmiiEaGBCXNSAND9qbL+nXwC3ASuCHEdHS/9QN/iy/ApwlaR7FLwifiohHWpkT\nQNK5wGRgM0n3U6y1MKYnZ1RgEORQGYHPAS8ETk2/jT8bEbtXMGclNPBvnv0caiQnFTmHKJoujwBu\nkdSzns1JwFY9WYd7HnmAm5mZ9VGlpiQzM6sAFwYzM+vDhcHMzPpwYTAzsz5cGMzMrA8XBjMz68OF\nwczM+nBhMDOzPlwYzMysDxcGsxGS9BZJl0hakqaIPi53JrNmqMxcSWZ1kibOOyAiDpB0EMX0Mhfm\nzmXWDL5iMBuZI4Fvp+cvAh5fw7ZmteLCYDYyLwTuS89fD1ybMYtZU3l2VbMRkPQqiumORbGw+q2Z\nI5k1jQuDmZn14aYkMzPrw4XBzMz6cGEwM7M+XBjMzKwPFwYzM+vDhcHMzPpwYTAzsz5cGMzMrI//\nD93U8pTbbuu4AAAAAElFTkSuQmCC\n",
"text": [
"<matplotlib.figure.Figure at 0x1095b8c10>"
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we'd hoped, these are pretty close to the true underlying values ($a=1.0$ amd $\\sigma=1.5$).\n",
"\n",
"Remember that these are slices in likelihood - no more, no less - they are not constraints we have on the parameters. They are also not the likelihood for a single parameter - they are likelihoods *after adding an additional, very strong piece of information* (the exact value of the other parameters). When you look at something like this you have to be very careful about what you're question you're asking, and usally when you want something more useful result is that you have to marginalize over the other parameters in your space.\n",
"\n",
"In particular, if our slice had been taken through points in the other parameters that were not the truth then we might not expect particularly good answers.\n",
"\n",
"Let's try to visualize this and do something slightly more sophisticated - a 2D plot of the likelihood, on a grid, in the $a$ and $b$ parameters. This is still a slice, but a 2D slice instead "
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#This is a handy numpy function to build two 2D grids of the parameter values\n",
"a, b = mgrid[0.0:2.0:0.01, 2.0:4.0:0.01]\n",
"#The T means transpose - numpy very irritatingly swaps the order around. I'm sure there are reasons.\n",
"imshow(exp(logP(a, b, sigma, s)).T, extent=(a.min(), a.max(), b.min(), b.max()), interpolation='nearest')\n",
"colorbar()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 6,
"text": [
"<matplotlib.colorbar.Colorbar instance at 0x109821b90>"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAATAAAAEGCAYAAADi2MlDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+MY9d13z9nSQ7n91BDZma91qYbxHJqN22jIJBV24U3\nbQNEaqEWhds6gJvUKFDDqBAjQIqihlFXQICgfwQ13MqCgtqGnAA2UDtV1ESCEaRexW4QpbFW8g+5\nqJVEhSV5Z0zOcnZ+/+De/nHvfe++Rw7J4QyH75HnAzzx1+PbS5D6zjnnnh9ijEFRFCWPXBr1AhRF\nUQZFBUxRlNyiAqYoSm5RAVMUJbeogCmKkltUwBRFyS0qYIqi9EREPisiayLyrXO41k+JyB+LyLdF\n5GUR+afBaz8mIi+IyPdE5IsiUup2LRUwRVH64XPAz5/TtXaAf26M+Ul3zU+KyKJ77T8Cv2GMuQ+4\nDfzLbhdSAVMUpSfGmK9hBSVCRH5cRJ4TkT8TkT8SkZ/o81rfM8b8ubv/A2Ad+BEREeBngS+5U58C\n/lG3axVP+TkURVE8vwl82Bjzqoi8C/g08HdPcwEReQCYMsb8uYjUgKYx5q57+Q3grd3erwKmKMqp\nEZF54G8B/80aTgBMudf+MfBYh7e9box5KLjGW4DPA7846DpUwBRFGYRLWGvp/vQLxpjfAX6n25td\nzOv3gI8ZY/7UPd0AKiJyyVlh92KtsK6LUBRFORXGmDvAX4rI+wHE8jf6ea+ITAH/Hfi8Ezt/TQN8\nFfgn7qlfAp7uei3tRqEoSi9E5AvA+4AasAb8e6zYPAG8BSgBXzDG/Fof1/og8FngO8HTv2SM+aaI\n/BjwRWAZeBH4oDHm6MRrqYApipJX1IVUFCW3qIApipJbLmwXUkTUV1WUEWGMkd5nncxp//8967/X\nLxecRvGJi/3nLpQbwPURr2GY3KD35ythf1IlbAx21T3/V4B3wl8FHgTea4/LP/EX1GgAUKBFgyqv\n/79r8Ccl+Dr2eAngL4FX3bXeAO4Ae8AxcGJ895TcYHy/v04pWaenZ3Te8fFz+df6Q/PAlCFQTN3O\n2Jt5oGKPS5UdZtljll0AjilQoAXHRatL++4AYrFSRknXquoRoQKmDIHQEgOYhWliAavB7Pwus+wy\nxaE7Z8reHEssXscAhqSldeQen6f1pfRDFsUii2vKKddGvYAhc23UCxgy10a9gMwzM+oFdEAF7Ny4\nNuoFDJlrfZ7nra4Z4p/8TMJ9pALzc1sssMUUBwC0KNhTj4Nj3z9xxPCtrWtDvn7+URdSmRC8++gF\nbCElYEcssMUMu5SdC3lI2YpYGP/ahzj+5WNgGgsbFVkUiyyuScklPu7lxSsUMLHC5URsumKtrzKH\nkQVWYNae2iZgnawvjX2NArXAlAnBC9lM/DCwwBaWthM7kJ5jCknxOob2dAkVr1GRRbHI4pqUXBNa\nX6n0iZo9FthyLuReHPsCWhST8a+EBaau46jJogWmpUTKOeJ/4t76ckci/rXPAlvMskuZAwq0bP4X\n0GqlLDAgGf8iuK+W2EVT6vM4CREpiMhNEfkfJ7z+KTfM42URaesz1gm1wJRzJoyBubjWPJGILdWa\nVGgyzxZTHFq30dE67iRge7THwdQaGwXnkEbxUeAVYCH9gog8DLzNGHOfa0/9BLZuoytqgSnnQLcA\nPkn3sbzlLLA9yi6ADzaN4mB/ygrXNin3Ua2tLFDs8+iEiNwLPAz8V6BTneQj2CEeGGNewHZmXe1w\nXtuaFEVRenLGGNh/Av4NsHjC628Fvh88fh3bUnqt20VVwJRzJoh/TbunghhYhSaz2DKiAi2KPv5F\ngdZxsUMKRZjZ6p9TRsFJYvENbOvUkxCRfwCsG2Nuisj1bqemHvfsgKECppwj3n0sAotWuCAu4K7t\nRDuQcQ2kpUWBI+9CJlIoVLCywkkW2IMkg1WfaT/l3cAjLs41DSyKyOeNMeE0ojeAq8HjngM9QGNg\nyrnRYQcyCN5Tg4oL4PtdyEIQjLc5YOWUBRYWbWssbNQMGgMzxnzMGHPVGPNjwAeA/5kSL4BncOPV\nRORB7MSjru6jX5OinBNhCVHJChfECawFa33NuxrIFoVoF7JFEfYlJWC7aA5YdjjHPDADICIfBjDG\nPGmMeVZEHhaRV4Ed4EP9XEgFTDkjXXYgQxeyRmB97VHmkEPK0VUOSe1AHkPnGkhtozMqzqMbhTHm\neeB5d//J1GuPnvZ6KmDKOeLdx9nYdQSo2fhX6D5OcUCBWWt5AQcE7uO2v16nHDBlVGQxE18FTDkn\nUiVEKQGz8a/b0S5kIdh9jG699dWWA6YuZBbIolhkcU1K7ggD+EVgIS4dwt5WCk0W2G6rgfS3u8x2\nSGLV9IksUepXLS7w740KmKIofVFUAVPGl3AHUqwLWXMv1XDu421mXBLrLrO2fMj1wj88mEq5kGH8\nS62vLFAq9D7nolEBU85IuAPp4l9Fki5kbZ8KPgdsmykOYuFyO5HtdZDpLhTKqOnbArtAMrgkJR+k\n96S8iM0m+39hO1Dc43YgfSfWPVqukY79s76/PZuywHZpbycNao2NjlK59zkXjQqYcg6EO5BBAN8J\nWKVsra973A6kLyM6pmDTJ8Bm4bflgIXuo1pjIyeDapHBJSn5Ir0DOdNmgVWpRz3AZoI20oeU2fPp\nkdvSIQdMyRQZVIsMLknJH6kAftT/y1pPcfzLZuEfuvhXK7TAtjmhD5hm3meGDKpFBpek5INicBsE\n8P0E7hrM15oA3BME8MtBDeQBU+x6C6xrDpiSCXQXUhkPTgrgLyZ7f81ZAavSiCywGXY5oEyLIoeU\n43pIL17bYN1HLSPKHBlUiwwuSckXHUqIXPvoKg3Au5C320qIrAXm+uYnXEgvXmqBZQrdhVQUJbdk\nUC0yuCQlH4RuZNCFwqdQXLaWFxDFvxbYSlhge8yyt+N3IekQA4NkU0NlpGRQLTK4JCVfhEF8oh3I\nUu0ONeoA1KhHZUSeA6asC7kduJBtraRVtDKFBvGV8SO1A+mbF1abKQtsi/LBIa1igVah4HLAZrm7\nnYqBAbHFpTGwTJFBtdCe+MopCTuwpuoggwz8GnWqNIIdyG1md+5GV4kC+NviDgL30ZcQ+cdKJhiw\nKb6ITIvICyLykoi8IiK/3uGc6yKy6SZ33xSRj/e7JEUZgFJwm2pi6HYgq86FrNBkYWcb2YfCXIuD\nwhR7zNocsKa7TGSBaReKzDKgWhhj9kXkZ40xuyJSBL4uIu81xnw9derzxphHLmBJiuJJBfBdBr6v\nfQQrYNM74Adxtyhy4FzIyG1MpFCo65hJzpBGYYzxAdApbDRto8NpnSZ2d0VdSGUA0oM8giEeFVi6\nbK2v0IVkh0iXDpliz3UGiyyvtjKicKCtWmKZYNC5aoCIXBKRl7CTtr9qjHkldYoB3i0iL4vIsyLy\nzn6XpCgD4F3IdvexUm5S88IFLOxsWwEDWsU4gXWrtZC0wCIXUkUrk5ywC3ljHW78sPtbjTF3gZ8S\nkSXgKyJy3RhzIzjlReCqczMfAp4G3t5rSSpgiqL0xwlqcf2KPTyPpW2rAGPMpoj8PvAzwI3g+a3g\n/nMi8mkRWTbGdHI1I9SFVE5Blx3IRAlRPQriV6kzvUlkgR27FIpdZtnddkH8Js59BA3eZ5jBdyFr\nIlJx92eAnwNups5ZFRFx9x8ApJd4+SUpyinx9Y/Q5kJehhoNl7zqthg3sAJVtgH8XWbZY4b9ZsqF\nxJCsg1QhyxSDJ7K+BXhKRC5hjabfMsb8YTiZG3g/8BEROca24/1APxdWAVMGxP90kjuQl2o7UfD+\nnpYTsGAH8pApDpliiwXYLqUELGwhreKVOQZPo/gW8NMdnn8yuP848PgFLUmZXLwbGVpgxPlfq/EO\n5GLDidAmkYvoA/i7zFrX0QtYlELhLTAlc0yPegHtqIApAxCkTqR2IGvU3Q7k7TjTx8W/KNo20ltu\nvEeifGgb2kuI1BLLFFoLqeSXdAA/sMDaMvAb1Gjg2oFZASvbt/kurFssdLDAdtE6yAyTQbXI4JKU\n7JOywKIMfKKdx+rOhnUdwYpTwR57zLLNgi0jCi2wRCNDJZNkUC0yuCQl24SWGIRzIKcvb7gdyAbT\nG8QCtgPM2bf4+Nf2TsoCA5LWlwpZ5sigC9k1D2yYVeSKouSMM5QSDXNJJzLMKnIlj4Q7kDPx0z6A\nv2TjXyus2fjXHfe6D2eVieJf281UEL+tjY6SOTLor/Vc0rCqyJW8Eg7xIJHAWqXBKms2gL9BLE4H\nQBGO5mwMzAbwp1MupLbRyTwZHOrRs5RoWFXkSt4I/YMZYJFojFq6hGhnw1pgOyS7UExfilMovHjt\nk2qjE1pgKmSZIm8uJAyvilzJE+kBHr4HGMEQjyNqrv4xCuD7/K9prPtYnnUCNh8LWJsLCTrUNqPk\n0YX0nE8V+Y3g/jV3KPkh6L7q5zkGU7itC7lubfUNovIhpnE7kDPtSaxREbcO8jg/XnPHOZPBXciu\nAiYiNeDYGNMMqsgfS52zCqwbY0zvKvLr57FmZSR0D+BX52z8q0rdildQPsQSMB3HvyIXsklgaGkR\n9/lxjaRx8Pz5XDaHFtjQqsiVvBK0z4EgA9/Gv2qtBqyTKOCmAMwRiVcyB8y4kzT7PvPkTcCGWUWu\n5A1vfRWJklfBBvBdC51V1lhcO7IB/E1iPSoC5UDA6pWgB5iPe6kLmXny5kIqSixaEAfvF1IW2BEr\nrFPFWV+bxNaXf5uzwJpUoFlKtZAG3YHMAdqNQlGU3JJBtcjgkpRsEiawugJugJqdQlSjzipr1gLz\nSazeSisDS7QH8BMWmLe4NIUis2TQhdSe+Mop8C6kJKZwV8uN7i7kNOzPwRbzyRywRBND7USReYY4\nmdud9ykR+Z5Lir+/3yUpSg+C4D0zNhYSWGC2iWGd1Y3NOIB/QPwXuwxbc/Nu1O09KQHzoqUpFJln\niJO5ReRh4G3GmPtE5F3AE8CDQ1qSMnmkdiC9gLkayBXWEe8++vIh/+uKrK+F9kaGiSC+kmnO4EL2\nUVP9CPCUO/cFEamIyKoxZq3bdVXAlC6ku7AG8a+aO6V2xAprNv61RlwD2SZgdgey2agEAhbmfmkK\nReY5wy6kyyV9Efhx4IkONdVvBb4fPH4duBf7qzoRFTClD0rBsRAXcOMD+K6FTmiBlYm7F8zh3McK\nR83QArtDbIGFqJBlkpMmc78IN252fs3TR001tHe1MfRABUzpk1QA3wmYDeCvxRn4voSoTPwXeyl2\nIalLKgcsvfuormRmOWky9wP28Dz2uZMvcVJNNfAGcDV4fK97riu6C6n0IFUD6QP4roWOTZ9YZ/HN\no/YdSG+FLcUWWDKFwruQGrzPBUOczA08A/yiO+dBoNkr/uWXpCiK0pvB1aJnTbUx5lkReVhEXsUG\nIT403CUpE0SQwBo0MARYYT0ZwPcWmCsfAjBL2AA+lQ45YKELqVZYphlwF7Kfmmr3+NHTXlsFTOlC\neg7kYtBC2oqN70LBm8RdKPwOpBOw28vT7QIGxAmsSi7IoFpkcElKNiil7idTKJYu26m1q6zHO5Bh\nF1YX+4IghWIniIG1tc9R6yvzZLAnvgqY0gMfxJ+x950LWS17AVtjdeeHUQrF7g7MzmF/7M4C89bX\ndiKF4qTyIRWyzJJBtcjgkpTskNqBLJLYfQSbhT/t3cdN2Nt3AjYPVO1Voh3I+nRKwEIrTFMoMk8G\n1SKDS1KyR9CF1Y1QW2EdIO5A4QL4xy3i+NeifXeTCrfDFIp9sM17NXCfKzKoFhlckpIdwh74gYDV\njO19jxMyZ4Ht7jg5ci2kYwvMBfDrdEhiDUuJlCxjMthORwVMSeGD92EXVhcDcwH8+cs2eRVI9ADb\n8gI2DczZ9AmABtW4C0XT/zvaPidvtDKoFhlckqIoWUQFTMkZYRfW2XgH0o1QA1jZ3IgC+HdwPyiX\nQnF72RZD3qZCs9WpC4XGwPLEQXmqzzMPh7qOEBUwpQepKdzBCDWA0g+ANTAN6xQugI1/LUPDpes3\nqdCsBzEwTaHIJa1C9oJgKmDKCaS6sAYlRKu+fAiiAP7WjpWkZbAC5sqHwArZ3frcCSkUYRqFkmVa\nGWyKrwKmdKFDF9agAytgayA3oHFgJWkm2IH0ApYoIdoGTaHIJ8cqYEr26bD7GFhg05c3bPZ9aIGt\nwZZ7V6mILSFahrrLo4hSKDpaYCpieaGVQbnI3oqUEVLq8Nh1YXVTiGpLtoHhyqZrab4O3LEB/GNg\nZhorYIEFVqcaCxigRdz5RF1IJUeENZCScB9rNGzwHmAd7rgA/gxQcvGvnZVLQRD/ng47kNoDP2+o\ngCk5IExkbe/CWqVuu0+86U5bt/EvL2B+B7JZrtDwLmRikMdJ1peKWdY5oN80iotDBUzpQqoGMtqB\nXI9nxTTsHI8jAgGrWrfxtnMhj+qLXXYgdfcxL2QxBqY98RVF6YsWhb6ONCJyVUS+KiLfEZFvi8gv\ndzjnuohsishNd3y8nzVlT1KVEeN/Eu1tpC9d3unoQm5h7agSRDuQNlLm+k4ndiA1hSKvnCEGdgT8\nijHmJRGZB74hIn9gjPlu6rznjTGPnObCKmCKwwftPT6NYjGOf602rAvpGxgCxpUQgXMhF4l2IH0M\nrD2FQkUsjwyaB2aMuQXccve3ReS7wBUgLWDpuZA9UQFTUpSC22QXVjtCbS1uYAisb8Qh+UWwLXRW\nfAcKGwNr70Khca88ch4xMBG5BtwPvJB6yQDvFpGXsfMgf7XD9O42VMCUE3DuIwQ7kI1E/y+wfQx9\nDeSCC+AfrUCdGvXQhdQUitxzkgv54o0tbt7Y6vhaiHMfvwR81BiznXr5ReCqMWZXRB4Cngbe3uua\nKmBKQOhGBjuQUQ2ky8APBMzvQJaAWR//WlqyFlgjsMASKRRh/aMKWV44PCGN4ievV/nJ69Xo8ece\n+0HbOSJSAr4M/LYx5un068aYreD+cyLyaRFZNsZsdFuTCpiSIswDm40y8Ll85CywoIU0cQA/TKFo\nuH4VR/Wop7RzIdNzIJU8MWgMTEQE+AzwijHmkyecswqsG2OMiDwASC/xAhUwJSI9G34GWIgssKXL\nVrxqrQasw5H7aSUC+Mu4Gki3A1l3L9YJ+uBrJ9a8coYY2HuADwLfFJGb7rmPAT8K0YDb9wMfEZFj\n7A/lA/1cWAVMoX0GpL+diSywatkmRiy+eQRvwtqmPSshYEtEAfyGr3+EE/rgq4jljUHTKIwxX6dH\nzqkx5nHg8dNeWwVMUZS+0FpIJeOEQXxXxB2kUKwEAXwfnNhz74pSKFaDFApvgWkKxVig/cCUDBMO\nsYVEEbfrQrHqUygasYD5AP5yAVgGE6ZQeOFqgnUZvQup7mMeOaQ86iW0oQKmpAgELEihsALmdiDX\n49iXL+JedvGv+vJ8nEIRxsCijmHhoUKWJ9SFVDKOz76HsAtF6fIdVlmLkliPNmIBc2cironhOis2\nheLWYmyBdR3koeQFdSGVDOJdx3QaRVADWbVTiH5kfRvWYGMzbiHtzox2IBMpFG0CpuKVZ7LYTid7\nK1JGRDgD0j2OAvhJ93GNZEvCRbA5YOkUCi9g+6ApFPlHXUgl4wQzICFVA5ncgfQC5tJdYdUea86F\nTPbABxWv/KMCpmSU1AxISHRhtV0o1qMSojvEMrSMzZ5gGY7e4vqAtWpB/SN0zsBXIcsbuRMwEZkG\nnscOi58CftcY8+86nPcp4CHsL/VfGGNups9RFCXfHOQtjcIYsy8iP+taXBSBr4vIe11pAAAi8jDw\nNmPMfSLyLuAJ4MHhLls5f0IrjDiFwhVxV6lbF7LDDqRPoYi6UNQrJ4xR0z74eSZ3FhiAMWbX3Z0C\nCsQ5jJ5HgKfcuS+ISEVEVo0xayg5INx9DIL4rgZyvtZkhTVWD6wL2WgkdyBngNIyUfxrjVXu3ppL\nuZBevNRtzDO5FDARuYRtNvbjwBMduiS+Ffh+8Ph14F7iuTVKZkkXcUch+cgCq87ZIu65N+/aFIpW\nhx3IqAtrLbkDmRAwFa+8k8s8MGPMXeCnRGQJ+IqIXDfG3Eidlu5lbc5pfcqFkBpiC+1zIF0Khd+B\n9NK37P+z0mEHMvIUNYViHMh1HpgxZlNEfh/4GeBG8NIbwNXg8b3uuQ6Eb7vmDiUbOAvM/yKiHcig\nBnIj7sDqWhVaAVsBrgSTiKL4l/87pi7kxfKaO86X3LmQIlIDjo0xTRGZAX4OeCx12jPAo8AXReRB\noHly/Ov6WdernDthEbcrH4JEDliVeiKFAuJsseUCsAr7zgJrHlSCHvg+WqYpFBfLNZLGwfPnctXc\nCRjwFuApFwe7BPyWMeYPReTDYDspGmOeFZGHReRVYAf40HCXrJw/gQsZCNilyzuuD75Nvz/aSCaw\nQrwDuTb3IzSosXmrmmohDTqJezw4OKEn/ijplUbxLeCnOzz/ZOrxo+e8LuVCCHcg3RQiN4eDGlRq\nTWuBbW7a+FcwA9KF+hEXwF9nlXVWoF5KDbEFdR/Hg0FjYCJyFfg8NthggN80xnyqw3mnzifNXlRO\nuSA67UDOJiywaqFOlTqlH2BTKIgncEczaFaIaiDr1Oz40sQUItBBtuPBMCdzD5pPqgKmKEpfnKEn\nfj+TuQfKJ1UBm3hSnVgDFzLagVwj6sK65872LiTLwBVYY7W9C0VkgWncaxw4jzywLpO5B8onVQGb\naFI7kL6FNAQ7kHaMGhvxnuIMtvkEuDtXOnWh8IF7f1/dx7xzUgzsBzf+Lz+48b2e7+8xmRsGyCdV\nAVNIpFAkBKwe9wELUigWgeU592AFNlambQ5Yo5bagUxbYBoHyzMnuZAr19/ByvV3RI9vPvZc2zm9\nJnNzqnzSGBWwiccH8JMCVrp8xxUGuUG2d5IpFLPL7sEVuwO5xoqdxB3lgIVNd0DdyPxzOGAaRT+T\nuTlVPmmMCthE0qmN9GLcgQKoVG0KRW0jbiPt41/LEG9Drvj4V6cdyDAPTMk7Z4iB9ZzMPWg+qQrY\nRBO2kS5FHSgANxitjrj4VwMrQ9EMyBV3iSt2kEcigJ9oIQ3qOo4Hg+aB9TOZ25136nxSFbCJJ2gj\n7cqHwLeRjncgt4hnQK5ALGArUKfKGqupGshwiK2K1ziQx1IiZWzpUAOZELB6XAO5kayBrJaJtiH3\nXQyscRDuQJ40wFaFLM+ogCmKklty2Q9MGWdSJUSBBVZzzXFYBzbjcPwCsOjqH8EWca+xyma90iGA\nH6ZP6C5k3sl1PzBl3PBffdKFnL5sO4ZXsLuQNMAEO5DRDMgr9t3rPgP/1nQqhUID9+PGoGkUw0QF\nbKLwBdzF4DaZA1ZZsnVAUQrFOmztxDuQUQNDb4G5PvjJNtJhB1ZQIRsP1IVUMoRPoXCBfOc+3uMK\nGRMpFAfxDuQyROVDEFpgBG2ktQf+OKIupJIxghQKP8SDBkCcQrGZzMCvloEVME7A1lg9oQZSRWzc\n0F1IJSOkJnH7AH7Fpk+AjYGxQZQDBjb+5QP4a8tLgO0DFtVAdgzgg8bDxgMVMCVDBCIW7EBWAheS\ndWAn1UInap9jg2BrrHJ0azGwwMIaSN15HCdUwJQREwbv25NYpy9v2NQJiHYg2ewcwF93maztcyC9\nBaYW17hxQHnUS2hDBUxRlL5QC0wZIaUOj9tTKHwQv7q5CRtgdqwjWMQ1oIgaGFoLbM0WQ3Yo4ta4\n17ihAqZkhLCdTjKFouaC+KUNYNPmgIEN9/sRar4DBZBso4NBxWt80TwwJUN4C4xECoUP4kclRAf2\n4QJQWgZW4Icr83YCEdDYSc+B7BT/UjEbBzQPTBkxYSZ+qgtFJSgfgiiA7yc7LoL1IVes9eUtsO1b\ntdQkbp9CofWP44a6kEpGODmFwueB4VzIPYIM/GACkY+BUZcghULnP44zKmDKCPExL38/aYGVaneo\nUY9dyAa2sa87e9Fl4PscsIbvKd3WRlpjYOPKweHAPfE/C/x9YN0Y89c7vH4d+F3gL9xTXzbG/Fo/\n11YBmwj62IGsNqnQ5J6WE7BNIKiBXKxidyBXbOC+7gWsbQ7kHsp40joeWC4+B/xn4PNdznneGPPI\naS+sAjZRhDGw5A5khdtUabC45iynO0QW2KL/zwpsXJm2FljLNQ5rq4HUGZDjSut44MncX3MDbbuR\nngnZFypgE0lqB7JiUyii+keATTD7Qbl3FMBfZZ1VNm4FFljUAyzdxFCFbJwYVMD6wADvFpGXsbMg\nf9UY80o/b1QBUxSlL46POguY+V9/hPnjr53l0i8CV40xuyLyEPA08PZ+3qgCNjH4nUd/PznIo+oi\nW2y6U3bguOVqIOdI7EBGHVjhhDmQmj4xjtxtnSAXD/4de3h+49dPdV1jzFZw/zkR+bSILBtjNrq9\nD1TAJowwBuZSKHwQ37uQLg2MbTg6dmfOEQXw18MOrNDHJCJlbBiSCykiq9gdSiMiDwDSj3iBCtiE\nEE7ghrQFdqm2EwuYn5+2A8fHtoQIl4GftMDceW1dWDX2NbbsDyYXIvIF4H1ATUS+D3wC99fUGPMk\n8H7gIyJyjE0m/EC/11YBG3tKHe6nUihqzVjAvAvpSohmyoCrgfzhyjxrrNgyIm+BbUNykK2K19gy\nYGTAGPMLPV5/HHh8kGurgE0MYQzMpVC4GshKoUmF21R2NhMCViw6AQtKiBrUaKxVYwFrgi0hSruP\nKmRjRwZDmypgE0EpdetSKFwN5AJb3EOT6U2i3C/2oVQE8QH8Ves+rrPC3VtzcfJqoohbayDHmgx+\nrSpgE4MvH4LIfUwH8DeIBewYigWs+1gFrtjuq5H7mHAhtQZyIsjg16sCNjGELuQsTJPIwo/iXzvx\nO2SaKP519JbYAov7fxE0MVQRG3tao15AOypgiqL0h7qQysUTplB0ciGPkjuQ+8Fbgx3IxtJScgcy\nKuAOdyCVsWa/9ykXjQrYWJNOoXDB+yCFYr7WjOsgd4jSJwArYC4HbI0VO4Xb70Bu+5PCLqzqPo41\nGfwbpQI2EYQ98N1jL2BzW1RossCWFTD/Iy1i42RLwKot4m5Q5W59LmghDZ27sKqQjSUqYMrFE6ZQ\nzMRPuxSHGGFQAAAOk0lEQVQKb30tHHQQMNdChytQ9zuQt0hZYBq8nxhUwJTREHRgBWtZOQtsgS0W\n2GJu864VML/T5ONfVdhZueTa6KykhtiCDrGdIDL4NauATQTehXTWmA/g1zg5gB+kUKyVbQF3ewAf\ntAZygtA0CuViSe9ABhZYMIkoin+lA/h+B5IaDapxAD8a4AEnj1JTxg51IZWLI90H3xVwQ9yForIT\nuZBR/Mt3TPE7kKt2B3KNFRvAj+JfYe+v8JetQja2ZDCN4lK3F0Xkqoh8VUS+IyLfFpFf7nDOdRHZ\nFJGb7vj48JarKMrIOO7zuEB6WWBHwK8YY14SkXngGyLyB8aY76bOG2iiiHIRlEjMgYQoBrZQSaVQ\ntIh/EXNAFUzkQtY6jFALb7WIe+zJ4NfbVcCMMbdwreuMMdsi8l3gCpAWsIEmiigXRTSaw+KLuMtW\nvOZDF9L/IpaAZagvz8dzIH0Rd9SBAjSBdYLIoIB1dSFD3Fik+4EXUi9FE0VE5FkReef5LU85G+Hu\nYymO5wdtdBbYYpa9OL5RJg7gr8Y9wNp3IMPgve5ATgRHfR4XSF8C5tzHLwEfNcZsp172E0X+JnZ4\n5dPnu0TlbKSG2AZtdLyALexsx1vk08QpFMtQp2YtsEQJ0RHJHmAqXhNBq88jhYh8VkTWRORbJ11a\nRD4lIt9zhtD9/S6p5y6kiJSALwO/bYxpE6fTTRS5Edy/5g5lOHToge8GCYUCNs8W0zvEFticu41S\nKOy8oru35gILzM+ABBWvLPKaO86ZwXchu07mFpGHgbcZY+4TkXcBTwAP9nPhrgImIgJ8BnjFGPPJ\nE845xUSR6/2sSTkzYQqFFzJngYFzIU0cwD8gDuCX3TnLcGe1lIx/JQL4Xri0BjJ7XCNpHDx/Ppcd\nvCd+r8ncjwBPuXNfEJGKiKwaY9Z6XbuXBfYe4IPAN0XkpnvuY8CPun/sTBNFlGHjdyADFxJgHkoV\n5z6yHQfwfewLoAr1gp8WWe2wA6kW2MQxvK/6rcD3g8evA/cCZxMwY8zX6REnO8tEEeUi8EF8SVhg\nCxUfwN+NXYMisQu5bNMn1lih0QoC+IkRaqCW1wQx3FKidCaD6edNmok/toQ7kEEGvrtdKNj41yy7\ncQnRHLEF5nqANaixcasatNDxAfwM7qkrw+Wkr/z1G/DGjbNc+Q3gavD4XvdcT1TAFEXpj5ME7PJ1\ne3j+9LHTXvkZ4FHgiyLyINDsJ/4FKmBjTpCF79MnINqBnGWPGXbjBFZf/wjsVC9FRdzcmk4lsIYu\npLqPE8OAX3WvydzGmGdF5GEReRUbkf1Qv9dWARtLwhQKV0Lke4BBIgeszGFcxB24kPVy1Y1RS+9A\n3iFOXgXdgZwgDnqf0olek7ndOY8Ocm0VsLEjnUKRSmLF3s6wyyy71gILUyiq9hSffd8I50BGO5Aq\nWBNJBsOeKmBjS4lEGoXvAQZQMSywxQy7lA8O7XNFoux7sAms66zQ2Kl22IH0baSViSKDf7NUwMaa\noJlhYIHZHLBtZtmjcHzXPlkG5uBO1Vpwdec+bt+qdSgh0v5fE4l2ZFUuhhNSKJyAzc5b93GWXcqp\nFIrbBWumRS106pIK4I+oalcZPRk0ulXAxhonYEWi+keAhbIN4E9xiPgA/jRQtcIFvgayUwC/Uw6Y\nitlEoAKmXByBFRZ2ocAG8GfYpey3ldwMyKNFK1xgXcjbVDrUQPoAfrgLqUwEGfw7pQKmKEp/DJhG\nMUxUwMYGnz7RYRKRzwFzFphPYp3yOWBl+1pzaZ4mQQys0akGUncgJ5YMfu0qYGOLdyFnEoNsAWbZ\no8xB7EJOA3PQpGITV7Gu5FF9MQjgp3cgM+hPKMMlg1+5CthY4q2xDp1YIUpgLdCKEliPFq2A+SB+\n5xY6ugM50WgahTJ8vOWVaqMzD8zbvjk2gH9oBQxgGraWpt2MbmumNbknlYGvO5ATj7qQysWRmsQ9\nD/MV2/3bx78KHGOc1m2xECWvAjQOOnVhDXcgM/hrVoZLBr9yFbCxJOhCEQTwZ+ZsF9XZIIXioAyt\n4qXI+vIu5Ga90qWJoTKRZPDrVwEbC7zbSHA7Q6ILxbwVLvvKrt2BxIrXbnmWLRZock+UB8at6cD6\n0ilECppGoVwUwSAPPyYtELCycx8BDspTbDHvrK9qFANrH2KrAfyJR11IRVFySwb/dqmAjR1hG51i\nIoVixllgUxxSdDuQh5Qj99HmgdkYWH81kMpEkcE0ir4mcyt5wIsWtLmQTsDKHLrjgAItWhTZZcYJ\nmE1ibTYqNBu9aiD9/Qz+SVaGx3GfRwdE5OdF5P+46dv/tsPr10VkU0RuuuPj/SxJLbCxI2yjE+aA\nmYQFBtCiwCFTbDsBa3KPzb6HQLxAdyAVYGADXEQKwH8B/h522tD/FpFnjDHfTZ36vDHmkdNcWwVs\nLAmsMWeBXZq3yatAHMBnKnAhK3H7HOhQQqQ7kBPP4F/9A8CrxpjXAETki8A/BNIClp4N2RN1IceC\ndAF3yn2ctk0MpzhgioNIyA4ps8ts7D769jkdh3ioyzjxDO5Cdpq8/dbUOQZ4t4i8LCLPisg7+1mS\nWmC5p5R67EWslEihmJo+jIQLrPt4gE2h8BZYc6cSW2CJ+JcG8JUz0c+U7ReBq8aYXRF5CHgaeHuv\nN6mAjR0dLLB5mC3EyauFYAdyj1luu/hX1P8e+gjgK4rnhjtOJD15+yrWCoswxmwF958TkU+LyLIx\nZqPbhVXAxoZS6v5M5D4yj4t22VTqAsccU+A42IG07qO4uBepAP7uBX0GJZ9cd4enbTL3nwH3icg1\n4E3gnwGJWZEisgqsG2OMiDwASC/xAhWwMSFdSjRDog5ymih1wnNImRaFRA4YTZIWGLvEFphaXspg\n378x5lhEHgW+gp3A8BljzHdF5MPu9SeB9wMfEZFj7A/vA/1cWwVMUZQ+GTwOaox5Dngu9dyTwf3H\ngcdPe10VsLHDW2OSDOKH/b+wQXy/A9mkQvMgSF4FbWKodCB7vwEVsLGhFNwmc8BiF9L+BW1R5IBy\nIgs/0T4HdAdS6cDeqBfQhgpYrgm7r3rac8CYPkrUP/oM/D1sG50tFqA+ncy+3wfdgVSSZO+7VwEb\nG/xX2Z4Ddmk66T4C7DKb2oHECpi3wDDoDqSSJHuWuArYWJAu5A76gE3D1LR1H72IHTBFkVZkfUXx\nr6j/PSTdR7W8FMjib0AFLPeEwgXRJKJAwMouC79FAbAxsLAGcrNeia2vSMC2SAbxFUUtMGUo+PIh\naOuFPw1T5cPE2YeuItILGPXp2PqKBGyXOPalKJDFP2QqYGNHKoWiQxLrAWX23A7kFgux9RUN8IDY\nhdQAvuLRXUhFUXJL9qxxFbDckh5gG7qQJCwwb30duxjYIVNxAH8n3X3VW1naxFBJk73fgwrYWBDm\ngs3YuykBK3BMy51zgE+jmGe7uZAK4Hs3IWzwlL0frjIK1AJThkZggQXiZQXM/vAOmXLnTCUTWNv6\n3/tb3YFUQrL3W1AByz3BBCIgnUJBEYpugEfsQsaTiBLW1zbEiau6A6mkyd7vQQUs16SSV/39hAV2\nlGhgaG+n7C7k5nyHDPy0BaZupOLJ3m9ABSz3pIP4aQssWf8IcRnRflv8yw/w8PcVJUTTKJSh4b9K\niYP4RbgUCNiBs8Di+Fcp5T76AbaglpfSTvZ+CypguaRTCkXQUjoM4BfDFApvgc2wdZBKYE30vwcN\n4CvtZC8G1nWsmohcFZGvish3ROTbIvLLJ5z3KTdx92URuX84S1UUZbQc9Xm002sytzvn1DrSay7k\nEfArxpi/BjwI/GsReUfqH30YeJsx5j7gXwFP9PMPjx+vjejfLaYOIveRIhSKcRPDXWbZZdb2AUvH\nv44hLh/q1MTwtWF/kBHz2qgXkAMGGwwZTOb+eeCdwC+cl450FTBjzC1jzEvu/jZ2ku6V1GmPAE+5\nc14AKm7CyITx2gj/7RIwa49UDlix6Hcgp6JjiwXuNueSLmTU/yv8SxrGwV670E908bw26gXkgIEt\nsGgytzHmCPCTuUMG0pG+J3O7kUj3Ay+kXuo0dffefq+rnIUwB8wd3vryt8RDbL0Floh/RQH8dPuc\n7MU7lFEz8GjufiZzD6QjfQXxRWQe+BLwUWeJtZ2SetzPJF7lTIQ5YEEQPxQv9+0eU+CQMrvMArC7\nPXtCAD8tXBrEV0IGTqPoVw9OrSM9BUxESsCXgd82xjzd4ZT01N173XMdaBt4OWY8P+oF2MTU33MH\nVpu+d24Xz8DnGyrj/vnOyn8Y9I09J3N3OKeLjsR0FTAREeAzwCvGmE+ecNozwKPAF0XkQaBpjFlL\nn2SMSauroig54Yz///aczE2fOpKmlwX2HuCDwDdF5KZ77mPAj4IdTGmMeVZEHhaRV4Ed4EN9fSRF\nUSaCfiZzD6ojYoyGqxRFySd970L2w7CS1bJEr88oItdFZFNEbrrj46NY5yCIyGdFZE1EvtXlnNx+\nf70+X56/u4nFGHMuB9Y0fBW4ht0Sewl4R+qch4Fn3f13AX9yXv/+RRx9fsbrwDOjXuuAn+9vY1Nl\nvnXC63n//np9vtx+d5N6nKcFNrRktQzRz2eE9u3gXGCM+Rpwu8spuf7++vh8kNPvblI5TwEbWrJa\nhujnMxrg3c7FelZE3nlhqxs+ef/+ejHO391Ycp7dKIaWrJYh+lnri8BVY8yuiDwEPA28fbjLulDy\n/P31Yty/u7HjPC2woSWrZYien9EYs2WM2XX3nwNKIrJ8cUscKnn//roy5t/dWHKeAhYlq4nIFDZZ\n7ZnUOc8AvwhwmmS1DNHzM4rIqksARkQewKaqbFz8UodC3r+/roz5dzeWnJsLaYaYrJYV+vmMwPuB\nj4jIMXZCxgdGtuBTIiJfAN4H1ETk+8AncEWW4/D99fp85Pi7m1Q0kVVRlNxyromsiqIoF4kKmKIo\nuUUFTFGU3KICpihKblEBUxQlt6iAKYqSW1TAFEXJLSpgiqLklv8PDrsHxYLFZLEAAAAASUVORK5C\nYII=\n",
"text": [
"<matplotlib.figure.Figure at 0x1095b8590>"
]
}
],
"prompt_number": 6
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There is a clear degeneracy here - why is it in the direction you would expect?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A few more exercises for the enthusiastic:\n",
" - Do the derivation of the analytic form of the likelihood\n",
" - Find the maximum likelihood for the three parameters"
]
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment