ex2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn.linear_model\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Practical exercise for Example Sheet 2\n",
    "\n",
    "## Confidence interval for climate change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4EAAAB7CAYAAADKS4UuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3LElEQVR4nO2de5DcV3Xnv6dfMz09mtGMRpb1siUZFMmWsbGFYyVUILFDZGqJzYYQTMo4hCqHqmTXkA0bEvLYdYosy2YTnEcVURFYm8o6BEPA3rLNYieY7CLHlliZkTxGyJIVjTzWY2bUo+npmX7d/aN/d/xTq6f79zi3f/fXfT5VXRr1TN++v3vuOfdxzj2XlFIQBEEQBEEQBEEQeoNE1BUQBEEQBEEQBEEQOocsAgVBEARBEARBEHoIWQQKgiAIgiAIgiD0ELIIFARBEARBEARB6CFkESgIgiAIgiAIgtBDpKKugBfGxsbUli1boq6GIAiCIAiCIAhCJBw8ePC8UmotR1mxWARu2bIFBw4ciLoagiAIgiAIgiAIkUBEJ7nKknBQQRAEQRAEQRCEHiIWnkBBEARBEJqTL5bx4NNHMTE1h53rh3D/7dsxnE1HXS1BEATBYsQTKAiCIAgx5sGnj+L5EzNIJxN4/sQMHnz6aNRVEgRBECxHFoGCIAiCEGMmpuYwmssgnUxgNJfBxNRc1FUSBEEQLMdYOCgRbQbwMIB1ABSAfUqpB4loFMBXAGwB8CqA9yulZk3VQxAEQSNhc0I3snP9EJ4/MYPRXAYzhRJu2ToadZUEQQiBjFVCJzDpCawA+A9KqWsB3Arg14noWgCfBPCMUurNAJ5x/i8IgmAcCZsTupH7b9+OW7aOolyt4Zato7j/9u1RV0kQhBDIWCV0AmOeQKXUFIAp5+eLRDQBYCOAOwG80/mzhwB8B8Bvm6qH0H3IDpl9xEUmEjYndCPD2TT+4D3XRV0NQRCYkLFK6AQdORNIRFsAvBXAvwBY5ywQAeB11MNFm33mPiI6QEQHzp0714lqCjFBdsjCkS+W8cDjR3D3vv144PEjyBfLoct88Omj2P/KNE5OL+ArL5zC3fv2s5TLzc71Q5gplFCu1jBTKGHn+qGoqyQIgiAIlyBjldAJjC8CiWgQwNcAfEwpdclWhlJKoX5e8DKUUvuUUruVUrvXrl1ruppCjJAdsnCYWERPTM0hXyxjbrGCdIpwcnrBysW5hM31DiY2O4TuRvqMYAtxGqtEb+KL0XsCiSiN+gLwb5VSX3fePkNE65VSU0S0HsBZk3UQosNUiKAkQQiHiUX0zvVD+MFkHukUoVqth6fZuDiXsLneQW92jOYyy5sdYWQfl5BnITjcfSYOcPdr0RMe4jRW9aLedAvGPIFERAD+BsCEUupPXb96DMC9zs/3AvimqToI0WIqbDNOO2Q2YiLM5P7bt+PqNQMoVxRW9acwPJCW8BUhUrg3OyQMvfvpxSgT7n4temIvpjx2vag33YJJT+BPArgHwDgRHXLe+10AnwHw90T0EQAnAbzfYB2ECDFlGEzskPXS7uX9t29fflauRfRwNo1H7ttzWRsK3umlPmgS3Y4nzhdQripsv3IQF4uV0BEDMtGxE0696cUoE92vAWCmUMKjBycBIHA7xklPes3m6rP7+WIZP5jM47nj03jkvj2hn7kX9aZboPqxPLvZvXu3OnDgQNTVEHzywONH8PyJGQz1p/DymYvoSyVwx671VhlaPQg8MT512YRRwhmETqL1xT2Q2tYHtb6MT15AoVRFri+F6zcOWxVCZsruxEE+ml6a3HLKJS7tZkJfZgolzBRKGM1lMJrLBG7HOOhJnMZ9TlnfvW8/Tk4vYG6xglQSKFcUfultm0M/c1z0plsgooNKqd0cZXUkO6jQm+iwzePnCyAQtozlrAsP0aErxXINS+UqJmeK1u9eCt2JqR10zhAgrS+nLyzi+LkCTs8WrQsh0+2YzaSwa8MwtqzJ4Q/ec13oSUmcwtBtDsnjDknj1BsdZfLIfXtY+owpOOWr+/X8UgWjuQzetG4wVDvGQU/iNO5zynrn+iHki2WkkmA9ux8XvREux2hiGKG30YZhYmoO6WSi/kokQhkd7h0nPYHIF8uYLZQwt1gOFc4Qpx0xm+tqyuNkM6ZCajgP7Wt9OX2hiP5MAsVylWXyxDmRN9WOcUrUYGNIXjPvC0cSiV4MRTOx8AWA50/MAAqh2tGknnCNWdzjPmfdVqorh6zvv307njtev8ppOJuWs/uCeAIF83AmIuHe4dZ12zySRSadwEAmGWr30uYd+EZsrqspj5PNmNpB515gzRRKyKaTWCzVkE0nWZILcdqIOHgiTKdUt/GOMW7vi27D8ckLqCmFYrlqrby5MZXcy1a90bL+uT97Ft889BpACDUeNI77/ekEakph/HTetz5y122lunLIWp/d/6W3bcbVawawZ9saq+QsdB7xBDbBZg9JHOFMRMK9w+2u2103bmTzLNq0A78SNteV2+MUB53WO+i6rh/98gGWunJ6SrS+jE9ewOqB9CUe2iC4Pb5cE/k4eOxMJWjQmEj+FBZu74vbw10s13D9xmHr5c6FqeRetrafewOhWq1hcqaIrWO5UCGr7nG/VKnh0KkLyKaTvr3T3HVrVVcOWXPKOQ7jqtAaWQQ2wcSdJ72sLJxGhzv0h3vg0/VzJ6V44PEjVsrb5jAqXTftaQpbR5P3GHHrNnddOScR3PrCOZE3ZWNNlDsxNYd8sYy5xQrSKcLJ6YVAcm6s24ffvhVf+j8nlv//+Xt2W2N3tE5vHsmiUKqgL5UI1R9t3sRqhLsP2bxgMwH3BkJj+929b3/gvuSu28z8Es7MLWKmUMLVawaQL5ZD65/Nspb7AeNP23BQIuonot8koq8T0deI6ONE1N+JykWFicHF5tA7N6bDlMJic8gKEI9kOBqb21LXbePqfmxbm8PGkay1E0Zu3W6Wsj2MLtp8aJ9TLqZsLGe52r6eOF/A2YtLSCbCJWhorNuvPXzA2nFG67QCcNeNG/HUx94Rqj+aDHnlHgdN9CFbx2gNZz25j22sVH6QvuSuW0Up1FRdnxUQSM5xkS8Qr40YoTlePIEPA7gI4C+c/38QwJcB/KKpSkWNCQ8Jt7KYSpxhu8fE5l0x9/MRAduvHMRAOhU6GY4pbG5LUx5aTp3W8n704CQG+1IYHcyw6LauqztlO4cu2hiNwCkXExOSfLGMJw9PYaFURb5YxqbRLMtCddtYDrMLJRSWqhgbzKBYruLV6YLvqIHGZz506gJu3LyapQ3yxTI++9QEnpk4CyLgth3r8Im9OwL3GW6dNhnyyj0Oct7FFxfvC2c9uY9ttCrfb19yf3Y0l8GWsRwG0imUq7VQGzu2yxewO5pI8IaXReAupdS1rv//ExG9ZKpCNmBicOFWFm0o3BPFYqnKNliZ9JjEwbgFwf18S5Uajr4+j+s2DFmbdayXMKHTWt6DfSnMFEo4dmZ++X4tjro+enByOWU7FHzpYrNNosJSBQrA2sE+a/SPUy4mJiQPPn0US5UaqlWF2UIJC6Uq9u5ahwcePxJIH9329earRnD8fAEAkKBLowa8yqXxmTeNZFnCp/WzP3n4dZTKNQDAk4dfRyaViLzPaEyebRo/nWcZB3W5J84XUK4qpFOEC4VyqDEwLt4XUxlMTRCmfPdn9R2J6UQisP7FQb4mznIL0eAlO+j3iehW/R8i+nEAXX1zu1bqz99Tv4vxo18+ENotzx16pw1FsVxlS9WeL5aRL5Zx8OQsJqbmcG5+iTXEJg7GLQzu59uxbhXSSWKRt04icXJ6AV954RTu3rff6hARGzEREqnl/aYrBrF6II1z80uYLiyhVKmFko+u6/tu3lT3Hjgp2/3oYrPsqienF5BfKFulf5xyMRHePDE1hx3rVmEkl0EymUA6SSBQ4NA+d9jZ3GIF775+PbaO5XDdhiEMpFO+5dL4zH/9od2h20CHoz16cBLzixWkUoRMOoFKLZhnwwSmwzULSxWWUFO351dBYWa+FPouPhszv2rccskXyzg3v4RiqYLDr+WXPd1Rj10mwy05bJDN8tXofp3NpJAgWj7LLZvT8cOLJ/BmAN8jon91/n8VgB8S0TgApZR6i7HaRYzu6EP9KXzj0Gk8eXgKd+xab0UYI3fiDKD+vMqJZ88Xy1g9kA48iWjmuer20AH38+kJHkf43qMHJ1GtKRCBNYmEjR7FONRR45Z3qVrDUH8ab75iFQ6dusDiZQvjJWuWXVXrtZ5cdJv+mfAWaBlvHcstt1mYzaxmMnVHEPiVS7NnDtsGbg93YamC+cUq+lMJ9KWT1kxIm0WVuNs2jIe2HlVTwfWbVof2ULvL3bVhGEfPXLxkYyeIDtqY+VXjlotSABGano+P0ptsMiKJwwaZkK+p+5Vt2lAUguFlEbjXeC0sRXf0E+cLKJVrqNaUb6NhalLLnaodqD/v2lV92LA6i3K1hnK1FqiuKxlZmwcvILysuJ/PPRmbyi86Xohk6CQSNofjdqqOHHrplndfKoEtYznWQTHMhKLZJtFwNo3VA2lrEgF1YsFvQqe5F2222UU97o3mMqgphelCCauyKdy2Yx3rpiBHHd36FsZ2NNug5MgIni+WL7mY+7adV6AvlQwla5vPcrvlsnZV3/IcIp1M1F8WnI+3fQFjQr7cmybdvqHfS7RdBCqlThLRCIDN7r9XSn3fZMVsQHf0ucUKAGCoP+3baDQq32efmlgeBMIMiCZ3vcMq9kpG1ubBCwi/AOF+PvdkLF8so1iuYqw/heGBdKAdedsHP8Bcco/GwS6MrBvL+/w9u98IJwtxFoSTVptEQeyNiYk898TE63eE1Wmb7+3iwD0OjA324d3Xr19+5iB3V5rY2Gk2VnF7aMP2+WaRNf9x705rIxs4WGkOEXReYcLu9OIChnvTxLaNKyE4bReBRPRHAH4FwCsAlPO2AvAz5qplB7qjn5+fQjlB2DSa9W00GpXvmYmzGBvsY/d0cHs2wii2yUyMJr0Gti2S3O141egAiLAcVhv0rIHtg5+p5B6Ng10YWbdbvATRHdvvETMxkeeemHj9jrDYtmjjpp33069cTMigEx5aneQjaF/kiqzRcNuIZuUBMBINE9Q2mrA7cV3AhJE/96ZJt9vAXsJLOOj7AVyjlCqZroxt6I4exmg0Kh8RjCw0OIwll2KbzMRoMkzQtsWrux33XLMm9KDfqcGP65m5duSbDXZhZN2svLC6Y3uoromJPPfExOt3CK1p1pfDyMWEDDrhoQ3bF01lBOeyEc3KAxDqO1ayg0HrKZs4b8DtuQuzaSJ0D14WgYcBrAZw1mxV7CWM0WhUvlKlhkOnLrArXhhjyb0jaMLIdsJLF3QS0WqREsZwc7ejCblwh1qa2JFvNhnj3Nix8R5Rblo9c9BFeicmJnHa9e+EpyfMxexB5WJiY6cZtt0t2olFaZh2XMnm2GSH4rqJ06lzsF6Jw1lkIRpIKdX6D4h2A/gm6ovBJf2+UurnzVbtDXbv3q0OHOiOWylMhTW6J8p+D7Y3+yxweRx/lLtnzerIfX6Is266re7et3/5ULwOCXrkvj0dr6Mpmj37xNQc6zOHbUObJ9eaMPrbCVptFD0xPoVyVWH7lYO4WKyEqnucssNyw90HOMvjlovt/R2wry9yj9NxGPdtk4GbVnUz0b/joDNCZyCig0qp3SxleVgEHgHw1wDGAdT0+0qpZzkq4IVuWgSaIoyxbDbJBhBo4m3KaLfzOEVpFFstUrrdcDd7dr17O9SfwstnLqIvlQh8tQrQ/W0I2D3ZWQktl1OzRVSrNYzkMtg6luu6jY5Owb1hZMMG1Er92oa6RUVQXW/2uY9++UDgdjRxJtDv99lu41rR6c3fbms/ITici0Av4aALSqk/5/gygYeVjEHYdPIcGb1MnW3iPqfCSauQFe6QCxMDAfeBc/3MT4xPsdwP1c1hK80yjcZlYNf6ly+WMVsoYW6xHKuQLdvgDn2zIZRupfHAhrqZop09DTpGNhsDw7Qj5/k9L2MI19zAlsVQq/lHULm0era4nmUU7Cbh4W/+mYj+CxHtIaKb9Mt4zSIiXyzjgceP4O59+/HA40eQL5ajrtJlaGOaTiYuOdAdlPtv345bto5ecn9Ys/e8wL0wayWPneuHMFMooViq4PBrebw6XYhEZq3aShvuR+7bgz94z3VsafW5ZB+2zGbPrp9561gO120YwkA6FaovcLehTZiQJzcr6aDWv80jWWTSCQxkkl23SO8kfm1uu7EqqA3nZKXxIGjdumF85hwjTcnYbzt7sWNcz22LzdT2r1ytYaZQuuTapqByseXZWhEHHRS848UT+Fbn31td73XtFRG2Z+kD+BdanDuC7XbA/O7itZJHGI8T525iJ3foTHg/TaWK7rbdftsO+3eqjivpoNtDe9eNGyU8KSR+7Ui7scpveZ28ky2ozbR5fNbt9+jBSQz2pTA6mGmq05x20dTY47edvdgx0/cQt4O7f7eKUNFy0d/p9X5NW6KbWmGzDgr+aesJVEr9dJNXVy4AgXgoYasdqKhptwPmd6erlTzCeJxs3nHz4v3klL2p/mSDJ4ITE32Gu+1N1HElHbTJQ8u9Ox2H3W7uscpE3+G2ATaPz7r9BvtSmCmUcOzMfFOd5vb4msBvO3uxY1x9IajN5O7fXuyf3++0eW6nsVkHBf94uSx+HYA/BrBBKXUHEV0LYI9S6m+M1y4C4uC94DojFUWa7iCDSzt5BJGZzYbMi/eT83xckDK99B3uXWq//ZW7f5voM7bdbdYMbpvYSW+lLeWZgFsuJvpOFNc2RHVmTLffaC6DY2fnMb9Uwd5dV16m09weXxP47Vte7Jif524lw6A2M4ox3+93xuH8exzmyIJ3vJwJ/B8AvgVgg/P/owA+Zqg+kcO1W9Vq9y7szh7XDnwU3jC/O11e5BFEZjbvuHnxfnJ6X4KUGUXf8fud3HU00We45WmijtzenE56K20pz4Q3h1suNttEjZdnjirKQ7cfUL9r7303b2LR6SgWL377VifPvgf9rij6t9/v5G7HONgdIVq8nAkcU0r9PRH9DgAopSpEVDVcr8jg2rlstXtnyy5zVIOLn50uL/IIIjObd9zisNMWh11V7jra3Gc0JurI7c2Jg7eSuzwTNp9LLtrrMj55ATWlUCxXre3fXp45qigPU/YhivEg6kyUcYi6sPU73dhsdwQ78LIILBDRGtSTwYCIbgWQN1qrCPASQuInzKSVEbMlFDGK0BpbDIgt9WhG1AOHF6KYmPj9Tu462txnNHGoo4m+w60z7crzaxdtsfnNcE8Ui+Uart84vGIfsiU9fyv89C/Ocd+U7sVhPODGxNwkCtsY9ZEIm+2OYAdeLou/CcBfANgF4DCAtQB+USn1ovnq1enEZfFeLqT2c2l1q7/Vv+O6TDsoXgxKL1zU3QtwT96imAxGfSbQRB07XZ6pMuP0/Rz4tYs221E/F1vb/BwaP/2Le9yPC1w6aEqXo5ibxMEudZPdEYLDeVm8l0VgH4AqgB8DQAB+CCChlFriqIAXOrEI9DIQ+hksWxkU/bsnxqdQripsv3IQF4sVKxXUzzPHhTgYe27aDQa92CYm4B50TQziMjEIj1+7aLN++ekP3TYecI/7cYHLBkRpS7jlEge7GKXdsdmG9Rqci0Av4aD7lVI3ATjiqsD3AXTVhfHcWShbhQHo301MzS0rdDqRWNFVH6XyxeF8ml/axcl3o7FrFxZiyznVuGN7ghJTZXJjuw76tYs2h+n6CTeMajww1R9MZZ+2HS4bEKUtccvl3MUlENUXSUH7R5TP4rV/m7A7Xr9b5gjdyYrZQYnoSiK6GUCWiN5KRDc5r3cCGGhXMBF9kYjOEtFh13ujRPRtIvqR8+8Ix0NwYCoLZSu8Zo7SykcAvnHoNPZ+7tmmmZ4kE5Q3vC6IbLxDMCjt+locFgYm4NYZ7gx0cboXkhPbdbCb7KKfjITcz+1V/0z1hyjGfRP4tWNcNiBKW+KWC1E9aUWY/hHls3jt3yb6otfv7tU5QrezYjgoEd0L4FcA7AbwAuqhoAAwB+AhpdTXWxZM9FMA5gE8rJTa5bz3WQAzSqnPENEnAYwopX67XSU7EQ4aBV53YHQIwInzBcwWSkgmCZtHBi4LV4hDOIMNtGunKMN/ojpj0at9x/ZzJd14JtALUelgHNqmm/Cqf90YksnZ1/zaMdvPBPqtw97PPYuFUhVD/WlsGs0CCr77R5TPEmX/9vrdvTpHsJGOhIMqpR4C8BAR/YJS6mt+C1ZKfZeItjS8fSeAdzo/PwTgOwDaLgK7Fa8hQjoEYG6xAgAY6k833YmRnRpvtAt/ijL8x1TIRbu+pttEp4kfP53HA48f6fpJMLfOeNXpKDMO2hyaqPGqg9wTNwl56ixe9a8bQzI5+5pfO8ZlA2ywJQ8+fRRLlRqqVYXZQgkLpSruvHFD+w82wHWZfRCi7N9ev7sXs9T2Am0viw+yAGzBOqXUlPPz6wDWrfSHRHQfER0gogPnzp1jrEL80CEA2XQCfekkNo1mm4YrxCHMywbahT9FGf4T1UJet8n1m1YjQYRsOmllGB43UemM7eGOJvATsuZVB722o9fvjsNGmomw/6jwqn8mbHLU7cjZ17pp7Pcrl4mpOexYtwojuQySyQTSSTI2Zuu6/dyfPYtvHnoNILDY7yjnHF6/m/sie8EO2mYHDVV43RP4v1zhoBeUUqtdv59VSrU9F9it4aB+abf7ZENoRi9hor29hlyYknW70BCTfawbrnTwSjeGt7XDRDgRdyiTiTpy97FuCsvqpisG/ML5/d009tt8DYL+rlOzRVSrNYzkMtg6lusJ+y3YQ6ezg3JyhojWK6WmiGg9gLMd/v7LiJPxbBeuYENoBgdxkYmJ0DGvIRemwtbahYaY+N5mV6YELTtfLOOzT03gmYmzIAJu27EOn9i7Y8X+E5XOdGN4WztMeNm8tqP+bgCYKZTw6MFJALjMtpgIeeLWGb/taLM9NaV/Xto8aq8vZ1+Ly9jvpS/6lUsnwxR13fLFMmYLJcwtlq203zbrvGAXbcNBAYCIfoKIPkhEH9KvgN/3GIB7nZ/vBfDNgOWw0YthWbYTF5mYmER4DbkwNYFpFxpi4nu1vIvlGpbKVUzOFAOX/eDTR/Hk4dcxv1jBxWIFTx5+naX/cIeOxSHjIDd+Q9a8tLnXdtTffezsPGYKJQz2pZraFhMhT9w647cd42JPOfHS5lGHUPZieJ2XvuhXLp1sR123zSNZZNIJDGSSVtrvXtR5IRhtPYFE9GUA1wA4hPql8UA9G+/DbT73COpJYMaIaBLAHwL4DIC/J6KPADgJ4P1BK85F1LuBUWFzdrCJqTkM9adw4nwBc4sVnJ+fYrnkdHzyAgqlKnJ9KVy/cbhjh7lNtJGp7263o2zCg8W5uzoxNYdqDcikEgABlVqNxVPC7c0xsXNv++6v3x17L23utR31dz96cBKjuQzetG4QUOiIvefWGb/tGMUYF3Vf9NLmvZjoIurLw730RZvl4q7bXTdujMzGtmr7fLGMJw9PYaFURb5YxqbRbM/MawX/tD0TSEQTAK5VJg8PtsHkmcCozwVEBddzmzjX8MT4FGYXykglCASgL53EnTduCFyuruNMoYSZQgmjuQxGc5mOXQcQxTmjZiGWF4uVjl+B4OdszlB/Ci+fuYi+VAJ37FofaIB94PEj+Mah0yiVawC89Z1W8tH1f/TgJAb7UsuLBxvPgNh+3YVfTJybjKKNom7HKMa4qMfVqNs8LKY2LjnlEqSsKPpF3PtCM1q1Y5AxUIgXnT4TeBjAlQCm2v1hHLFh1ykKI8W1O8y5y6x3/reN5fD8q7NQSmHdUH/onSxdx9MXiujPJFAsV0PX1Y/MTIaNroQ7xLJarWFypoitY7mOXYHQWI9W3hzO3dX7b9+OpUr1kjOBYTwluv6Dfal6OOGZ+eUNhCAE1XUTZ2naEfV1CSa8ztz2ntNbaYooxrioI2yibvOw6H7l3rgslqodP0/KXVYUfTFqO2aCVm2vs6Wemi1ibrFiNFuqEH+8LALHALxERM8DWNJvKqV+3litOogNg0UURoprgsU5UXMbtitW9SFfLGPrWC50ubqO2XRyeUANW6Yfmenvd3u6TN/BZ8sBdi8TBU4dHM6m8en3vgWffq/3z7Tqw7r+o7kMjp2dx/xSBXt3XRl4UA2q660+pxeIJ84XLvP6hoFzwhhk8dtqwhh0Mc1t76Ne7Hih1TMHaUcvn/E7LnSjtyYMJjYuAd7xOkhZUcy34qCjfmnV9vp37rlTL+uS0BoviWH+E4C7APwxgP/ueglMRGGkgiSmaJaogTPBhftA+HA2javXDHgut1USCV3Hjav7sW1tDhtHsqHr6kdm+vuPny+AQNgyljN+WNuWA+xRJ1/wQrM+rPvTifMFHHltDuVaDaO5DN5386ZQyQeC6roXb+W2sRwUFF49X2CRNUcilzD3arVK+GBL4oNWbRT1PXReCNKOXj7jd1xoVabt7WiifrpfZdNJLJZqyxuYYe0n53gdl+RWcRiD/NKq7eMiF8EOjN4TyIWN9wRy7lxGfX7CK17OTgVpD47zD50+dxfk+1qdceLeCbdlZ71T9TB1FxvHOcVm5frtp60+F+TsnImzbM3qCMDIvVomzgtye8VstOuN9R0/nUc2nfTVjo1tXyxVcP2m1aF0r5U8bT/ranLs8TMm2mLzgxLlnZGm6eZnEzoP55lAL4lhbgXwFwB2AsgASAIoKKU6tp1i4yKwFy96NTVQc7Rlp5NIBJFZu8Pctk0Y4wR3+5m6zN3EmUBbEjQ0azMASCcTOHG+gNlCCckkYfPIQOjvM1H/uPShMDQ+Y00pJIhC9Z0gZbQrM+wmR9Dv8ku+WMbezz2LhVIVQ/1pbBrNAgqRyLlTG6EAjMxXunkMNPVs7nLPXVwCUT2CgnMjmTs5kRCeTieG+UsAHwDwVQC7AXwIQM/7lzlDOLnj5NulDw5qwL2cnQrSHhxtafLqgmb1CiKzVmecbD63EIdNCu72M3WZe1Bdb/Y59yBdUwrFctVz+I+J/rZSmz1/YgabR7IolCroSyVYQpRMJJiwpQ+Z9Bo0po6vVGrLXjyv7djY9m5v4lB/Ck+MT/mueyt5cusidzKzpUoN1arCbKGEhVIVd964IVT9gmJCp3WY7lB/Ct84dBpPHp7CyEAGCsDawT7WHAbdPAaaejZ3uTr8/uarR3zLpfH5Pvz2rfi1hw/g5PQCAKBaU2zJiQS78LIIhFLqGBEllVJVAF8iov8H4HfMVs1uWg1MUU+aWyWRCJOEpnGg/vDbt+KBx49gYmoO+WJ5eWDwmwjA72e91I1jUsg9+WicyOuzJFxtwMlK8rE1u1rUd7FFgVuXi+Uart847Fkuncy8GTbz60r2lLsPcrSJu67b1uZw4+bVeOXcvK8+ZCpR2EoLFr9lN7a92xvx8pmLINAlZ/u8lN9Knty6yJ3MLGwmRq75gsmN0BPnC1gqVTG/WMGZuSVk0wmsXdXHuqAxtfHGQVidNPVs7nLzxTKGs+lAC83G53vu+DROTi8gnSJcLFaQTBBbciLBLryEg34XwO0AvgDgddSvivgVpdQN5qtXx8ZwUO7wLM6FY6vwGY7QmmZ30F0olH2FIpgMY+DC9GLe5jZw1+3gydnletkQ1tbJECWbCaPLcQr16VSYGIe+2xrarsslAKdmi8gXy6hUa7hu43Boubvb7dXpAraM5TCQTllhK5ph23n+Bx4/gv2vTC97cq5eM4BH7tvDtlkSBv18p2aLKCyVQVS/u7dcVVg/3I/RXPj7dk3VP2zovds2Ts4uhOrXnTgT2Gwj2atctM0BgGNn53FmbhH96SSIgFKlZkTeQnA6HQ56D+pZRH8DwMcBbAbwCxxfHmda7VwGcf1z7v62upKAY0dK17XxDjo/hlG3EVA3ZPNLFdy6bY1Vk9BWIXgcxtzdT9au6gs8YTIxwLjrNpxNI18so1iqdOyKi1aspCsmEwGFKc/UQiuMLuu+rSd52XRyuS3dnpeg7WBKT0zsRDfW9fP37I60rqa9BlvHcjj8Wh7pZOoSuQfVH7ed1P0pnUhY583RBPUiN+vTHF5KPYGfW6wgnSKcnF4IJA8T3nH9fOfnp1BYIgxkkqhWFdJJhLoux5R3nyOCRY8v7jsalyo1HH19HtdtGArUr01djeEu1/3sfvuitg36mfvTSVSqNaSSCaQShMG+FDaOZJfHL6F7aHtFhFLqJAACsF4p9Z+VUr+plDpmvmrxozGt/EK54jklMedEx30lgVKAAvCVF07h7n378eG3bw2dPljXdai/voeg76Dzk3pZp20+dnYeM4USBvtSkaZ69wpnanqu1NW6TgTgG4dOY+/nng2dqrzZdR2dvOKiFaYWBZztmC+Wcfe+/fjKC6fw0tRFHD9XwOnZImu7caQCb9aWHH3cRj1pJMz1FSuVFcT2N7KSXMNeReAuty+VwPYrB9l0SNctyPnUuNCsT7e6xsQrO9cPIV8sI5UEqtX6xN6WkDv9fN/6+Duwfd0gKlWFoWwaV60ZCHVdDqd9aFYuATh6Zh6vnJ3HifMFrMqmfF/HUyxXl+9o3LFuFdJJYrl2wa3Hn/qHH+D3/mGc5XqRMH1R24b5pQpGcxncsHkYawb7kEwQPnDLVfinT/w0vvrRnwh1PZJgJ209gUT0HgB/gnpm0K1EdCOAB7rlsnhOtAHaNpbDy2cu4tXzheW08u3g3P3VxmBiag4npxcwt1hBMgkcPXMR7//893DHrvW+drsbd+2uWTuIQ6cuhEr2oHcYHz04idFcBm9aNwgoWDP4rQT3Yj3sLrI72UO1pqBU/RB32N19d932XFP30H70yweWw9TSiURksjLlKXGffymVa6Ha8cGnj152pqJQqgAAHj04CQCBvWOcnqtmbcnRx23TEzeN4ezlag21mlqOaAhS1zC2v5GVvAZho0W4PXYreV38nk+NC6Y2n+6/ffvyGazhbBrDA+lQGx0mIkOGs2k8ct+epmH4QTCdKOXE+QIAhUoNvhP2aJuo72YczWUwt1jBu69fz3Y2V+vxk4dfB4Fw3YahSM/bu23D8ydmkErU5bJ315Vdp8fCpXi9LP4WABcAQCl1CMBWYzXqINyXvGoDlM2ksGvDMLasyXneOTFxwad7h7G4VANAWChVfe+8Ne7aKSjcsnUUCsBdN27EUx97h+8dIm103nfzpnpYqEIsLnLl9Epw7CK7kz0slKpYqtRTlXNkrNXhw9o7dM3aQSsu3TV1Ga6W7dxifbEWph0npuYwnE2jUgWSCUK5qlCpqtBeb7eH8eT0AvYfnw61i96sLTn6uG164sYdzr5UrqLmnIsPEtGgCWP7m9FsbDIRLRJGh9zjwsnpBeQXylZmduTClEdaL7B+6W2bcfWaAexxjkUExZSXjVMPTbWl24b3JRMYyCSRTCZ8JezRurFxdT+2rc1h40iWbfPpgceP4NGDk5gplAACqjWgUqtZozemxlbBXrycCSwrpfJE5H7P/hvmPcCdhY3jjA4n7h1GIiCTpEAT28bJx/FzBbbD/qZ2+U0lCLGtvu7sdKW5RdQUsGk0y+Iha9SPGzYPL3uKohwgTJ2vcJ9/KScoVDvuXD+ECwv1TaULCzUM9qVQrioM59KhvN6XeBidxWrYxX5jxtqlShXn55cwXVjCbTvWBZKzzdlVtT3LF8uYLZSQICCVSoS6vsJt+3WSp7v37Q9sg5qNTSaiRcLQ7NywntRzeOejzrLdiLtP37B5GKVKLZSM3XDaNJuvWtDottThw+On8yznzBtt+PYrB3GxWMEtW0c9l2tqfNE6PdiXqh+FOTOPZAIgJFj1Jgymnl2wFy+LwCNE9EEASSJ6M4B/D+B7ZqvVGbiNpW0TH3cIhw59CjKxNZm6mdvocC/sm01EbKqvO9nDYF8KRAAUWPqfycW/jei+yKHH9c/Uy9DJE7Sstdc7iB5pD+PcYgWpZL1/cnpkH3z6KF48lcf2daswUyghk0oEmpSZmEzki2V89qkJPDNxFkTAbTvW4RN7d/iun9YZdzi7Dt0MOgF19xmi+i6p32sS3LgTZ80USvjqgVPYsDqL2YVyqMU5J+5xYTibxuqBNKsHwdRVGUFpFk5rS93c2HzVgsadmGr/K9M4PVvES6/N4bnj04EyozaWa9tcDHhDp0dzGRw7O4/5pQrec8N6EMj3NTKCwIWXReC/A/ApAEsAHgHwLQB/ZLJSnSKIsWy1O2njLgqHUbTRoK4E98JeT0Tcl+WGnTBy1rfZ2T2u3XLTkwnbdvo1HHrcrAwOPXJ7GHU6eU59tNmL8ODTR/Hk4ddRKtcAAE8efh2ZVMK3rNxyCHp3YSNuebuveAibIVRn60sm6lkj9SQy6OKcE27b02gPxk/nre2LNutJ3MZrjsyojdg4F3OPp3LeTrCFtotApdQC6ovAT5mvTmcJYiw7uTvJOUkOYxRtNKgr0ep6jCBtx5kspFV9gy60TMrGxGQibpfQc8IhK7eH0cTC2WYvwsTUHKo1IJNKAFQ/SxNk8m3annG0YWPirEKpinRKWXVhs+kojppSKJaqVvZFm/UkbuP1DybzSKfIusyoGq55mKnFua2bqUI88HJZ/G4AvwtgC1yLRqXUW4zWzIVNl8U3XuRbLFVw/abV7AqoE0C4M4bt2bYmNsY9KppdZK/PBARpO/dludVqDSO5jO87EVeq48TUHLatzS2Hg3D1H5sHBZsvoRfs7zvfOHR62RPYl07izhs3WGcTTVxG7r6zrFsvbDYxtnbiku44lGsrcZjnuMcsv5ewdwLb6yfw0+nL4v8WwCcAjAOocXxpnGncAawpZcQzyJ0AoldwX4/BcZ0BZ7IQjXvH+8VTedyydZR18WPbWRo3ppNJCMHgvHbC1ET2/tu3Y6lSveRMoI2hbpyeGHcCjdUDaeT6Ul17YXMz71rYdnTbwv2vTOO54/uXN56iiqxpBaft1no4PnkBc4sVzC6UkUwEP0trAu6rJ0xgc+gvYH/9BLvxsgg8p5R6zHhNYkKjS3/8dB7ZdNLIfTcmE0B0O1zhOiYOmps22jYPCqaTSQjB4Jx8mtqEGM6m8en3vgWffm/oomJDnEL7wmIiY6TbFuaLZeSLZdx89Yh1m2MaTtut9XCmUMKZuUUQEfpTicBnaU1hex83FfrLtVlmc2iyYD9e7gn8QyL6AhHdTUT/Vr+M18xSGu/KuX7jsLH7boazaQz1p1CuKPYEEN0O9303cbgjqVPlh8Etlz3XrMEj9+0J3Kbc93z2MpyTT5s3IQR70Tb2+k2rkSBCNp0Mfc+d2xbmi2UMZ9NW90tO2631sFiuAgRAKWTSicBnaXsVU3fncd3nKHf7CWHw4gn8MIAdANJ4IxxUAfi6qUrFCVOHfU0ngOh2bN5dNJ29zebscJxysTnsNW5w7ibLzrQQBs5NBLctvHrNABRgdfg5p+3WephNJzFXLANEKJVr6EsnrdoYtB1Tcwmufm7zXEewHy+JYX6olPqxDtWnKTYlhhEEITwcoTCNiSQksUxwOM/x9VpyC4EXU4kueq1fxuFMYC8jCV2EoHAmhvGyCPwSgP+mlHqJ4wuDIItAQYg3jROwUqWGQ6cuhBoA9UXD+qzP1WsGQl00LAhC9JhYrPXaAlCwH+mTQlA6vQicAHANgBOoXxhfjy6P8RURonyC0Fkadz2nC0t48xWrQnnx4pBeXBCE6BGviyAI3UKnr4jYy/FFNiFniXoDWezbQ+P5h/PzS8v3ngU9nzOcTWM4m8bNV48sLyYl4YEgCI1IsiJBEITLaZsdVCl1stmrE5UzhQwIvQFX9i0hPI1Z727beQVLRjObM6EKgmAHYicEQRAux4snsGvQnqET5wsoVxW2XzmIi8WKlVnChPDIYt8emmW947pA3NZMqIIg2IHYCUEQhMtpeybQBrjOBOpzAUP9Kbx85iL6UgncsWu9VWGCEsLIh5wDEQRBEARBELoFzjOBXi6LZ4eI9hLRD4noGBF9slPfqz1D2UwKuzYMY8uaXOiLv7mREEY+5BJVQRAEQRAEQbicjoeDElESwF8B+FkAkwBeIKLHOnEFRRwuMZYQRj7kElVBEASBC4nUEQShm4jCE3gLgGNKqeNKqRKAvwNwZye+OA6eITnALgiCIAj2IZE6giB0E1EkhtkI4JTr/5MAfrzxj4joPgD3AcBVV13F8sVx8AzJAXZBEARBsA+J1BEEoZuwNjuoUmofgH1APTFMxNXpGHFYqAqCIAhCrxGHIyWCIAheiSIc9DSAza7/b3LeEwRBEARBsJI4HCkRBEHwShSewBcAvJmItqK++PsAgA9GUA9BEARBEARPSKSOIAjdRMcXgUqpChH9BoBvAUgC+KJS6kin68GJZAwTBEEQBEEQBCEuRHJPoFLqCaXUdqXUNUqpT0dRB04kY5ggCIIgCIIgCHEhkkVgtyEZwwRBEARBEARBiAuyCGRA7vYTBEEQBEEQBCEuyCKQAckYJgiCIAiCIAhCXLD2nsA4IRnDBEEQBEEQBEGIC6SU/fewE9E5ACcZixwDcJ6xPIEHkYudiFzsRORiLyIbOxG52InIxU5ELnbyY0qpVRwFxcITqJRay1keER1QSu3mLFMIj8jFTkQudiJysReRjZ2IXOxE5GInIhc7IaIDXGXJmUBBEARBEARBEIQeQhaBgiAIgiAIgiAIPUSvLgL3RV0BoSkiFzsRudiJyMVeRDZ2InKxE5GLnYhc7IRNLrFIDCMIgiAIgiAIgiDw0KueQEEQBEEQBEEQhJ5EFoGCIAiCIAiCIAg9RNcsAonoi0R0logOu967gYj2E9E4ET1OREPO+2kiesh5f4KIfsf1mVed9w9xpmHtVXzKJUNEX3Lef5GI3un6zM3O+8eI6M+JiDr/NN0Do1y+Q0Q/dPTlEBFd0fmn6R6IaDMR/RMRvURER4jofuf9USL6NhH9yPl3xHmfHH04RkQ/IKKbXGXd6/z9j4jo3qieqRtglkvVpS+PRfVM3UAAuexwbNwSEf1WQ1l7HVt2jIg+GcXzdAvMcpE5GRMB5PLLjv0aJ6LvEdENrrJEX5hglos/fVFKdcULwE8BuAnAYdd7LwB4h/PzrwL4I+fnDwL4O+fnAQCvAtji/P9VAGNRP0+3vHzK5dcBfMn5+QoABwEknP8/D+BWAATgSQB3RP1scX4xyuU7AHZH/Tzd8gKwHsBNzs+rABwFcC2AzwL4pPP+JwH8V+fndzv6QI5+/Ivz/iiA486/I87PI1E/X1xfXHJxfjcf9fN0yyuAXK4A8DYAnwbwW65ykgBeAbANQAbAiwCujfr54vrikovzu1chc7Ko5PITetwAcIdrfBF9sVAuzv996UvXeAKVUt8FMNPw9nYA33V+/jaAX9B/DiBHRCkAWQAlAHOdqGev4VMu1wL4R+dzZwFcALCbiNYDGFJKPafqvfxhAHeZrXl3wyEX87XsPZRSU0qp7zs/XwQwAWAjgDsBPOT82UN4o//fCeBhVec5AKsdffk5AN9WSs0opWZRl+fezj1Jd8EoF4ERv3JRSp1VSr0AoNxQ1C0AjimljiulSgD+zilDCACjXARGAsjle874AQDPAdjk/Cz6wgijXHzTNYvAFTiCNzrmLwLY7Pz8KIACgCkA/wrgT5RSekKsAPxvIjpIRPd1srI9xEpyeRHAzxNRioi2ArjZ+d1GAJOuz0867wm8+JWL5ktO6MHvE0mYLhdEtAXAWwH8C4B1Sqkp51evA1jn/LwRwCnXx7RurPS+EJKQcgGAfiI6QETPEdFd5mvcG3iUy0qIvhgipFwAmZMZIYBcPoJ6dAMg+mKMkHIBfOpLKlx1redXAfw5Ef0+gMdQ9/gB9V2MKoANqIdK/TMRPa2UOg7g7Uqp01Q/2/RtInrZ8ZoIfKwkly8C2AngAICTAL6HupyEzhBELr/s6MsqAF8DcA/qnlohBEQ0iHp7fkwpNedeWyulFBHJ3T4RwCSXqx2d2QbgH4loXCn1iqEq9wSiL3bCJBeZkzHjVy5E9NOoLzbe3tGK9hhMcvGlL13tCVRKvayUepdS6mYAj6AewwzUzwQ+pZQqO+Ft/xdOeJtS6rTz71kA/4D6glFgZCW5KKUqSqmPK6VuVErdCWA16rHRp3Gpu3uT857ASAC5uPXlIoD/CdGX0BBRGvWB4G+VUl933j6jwwmdf88675/GpV5ZrRsrvS8EhEkubp05jvqZ2rcar3wX41MuKyH6wgyTXGROxoxfuRDRWwB8AcCdSqlp523RF2aY5OJbX7p6EeishEFECQC/B+Dzzq/+FcDPOL/LoX5w/2UiyjkeDf3+uwAcbixXCMdKciGiAafdQUQ/C6CilHrJcYfPEdGtTrjhhwB8M5rady9+5eKEh44576cB/BuIvoTC6d9/A2BCKfWnrl89BkBn+LwXb/T/xwB8iOrcCiDv6Mu3ALyLiEacjGLvct4TAsAlF0cefU6ZYwB+EsBLHXmILiSAXFbiBQBvJqKtRJQB8AGnDCEAXHKRORkvfuVCRFcB+DqAe5RSR11/L/rCCJdcAumLsiAzDscLdc/FFOoHiydRd5Hej7rH4iiAzwAg528HAXwV9TNQLwH4hPP+NtTPP73o/O5TUT9X3F8+5bIFwA9RPxT7NOphU7qc3U5nfgXAX+rPyCs6uQDIoZ4p9AeOvjwIIBn1s8X5hXpYh3La9JDzejeANQCeAfAjRwajzt8TgL9y9GIcrkytqIf3HnNeH4762eL84pIL6lndxp0xZhzAR6J+tji/AsjlSsfezaGe4GoS9aRjcD531JGZjP0WyAUyJ4taLl8AMOv62wOuskRfLJNLEH3RkzxBEARBEARBEAShB+jqcFBBEARBEARBEAThUmQRKAiCIAiCIAiC0EPIIlAQBEEQBEEQBKGHkEWgIAiCIAiCIAhCDyGLQEEQBEEQBEEQhB5CFoGCIAiCIAiCIAg9hCwCBUEQBEEQBEEQeoj/D7RG+2jZ91g8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x122.4 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load the dataset\n",
    "url = 'https://www.cl.cam.ac.uk/teaching/2021/DataSci/data/climate.csv'\n",
    "climate = pandas.read_csv(url)\n",
    "climate = climate.loc[(climate.station=='Cambridge') & (climate.yyyy>=1985)].copy()\n",
    "t = climate.yyyy + (climate.mm-1)/12\n",
    "temp = (climate.tmin + climate.tmax)/2\n",
    "\n",
    "fig,ax = plt.subplots(figsize=(15,1.7))\n",
    "ax.scatter(t, temp, s=15, alpha=0.7)\n",
    "ax.set_xlim([1984, 2025])\n",
    "ax.set_ylabel('mean temp')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* **Step 0.** Write out the likelihood function for the observed data, given all the unknown parameters. I suggest\n",
    "the probabilistic model\n",
    "$$\n",
    "\\texttt{temp} = \\alpha + \\beta_1 \\sin(2\\pi\\texttt{t}) + \\beta_2 \\cos(2\\pi\\texttt{t}) + \\gamma(\\texttt{t}-2000) + N(0,\\sigma^2).\n",
    "$$\n",
    "This is like the linear model from section 2.2 of the notes, but with an explicit probabilistic noise term in the style of section 2.4.\n",
    "\n",
    "* **Step 1.** Invent a prior distribution for the unknown parameters. Take $m$ samples from your prior distribution, to create a list\n",
    "of samples $[(\\alpha_1,\\beta1_1,\\beta2_1,\\gamma_1),\\dots,(\\alpha_m,\\beta1_m,\\beta2_m,\\gamma_m)]$. I'll assume you know $\\sigma$; if you are uncertain about $\\sigma$ then you should include it in the unknown parameters.\n",
    "\n",
    "* **Step 2.** Define a weight $w_i$ for each sample $(\\alpha_i,\\beta1_i,\\beta2_i,\\gamma_i)$, equal to the likelihood function for the observed data. Then rescale your weights so they sum to one.\n",
    "\n",
    "* **Step 3.** Plot a posterior histogram for γ, and then find a confidence interval for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0.\n",
    "# Define a function logpr(temp, α,β1,β2,γ) which returns the log likelihood of the observed data temp.\n",
    "# Use the probability model in the equation above. Just pluck a value for σ out of thin air.\n",
    "# Also, use scipy.stats.norm.logpdf, rather than typing out the formula yourself.\n",
    "\n",
    "def logpr(temp, α,β1,β2,γ):\n",
    "    # temp is a vector, the parameters are scalars, and the function must return a scalar\n",
    "    t = climate['t']\n",
    "    return ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1.\n",
    "# Instead of creating a list of 4-tuples, let's create 4 lists of values.\n",
    "\n",
    "m = 20000\n",
    "α_ = ???\n",
    "β1_ = ???\n",
    "β2_ = ???\n",
    "γ_ = ???"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2.\n",
    "# This code uses a trick to avoid underflow. See section 7.1 of notes.\n",
    "\n",
    "logw = np.array([logpr(climate['temp'], α,β1,β2,γ) for (α,β1,β2,γ) in zip(α_,β1_,β2_,γ_)])\n",
    "logw = logw - max(logw)\n",
    "w = np.exp(logw)\n",
    "w = w / np.sum(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 3.\n",
    "\n",
    "# Plot a histogram of the posterior density.\n",
    "fig,ax = plt.subplots()\n",
    "ax.hist(γ_*100, weights=w, density=True, bins=np.linspace(0,8,60))\n",
    "\n",
    "# Find a posterior 95% confidence interval using γ_ and w\n",
    "# (see section 7.2 of notes)\n",
    "lo,hi = ???,???\n",
    "\n",
    "ax.axvline(x=lo*100, linestyle='--', color='black')\n",
    "ax.axvline(x=hi*100, linestyle='--', color='black')\n",
    "ax.set_title(f\"Temperature increase ≈ [{lo*100:.3},{hi*100:.3}] °C/century\")\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}