fermi_likelihood_lecture-solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likelihood lecture: Solutions to exercises\n",
    "============"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import seaborn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1\n",
    "\n",
    "This is one example of what the Mock data could be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFJCAYAAABDx/6zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG5hJREFUeJzt3X9sVfX9x/HX7b0D1tteqXJlkk7gaklWyLLN7tYyrZrv\nWDVMkxk2hQVdYAs/mmgddlSkLQIBWbdOZtYv3WAxaa1IkC1Es+S71WnVkhviV7fZ4pbxBRJKaUqt\nob1CuW3P9w/WIvTHvZfeez/33vN8/MU99/ae933z48Xn3M/5fByWZVkCAAAJlWG6AAAA7IgABgDA\nAAIYAAADCGAAAAwggAEAMIAABgDAAFciTtLd3RfV63NyMtXb+1mcqrEnehof9DX26Gl80NfYC9dT\nrzd70p9PyhGwy+U0XULaoafxQV9jj57GB32Nvan2NCkDGACAdEcAAwBgAAEMAIABBDAAAAYQwAAA\nGEAAAwBgAAEMAIABBDAAAAYQwAAAGEAAAwBgQNi1oEOhkCoqKtTR0aGMjAxt27ZNLpdLFRUVcjgc\nysvLU3V1tTIyyHIAACIVNoDffvttDQ4Oav/+/Xrvvff0wgsvKBQKqaysTIWFhaqqqlJzc7OWLFmS\niHoBAIi58rpWSVLN+sUJO2fYYev8+fM1NDSk4eFh9ff3y+Vyqa2tTX6/X5JUXFys1tbWuBcKAEA6\nCTsCzszMVEdHhx544AH19vZqz549Onr0qBwOhyTJ7Xarr2/y7QZzcjKj3jUi3DZOiB49jQ/6Gnv0\nND7o68SczsuZFm2PptLTsAH80ksv6a677tKGDRvU2dmpxx9/XKFQaPT5YDAoj8cz6XtEuwel15sd\n9R7CmBw9jQ/6Gnv0ND7o6+SGhixJ0e1fH66nU94P2OPxKDv78pvccMMNGhwcVH5+vgKBgCSppaVF\nBQUFERcMAAAiGAH/6Ec/0qZNm7RixQqFQiE99dRTWrRokSorK1VbWyufz6eSkpJE1AoAQNoIG8Bu\nt1u7d+8ec7yxsTEuBQEAEEsmZjhHgpt3AQAwgAAGAMAAAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgA\nAAMIYAAADCCAAQAwgAAGAMAAAhgAAAMIYACArQXau/Rp/4B6zl9U1b6AAu1dCTlv2M0YAABIV4H2\nLtUfbht9fLo7OPq4MH92XM/NCBgAYFtvHDk5wfFTcT83AQwAsK0z5z4b93hnTzDu5yaAAQC2NWdW\n5rjHb7nJHfdzE8AAANtaWjRvguNz435uAhgAkLbCzXAuzJ+tNQ8tlDPDIUnK9WZpzUML4z4BS2IW\nNAAgTUU6w7kwf7YOvnVckrR1tT9h9TECBgCkJZMznCNBAAMA0pLJGc6RIIABAGnJ5AznSBDAAIC0\nZHKGcySYhAUASEsjE632vt6uoWFLud4sLS2am5AZzpEggAEAacvUDOdIcAkaAAADCGAAAAwggAEA\nMIAABgDAAAIYAAADCGAAAAwggAEAMIAABgDAAAIYAAADwq6EdejQIf3hD3+QJA0MDOjYsWNqamrS\njh075HA4lJeXp+rqamVkkOUAAEQqbGo+/PDDamhoUENDgxYuXKjNmzfrN7/5jcrKytTU1CTLstTc\n3JyIWgEAiIua9YtVs35xQs8Z8bD1H//4h/7973/rkUceUVtbm/z+y2tqFhcXq7W1NW4FAgCQjiLe\njKG+vl6lpaWSJMuy5HA4JElut1t9fX2T/mxOTqZcLmdUhXm92VG9HuHR0/igr7FHT+PDrn11Oi/n\nVTw+/1TeM6IAPn/+vE6cOKE777xTkq76vjcYDMrj8Uz68729n0VVlNebre7uyUMd0aGn8UFfY4+e\nxoed+zo0ZElSzD9/uJ6GC+eILkEfPXpURUVFo4/z8/MVCAQkSS0tLSooKIjkbQAAwH9ENAI+ceKE\ncnNzRx9v3LhRlZWVqq2tlc/nU0lJSdwKBABgKhI9uSpSEQXwj3/846sez58/X42NjXEpCAAAO+Dm\nXQDAuMrrWlVex10u8UIAAwBgAAEMAIABBDAAAAYQwAAAGEAAAwBgAAEMAIABBDAAAAYQwAAAGEAA\nAwBgAAEMAIABBDAAIKZYwjIyBDAAxAChg2gRwAAAGEAAAwBgAAEMAIABBDAAAAYQwAAAGEAAAwDG\nCLR36dP+AfWcv6iqfQEF2rtMl5R2XKYLAAAkl0B7l+oPt40+Pt0dHH1cmD/bVFlphxEwAOAqbxw5\nOcHxUwmtI90RwACAq5w599m4xzt7ggmuJL0RwACAq8yZlTnu8Vtucie4kvRGAAMArrK0aN4Ex+cm\ntpA0xyQsAMBVRiZa7X29XUPDlnK9WVpaNJcJWDFGAAMAxijMn62Dbx2XJG1d7TdcTXriEjQAAAYQ\nwAAAGEAAA0ACsW8wRhDAADBFLNt4Bb2IHJOwAGAKWLbxCnoRHUbAADAFLNt4Bb2IDgEMAFPAso1X\n0IvoRHQJur6+Xm+++aZCoZCWL18uv9+viooKORwO5eXlqbq6WhkZZDkA+5kzK1Onu8cGjB2XbaQX\n0QmbmoFAQB988IFeeeUVNTQ06OzZs9q5c6fKysrU1NQky7LU3NyciFoBIOmwbOMV9CI6YQP43Xff\n1YIFC1RaWqq1a9fq3nvvVVtbm/z+yyujFBcXq7WVKfUA7Kkwf7bWPLRQzgyHJCnXm6U1Dy205aQj\nehGdsJege3t7debMGe3Zs0enT5/WunXrZFmWHI7LDXa73err65v0PXJyMuVyOaMqzOvNjur1CI+e\nxgd9jb1U6+l378nWH975P0nSf1f816SvdTov/9tp4jNGe87rqTWaXqSDqfw+hg3gmTNnyufzadq0\nafL5fJo+fbrOnj07+nwwGJTH45n0PXp7x/9ifiJeb7a6uycPdUSHnsYHfY29VO3p0JAlSWFrj/R1\nsXY9fb3eWk19xkQL19Nw4Rz2EvQdd9yhd955R5ZlqaurSxcuXFBRUZECgYAkqaWlRQUFBVGWDQCA\nvYUdAd933306evSoli1bJsuyVFVVpdzcXFVWVqq2tlY+n08lJSWJqBUAgLQR0W1IP/vZz8Yca2xs\njHkxAADYBTfvAgBgAGtBAwDGVbN+sekS0hojYAAADCCAAQAwgAAGAMAAAhgAEoTN6vF5TMICgARg\ns3pcixEwgKSyevv/qLwu/TZ4YbN6XIsABoAEYLN6XIsABoAEmDMrc9zjbFZvXwQwACQAm9XjWkzC\nAoAEGJlotff1dg0NW8r1Zmlp0dy0nIDFClqRIYABIAYiCZ3C/Nk6+NZxSdLW1f54l4QkxyVoIMHK\n61rTcpYvgOgQwAAAGEAAAwBgAAEMAIABBDAAAAYQwAAAGEAAAwBgAAEMAIABBDAAAAYQwAAAGEAA\nAwBgAAEMIGkE2rvUc/6ies5fVNW+gALtXaZLAuKGzRgAJIVAe5fqD7eNPj7dHRx9nI47BgGMgAEk\nhTeOnJzg+KmE1gEkCgEMICmcOffZuMc7e4IJrgRIDAIYQFKYMytz3OO33OROcCVAYvAdMICksLRo\n3lXfAV85PtdANfFTs36x6RKQJBgBA0gKhfmzlfXFL4w+zvVmac1DC5mAhbRFAAMJFGjv0qf9A9xm\nM4HpX3AqI8OhmzwztHW1n/BFWuMSNJAg3GYD4PMYAQMJkgq32ZTXtaq8rtV0GYAtRDQC/t73vqes\nrCxJUm5urtauXauKigo5HA7l5eWpurpaGRlkOTAZbrMB8HlhA3hgYECWZamhoWH02Nq1a1VWVqbC\nwkJVVVWpublZS5YsiWuhQKqbMytTp7vHhi232QD2FHbY+vHHH+vChQtatWqVHnvsMX344Ydqa2uT\n3++XJBUXF6u1lUtWQDhLi+ZNcDy9brMBEJmwI+AZM2Zo9erV+v73v6+TJ0/qJz/5iSzLksPhkCS5\n3W719fVN+h45OZlyuZxRFeb1Zkf1eoRHT+Mj0r5+955seTwz9KtX/leDQ5bm3eLR9/8rT8Vfz41z\nhZFzOi//vTb1Z2Xk/E6ngz+vcUBPY28qPQ0bwPPnz9fcuXPlcDg0f/58zZw5U21tV2ZyBoNBeTye\nSd+jt3f8774m4vVmq7t78lBHdOhpfETb16/k3qAb3NMlSVWPF0hSUv2+DA1ZkszVNHL+oSErqfqS\nDvg3IPbC9TRcOIe9BH3w4EE9//zzkqSuri719/frW9/6lgKBgCSppaVFBQUF0dQMAIDthR0BL1u2\nTM8884yWL18uh8OhHTt2KCcnR5WVlaqtrZXP51NJSUkiagUAIG2EDeBp06bpl7/85ZjjjY2NcSkI\nAAA74OZdAAAMIIABJIWRdbKHhy192j/AOtlIe6wFnSZGlg9kqzOkomvXyR4atlgnG2mPETAA41Jh\nnWwg1ghgAMaxTjbsiAAGYNycWZnjHmedbKQzAhiAcayTDTtiEhYA40YmWu19vV1Dw5ZyvVlaWjSX\nCVhIawQwgKRQmD9bB986LqfToa2r/abLAeKOS9AAJF25D7fn/EVV7QtwHy4QZ4yAAYy5D/d0d5D7\ncIE4YwQMgPtwAQMIYADchwsYwCVoIMGScbnQObMydbp7bNhyHy4QP4yAAXAfLmAAI2AkJTaXSCzu\nwwUSjwBOYoQQEmnkPlxJ3IcLJACXoAEAMIAABgDAAAIYAAADCOA0wBKCAJB6mISV4lhCEABSEyPg\nFMcSgrFVXtc6OvscAOKJAE5xLCEIAKmJAE5xc2ZljnucJQQBILkRwCmOJQSRTmrWL9a+zd8xXQaQ\nEEzCSnEsIQgAqYkATgMsIQgAqYdL0AAAGEAAAwBgAAEMAIABBDAAAAYQwAAAGMAs6CQ1ssHC0LCl\nqn0BLS2aZ5tbi+z82U2rWb/YdAmAbUQ0Au7p6dE999yj48eP69SpU1q+fLlWrFih6upqDQ8Px7tG\n2xnZYGFo2JJ0ZYMFO+xyZOfPDsBewgZwKBRSVVWVZsyYIUnauXOnysrK1NTUJMuy1NzcHPci7cbO\nGyzY+bMDsJewAbxr1y49+uijuvnmmyVJbW1t8vsvL/ZQXFys1lZ2jok1O2+wYOfPDsBeJv0O+NCh\nQ7rxxht1991367e//a0kybIsORwOSZLb7VZfX1/Yk+TkZMrlckZVmNebHdXr08mtX8rWyc7zY45/\neXb2hH1xOi//nkzWt1To6fV89liKpI/XSoW+php6Gh/0Nfam0tNJA/i1116Tw+HQkSNHdOzYMW3c\nuFGffPLJ6PPBYFAejyfsSXp7xx/VTMTrzVZ3d/hgT1cl3/yy6g+3jXt8or4MDV3+znSi51Olp9fz\n2WMpXB+vlSp9TSX0ND7oa+yF62m4cJ40gF9++eXRX69cuVJbtmxRTU2NAoGACgsL1dLSojvvvDPK\nkhGOnTdYsPNnB2AvUd8HvHHjRr344ot65JFHFAqFVFJSEo+6bK8wf7ZmZk3XTZ4Z2rraH5cAKq9r\nVXld8n2Hn4jPPp6R2596zl9U1b4AM68BxFXE9wE3NDSM/rqxsTEuxQCmjNz+NGLk9idJjL4BxAUr\nYQHi9icAiUcAw6hkuQzO7U8AEo2lKNMESwhOzZxZmTrdPTZsb7nJbaAaAHbACBiQtLRo3gTH5ya2\nEAC2wQgYELc/AUg8Ahj4j8L82Tr41nFJ0tbVfsPVAEh3XIIGAMAAAhgAAAMIYAAADOA74CTW23fR\ndAkAgDhhBAwAgAEEcJpIlhWlAACRIYABADCAAAYAwAACGAAAA5gFncRysmeYLsEYNpcAkO4YAdtU\noL1Ln/YPqOf8RVXtCyjQ3mW6JACwFUbANhRo71L94bbRx6e7g6OP2XwAABKDEbANvXHk5ATHTyW0\nDgCwMwLYhs6c+2zc4509YzekjycugwOwMwLYhubMyhz3+C03uRNWw8hl8KFhS9KVy+CEMAC7IIBt\naGnRvAmOz01YDcl6Gbxm/WJmYANICCZh2dDIRKu9r7draNhSrjdLS4vmJnQCVrJcBgcAU2w5Ambd\n5MshPDNrum7yzNDW1f6Ez35OhsvgAGCSLQMY5iXDZXAAMIlL0DAiGS6DA4BJBDCMKcyfrYNvHZck\nbV3tN1wNACQWl6CTFPfIAkB6YwSchFgqEgDSHyPgJBTtPbKMlgEg9TACTkLR3CPLaBkAUhMj4CQU\nzT2yybqiFABgcgRwEormHllWlAKA1MQl6CQUzT2yc2Zl6nT32LBlRSkASG5hA3hoaEibN2/WiRMn\n5HA49Nxzz2n69OmqqKiQw+FQXl6eqqurlZHBYDqWIr1HdmnRvKu+A75ynBWlACCZhQ3gv/71r5Kk\n/fv3KxAI6Fe/+pUsy1JZWZkKCwtVVVWl5uZmLVmyJO7FYixWlAKA1BR22Prtb39b27ZtkySdOXNG\nHo9HbW1t8vsvj8qKi4vV2mrvjQ1MM72xAgAgehF9B+xyubRx40b9+c9/1q9//Wu99957cjgckiS3\n262+vr5Jfz4nJ1MulzOqwrze7KheHw2n0xH3c8RCNHVG8tprn0uGPiRDDVOVyrUnK3oaH/Q19qbS\n04gnYe3atUtPP/20fvCDH2hgYGD0eDAYlMfjmfRne3vHn6k7Ea83W93dk4f6VAwNWZIU13PEQjR1\nhnvteD1Nhj4kQw1TEe8/q3ZET+ODvsZeuJ6GC+ewl6D/+Mc/qr6+XpL0xS9+UQ6HQ4sWLVIgEJAk\ntbS0qKCgIJqaAQCwvbAj4O985zt65pln9MMf/lCDg4PatGmTbrvtNlVWVqq2tlY+n08lJSWJqBUx\nVrN+sekSAMC2wgZwZmamdu/ePeZ4Y2NjXAqKt5F1k4eGLVXtC2hp0TwmLQEAEs5WC3GwbjIAIFnY\nKoAnWzeZADaDy+AA7MpWy1exbjIAIFnYKoCj2WUIAIB4slUAR7PLEAAA8WSr74BZNxkAkCxsFcBS\n5LsMAQAQT7a6BA0AQLIggAEAMMB2l6BTCffIAkD6YgQMAIABjIDTBKNlAEgtjIABADCAAAYAwICU\nDODyulaV17WaLgMAgOuWkgEMAECqI4ABADCAAAYAwAACGAAAAwhgAAAMIIABADCAAAYAwAACGAAA\nA2y5FjTrJgMATGMEDACAAQQwAAAGEMAAABhAAAMAYAABDACAAQQwAAAGEMAAABhAAAMAYEDKBXCg\nvUuf9g+o5/xFVe0LKNDeZbokAACiNulKWKFQSJs2bVJHR4cuXbqkdevW6fbbb1dFRYUcDofy8vJU\nXV2tjIzE5HigvUv1h9tGH5/uDo4+LsyfnZAaAACIhUmT8/Dhw5o5c6aampq0d+9ebdu2TTt37lRZ\nWZmamppkWZaam5sTVaveOHJyguOnElYDAACxMGkA33///XryySclSZZlyel0qq2tTX6/X5JUXFys\n1tbW+Ff5H2fOfTbu8c6eYMJqAAAgFia9BO12uyVJ/f39euKJJ1RWVqZdu3bJ4XCMPt/X1xf2JDk5\nmXK5nFEV5vVmjzl265eydbLz/JjjX56dPe7rcTV6FB/0NfboaXzQ19ibSk/D7obU2dmp0tJSrVix\nQg8++KBqampGnwsGg/J4PGFP0ts7/sh1Il5vtrq7xwZ7yTe/fNV3wJ8/Pt7rccVEPcXU0NfYo6fx\nQV9jL1xPw4XzpJegz507p1WrVqm8vFzLli2TJOXn5ysQCEiSWlpaVFBQEG3N160wf7bWPLRQzozL\nI/Bcb5bWPLSQCVgAgJQz6Qh4z549On/+vOrq6lRXVydJevbZZ7V9+3bV1tbK5/OppKQkIYWOKMyf\nrYNvHZckbV3tT+i5AQCIlUkDePPmzdq8efOY442NjXErCAAAO0i5hTgAAEgHBDAAAAYQwAAAGEAA\nAwBgAAEMAIABBDAAAAYQwAAAGEAAAwBgAAEMAIABBDAAAAYQwAAAGEAAAwBgAAEMAIABBDAAAAYQ\nwAAAGDDpfsDJqmb9YtMlAAAwJYyAAQAwgAAGAMAAAhgAAAMIYAAADCCAAQAwgAAGAMAAAhgAAAMI\nYAAADCCAAQAwgAAGAMAAAhgAAAMIYAAADCCAAQAwwGFZlmW6CAAA7IYRMAAABhDAAAAYQAADAGAA\nAQwAgAEEMAAABhDAAAAY4DJdwOcNDw9ry5Yt+uc//6lp06Zp+/btmjt3rumyUk4oFNKmTZvU0dGh\nS5cuad26dbr99ttVUVEhh8OhvLw8VVdXKyOD/39Fq6enRw8//LB+//vfy+Vy0dMYqK+v15tvvqlQ\nKKTly5fL7/fT1ykKhUKqqKhQR0eHMjIytG3bNv68TsHf/vY3/eIXv1BDQ4NOnTo1bh8PHDig/fv3\ny+Vyad26dbrvvvvCvm9Sdf8vf/mLLl26pFdffVUbNmzQ888/b7qklHT48GHNnDlTTU1N2rt3r7Zt\n26adO3eqrKxMTU1NsixLzc3NpstMOaFQSFVVVZoxY4Yk0dMYCAQC+uCDD/TKK6+ooaFBZ8+epa8x\n8Pbbb2twcFD79+9XaWmpXnjhBfp6nX73u99p8+bNGhgYkDT+3/vu7m41NDRo//792rdvn2pra3Xp\n0qWw751UAfz+++/r7rvvliR97Wtf00cffWS4otR0//3368knn5QkWZYlp9OptrY2+f1+SVJxcbFa\nW1tNlpiSdu3apUcffVQ333yzJNHTGHj33Xe1YMEClZaWau3atbr33nvpawzMnz9fQ0NDGh4eVn9/\nv1wuF329TrfeeqtefPHF0cfj9fHvf/+7vv71r2vatGnKzs7Wrbfeqo8//jjseydVAPf39ysrK2v0\nsdPp1ODgoMGKUpPb7VZWVpb6+/v1xBNPqKysTJZlyeFwjD7f19dnuMrUcujQId14442j/0GURE9j\noLe3Vx999JF2796t5557Tk8//TR9jYHMzEx1dHTogQceUGVlpVauXElfr1NJSYlcrivf1o7Xx/7+\nfmVnZ4++xu12q7+/P+x7J9V3wFlZWQoGg6OPh4eHr/rgiFxnZ6dKS0u1YsUKPfjgg6qpqRl9LhgM\nyuPxGKwu9bz22mtyOBw6cuSIjh07po0bN+qTTz4ZfZ6eXp+ZM2fK5/Np2rRp8vl8mj59us6ePTv6\nPH29Pi+99JLuuusubdiwQZ2dnXr88ccVCoVGn6ev1+/z35uP9PHa7AoGg1cF8oTvFZcKr9M3vvEN\ntbS0SJI+/PBDLViwwHBFqencuXNatWqVysvLtWzZMklSfn6+AoGAJKmlpUUFBQUmS0w5L7/8shob\nG9XQ0KCvfOUr2rVrl4qLi+npFN1xxx165513ZFmWurq6dOHCBRUVFdHXKfJ4PKMBcMMNN2hwcJB/\nA2JkvD5+9atf1fvvv6+BgQH19fXp+PHjEeVXUm3GMDIL+l//+pcsy9KOHTt02223mS4r5Wzfvl1/\n+tOf5PP5Ro89++yz2r59u0KhkHw+n7Zv3y6n02mwytS1cuVKbdmyRRkZGaqsrKSnU/Tzn/9cgUBA\nlmXpqaeeUm5uLn2domAwqE2bNqm7u1uhUEiPPfaYFi1aRF+v0+nTp/XTn/5UBw4c0IkTJ8bt44ED\nB/Tqq6/KsiytWbNGJSUlYd83qQIYAAC7SKpL0AAA2AUBDACAAQQwAAAGEMAAABhAAAMAYAABDACA\nAQQwAAAGEMAAABjw/4Xjz4G6YmymAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f9e24e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def make_data(intercept, slope, N=20, dy=5, rseed=42):\n",
    "    rand = numpy.random.RandomState(rseed)\n",
    "    x = 100 * rand.rand(N)\n",
    "    y = intercept + slope * x\n",
    "    y += dy * rand.randn(N)\n",
    "    return x, y, dy * ones_like(x)\n",
    "\n",
    "theta_true = [25, 0.5]\n",
    "x, y, dy = make_data(*theta_true)\n",
    "errorbar(x, y, dy, fmt='o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# x, y, dy are all assumed to be arrays\n",
    "# theta is a list of model parameters, [intercept, slope]\n",
    "\n",
    "def model(theta,x):\n",
    "    \"\"\"\n",
    "    Computes model value given parameters and data.\n",
    "    Returns ym array.\n",
    "    \"\"\"\n",
    "    return theta[0] + theta[1]*x\n",
    "\n",
    "def loglike(theta,x,y,dy):\n",
    "    \"\"\"\n",
    "    Computes log-likelihood given data, uncertainties, model and model pars.\n",
    "    \"\"\"\n",
    "    ym=model(theta,x)\n",
    "    like=-0.5*(log(2.*pi*dy**2)+(y-ym)**2/dy**2)\n",
    "    return like.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "theta=[20,-0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ym=model(theta,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAFJCAYAAABDx/6zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1AVHXe9/HPAIIColiopPlA0pa6XVQEWoniptTdel27\n5rbpXm2b9qCXexVmKpqI3VqmJj1tlpVum2bmbq21u3dbLZmUFKtutkmaralbCGhCCqPiMJz7j5Ki\nMzAMzMyZh/frL+d7hjnf+frwcR7O72czDMMQAADwqwirGwAAIBwRwAAAWIAABgDAAgQwAAAWIIAB\nALAAAQwAgAWi/HGSI0dqPbp/YmKsampO+Kib8MRMfYO5eh8z9Q3m6n3uZpqU1LXVnw/IV8BRUZFW\ntxBymKlvMFfvY6a+wVy9r6MzDcgABgAg1BHAAABYgAAGAMACBDAAABYggAEAsAABDACABQhgAAAs\nQAADAGABAhgAAAsQwAAAWMDtWtAOh0N5eXkqLy9XRESEFi1apKioKOXl5clmsyk1NVUFBQWKiCDL\nAQBoK7cBvGXLFjU0NGjDhg3aunWrHn74YTkcDuXm5iozM1MLFixQUVGRxowZ449+AQDwuvySJZKk\nRZfP9ds53b5sHThwoJxOpxobG1VXV6eoqCiVlZUpIyNDkpSVlaWSkhKfNwoAQChx+wo4NjZW5eXl\nuuaaa1RTU6Mnn3xS27Ztk81mkyTFxcWptrb17QYTE2M93jXC3TZO8Bwz9Q3m6n3M1DeYa8siI77O\nNE9n1JGZug3gZ599VldeeaVmzpypiooK3XTTTXI4HE3H7Xa7EhISWn0MT/egTErq6vEewmgdM/UN\n5up9zNQ3mGvrnI2GJM/2r3c30w7vB5yQkKCuXb9+kG7duqmhoUGDBw9WaWmpJKm4uFjp6eltbhgA\nALThFfCvfvUrzZs3T5MmTZLD4dCMGTM0dOhQ5efnq7CwUCkpKcrJyfFHrwAAhAy3ARwXF6dHHnnE\nVF+3bp1PGgIAwJus+IZzW3DxLgAAFiCAAQCwAAEMAIAFCGAAACxAAAMAYAECGAAACxDAAABYgAAG\nAMACBDAAABYggAEAsAABDACABQhgAEBY2161U8fqj6v6VI3uKy3U9qqdfjmv280YAAAIVdurduq3\nZeubbh+yVzbdTu+V5tNz8woYABC2Xj/wlsv6Gwc3+/zcBDAAIGxVnjjssl5hr/L5uQlgAEDY6h3b\n02U9Oa6Xz89NAAMAwlbOgNEu62P7Z/v83AQwACBkufuGc3qvNN08ZJIibZGSpD7xybp5yCSffwFL\n4lvQAIAQ1dZvOKf3StMr+16TJM3LmOG3/ngFDAAISVZ+w7ktCGAAQEiy8hvObUEAAwBCUlu+4exs\nbNQTm3apsjRNDSdj/NWaJD4DBgCEqJwBo5t9BnzG2P7ZcjQ06vE/fqR/7jv6TTVGzvpov/ZHAAMA\nQtKZL1o99/GLchpO9YlP1ug+I/W3vzVq5edvN7vvWUP3KKZ7rV/74y1oAEDISu+Vpm4xCeoedZYc\nZVfqqbVfae/nXzUdnz3xYq3JG63OPY75vTdeAQMAQtbhr06qvDjzm1vfvsK958ZLdV6fbtY09Q0C\nGAAQcr44UqcFq/9uqi+8+TL169XVgo7MCGAAQMj47NBxLX5uu6ne48K9evC/plrQUcsIYABA0Nt9\nsEbLX/jAVJ/58zRtqFhlQUfuEcAAgKD1wadH9NhLH5nq8268VIPOfMZb4eem2ogABgAEnZJdFXrm\nz7tN9UD6jNcdAhgAEDSKdnyh59/ca6ovuW2YevWItaCj9iOAAQAB708lB/TH4s9M9Qf/53L1SOhs\nQUcdRwADAALWhqJP9ca2z031h++4Ugmx/l060tvcBvDLL7+sP/7xj5Kk+vp67d69W+vXr9f9998v\nm82m1NRUFRQUKCKCRbUAAN6x+i8fa+tHlab6b3KzFNs5NF47un0W48eP1/jx4yVJ9957r6677jo9\n/vjjys3NVWZmphYsWKCioiKNGTPG580CAELbI7//UB82bZDwrSdmjlRMp0ifnXfR5XN99tgtafN/\nIz766CP961//UkFBgX7zm98oIyNDkpSVlaWtW7cSwACAdlv83HZ9dui4qf7UrFGKigzNd1jbHMCr\nVq3S9OnTJUmGYchms0mS4uLiVFvb+g4SiYmxiory7H8uSUnB8TXyYMJMfYO5eh8z9Y1AnOst972p\nquoTpvqm5f+pyAibV85x5nF88fw78phtCuDjx49r//79GjZsmCQ1+7zXbrcrISGh1Z+vqTEPtzVJ\nSV115Ih/t4UKdczUN5ir9zFT3wi0uU5bsUX1DqepvnpOtmw2m6qP1nntXM5GQ5K8/vzdzdRdOLcp\ngLdt26bhw4c33R48eLBKS0uVmZmp4uLipmAGAKAlhmFoytLNpnpkhE1Pz862oCNrtSmA9+/fr759\n+zbdnjNnjvLz81VYWKiUlBTl5OT4rEEAQHBrNAzd4iJ4E7vGaMX0K3x+fiu+YNUWbQrgW265pdnt\ngQMHat26dT5pCAAQGhqcjbpt+dumer9e8Vp4c4b/GwowoXExFQDA6/JLlkjy/BXkaYdTU1dsMdWH\nDuyhu36e5pXeQgEBDADwitoTp3Xno++a6sOG9NJt44ZY0FFgI4ABAB1SffyU7l5ZYqqPvqSP/nvs\nDyzoKDgQwACAdqk4atc9T5ea6j84t7vm/OISCzoKLgQwAMAj+yuOa9HvtpvqGRf21NT/GmpBR8GJ\nAAYAtMnuA9VavmGnqX5Vel9Nuup8CzoKbgQwAKBV2/cc1spNu0z1n2alaNzlA/zfUIgggAEALtkr\nkvTVpylaqebh+8ucH2jUxX0s6ip0EMAAgGb+3/sH9Ye390lKaVaf9pOhuuyCntY0FYIIYACAJOnF\ntz7V63//3FSf+fM0DRnYw4KOQhsBDABh7uk/lem9sipTPSltl6IT7BoycLQFXYU+AhgAwtSKDR+o\n7ECNqb7olkz1OTtO+SXma3zbor1LWIYbAhgAvCCYQmf+M6U69KXdVF82bbjO7tbFgo7CEwEMAGFi\n+kNbdLLeaao/fMeVSoiNtqCj8EYAA0CIm/zAWy7rj8/IUpcYYsAqTB4AQlRLwbvq7pHqFBXp527w\nfQQwAISYloL36dmjFBkR4edu0BICGABCREvBu3pOtnr2TNCRI7V+7gitIYABIIgZhqEpSze7PLYm\nr/3X726v2qlj9cflNJy6r7RQOQNGK71XWrsfD2YEMAAEocZGQ7cs837wSl+H72/L1jfdPmSvbLpN\nCHsPAQwAQcTR0KjbH3zb5bGOBu8Zrx9w/Vb2Gwc3E8BeRAADQBA4Wd+g6Q8VuzzmreA9o/LEYZf1\nCrt5uUq0HwEMAAGs9sRp3fnou6Z6dKcIPTlzlE/O2Tu2pw7ZK0315LhePjlfuCKAASAAVR8/pbtX\nlpjqvXrEasltw3x67pwBo5t9BnzG2P7ZPj1vuCGAASCAVBy1656nzZsgXNCvu2ZPusQvPZz5nPe5\nj1+U03CqT3yyxvbP5vNfLyOAASAA7K84rkW/226qZw7updv/c4jf+0nvlaZX9r0mSZqXMcPv5w8H\nBDAAWOjjA9V6cMNOU31M+rmaeFWqBR3BXwhgALDA9j2HtXLTLlN9fFaKfnz5AP83BL8jgAHAj+56\n6Rl99WmKqf7Lq3+gUWl9LOgIViGAAaCD2rJs41/eO6CXtnwmqXn4TvvJUF12QU//NetjLGHZdgQw\nAHSAu2UbNxR9qje2fW76uZk/T9OQgT381qc/sISlZwhgAOiAlpZtfOGfr2nljmpTPSltl6IT7Boy\n0LurVwUClrD0DAEMAB3Q0rKNJ21fNbu9+JZMnXN2nPJLzNf4hgqWsPRMmwJ41apVeuutt+RwODRx\n4kRlZGQoLy9PNptNqampKigoUASbPAMIQy0t22icipckLZs2XGd36+LvtizBEpaecZuapaWl+uCD\nD/TCCy9o7dq1qqys1JIlS5Sbm6v169fLMAwVFRX5o1cACDgHPnQdLpMuukZr8kaHTfhKXy9h6QpL\nWLrm9hXwu+++q/PPP1/Tp09XXV2dZs+erY0bNyojI0OSlJWVpa1bt2rMmDE+bxYAAsXkB8583pms\n0/+Sos75TLbOdTonvpeuHhie3/xlCUvPuA3gmpoaHTp0SE8++aS++OILTZs2TYZhyGazSZLi4uJU\nW1vb6mMkJsYqKirSo8aSkrp6dH+4x0x9g7l6XyDPdNzMV0w1Z3WykoeUKyKykx4ZV9Dqz0dGfP1v\npxXP0dNztqfXa5JG6M/7/ypJeujaBR6dLxh15PfRbQB3795dKSkpio6OVkpKimJiYlRZ+e17/Ha7\nXQkJCa0+Rk3NCY+aSkrqqiNHWg91eIaZ+gZz9b5Anem3r3ibe2rWKEVFRii/pFTORrnt3dloSHJ/\nP29rz1zb26tVz9Hf3M3UXTi7DeBLL71Uzz33nG6++WYdPnxYJ0+e1PDhw1VaWqrMzEwVFxdr2DDf\nbo0FAFZpKXhXz8lueicQaA+3AZydna1t27ZpwoQJMgxDCxYsUN++fZWfn6/CwkKlpKQoJyfHH70C\ngF8YhqEpSze7PLYmL/Su34U12nQZ0uzZs021devWeb0ZALCSs7FRty572+UxghfexkIcAMJevcOp\naSu2uDxG8MJXCGAAYct+yqH/ffgdl8cIXmnR5XOtbiGkEcAAwk5Nbb1mPr7V5TGCF/5CAAMIGxVH\n7brnaddrMRO88DcCGEDI++zQcS1+brupHhlh09OzWSYR1iCAAYSsXfuPqvDFD031nold9MDtw/3e\nD5vV47sIYAAh5/2ySj31p49N9Qv7J2rWxIst6IjN6mFGAAMIKNP/dI+cjUa7voH75rbP9ULRp6b6\nFUN7a8qPB3ujvXZjs3p8HwEMIOi9tGWf/vLeQVP92uH9dd3I8yzoyIzN6vF9BDCAoLX6Lx9r60fm\nDeAn/ihVYy4714KOWsZm9fg+AhhA0Fn+wgfafbDGVL9t3GANG9Lbgo7cyxkwutlnwGewWX34IoAB\nBI28Ve/pcM1JU/2u6/9DQ1POsqCjtgunzepZQattCGAAAe/WZZub9pj9rvm/TFfKOa3vR+4vbQmd\n9F5pemXfa5KkeRkzfN0SAhwBDPhZfskSSbxKaIuW9uK979ZMJZ8V5+duAO8igAEEnPLiTE0uNofv\niulXKLFrjAUdAd5HAAMIGF+/4s0w1R/LHaG4zp383xDgQwQwAMu19FbzEzNHKqZTpJ+7AfyDAAZg\nmZaC95wRpVp8BZ+RI7QRwAD8rqXgXT0nW/eWLpWz0c8NARYggAH4hWEYmrJ0s8tj7MWLcEQAA/Ap\nZ2Ojbl32tstjBC/CGQEMwCccDU7d/uAWl8cIXoAABuBlJ+sbNP2hYpfHCF7gWwQwAK84fuK0ch99\n11TvHB2plXeNtKAjILARwAA65OixU5r1RImpnnxWrO67dZhHj7W9aqeqTx6T03DqvtJC5QwYHZKb\nFQASAQygnQ59adf8Z0pN9Qv7J2rWxIs9frztVTubbdd3yF7ZdJsQRigigAF45LNDx7X4ue2m+rAh\nvXTbuCHtftzXD7i+NviNg5sJYIQkAhhAm5QdqNaKDTtN9bGXnasbfpTa4cevPHHYZb3CXtXhxwYC\nEQEMoFXb9hzWE5t2merXjUzRtcMHeO08vWN76pC90lRPjuvltXMAgYQABuDS2zvL9dxfPzHVf3XN\nBcr6j3O8fr6cAaObfQZ8xtj+2V4/l5XYBxpnEMAAmvnLewf00pbPTPX/+clQpV/Q02fnTe+Vpo17\nX5HdYZck9YlP1tj+2Xz+i5BFAAN+tL1qp47VHw/Iy2yWrf+H9vz7K1P97hvSNHhAD7/0EBMZrZMN\nJ9U9ppvmZczwyzkBqxDAgJ8E6mU2s58o0ZfHTpnq+Tela2ByggUdAeGBAAb8JNAus5nywFsyXNST\n0nZp6dV3+L0fINy0KYB/+tOfKj4+XpLUt29fTZ06VXl5ebLZbEpNTVVBQYEiIiJ82igQ7ALlMpuW\n9uK979ZMrfzkUb/2AoQztwFcX18vwzC0du3aptrUqVOVm5urzMxMLViwQEVFRRozZoxPGwWCndWX\n2bQUvCumX6HErjF+6QHAt9wG8J49e3Ty5ElNnjxZDQ0Nuuuuu1RWVqaMjAxJUlZWlrZu3UoAA25Y\ndZlNS8H7WO4IxXXu5NNzA2iZ2wDu3LmzpkyZop/97Gc6cOCAbr31VhmGIZvNJkmKi4tTbW1tq4+R\nmBirqKhIjxpLSurq0f3hHjP1jbbO9ZqkEUpI6KzfvP87OQ2n+nfro58MztEV/S7zSV/jZr7isv77\nJdeqc7Trv/qREV//vbbqz8qZ80dG2Pjz6gPM1Ps6MlO3ATxw4ED1799fNptNAwcOVPfu3VVWVtZ0\n3G63KyGh9W9K1tSc8KippKSuOnKk9VCHZ5ipb3g61/O7XKBuMV//fZl96Z2S5PXfl5Ze8T49e5Qi\nIyJUe+ykWjqjs9HwSU9tdeb8zkaDP69exr8B3udupu7C2W0A/+EPf9DevXu1cOFCVVVVqa6uTldc\ncYVKS0uVmZmp4uJiDRvm2ZZjALyvpeBdPSe76R0rAIHDbQBPmDBBc+fO1cSJE2Wz2XT//fcrMTFR\n+fn5KiwsVEpKinJycvzRK4DvaTQM3bJ0s8tja/JG+7kbAJ5wG8DR0dFasWKFqb5u3TqfNATAvQZn\no25b/rbLYwQvEBxYiAMIIvWnnZpWuMXlMYIXCC4EMBAE6k46dMcj77g8FirBe2ad7EajUcfqj2t7\n1c6AWScb8AUCOETklyyRxFZnoab6+CndvbLE5bFQCV7JvE6203AGxDrZgC8RwEAAOvSlXfOfKXV5\nLJSC94xAWycb8AcCGAgg+8qP6b61O0z1qEibnpoVWhvTf1egrJMN+BMBDASAf+47qod//6Gp3qtH\nrJbcFvrX2Vu9TjZgBQIYsNB7uyr19J8/NtWHDEjUzBsutqAja1i1TjZgJQIYsEDtF71drlx15Q+T\nNfnaCy3oyFpnPud97uMX5TSc6hOfrLH9s/n8FyGNAAb86Peb/6Xy0kxT/drh/XXdyPMs6ChwpPdK\n0yv7XlNkhE3zMmZY3Q7gcwQw4AdP/alM75eZv1A06apUXZV+rgUdmZ25DtdpOHVfaaFyBozmFSjg\nQwQw4EMPPP8P7f38K1M98YJPteInt1vQkWvfvw73kL2S63ABHyOAAR+YtbJER4+fMtVn/jxNGypW\nWdBR67gOF/A/Ahjwopa2BMy/KV0Dk7/ZN7vCjw21EdfhAv5HAANe0FLw3n/bMPXuEdusFojLhXId\nLuB/BDDQAeNmvuKyvmL6FUrsGuPnbtqP63AB/yOAEZACfXOJll7xPpY7QnGdO/m5m47jOlzA/wjg\nABboIRSOWgreJ2eOVHSnSD93411nrsOVxHW4gB8QwEAbtBS8m5aNU3W13c/dAAgFBDDQipaCd/Wc\nbNlsNkVGRvi5IwChggAGvqfRMHTL0s0uj4XiXrwArEEAhwCWEPSOBmejblv+tstjBC8AbyOAgxxL\nCHZcvcOpaSu2uDxG8ALwFQI4yLGEYPudOOXQrx9+x+UxgheArxHAQY4lBD1XU1uvmY9vdXmsT1bp\nN78igAH4FgEc5FhCsO2qqk9o7lPvm+qdoiK06u5RkqT8klLTcQDwBQI4yLGEoHsHK2t177PbTPXe\nPWJ1/23DLOgILVl0+VwlJXXVkSO1VrcC+BwBHORYQrBluw/WaPkLH5jqF/ZP1KyJF1vQEQB8iwAO\nASwh2NyOTw7r8T/uMtWHD+mtW8cNtqAjADAjgBEytuws1+/++ompfk1mP/0se5AFHQFAywhgBL0/\nlxzQy8WfmerXZw/S1Zn9LOgIANwjgBG01r+5V3/b8YWpPuXaC3XFD5Mt6AgA2o4ARtBZuWmXtu8x\nX/98x3UXKS31bAs6AgDPEcAIGvev3aF/lR8z1ef+9yVK7dvdgo4AoP0I4AAVzhssfP+5H9nbR8fL\nk0z3u3dyhs7tGW9Bh6Fr0eVzrW4BCBttCuCjR49q/PjxWrNmjaKiopSXlyebzabU1FQVFBQoIoI9\nUb0pnDdYcPXc1adSkSf/Q87qrz/XXTp1uJK6d7GqRQDwCrfJ6XA4tGDBAnXu3FmStGTJEuXm5mr9\n+vUyDENFRUU+bzLctLbBQqhr6blHnfOZHvrfK7UmbzThCyAkuA3gpUuX6oYbblDPnj0lSWVlZcrI\nyJAkZWVlqaSkxLcdhqFw3WBh8gNvqbzO9XOMirOrW1y0nzsCAN9p9S3ol19+WT169NCIESP01FNP\nSZIMw5DNZpMkxcXFqbbW/ZqtiYmxioqK9KixpKSuHt0/lPRNSNa/j5Wb6ucmJLc4l8iIr39PWptb\noM503MxXmn5tnIyTLbbOdJ/Wnrs3tWWO3xeocw1mzNQ3mKv3dWSmrQbwSy+9JJvNpvfee0+7d+/W\nnDlzVF1d3XTcbrcrISHB7Ulqak541FS4L8b+o74j9dtj5g0WRvcd2eJcnI2GJLV4PBBnOvkB89vN\nDYfOU/SgD0311p67N7mb4/cF4lyDHTP1Debqfe5m6i6cWw3g559/vunXN954oxYuXKjly5ertLRU\nmZmZKi4u1rBh7CbjbaG8wYJhGJqy1PVn2avnZMtms2l71ZCQfO4A8F0eX4Y0Z84c5efnq7CwUCkp\nKcrJyfFFX2HPHxss5JcskeSfS08aDUO3tBC8a/JGN7tt1eYS4XzpFwD/a3MAr127tunX69at80kz\nCD3Oxkbduuxtl8e+H7xWCudLvwBYg4U44BOnHU5NXbHF5bFACt4zWrv0iwAG4AsEMLzqxKkG/frh\nYpfHXAWvP98Gb024XvoFwDoEcIiwOsCO2U9rxmPvmupxnaP0WG6WBR15pndsz69X3fqe5LheFnQD\nIBwQwOiQL786qdlPvmeq902K0/+dkmlBR+2TM2B0s8+AzxjbP9uCbgCEAwIY7fLFkTotWP13U33o\nwB666+fB95lpKF/6BSAwEcDwyL7yY7pv7Q5T/Yof9taUawdb0JH3WHX5E4DwRACjTXZ9dlSFG80r\nVF2T2U8/yx5kQUcAENwIYLSq9OMqrXq1zFS/PnuQrs7sZ0FHABAaCGC4tPkfX2jtG3tN9Zv/zwUa\ncdE5FnQEAKGFAA5gX5065vdzvrp1vza9s99U//X4H+qS85P83g8AhCoCGJKk59/Yq6J/fGGqz554\nsS7on2hBRwAQ2gjgENHeFaWqPx6kycXmZRgLfnWZ+vdm71AA8BUCOEwtWbdD5V+YF8q4/7Zh6t0j\n1oKOACC8EMBhZs6TJTry1SlT/cH/uVw9Ejpb0BEAhCcCOEz8dParanAapnrv4Tt0/8hZFnQEAOGN\nAA5g3Tt36/BjTH7A9TZ7K+/K0n3bl3f48X3F6s0lAMDXCOAQ1VLwrrp7lDpFRWh71U4dqz8up+HU\nfaWFyhkwmnWPAcCPCOAQ01Lwblr+n6o+WidJ2l61s9nOP4fslU23CWEA8A8COAQYhqHyYtdb/62e\nky2bzabICFtT7fUDrkP6jYObCWAA8BMCOIgZhqEpSze7PLYmb3SLP1d54rDLeoW9yit9tRVvgwMI\nZwRwEGpsNHTLMs+D94zesT11yF5pqifH9epwb23F2+AAwh0BHEQanI26bfnbpvp5fRJ06rw3v7nl\nPoBzBoxuFn5njO2f3cEO2y5Q3wbn29cA/IUADgL1DqemrdhiqqcNOlt3TLhIkpRf8qbpeEvOBNxz\nH78op+FUn/hkje2f7dfgC5S3wQHAKmEZwO1dN9nfGhsiVVGSrmnFzcP3yh8ma/K1F3bosdN7pemV\nfa9JkuZlzOjQY7VHILwNDgBWCssADnTH7aeV+9i7ktKb1a/O6KfrRw+ypikvC4S3wQHASgRwADl6\n7JRmPVFiqo/PStGPLx/g/4Z8KBDeBgcAKxHAAaCy+oTmPfW+qd5t0H7Fn3NYP77c/RergpHVb4MD\ngJUIYAv9u6pWC3+7zVS/ddxgdTq7Us99/KaqT3GNLACEIgLYAns//0oPPP8PU/2O6y5SWurZXCML\nAGGAAPajjz47qoc2fmiqz554sS7on9h029NrZFlRCgCCDwHsB9v2HNYTm3aZ6vk3pWtgcoKp7sk1\nsrxaBoDgRAD7UPGHh/Tsa3tM9UVTMtQnKb7Fn/PkGtlAXVEKANA6AtgH/lr6b23c/C9TfenU4Urq\n3sXtz3tyjSwrSgFAcCKAvejl4n36c8lBU73w11eoe3xMmx/Hk2tkWVEKAIKT2wB2Op2aP3++9u/f\nL5vNpnvvvVcxMTHKy8uTzWZTamqqCgoKFBER4Y9+A9La1z/R5g/KTfVH7xyh+C6d2vWYbb1GlhWl\nACA4uQ3gzZu/3vZuw4YNKi0t1UMPPSTDMJSbm6vMzEwtWLBARUVFGjNmjM+bDTRPbNqlbXvMbwE/\nPiNLXWL88+YCK0oBQHBymxJXXXWVRo0aJUk6dOiQEhISVFJSooyMDElSVlaWtm7dGlYBvPyFD7T7\nYI2pvurukeoUFen3flhRCgCCT5tepkVFRWnOnDl688039eijj2rr1q2y2WySpLi4ONXW1rb684mJ\nsYryMJiSkrp6dH9PREbY2nWOXy9/Swcrzc9107Jxioz0/lvwnvTZlvt+/1h75+BNgdBDRwVz74GK\nmfoGc/W+jsy0ze+TLl26VHfffbeuv/561dfXN9XtdrsSEszXsn5XTc0Jj5pKSuqqI0daD/WOcDYa\nktTmc+Q++o6On3CY6s/MyVaEzabqartX+zvDkz7d3dfVTD2dgy8EQg8d4es/q+GImfoGc/U+dzN1\nF85uA3jTpk2qqqrS7bffri5dushms2no0KEqLS1VZmamiouLNWzYMM87D3CGYWjK0s0uj62ek930\nDgAAAO3hNoDHjh2ruXPn6he/+IUaGho0b948nXfeecrPz1dhYaFSUlKUk5Pjj179oqXg7RwdqZV3\njbSgI99ZdPlcq1sAgLDlNoBjY2P1yCOPmOrr1q3zSUO+1tK6yY2Nhm5ZZg7eXoldtOT24RZ0CgAI\nZWG1EEdPWYJKAAAJPElEQVRL6yaveqVMzurkZvcd1Leb5v33pf5uEQAQJsIqgFtaNznqnM+aAjht\n0Nm6Y8JF/mwrrPE2OIBwFVYB3NK6ybbOdbryh8mafO2Ffu4IABCuwmb9yFOnGxRR7/or4X0TehO+\nAAC/CvlXwCdONWjJuh0q/9KuyB4DFD3oQ9N9WDcZAOBvIRvAtSdOa9HvtuvLY6eaas7qZP34nBS9\nVvEq6yYDACwVcgH8VV29Fqz+u+pONl+56t7JGTq3Z7wkqaT6bUmsmwwAsE7IBPDRY6c07+n35Who\nbFa/79ZMJZ8VZ1FXAAC4FvQBXFVzQnNXvd+sFhlh05Lbhuns7l0s6goAgNYFbQCXH6lT/uq/N6vF\nxkRp0S2ZSuwaY1FX3sU1sgAQuoIugA3D0Jwn32v25arErjEquPkyJcRGW9gZAABtF3QBfOq0syl8\nk8+K1bwbL1Vc504Wd2U9Xi0DQHAJugDuEhOlR+8coeioCEV3irS6HQAA2iXoAliS4rvwihcAENyC\ncinK/JIlyi9ZYnUbAAC0W1AGMAAAwY4ABgDAAgQwAAAWIIABALAAAQwAgAUIYAAALEAAAwBgAQIY\nAAALBOVKWB3FuskAAKvxChgAAAsQwAAAWIAABgDAAgQwAAAWIIABALAAAQwAgAUIYAAALEAAAwBg\ngaAL4O1VO3Ws/riqT9XovtJCba/aaXVLAAB4rNWVsBwOh+bNm6fy8nKdPn1a06ZN06BBg5SXlyeb\nzabU1FQVFBQoIsI/Ob69aqd+W7a+6fYhe2XT7fReaX7pAQAAb2g1OV999VV1795d69ev1zPPPKNF\nixZpyZIlys3N1fr162UYhoqKivzVq14/8JbL+hsHN/utBwAAvKHVAL766qt15513SpIMw1BkZKTK\nysqUkZEhScrKylJJSYnvu/xG5YnDLusV9iq/9QAAgDe0+hZ0XFycJKmurk533HGHcnNztXTpUtls\ntqbjtbW1bk+SmBirqKhIjxpLSupqqvVNSNa/j5Wb6ucmJLu8P5pjRr7BXL2PmfoGc/W+jszU7W5I\nFRUVmj59uiZNmqRx48Zp+fLlTcfsdrsSEhLcnqSm5oRHTSUlddWRI+Zg/1HfkfrtsfWm+ui+I13e\nH99qaaboGObqfczUN5ir97mbqbtwbvUt6C+//FKTJ0/WrFmzNGHCBEnS4MGDVVpaKkkqLi5Wenq6\npz23W3qvNN08ZJIibV+/mu4Tn6ybh0ziC1gAgKDT6ivgJ598UsePH9fKlSu1cuVKSdI999yjxYsX\nq7CwUCkpKcrJyfFLo2ek90rTK/tekyTNy5jh13MDAOAtrQbw/PnzNX/+fFN93bp1PmsIAIBwEHQL\ncQAAEAoIYAAALEAAAwBgAQIYAAALEMAAAFiAAAYAwAIEMAAAFiCAAQCwAAEMAIAFCGAAACxAAAMA\nYAECGAAACxDAAABYgAAGAMACBDAAABZodT/gQLXo8rlWtwAAQIfwChgAAAsQwAAAWIAABgDAAgQw\nAAAWIIABALAAAQwAgAUIYAAALEAAAwBgAQIYAAALEMAAAFiAAAYAwAIEMAAAFiCAAQCwgM0wDMPq\nJgAACDe8AgYAwAIEMAAAFiCAAQCwAAEMAIAFCGAAACxAAAMAYIEoqxv4rsbGRi1cuFCffPKJoqOj\ntXjxYvXv39/qtoKOw+HQvHnzVF5ertOnT2vatGkaNGiQ8vLyZLPZlJqaqoKCAkVE8P8vTx09elTj\nx4/XmjVrFBUVxUy9YNWqVXrrrbfkcDg0ceJEZWRkMNcOcjgcysvLU3l5uSIiIrRo0SL+vHbAhx9+\nqAcffFBr167VwYMHXc5x48aN2rBhg6KiojRt2jRlZ2e7fdyAmv7f/vY3nT59Wi+++KJmzpypBx54\nwOqWgtKrr76q7t27a/369XrmmWe0aNEiLVmyRLm5uVq/fr0Mw1BRUZHVbQYdh8OhBQsWqHPnzpLE\nTL2gtLRUH3zwgV544QWtXbtWlZWVzNULtmzZooaGBm3YsEHTp0/Xww8/zFzb6emnn9b8+fNVX18v\nyfXf+yNHjmjt2rXasGGDVq9ercLCQp0+fdrtYwdUAO/YsUMjRoyQJKWlpWnXrl0WdxScrr76at15\n552SJMMwFBkZqbKyMmVkZEiSsrKyVFJSYmWLQWnp0qW64YYb1LNnT0lipl7w7rvv6vzzz9f06dM1\ndepUjRo1irl6wcCBA+V0OtXY2Ki6ujpFRUUx13bq16+fHnvssabbrub4z3/+UxdffLGio6PVtWtX\n9evXT3v27HH72AEVwHV1dYqPj2+6HRkZqYaGBgs7Ck5xcXGKj49XXV2d7rjjDuXm5sowDNlstqbj\ntbW1FncZXF5++WX16NGj6T+IkpipF9TU1GjXrl165JFHdO+99+ruu+9mrl4QGxur8vJyXXPNNcrP\nz9eNN97IXNspJydHUVHfflrrao51dXXq2rVr033i4uJUV1fn9rED6jPg+Ph42e32ptuNjY3Nnjja\nrqKiQtOnT9ekSZM0btw4LV++vOmY3W5XQkKChd0Fn5deekk2m03vvfeedu/erTlz5qi6urrpODNt\nn+7duyslJUXR0dFKSUlRTEyMKisrm44z1/Z59tlndeWVV2rmzJmqqKjQTTfdJIfD0XScubbfdz83\nPzPH72eX3W5vFsgtPpZPOmynSy65RMXFxZKknTt36vzzz7e4o+D05ZdfavLkyZo1a5YmTJggSRo8\neLBKS0slScXFxUpPT7eyxaDz/PPPa926dVq7dq0uvPBCLV26VFlZWcy0gy699FK98847MgxDVVVV\nOnnypIYPH85cOyghIaEpALp166aGhgb+DfASV3O86KKLtGPHDtXX16u2tlb79u1rU34F1GYMZ74F\nvXfvXhmGofvvv1/nnXee1W0FncWLF+u1115TSkpKU+2ee+7R4sWL5XA4lJKSosWLFysyMtLCLoPX\njTfeqIULFyoiIkL5+fnMtIOWLVum0tJSGYahGTNmqG/fvsy1g+x2u+bNm6cjR47I4XDol7/8pYYO\nHcpc2+mLL77QXXfdpY0bN2r//v0u57hx40a9+OKLMgxDt99+u3Jyctw+bkAFMAAA4SKg3oIGACBc\nEMAAAFiAAAYAwAIEMAAAFiCAAQCwAAEMAIAFCGAAACxAAAMAYIH/D7cHUMBU1CCwAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112d81748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(x,model([25,0.5],x))\n",
    "errorbar(x, y, dy, fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1376.11719488 -59.3058585913\n"
     ]
    }
   ],
   "source": [
    "print(loglike(theta,x,y,dy), loglike([25,0.5],x,y,dy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 1.75304421234e-26\n"
     ]
    }
   ],
   "source": [
    "print(exp(loglike(theta,x,y,dy)), exp(loglike([25,0.5],x,y,dy)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What you really have to do is to turn a minimization—as provided by `scipy.optimize`—into a maximization problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import scipy.optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function conveniently designed such that the maximum of `loglike` corresponds to its minimum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fun= lambda theta : abs(loglike(theta,x,y,dy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1376.1171948829217"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fun(theta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 58.32162215747249\n",
       " hess_inv: array([[  4.16161962e+00,  -6.36138321e-02],\n",
       "       [ -6.36138321e-02,   1.38985175e-03]])\n",
       "      jac: array([  0.00000000e+00,  -4.76837158e-07])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 24\n",
       "      nit: 4\n",
       "     njev: 6\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([ 24.81802587,   0.47417324])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.optimize.minimize(fun,[0.,1.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 4\n",
    "\n",
    "The probability of detecting $n_i$ counts in this bin is \n",
    "\n",
    "$$P_i={m_i}^{n_i} e^{-m_i}/n_i!$$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}