{ "metadata": { "name": "", "signature": "sha256:133e791b47e184498d5bea51d0829c1dafa7ae5909abe09c4f170c9084951099" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Frequentist Model Fitting Breakout\n", "\n", "In this session, we're going to fit some **Generalized Linear Models** to astronomical data.\n", "\n", "As usual, we'll start with the standard imports and IPython notebook setup:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# use seaborn plotting defaults\n", "# If this causes an error, you can comment it out.\n", "import seaborn as sns\n", "sns.set()" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part I: Fourier Fit to RR Lyrae\n", "\n", "We'll start by doing a multi-term Fourier fit to an RR Lyrae star.\n", "Note that downloading the data will require installation of [astroML](http://astroml.org), which can be done easily by running\n", "\n", "```\n", "$ pip install astroML\n", "```\n", "\n", "at the command-line." ] }, { "cell_type": "code", "collapsed": false, "input": [ "from fig_code import sample_light_curve\n", "t, y, dy = sample_light_curve()" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[10003298 10004892 10013411 ..., 9984569 9987252 999528]\n" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize the data:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.errorbar(t, y, dy, fmt='o')\n", "plt.gca().invert_yaxis();" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAFVCAYAAAD7Sga4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+QHOV95/HP6AdahEbL3nqTtUNOOJH9FCKS2cOOQNgY\nZBfkEhwHFHOKkwtEQIwF0nHyGvmHwGfAdaK0yAQbg4MgUnwp76EEO5ekAq47ISugGOOrhfiC7zEy\nhourRCwrirQSrGC1c39Mz27vTM90z0zPdD/d71cVxezM7Op5pqf7+/z4Pk8XSqWSAACAu+YkXQAA\nANAegjkAAI4jmAMA4DiCOQAAjiOYAwDgOII5AACOCw3mxpiVxpgnq577qDFmf4Pf+TljzD8ZY94Z\nRyEBAEB98xq9aIy5VdLvSTrue25I0roGvzNf0lclnYipjAAAoIGwnvkBSVdJKkiSMaZf0hck3VJ5\nLsA2SQ9IOhhTGQEAQAMNg7m19jFJk5JkjJkj6WFJm+TrqfsZY66VdMha+y3vqXoBHwAAxKQQtp2r\nMeZsSV+XtFHSn0g6JKlH0jJJD1trN/ne+21JJe+/8yRZSR+21v5zvb9fKpVKhQIxHwCQK7EGvsjB\n3Fp7oe+5JZJG/c8F/N6Tkj5mrf1hSBlKhw6NRy+xYwYGiqJ+7spy/bJcN4n6uS4H9Ys1mEddmlYd\n8Qv+54wxu4wxvxhbqQAAQGQNs9klyVr7sqRVjZ6z1l4T8HuXtl88AAAQhk1jAABwHMEcAADHEcwB\nAHAcwRwAAMcRzAEAcBzBHAAAxxHMAQBwHMEcAADHEcwBAHAcwRwAAMcRzAEAcBzBHAAAxxHMAQBw\nHMEcAADHEcwBAHAcwRwAAMcRzAEAcBzBHAAAxxHMAQBwHMEcAADHEcwBAHAcwRwAAMcRzAEAcBzB\nHAAAxxHMAQBwHMEcAADHEcwBAHAcwRwAAMcRzAEAcBzBHAAAxxHMAQBwHMEcAADHEcwBAHAcwRwA\nAMcRzAEAcBzBHAAAxxHMAQBwHMEcAADHEcwBAHAcwRwAAMcRzAEAcBzBHEBmfPIr+/XJr+xPuhhA\n1xHMAQBw3LywNxhjVkraaq291PfcRyXdbK1dFfD+T0v6kKT5kr5srd0VY3kTMTI6phdePiJJWnZ2\nn4bXDiVcIgAAZjTsmRtjbpX0kKQFvueGJK2r8/5LJF3oBflLJP1SXAVNij+QS9ILLx/RJ+5/Wq+8\nOp5gqQBI0vV379H1d+9JuhhA4sKG2Q9IukpSQZKMMf2SviDplspzVS6T9H1jzDcl/ZWk/xFfUZPx\nA18grzgyflL3/cU/JFAaABUb7t2nqZI0VSo/BvKsYTC31j4maVKSjDFzJD0saZOk43V+ZUDS+ZJ+\nW9KNkv4stpICgGfDvft0YmJy+ucTE5Nat3WPDh+b0OFjExoZHUuwdED3hc6Z+5wvaamkByT1SFpm\njNlurd3ke8/PJP3AWjsp6YfGmAljzFustT9r9IcHBorNlrtr3vWOAT334qFZz/X39mjLupWRy53m\n+sWB+rnL1br5A3mQF14+omvveEJb1q3U0rPO7FKpus/V4xdV1usXp8jB3Fr7rKRfkSRjzBJJo1WB\nXJKekvSfJG03xrxN0hmSDof97UOH0jH/PDI6Nj2sfo6X6LZxzXJ94v6ndWT8pCSpr7hA2z5ezvuL\nUu6BgWJq6tcJ1M9dWa6bJB0+OqE7dnxH99x0UdJF6YisH7881C9OUYN5qerngv85Y8wuSZ+11v6N\nMeZiY8x3VR7CX2+trf7dVKqX6LZxzQptXLNieo5845oVSRURGcMqiVpRP5MzeuaF9s7TIKiDAHRC\noVRKPNaW0tD6um7rnpoWi1TuibfTss9D6zLv9WslKFc3HqXyd23jmhVaMtidocW0HbtmP5Pr7t6j\nRpev/t4e3Xzl8q59ntU6fYzTdvziloP6BSWRt4xNY4A2tLp0kVUStZr9TG64Ytn04z/80DL1FadX\n0KqvuEA7b788sUAu1a/PnbueZZc6xI5g7jnn7L6a5yqtaKAegnJyLjh3UI98arUe+dRqXXDuoDau\nWaE5BWlOId3TYVMlkXGP2BHMPcNrh2pa9vfcdFGiLXtkF43HWu1+JksGi+or9qiv2JOK8zaoPn5s\nQIU4Ecx9Nq5Zob7igtxfVBFdqwGIxmOtrH0m1fUJwigO4kIw91kyWNQ9N13k9AUE3dVOAKLxWMuV\nofKo/McY6KRmNo0BEKDVpYuVxiNmLBksasfm1UkXIzb+Y9woux1oF8EcaBNBOT22ra+5kWNqDK8d\nqtmAiu8N4sIwOwBnffIr+51a5sXUCjqFnjkAJ42MjunwsYnpxy7srsYoDjqFYI5I2Hq0Pj6b7mu0\n/TLJq8gjhtl9RkbHdN3WPbpu6x42dPBpdZezPOCzSQab9QCzEcw9lYtySeU7yHBRnsGFsz4+GzRj\nZHRM67bu0To6DIgZw+yeRhdl5riAdDnn7D4nlnn575q2sOpOb0wNIE70zBGKrUfra/ezcS0bOy1c\n2C2uerQv6JatjOIgLgRzDwGrPhcunElp57OpZGNz043WpH23uKDRPqBTCOaesIty3ntQ/vWxfYsW\nTCcK3vZgfj+TCn8giRpUSJxrX2W3uB2bVzt7ntJhQFwI5j71NnSgBzWzPvat/Qv10sFj00OHz714\nKPdBaMlgUf2Le9S/OPrdukicy76g0b5CYeYxI1yIE8HcJ+hGK/SgZiMI1aKxhyBBo323X/MedoBD\nRxDMQxC80Eirjb1O5Gi4NsTcLBfrVz3at3vvAR0ZP6kj4ye1e++BpIuHDCGY+7AGNByJgrO12tiL\nO6mQ0YF08o/27d57gFE+dAzB3FOvh/X2ty6ueW+eg1d1EOrv7WHer0Vx3XSDqaCytDdoGOVDJxHM\nPfVOtCPHT7Isq4o/CG1ZtzLp4iSqnZGKoByNVhAkghs0197xRO4aNMgvgnkE3LZwNn8QWnrWmUkX\nJ1Gswe+eI+MTOjI+EfhaUIPm8NGJVDVomKJCJ7Gdq6fR9pDcthCNbFyzQnfuenb6cVRx3W3Nla1N\n82547ZA+cf/TOjJ+UtJMww+IAz1zDz0stGrJYFF9xR71FaOvM49znjsP392R0TFNlaSpkgLnw4N6\nvf29Palr0DDKh04hmPtwoqFb4p7nTvvWpu2I0vAJatDsvP3y1DVo4sqTAKoxzO7DcDpatW39qkT/\n/crWpllUPYUglRs+2x99Tn+08X3Tz7U63QFkAT1zIAEkQ7Xv+Otvzvq5lekOICvomQMJ6EQyVGV3\ntKRHCbpl0enzky5C0+JKegSq0TMHEhJnjkbaN0xpx7I6oxibrj4vgdK0js190EkEcyAhcSVDZT1I\nDK8dUnHhTC+8uHC+kwlkbO6DTiKYA47LQ5Dw98Jd65ED3cCcOYDERJ3nr9wzvvK4njTnC7C5DzqJ\nnjngOFcz47M8zx8kD5v7IDkEc8BxLgaJrM/z15PlzX2QLIbZgRgkvSzMtQ1TGm0Ec9q8uZLSPWTe\nqixv7oNkEcwRi5HRselErHNytn62MlxceZxE3bMSJMZfe1NSeTOYpD5LwEUMs6NtlSHTkqSS8jNk\nKrU/XPzJr+yf7tVjturPctv6VZnsrQNxIJijbXlYGlVPO3XPWwKYX9BGMEHy8j0C2kUwBxKQ1wSw\niqCNYAoJlgdwHcEcbXN1aVQcWq17nkczKjZdfd50Zvemq8/L9fcIaBfBHG1zcWlUXPJc93ZV3+WM\nzxJoHcEckYyMjmnd1j1at3VP4PxunDcNcU0ra4fphZZVJ7Xl+XsEtKNQKpWSLkPp0KHszhMODBTl\nev2q53elmcDz7uVvc75+jXTy+MV9C9RmZeG72Qj1c1sO6hdrmkhoz9wYs9IY82TVcx81xtSspzHG\nzDHGPGKMecoYs88YY+IsLJLB/G5nxN0LZZkbkF8Ng7kx5lZJD0la4HtuSNK6Or9ymaQzrLXvlXSH\npC/EVE4gc+K6BaqU72VuAMJ75gckXSWVV40YY/pVDtC3VJ6r8rqkXmNMQVKvpDfiKyqSEmV+N2xO\nHZ3j8jI3RhOAeITOmRtjzpb0dUkXSXpM0qckTUj6urX2wqr3zpP0PyW9VVK/pA9Za/8+pAyJT9oj\n3LV3PKHDR8tblvb39mjn7ZdPv3bbg/v13IuHZr2/v7dHW9at1NKzzuxqOfPoN4f/UkGncfVxShv/\n9+a8dwzozhvZ3Q25EuuceTN7s58vaamkByT1SFpmjNlurd3ke8+tkp621n7WGHOWpD3GmF+x1jbs\noWc8ySET9bv5yuXTc+Q3X7l8uk4DA0U9XxXIJenw0QndseM7XU/qipsTx69Oc3hqqtSw7EnWrXo0\n4bkXD+n3/8vj2rhmRWxL0Zw4dm2gfm4bGIh3yWXkYG6tfVbSr0iSMWaJpNGqQC5JZ0g65j0+Imm+\npLkxlBMJq8zvIn0W9szTiYnJWc8VCtJHLvnlhEoUrlFSJd8zoHlRg3l127/gf84Ys0vSZyVtk/Qn\nxpi/UzmQf9pa+3ocBUV6nXN2X92la+i816oCuSSVStLuvT/SBecOJlAiVOT5boLortBgbq19WdKq\nRs9Za6/xvXxlTGVLlaTvV51mw2uHEl8zDbfkoQFYLzExzqkEoIId4CJg2U84/0U4SxfkdnUjW9vF\n3eTysHUr+zOgmwjmIVxe9tNNu/ceCHycZ91qBLoaGFvZBhdAMIJ5CFrX4Wjw1Or2Z+JiYFwyWNSO\nzau1Y/Pq1Dc8WuHiiAncRTBH22jw1Or2Z5L1wOgiV0dM4CaCeQha10DnZH0HOO4C1xlZ/960gmAe\ngtZ1OBo8tfhMwuUhsTTO/feBRgjmEdC6bowGT63htUOaN3dmt8Z5cwu5/0z8yLMA4kUwj4DWdTga\nPLONjI5p8tTMXkuTp0oEKx/yLIB4EcwRiyWDRc0pFDSnUKDBI4IV0Cl5mJ5pBcEcseAEQzPIKUAr\nmJ6pj2COtnGC1SJYNUaeBVrBiFd9BHO0jROsFsEqnIsb3QBp1cz9zDPP38Ncxh2O0KaNa1bozl3P\nTj/uNNduBuRv2NDIQRR5uEFPq+iZexgqbh1DysG6uSubizkLI6NjmipJUyU5U2YkixGv+gjmnjiH\nivO2OxEnWLJcbIh2qsx5O/fyiDs0BiOYIxasM0+OizkLLpYZ6cAdGoMxZ+5hLqY9lY11kG+uzdvD\nLfVGdDauWZH7kUB65h6GiuGqtOQsNDNvn5Yywy2M6NRHMPdhqBguSkNDtNk58DSUGcgSgrkPe7C3\nh+Sj5CS9ZruVHlPSZYZ7GNGpj2AeMxeXCMUhr/UO060GzpLBovqKPeor9jjTEK1X5lY/M76DyLNC\nqVQKf1dnlQ4dSu8SmmZUDzVKUn9vj26+crkzF9hmDQwUtflL++omD7pe74GBolr9fqZ9E6J26lYt\n6Lvfyneg1c8sr+deVq6dQYLqd93WPQqKWJVpGpcMDBQL4e+Kjp55jIKGGg8fnXAyOaOZ3hFJKbVc\nXPvdjjjmwNv5zLJ07gGtIJgDHZDHBk7fogWBj6PK42eG5jBnXh/BPEZBX7T+3p7Mf9E4wWrVm7ya\nPDXV1XJ0y8jomF46eGz655cOHuvqSERez728YRVEfQTzGAV90Xbefnnmv2icYKier5bKvertjz4X\n+W+00yjM67mXR2znGoxgHrO8ftFYZjRbvcyWeXPzdcodf/3NyO9tt1HIdzAflgwW1b+4R/2L3Vm5\n0Q35urJ0gX9JjP+xS2uwW1ni4+LSqE5i6qFs0enzm3p/OwGZ7yDyjGAeow337tOJicnpn09MTOrD\nw3+p7/zjqwmWqjntZBRvW7+KPbk9eZt6WFan8bLp6vOa+jsEZKA1BPMY+QN5xVRJ+uO/esGZzSzI\nKI5PnoZ9h9cOqbhwphdeXDi/5cYLjUKE4TtSi2DeZVlfb4wZeetlbrr6vOnGS7M9cqAZLk1bdgvB\nvE3+L9UZPdHuKJvmni5zvfHJ2/aiSwaL2rF5tXZsXp2LxguSc2R8QkfGJ5IuRqoQzGP0pVsuVsGX\nxlwoaNbPLsjbXG+n5G0HuLRg+BV5RTBvQ1DP64Yrlk2/fsMVy/SupQM1v5f2nm5el9fFidwDAN1E\nMG9RvZ7XW/vPmJ43vODcQd154yrnerq79x4IfAwASCeCeYua6Xm5lNXM8HA8yD0AOmPDvfs0VSqv\nFNpw776ki5MaBPMucCkxqJXhYTJLa5F7AMQvaC+P6+7e49ReHp1CMI9gZHRM67bu0bqte6bnxul5\nld324P5cZWw3w6URGcAFQXt5lErSQ3/9QgKlSReCeYh6w84fuWRp3Z5XZQjIRc00UkZGx/Tci4em\nf2ZIfjaXRmTiwigNkAyCeYhGw871sr4rNwFwUaPh4eoLNRnbQHxoCIUL2sujUJi9iiivCOZtWDJY\nnB5GzVLPK2h4OG8boABIn6C9PB7evFoXnDuYXKFSgmAeotGw88jo2PSQepYCXPU2pPWmGnpOm1vz\nu4WC9JFLfrmbxUVK0OBrHZ9ddNV7eaCMYB6i3rDz7r0HcrOEq95w+utvnKp5vlSSdu/9UTeKhRRh\nSWPr+Oyac8G5g7P28kBZaDA3xqw0xjzpPR4yxvzEGPOk99/VVe+dY4x50Biz33s9E120oGFn5ouB\nGZwPreOzQxwa3hnEGHOrpN+TdNx76nxJ26212+v8ym9JOs1au8oYs1LSPd5zTqtkJUcx/tobmjxV\nTmUfGR3T3Rsu7mTRuuKcs/tm9Ryk8ghF36IFeungsZrnWYYFAN0V1jM/IOkqSZWUg/Ml/YYx5tvG\nmB3GmEVV779I0uOSZK19RtK74yxsUqKuM583tzAdyKXycNm1dzzh5HCZ/4YV9aYatlzzbvX39tQ8\nn6VkQETDvgut47NrXiWnBzMKpVLjBdHGmLMlfd1ae6Ex5lpJz1trx4wxn5HUZ639pO+9D0n6C2vt\n497Pr0h6u7V2qsE/keoV2bc9uH/WWmpJ6u/t0ZZ1K3XXI8/o8NGJ6ef+5diEgj7O/t4e7bz98m4U\nt2MO/ORfddcjz0iStqxbqaVnnTn9/Cfu/bYk6Z5b3j/9PPLn2juemHU+uP6d7yY+u1yK9Z6a0W7A\nPeMb1tqj3uNvSrqv6vVjkvzdsjkhgVySdOhQenuuz1cFckk6fHRCd+z4jjauWaE7dz0rSbr5yuW6\nY+ezgX9jaqqU6jpG0btgrrZ9fObWkpX6LD3rzFlTEK7Xs9rAQDFzdaqIu269C0+bDki9C09L/HNz\n6djdfOXyWdeSKOV2qX6tyEP94tRsNvvjxpj3eI8/IOl7Va8/LenXJckYc4Ek5zM46g0bTJ6aqlnC\nFTRc1t/bw3AZMm9kdGxW/sRLB4+Rka3gKbog1dcSoFlRg3klpt0o6YtedvuFku6SJGPMLmPMWZK+\nIWnCGPO0yslv/znm8qZO2Nzyztsvz/zJyc5VICO7VrNLzvzXEqBZocPs1tqXJa3yHj8v6b0B77nG\n9+PH4ypcGhQU3DufN7fcDqoEscpJuHHNCn3eG27PQ4+8cqMVqXzxGl47lHCJgHRo1MC556aLal6r\nvpYAzWDTmBAL6+wFXG+Xs6xs8Rqlt82NVlBBRjaQLIJ5iNfq3HJv994f5X4LRoZWUcH922s1ewfC\nPF9L0D6CeYvGX3uDLRgBH+7fPlvUBg7buSIOBPMQ9VrX/s1hKo6Mn9QXvva9TN58JUizUxDItjze\nvz3MxjUryrslNphyYIQLcSCYhxheO6R5c2fW9s+bW9A9N11Ud7V/VnaAi6LRFASAcgPnnpsuyv2U\nAzqPYB5iZHRsVoCePFXSJ+5/Wm9/6+JIv3/46IRzLWzm74DuIXkQcSCYh6g3BHbk+Mma+bBY9+ZL\nSDPzd1yEgPaRPIg4EMzb4A9aG9esyMQOcM3M3w2vHeJGK0AMzlx0WuBjICqCeYhGvc/qNeV53AFu\ny7qVZDADbRgZHdOPD86MfP344DjZ7GgawTxEoyGwkdGxmsx115fnNDt0XrnRChnMQGteqDMatv3R\n5xIoTbqxdXR9BPMI+hYtqHlcb25ZktM3TGD+DkiH46+/mXQRUufI+ISOjE8kXYxUIpiHqHc3qHqt\nadcy14M0O7pAaxmI36LT5yddBDiEYB6iXkJYlvl74WE98sqNVljGBrRmWZ2prU1Xn5dAaeAqgnmL\n/BvJVORtWRY3WgHaN7x2SMWFM73w4sL5TG2haQTzEPUSwj77H9+d+7lltqEE4uHvhdMjDxaUcIwZ\nBPMQjRLCghLjsoCTBuiurNw6uVO4GU04gnkEQQlh9RLjXnl1XNvWr9K29auSKm5b2AEOSEZliSdq\nsXwvHME8gqC7QWV1iJkd4BC3Zlc75HV1RF7r3Q6W780gmKMtW9atDL3FI/KLm/agk1i+N4Ng3qKs\nDjG3sgMct3hEkFbmOQn+CMLyvXAE8xZldae0rNYL3dfsPGeek5xoxDQ2vHZIc+fMLAeeO6fAdakK\nwbwNru/DXk9W64V0qDfPmdU8lDB5bsRENTI6plNTpemfT02V+IyqEMzbsGSw6PQ+7PVktV5IB+Y5\nZ8trI6YZfEbh5iVdANe5ugQN6LRlZ/fVDLU3yr84p8n3A5hBzxyBXF4rj3RodpvSvOZrZDWZNk58\nRuEI5gA6ZtPV503nX0TJPM5jvkZeGzHN4DMKRzDvEjaEQB41cwe+ynvymK+Rx0ZMs/yfC59RLebM\nAaRKHqd3Ko2YymPUquxfX3mM2eiZdwFrSJFXG+7dN33Tng337ku6OHAYN4BqrFAqlcLf1VmlQ4ey\nu1bwvr/4/qx7fksziRtZaF0ODBSV5eOX5fp1um4b7t2nExOTs54rFKQbrlimC84d7Ni/W5HlYyfl\nq37Va/El96+jAwPFQvi7oqNn3mHPHzhU8xzrI5EH1YFckkol6aG/fiGB0sBlrDMPRzAHAMBxBPMO\ne9fSgZrnWB+JPDijpza/tjLMDjSDdebhCOYddueNq1gfiVz60i0Xq+CbFSwUpIc3r+7KfDmyhXXm\n4QjmVTqxHpw1pMgrfy+cHjnawTrzxgjmPp1aQpbXjTCAC84dVP/iHvUv7qFHjrbs3nsg8DHKWJrm\n6dTShzwtH8miLNcvy3WTqJ/rWJrWHHrmHpY+AEA6cX0ORzAHAMBxBHMPSx+AbLnurm9xc6OM4Poc\njmDuYelDvLhLHIC4cH0ORzD3YQkZED8adogD1+fGCOY+SwaL2rF5tXZsXk2Lrw3cJQ5JGxkd00+P\nvM53MENY4ttYaDA3xqw0xjzpPR4yxvzEGPOk99/VVe+db4z5mjFmnzHmGWPMhzpVcKRT9RKSF14+\nok/c/7ReeTW7S2iQLnwHkUcNg7kx5lZJD0mqTFacL2m7tfZS779Hq37ldyUdstZeLOnXJH057gIj\n3VhCAr8kRmn4DmYHUzTRhfXMD0i6SlJlcfv5kn7DGPNtY8wOY8yiqvfvlnS772/X3gMRQC7QQwa6\np2Ewt9Y+ptkB+RlJw9ba90t6SdLnqt5/wlp73BhTVDmwfzbm8iLlWEKCiqR6yHwHs2vb+lXatn5V\n0sVIpdp7FDb2DWvtUe/xNyXdV/0GY8wvSnpM0v3W2tEof3RgINvJDHmq390bLta1dzyhw0cnJEn9\nvT3aefvlsf+b1931LUnSw1sui/1vV8vy8eto3QqSAnaLnjOn0NF/t1vfwTTI8ndTkubOLQ8KZ72e\ncWg2mD9ujNlorX1W0gckfc//ojHm5yV9S9J6a+2TUf9oXvYXdl1l7srfMg6q381XLtfndz47/bgT\n9T91qhwlOv3ZZun4Vet03c5Z0he4n3anvhN+vQtPmw7mvQtPy+QxzPJ3UyrXr1vneRLibqBEXZpW\naV/fKOmLXnb7hZLukiRjzC6vR/5pSb2SbvdlvPfEWmIkoplEpiWDxek7ZXViCQlL39yQ1EYfI6Nj\neungsemfXzp4jLl6B9324H7O8yZw17QOy0LrudEdi969/G019fO/f9nZfRpeO9SVsnQiSGTh+NXT\njbq98uq47txVHqW57Zr3dGV98HVb9wSN7k83JrIiy9/NLN4lrRp3TUPXNZPI1OkMZpYduYWNPtAK\nzvPmEcwRK05CVOt2BvLCntpUoEJB+sglv9y1MgDdRjBHqDQt9UlTWZBOr03Ubm9RKkm79/4ogdKg\nFZznzSOYI1QziUydPgm5exKQfcNrhzR/3kx4mje3wHkegmCOSNJ0xyL/v590WZA+9OrcNzI6pjcn\np6Z/njxVYkVCCII5IomayFSdgSqV58y3P/pcrGUJegxIjN5kAbk3zSOYI7J2EpmOv/5mbOXwrzll\n/SmCpGkkCegGgjm6YtHp82P5O9y8A1EsGSzqLWeezpI4RzFV0jyCOWK1rM5JuOnq82L5+90YxgcA\n1xDMEavhtUMqLpzphRcXzu/KfGWcw/idxj2au+NnRyd0ZHwi6WKgBcyZN49gjtj5e+Fx9cjDxDWM\nj+x4S285YRPIA4I5Yrd774HAx3Ho9DA+gOSxi1/zCOaIVacT1JIaxo8Ld3zrjpHRMf30yOt8zo5i\nF7/mEcwRq27Mdb2ltyfwcdqRid8dfM7II4I5nDIyOqYfH5y5KP/44LgzF2qSerqDz9l9LE1rHsG8\nBdffvUfX370n6WKkUtSTsNWMbpam5QMZ//k2vHZI/b5Rt7Bd/Pi+EMwRs6S20nRhaRq9jeiOjLe+\nrIzPORu2rFs5/ZhjF45g3qQN9+7TVEmaKpUfo1YSW2myNA3IlqVnnTl9HXElwTVJBPMmbLh3n074\nsixPTEzqurv36Dv/+GqCpUqfJYNF7di8Wjs2r479JHR5aRpzud3B54w8IphH9Mmv7J8VyCtKJemh\nv34hgRK5q5XlWZU5MdeXpiHcyOjY9OgXy8oQRTvTMllBMI9BqZR0CdwRx7KhJHaYi0OeNsJoNSEp\nju8Hc+bII4J5BJWeZD1nBFykESyOIdAlg0X1L+5R/2K37ohVb2Tn6//rxQRKk05xrFbgfubZ0VcM\n35KXkZwygnmI6p5CkDlzCg1fZ9kEGnEhE78ZndjlrtnPiPuZ5wMbBM2gSxkiqCeJ1i3smVfTQ21l\nqHnb+lUrcRdSAAAPN0lEQVRxFitRWcrEr3dx3bhmRVs949MXNHepqtzP/NSpEj1yh4Wd541G+u65\n6aJOFSuV6Jl3AckZM/K857LLmfhRdSqTfG7I6BeQdwTzEEHJNEAryMQPVy9kz5vLpQq18pRUGoYz\nJER1Mk2QKBeaqZKYN1drJ1+W7jTmaiZ+VO1u5xtnJvrDWy7L1HQMauV5pK8awTyCSjJNkLALzW0P\n7teUt3Tt2Ik3OlA6tzSb0Z21BBd/LzyLPfJ2M8nJRAdaQzCPYMlgUX3F8lKoZi40I6Njeu7FQ9M/\nv3lqKrWBKOmM+3rZylnbzcs/suD6KEM9/sZtKz3qdn8f+cGeAjMI5k1qZslL1gIR2pO1UYZ6du89\nEPg4qqyPXiA+w2uHNG/uzLDpvLmF3I7kEMwj2rZ+lbatXzXdS+8rhm9YUm9juH89fjJzF/AoGvVE\n6y3PylLLOw+NuzgaLHkYvUA8RkbHNHlq5ko7eaqUyQZyFATzFlQCe6tKJaXuAt7p5XONNt9ptDyL\nOVS31NvBbVvEoJyX0QvEIw8N5KgI5mi4HWJcc+mNNt8JC85Z2c0rz8to/FnHjVYnxLGdK5BHBPME\nFFIUlLrVE6o35RAU4Ko1M7WRZnnem73gTWu2+n3L2pa3iEeWpuHaRTDvoHobYJy5aEFqglLSw1Ts\n6xU9UCW94iCKoJsOFQrSDVcsk9T69y1LW94iPkzDzSCYd1BQq7G4cH6qWo31esyTp6a68u+/drK2\nt5o3ze47nmZfuuXi6V64VA7kD29erQvOHYz0+3nY8hbxYiljGcG8g4bXDqm/d+b2fXMK0h9tfJ8T\nrcapqVJXdl6jx5W9fccrvfDqx1L4sOjw2qFZjQG2vEUYljKWEcw7bMu6ldOPiwtPS7AkzTkxMRnr\nXHq7Pa4s3KymnX3HXdrS1t8Lr+6RRxkWXew7T7LeI3dh6gRuIJh32NKzzpzOxP7ihvcmXZy25HXJ\nR1xazWZ3cblW/+LyjolBwoZFK42bOYV897QQDfsSlBHMu6CSiZ1G9RKW4tbukqM0f4ZRtXpTiKST\nFJsVNorQaFi08ruSpu9pkGWVzwmtcbGh2ykE85yrF2D8WyRWdGLJB0uOsiXKxbVeTypoY6G8XpgR\njWsN3U4imHdBuzvGJaG48LSuLPnIUwJcq2tiXVpLG3ZxbRTsuTADrSOY51yjQNG3yBfMFzW+p3sY\nlhy1viY2S2tp2eFtBnO97XOpodtpocHcGLPSGPOk93jIGPMTY8yT3n9X1/mdnzPG/JMx5p1xFxjx\nqhcodu89oJcOHpt+/qWDx9oa8hxeO6TiwpleeF6XHLW6JtaVtbStXlyPv/5mri7MzPUibg2DuTHm\nVkkPSapc7c+XtN1ae6n336MBvzNf0lclnYi7sOiMoEDRiSFPfy88Tz1yv1bXxLZ7W9FuaXUUYdHp\n8zM1AhGGEYp4MDUzI6xnfkDSVZpZInu+pN8wxnzbGLPDGLMo4He2SXpA0sH4iolO6lagWDJYnF6y\n1MwF2qU11mFaGVp1rRfX6MY49VZP/M4H3hH6u3lAQiha1TCYW2sfk+RPd35G0rC19v2SXpL0Of/7\njTHXSjpkrf2W91S2trbKoHqB4u1vXVzz3jiGPJtNBnQtkDXSal1c6300ujFO2PK8rNxUB92Rp6mZ\nMM1uCv0Na+1R7/E3Jd1X9fofSCoZYz4o6TxJu4wxH7bW/nOjPzowkO2TNs31+8ErwYFizpyC+nt7\ndPhoeQ1sf2+Pdt5+eeDf6GT96pXvy9/4ft3yxC2u+rVcl4ICN9GfM6fQdtk6dex2fi64PvWWjk+V\nStNlmesti4yjbGk+94IsPuO0psrsWv2aFVa/uzdcrGvveCLSdSrrmg3mjxtjNlprn5X0AUnf87/o\n9dglSV7S3MfCArkkHTrkXi8rqoGBYrrrV+fqOjVV0sY1K/T5nc9Kkm6+cnlgPTpevwbl68bnGmv9\nWqzLOUv6auZY+4oL6h6TqNL03fR/Bls/dqGk9q8LaapftWVnBx/TjWtWRC5zmusXh6j1u/nK5aHX\nqTSKuyEWdWla5TJ0o6QveoH6Qkl3SZIxZpcx5hdjLRm6otEwVfVcehL7SGdpGK3VumQpMayd/emB\nIEsGi9N5Fi6eE3EJPYOstS9ba1d5j5+31r7Xy2T/qLX2uPf8Ndbaf6r6vUuttT/sTLE7J283Pmi0\nNK16fvfI+ETXbo0aVj4XT9p26uLK0rQwWWqctcO1PIi0y8J2z+2iOYzADOKgi81USTp64o2uZ5Rn\nKcO51bpkpfeRpcYZkCYEc58sLYFqRlAGcaN7XHQ7ozxLGc5ZqkursjLK0A5GKBA3grknS0uguoEh\nwe4aGR3TVKk8OuJ6Q7PVjXMA1Ecw9+R9Dqt6/TcbBKQHDc3syfv1Jm4u3swqbgTzLnEtsS5oGNCP\nIcHuyeKFvzL/DyAeBHMPc1izVScqFXwXXpKW2kMvAlxvEDeCuYcs21r+C8sNVyzLTEa5a7J44c/7\nUiKuN4gbwdyHLNvZ/JvGPPX9g4lmYee5N8uFP5uytOQSySOY+3Tq7mEuLnkLSrqaKpV081XLEyxV\nftHQzB6WKSJOBHNPpzKGb3twv5OZyFlMunJZVjaNqcjzSAvQCQRzT6eC1/MHDnXk7wJwH40axIVg\njkBZTLoCgKwimHs6FbzetXSgI3+300i6AgB3EMw9w2uHNG/uzGLqeXMLsQSvO29c5WxQJNs2XfK+\nnAtAfQRzz8jomCZPzdxeZPJUKbZENVczkcm2BQA3EMw9nczezlomMgAgXQjmAAA4jmDuIXsbAOAq\ngrmnk9nbLt+LmnWw6eDiLoIAuodg7tOJRDVXd4BDenA/cwBhCOY+nUhUYwc4tIutdQGEIZj7uDwc\nDgDIL4K5p1NDma7uAIf0IDkTQBiCuadTQ5ku7wCHdGBrXQBhCOZdwLaoaBffIQCNEMw9nRzKZFtU\ntIvvEIBGCOYehjIBAK4imPswlAkAcNG8pAuQJpWhzMpjAABcQDCv0qmtS9kSFQDQKQyzAwDgOII5\nAACOI5gDAOA4gjkAAI4jAQ5wBEmUAOqhZw4AgOMI5gAAOI5gDgCA4wjmAAA4jmAOAIDjCOYAADiO\nYA4AgOMI5gAAOC500xhjzEpJW621lxpjhiT9laQXvZcfsNY+WvX+T0v6kKT5kr5srd0Vc5kBAIBP\nw2BujLlV0u9JOu49db6k7dba7XXef4mkC621q4wxZ0i6NcayAgCAAGE98wOSrpL0Ne/n8yW90xjz\nYZV757dYa4/73n+ZpO8bY74pabGkT8ZcXgAAUKXhnLm19jFJk76nnpE0bK19v6SXJH2u6lcGVA74\nvy3pRkl/Fl9RAQBAkGZvtPINa+1R7/E3Jd1X9frPJP3AWjsp6YfGmAljzFustT9r8DcLAwPFJovh\nFurntizXL8t1k6if67Jevzg1m83+uDHmPd7jD0j6XtXrT0n6NUkyxrxN0hmSDrdVQgAA0FDUnnnJ\n+/+Nku43xrwp6aCkP5QkY8wuSZ+11v6NMeZiY8x3VW4orLfWlgL/IgAAiEWhVCLWAgDgMjaNAQDA\ncQRzAAAcRzAHAMBxBHMAABzX7Drzhowxp0naIWmppDclbZR0QtJOSVOS/o+km6y1JWPMH0m6SNK4\n9+u/qfIGNf9N5c1nxiVdY639mTHmAkn3eq9/y1p7R5zljiKobtba573Xvijp/1prv+r9fIPKmf6T\nku7ysvxPV0rr5pW5mfo5deykut/Ngsp7JZySdFLS71trf5qV46f69cvK8XtT0h97b3lR0vXW2lOu\nHb8m65aJY+e7tnxU0s3W2lXez04dO6/MzdSvY8cv7p75DZJe8wp+g6Q/kXSPpM9Yay9W+eLyYe+9\n/07SZdbaS73/xiV9XNLz3nv/VNIW770PSvoda+17Ja00xpwXc7mjqK7bI8aYtxhj/lblG8uUJMkY\nMyhpg6RVki6X9F+9g53mukkR6+dx7dhJwd/NL6p8ol0q6TFJm40xP69sHL/A+nnvzcrx+4KkT3ll\nk6QPOXr+Raqb9/8sHLtHJMm7cde6ypscPXZSxPp5Onb84g7myyQ9LknW2h9K+gVJq621+7zX/1bS\nB40xBUnvkPSQMeYpY8wfeK9fVPl97/8fNMYUJZ1mrf2x9/wTkj4Yc7mjCKrbL6i8pe3XVG6oSNKv\nSnraWvumtfaYyvvbr1C66yZFrJ8xZo7cO3ZSbf3eJmmttfYfvNfnS3pd2Tl+gfVz9NyTguu3zlr7\nlHfBH5T0r3Lz+EWqW4bOvV8wxvwblRsstyiD186g+nX6+MUdzJ+TdIVX8AtUHjZY6Hv9uKRelXeG\nu0/S76q8Y9x6Y8xylW/OUtkudtx772JJx3x/o/J8twXV7afW2u9Wva+omTpIwfVIW92k6PVbKPeO\nnRRcvznez6sk3aRyT9ZfD8nt4xdUPxfPPSm4fj3GmH8r6R8l9Uv6B7l5/kWtW1bOvZ+TNCppk2bu\nyCll59yrV7+OHr+4g/kjko4ZY/5O0m9JspL+xfd6UeXW82uS7rPWTtjyXdf2SHqXyoVfXPXeY97j\nisXe891WXbcfanbdKqrLG1SPtNVNil4/F4+dFFy/I8aY/yDpAUm/bq09rGwdv6D6Zer4WWv/n7X2\nHZK+Kmm73Dx+UeuWlWNXkvR2lb+XX5e0zBizXeWA5tqxk6LXr6PHL+5g/quS9lhr3yfpzyW9Kmm/\nMeb93uv/XtI+Se+U9JQxZo4xZr6k90r635KelvTr/vd6cwpvGGN+yRsivMz7G91WXbeD1tqTAe97\nVtL7jDELjDG9ks5ROfEvzXWTotfPyL1jJwV/N39b5R7rJdbal733fVfZOH716pel4/eoMWap9/px\nlRP9XDx+Uevm4nVTqq3ft6217/ByOdZKesFau0nZuXbWq19Hz71Ys9lV7on/d2PMZyRNSLpe5QbD\nQ97czwuS/tyWs9n/VNLfq5z9t9Na+wNjzMuSdnktnJOSPur93crtVOdKesJa+2zM5Y6ium43VL1e\nkiRr7avGmPsk/Z3Kdf+MtfakMeYBpbduUvT6/cDBYyfNrt/rkj6mch1ekfSYMUaS9lprP5+B4xdW\nvywcv+tVHs7caYx5Q+VVM9dba//ZwePXTN1cP3bV15aCsn3t9Nevo9dO9mYHAMBxbBoDAIDjCOYA\nADiOYA4AgOMI5gAAOI5gDgCA4wjmAAA4jmAOAIDj/j+L3T4EIjazSAAAAABJRU5ErkJggg==\n", "text": [ "" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This data has already been phased, so if we just take the fractional part of *t* then we can see the folded light curve" ] }, { "cell_type": "code", "collapsed": false, "input": [ "plt.errorbar(t % 1, y, dy, fmt='o')\n", "plt.gca().invert_yaxis();" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAesAAAFVCAYAAADPM8ekAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+UXGWd5/FPVXenO0lXkib2GGZxExV9DlEiOcACwR+I\nLs7OyDgS5WTcmYEBVCZARkNCRg3MCngOnARUBIQhOGTdWXNkQV1n9qAzE5CVzDJhNwJzcJ8hAh6d\njTNtpk13Ap2ku2r/qLqVW7fuz/p569b7dQ6H6qrqqpvb3fW9z/N8v98nVyqVBAAA0ivf7QMAAADh\nCNYAAKQcwRoAgJQjWAMAkHIEawAAUo5gDQBAykUGa2PMOcaYxz33fcwYsyfke37NGPMzY8xbW3GQ\nAAD0s8GwB40xN0j6PUmHXfetlnRFyPcMSbpf0pEWHSMAAH0tamS9X9IlknKSZIxZKukLkj7l3Odj\nm6SvSjrQomMEAKCvhQZra+2jkmYlyRiTl/SgpI1yjbTdjDGXS5qw1n6/cldQQAcAADHlotqNGmNW\nSPqGpA2S/lzShKQRSSslPWit3eh67g8klSr/nSHJSvqQtfafg16/VCqVcjliOgCgryQKfLGDtbX2\nPNd9yyXtct/n832PS/qktfYfI46hNDExHf+Ikdj4eEGc4/bjPLcf57j9OMedMT5eSBSs45ZueSN6\nzn2fMWanMeYNSd4YAADEEzmy7gBG1m3GlXJncJ7bj3PcfpzjzmjXyBoAAHQJwRoAgJQjWAMAkHIE\nawAAUo5gDQBAyhGsAQBIOYI1AAApR7AGACDlCNYAAKQcwRoAgJQjWAMAkHIEawAAUo5gDQBAyhGs\nAQBIOYI1AAApR7AGACDlCNYAAKQcwRoAgJQjWAMAkHIEawAAUo5gDQBAyhGsAQBIOYI1AAApR7AG\nACDlCNYAAKQcwRoAgJQjWAMAkHIEawAAUo5gDQBAyhGsAQBIOYI1AAApR7AGACDlCNYAAKQcwRoA\ngJQjWAMAkHIEawAAUo5gDQBAyhGsAQBIOYI1AAApR7AGACDlCNYAAKRcqoL15nv3aPO9e7p9GAAA\npEqqgjUAAKg3GPUEY8w5km6z1r7Xdd/HJF1rrV3j8/zPSLpY0pCku621O1t4vAAAdMX2Xfv0wiuT\nkqSVK8a0ad3qjr136MjaGHODpAckDbvuWy3pioDnXyDpvEoQv0DSmxo9MKbEAQBp4Q7UkvTCK5O6\n/p6n9NNfTCd6nUZjW9Q0+H5Jl0jKSZIxZqmkL0j6lHOfx0WSnjfGfFvSdyX998RHpPJJOTg1o4NT\nM9q+a18jLwEAQMv82BWoHZPTR3XXI8/Ffg13bLv4+u/8dZL3Dw3W1tpHJc1KkjEmL+lBSRslHQ74\nlnFJZ0r6iKSrJf1F3ANx/yNacfUCAEBaeEfmkt6f5Psj16xdzpR0qqSvShqRtNIYc6e1dqPrOb+U\n9GNr7aykfzTGzBhjXmet/WXYC9/1yPPef0SNyemjuvtbz+uhmz6Q4HDhNj5e6PYh9AXOc/txjtuP\nc1zvHW8Z149enKi5b+niEW294pxY5+vHPw2OcXHEDtbW2r2S3i5JxpjlknZ5ArUk/VDSH0u60xjz\n65IWSjoY9drPek6An2KxpIkJRteNGB8vcO46gPPcfpzj9uMc+9uw9nRdf89Tmpw+Wr3vXw/N6IFv\nPRcv0azU3PvHDdbet8m57zPG7JT0OWvtXxlj3m2M+XuVp9jXW2ubPMSysdHh6CcBANAmG9au0he+\n/oxm58phraQTS7Ub1q7S8mUnRtjezPHTVoyFziBHyZVKLYmlzSit+9xf6cjMbOQTxwrDdScE0bhS\n7gzOc/txjtuPcxzuytt2+w6SxwrDuuOa8yX5rk9rcCBXDfKO797xIb9EbV9db4py43176gJ1LuDw\nk2beAQDQaX6Z495AnVTXg/Wz++vXq7s/2AcAoN5pK8bq7nNmfdup68E6iU6cEAAAgmxat1pjhRM5\nVM70t3t51i+gDw7EnvH21fVg/Y5Tx2M9L5dT3QkBAKDTNqxdpbHCcOAA0i+g/9nm99bcJ+mfkrxn\n14P1LVev8f4DfI3OH+rA0QAAEG75soLuuOb80AHkhrWrlM9J+ZyqAd0d5CX9dpL3TEU2+DPP/z/d\n9chzOnT4qPK5nGaLtcdEFnhzyO7sDM5z+3GO249z3Bnj44VE8+JJOpi1jXOV4nAXnrvT4QEA6Edd\nnwb3E7UeILErFwCgf6RiZO3lHWl7OZt+OLc7uacoAACdlspg7eXuBrNwZLCmiUpQqzcAAFrF2z60\n04PEVE6Du3nbtvm1JaWzGQCgXbxxKMnWzdt37dOVt+3Wlbft1vZd+xo+htQHa7+2bQAAdIpfHIoz\nSHSCfEm1m37ECfJeqQ/WceRy0kcveHO3DwMAgKqgIP/5h/bq4uu/89dJXiv1wdqvbZtXqSQ9/MRP\nOnA0AIB+06Z+4O9P8uTUB2tv27Yghw4fpZQLANBycfqB+4kz2Iwr9cFaqm3b5tfyJZeTFgz3RGI7\nACClvP073MlhY6PDkf0/vOIONuPoiQi3fFlBO7ZcKEm64rbddY+XStKrR2c1PK8n/jkAgJTx9u+Q\nVJMB/tKBqYZaX29Yu0p3PfKcpl891tSe1j0xso6jWJIOTs00lRoPAOg/fqVZLzSYAe7lNPlqdtet\nnhuKDg3kdXyuGPj4C69M6qrbd2vFskV66cCUpO4UsAMAekOnSoSdUbYkTU4fTbTrVs8F60UL52ly\nekbFkNmEYknVQC3R5QwA0LxmM8DdrbTHxwv/J8n39tQ0uLOmUCxJgwO5RAv3k9NHdec3f9TGowMA\n9KKgrO2cK6M5bgZ4u/RMsPauKTgL9W88Of6JO/za8ZYfFwCgtwUtk5ZK5YDdgprqpvXMNHhQJxip\nfCKd2/mcAqfIR+cPte34AADZs2R0OHQXSDen7Gvb+jUtP46eGVmHcddhFxbM863FHisMa+OlZ3T8\n2AAA6beyPV3KWqZngnVYu7flywoaK4xorDCiV2dm5TewJrkMABCk0S5lDienql0lxD0TrOOcyKkj\nxwLLuthCEwAQZsPaVYm7lEnNbaEZV88Ea6l2utt7Il9/0vzQ+msAAMI4pVVJs74b3UIziZ5JMJNU\nne52bruFFbWnad0BAICkeipYNyKXk28mn3vaYmggr/s3X9DhIwMAZMFpK8bq2pO2epDYU9PgYfwS\n0HI56eMfXFl3v3d94fhcUZ/Y9nhL1xcAAP2h2eS0ODITrP1O1oNbLtS5b1tW91y/KfPZuZJu3rmX\ngA0ASCwsp6oVMhOspeZPVqlE1jgAIDl3CXE7yoR7bs06rDNMWAKadKK7jN/6AgAAzWhH5zJHpkbW\nYdwF61J5IxCvfE4qhm3nBQBAF/RFsPYrWJ8/XDupMDiQU7EkHTpyrC3dZwAAaFRfBGu/hLLpV4+r\nsGBI+coA29nFS2pP9xkAQLZsvndPdXm13TIdrLfv2qcrb9vt2ytckgYH8tqx5ULfx1rdfQYAkB3t\n7gXu1XMJZlGcBX7v1LeXsz7NlDcAIImgXuDt3DAqsyPrsPaj7vXpoICea1OtHACgt3WiF7hXZoN1\nGPf6dJAlo8NsqQkASIXMBuug/a/jKCwYYlQNAPAVFF/aGTcyG6yDerXWV1fXyuekL294F6NqAICv\nTvQC98pssJb8NxL3uyIqLBhSTuVAfeNlZ3f4KAEA3bZ91z5dcdtuXXHb7liJx37xpZ1ypVLXO3aV\nJiY6W898/T1PaXL6qKQTV0RZNj5eUKfPcT/iPLcf57j9+vEc+1UPOUG4XaPl8fFC1ERvjciRtTHm\nHGPM4577PmaMqasEN8bkjTFfM8b80BjzpDHGJDmYTkl6RdTJwncAQGd1I7s7qdBgbYy5QdIDkoZd\n962WdEXAt1wkaaG19p2Sbpb0hRYdZ0stX1bQHdecH2uNodOF7wAAeEWNrPdLukQq52UZY5aqHIA/\n5dzn8ZqkxcaYnKTFko617lBbL2qNwq/w/arbd9e1IWXkDQC9qxvZ3UlFrlkbY1ZI+oak8yU9KulP\nJM1I+oa19jzPcwcl/Y2kkyUtlXSxtfbvIo6hK4vmN963Rz96caLmvqWLR7T1inN06ilLJEm/vek7\n8js9Q4N5PXr7xdWvr7z1+5KkB7de1L4DBgC0zeU3f08HD5V3ZVy6eEQP3fSBdr9lojXrJMF6g6Q/\nlzQhaUTSSkkPWms3up77WZWnwT9njDlF0m5Jb7fWho2wO55gJimwZ7g74Sysr7hz1fXwE/uro++V\nK8a0ad3q9hxwE/oxYaQbOM/txzluv349xz/9xXR1jbqdiWWOpAlmsXuDW2v3Snq7JBljlkva5Q7U\nFQslTVVuT0oakjSQ5IDS5LQVY4HtSCenj+oLX3/Gd7euTvygAQCt4+QypVXcOmvvADPnvs8Ys7My\nkt4m6VxjzP+U9LeSPmOtfa0lR9picdYovIXvXn5tS9OWQQgAaC9nh8crY9ZoN6Iv66wdceqtf/qL\naX3+ob11948Vhqvf6/dYmq7Q+nVaq9M4z+3HOW6/rJ5jJwnY2ZmxVRqt0W55nXWWbVi7SvlcuXNZ\nUNbf3Y8+X32OY6wwrHwup6GB+tOXtgxCAOh37SzB7VSNdub2s05i+bKCdmy5MPBx5wcsSW88uVAt\n2VoyOk8vHyjfHhzIVafD0zaiBoB+1429p9uhr0fWYbw/4JcPTGvx6LBWLFtUDdTSiXVr9r8GgPRp\n98i3UzXaBOsAQT/glw5M+Tyb/a8BoB91agcugnWLbb53j666fXeijmZ0QAOA9ujEyLcTO3D19Zp1\nGL8a67HCsMZGh+tG184PyL3GPXUk1Z1WAaAvbFq3uu07LXaiRpuRdYCgqY2tl53le7+7k5kkHZ8r\n6hPbHq/rIw4A6KxO7z3dDgTrEEE/YL+SL7817tm5km7euZeADQBdlGSnxbQiWIcI+gEvX1bQWGFE\nY4URLV9W0OZ79wT2EC+VFJp1yBacAIAorFk3yOmC416nTior9X8AgPZiZN0EvzZzXmFrJJ3qfAMA\n6G2MrJvgF2zdOtHRrF39bgEA6UGwbhNvv3FvHfW29WsCy8N6NVsRANLGPQO6csWYNq1b3eUjagzT\n4E0IKrb/08vP1o4tF1bXnd1JZO5ksk51vgGAfhSUF9SLFToE6yaEBVunK1nQurbzS/PRC94cufMX\nHc4AILl25wV18rOZYN2kqGL7sHXtyemjeviJn9SUgSVB2RcA9AfWrJvUzjZzYVdslH0BQLgs5QUx\nsm4zv3VtR9gvjXvU7NdnnLIvAAiXpbwggnWbeX9ZHO5fmm3r19SUXnlHzcfnipqcnunJpAgA6KYs\n9AWXCNYdUSyWm5HmJOUikskk/1FzsSR94evPVL/u1IbnANDLstAXXCJYt4V7Cvu6Lz2pQ5Vp7MGB\nvE4KSSaLyiycnStVyw6yNL0DAGkT9Xnc6QRfgnWLeaewj8zMVm8709mzc0Xf73N+8GHc69J+u38B\nAJoTFYi7Ub9NsG6xqBakxZI0/eqxmqs27w8+aAcvL+/uXwCA5sQJxN1I8KV0qwuKJWnqyDEtWjhP\nUnSAd/OOoienw0fiWWm1BwCdEBaI273XQxhG1i0WVqrlFifD28lgdORzqhlFb9+1T8VSOfinZaoG\nALKuGwm+BOsW8yZ+5XLBzy2WpLseec73B++MoN3r0oUF86qPp3WqBgB6WZxA3I0EX4J1G7gD7Mc/\nuFL5kIAt1f/g3SNo97r04MCJHxeBGABaL24g7nSCL8G6DdwB9ty3LdNYYURDA/6nesloebTs/OCl\n8ojbm4U4deRY4jIBarEB9LNGN9qIE4g7neBLsO6AbevX6P7NF/h2Mnv5wLSuv+cpSdJAvvbH4Uxt\nHzp8VMdd5V4vvDKpgYH64XoapmoAoNelsdKGYN1BQVdozvS1X/315PRRzRbri7lm50o16+FpmaoB\ngDRotmmJtw10t1G61UHLlxWUU/w66iij84d05LXjkk4EYmfKZ9v6NdXb39n+IU1MkAEOoD9kcVdC\nRtYdFraOHJQV/qaTF/l+z8ZLz9COLRdqx5YLG/oF7OTG6QDQKVlMwGVk3SZB0yeb1q3W9fc8pcnp\no5JOTF/7PeZkhY8MD3TmoAEAsXVympyRdReErSP7ZYV7N0+Xev8qEQDaJYuVMATrLgjLNFy+rFCX\nFR5Xp3eBAYA0ymIlDME6hY77ZIV7ea8S/RIqnMB9433169IEdgBZlrVKGIJ1DxorDCufy+nuR5+v\n3he2GciPXpyoaUVKz3AAWbd8WaGpBNy0IVinkF+3M/ddS0bnJR4Vu9e4s5gpCQBZRrBOIWfrTLe5\nYnlE/aaTF+nlAydGwM6oeMXJzV85Hjp8tOnXAAC0HsG6S8K64wTdPzl9VC8dmPK9/1eHj/m2M3W4\n17iD6rndu3p5UZMNoJf1+mcYwTqlIjbq8uUu+3IbGszXZEL6ZUp6d/UCgKzIQkItn84pFVQn6NfN\nTCqvYzslYV7HZ4t1CWTuTMmoNfAs/KID6E9ZSaglWKdUUJ3g1svOCt29y28zEKk+gcwJ7AP5vO8a\nOJnjALIgKwm1kcHaGHOOMebxyu3VxpifG2Mer/x3qee5eWPMfcaYPZXH39yuA+8HQXWCYbt3Tb96\nLNF7+NV0kzkOAOkS2hvcGHODpN+TdLhy15mS7rTW3hnwLb8jaZ61do0x5hxJd1TuQwOcOkG/+4N2\n7/LZTVNS77faA4BGnLZirK5lcy9+HkaNrPdLukQn8p3OlPRbxpgfGGN2GGNGPc8/X9JjkmStfVrS\nWa08WJzgt6YdZPHovMBWe3413VGZ4734iw6gP2Wl9WiuVArfXdkYs0LSN6y15xljLpf0rLV2nzHm\ns5LGrLWbXc99QNIj1trHKl//VNIbrbVh/TNbtb1z37n85u/p4KGZyOctHp2n//L5/xDrdZYuHtFD\nN30g9PGBSsr5g1svavTQAaBj9v/8V7r1a09LkrZecY5OPWVJl49IUsKin6RbZH7LWnuocvvbku7y\nPD4lyX25ko8I1JKkiQmSlZJwagWvveR0ff6hvZHPHxzIa8tXnqxOBa1cMaZN61ZXH7/2w6frlp17\nq7e9Pw/344sWDFUT0rZ85cma1+l34+MFfpfbjHPcflk8x4uHB7Ttj070r0jDv298PNnIPmk2+GPG\nmLMrt98n6RnP409J+k1JMsacK4kspDaKM40zVhjW0kUjoRndYbuAuR+PyhwHgKzZvmufrrhtt664\nbXdXS1fjBmtnqvpqSV+sZIefJ+lWSTLG7DTGnCLpW5JmjDFPqZxc9ukWHy88VoasXTtrMy/+/Fd1\nj01OH9UtO/dWR+lhHdWcx/3KwsgMB5AG7ehQlqbS1chpcGvtK5LWVG4/K+mdPs+5zPXlH7Xq4BBt\n07rVuv6epzQ5Xe7rXVgwpCOvHZd0osQrKC2hWFK10Um7prPdFwMA0A5O4ybndqs+z8JKV++45vyW\nvEdcNEXJgA1rV1Vahg5r46Vn1E1pDw2G/5jjXi2SGQ4gbdI0+m0ngnUGLF9W0B3XnO9bjrB91z4d\nn43M8fOdzvZOKyUtgaBNKYB2a2fjpjQNUAjWGea94kz6vX6BNqirWtR7Z/VqF0B2palGm2DdY5KM\nVv2uOCVpcKC+vM+5Wtx87x59ctsTgYHW6aq2Y8uFob+wtCkF0AntHv26lxm7ueRHsO4hrRqtFhbM\nq9lK07lafPiJ/To4NRPZLxwA0qLdo9+wZcZOIlj3kLijVacMK+iKc2x0uKaH+Ia1q5qaMvcrmYhz\ntdvrm8EDSIe4y3O9jGCdYX5XnCcvXaCXDkzVPO+uR56LDNRJp4DStNYDINviLs/1MoJ1D2lkbcZ7\nxRk0Og8TFmjD1tD74WoXADqBYN1DGhmtLl9W0OuWzA9sJRolFxJoo9bQw9qYUtYFAPEl3cgDXbZh\n7arqphpxR6sPbr2o2rjeb2/XwYGcZudq25w5FwJ+nHXmf52q3/ErTncfvyB/1e27VVgwT1+8rq5B\nHgD0PYJ1j3FGq87tpLztSf0CtSR99II3+36/u61fHH5tRv2m4oslafrVY7Ff14u2pgCyjGnwPuRe\nS/YL1JL08BM/qbsvTsZ4t2sRASCLCNZ9yL2WHLX7ubu8KqjJiiNuxrdfolw+V67/BgDUI1j3uSQZ\n5gGbd1XNzhVjNWjxS5QbK4xocIBfRwDww6djn0uSYT4UEUynXz0eu8uZeyp+bHSYzHAACEGwRmA9\ntLu86rovPenbhrRRzlT8QD5f06QlqIVqWLczysAAZB3Bugc57URb9Rp+9dDeZLIjM7ORr5k0uWzb\n+jWabbIPObt7AegHBGv4ikomk8oNUxztaCcap3c4u3sBaESv7U1AnTUkNVafPDp/qJoU1mi5VlCT\nlg1rV+nuR5+vez711AD6ESNr+PLLEncbKwxr46VnNL113KZ1q+v2156dK+muR56rmyKfOnKsbm26\n3XvZAsieXsxzIVjDlzdLvJ1T3n6NWSanj+pQJThf96Un6/bZdtamP3rBqezuBSC2Xs1zIVgjkDtL\n/OMfXNm2HbSiGrMEJbc5a9Ps7gUgrl7NcyFYI5A7S/zcty0L3EGrWVFT7lHCdvdy9FoyCQC4EawR\nWytKxvxsWrda+ajhtQ/WpgEk1at5LmSDI1Snsq4LC+bp0JFjyuekFcsW1TRK8X/+UOg2nADgx7vz\nYNh2wGnCyBqpMDiQ19JFI9qx5UJtveyswOQ2qbw2vfHSMzp8hACyohfzXBhZIxW8I/gNa1fplp17\nJUlXfXCldvzlC5LKI/DBgXzd2rTz/dRhA/2lkb95J8/Fud0LCNZIJfcf07lvW6ZHfvCSpNYEYQI6\ngF77+2caHJkR1OjAuyFJrzVDAOCvF5ubNCpXKkXtUtx2pYmJdBej97rx8YJ6/RxHjYa9jQ6kyj7Z\no8OhyWpOFqjT2rSZq+0snOe04xy3X6+c46C/+Q1rV/XE1Pb4eCFRDQwja2RCUKODqKzyXmiGAKBe\nrzY3aRRr1ugJvba+BACtxMgamRDU6OBNJy8K/b5eaIYAoF6vNjdpFMEameDdeMRpdBBWs82mH0Dv\nCvqbz+rfM8EamRHU6KBTG5IA6KxisZwg3Q9/z6xZIzOWLytox5YLfe/3q9l2rsCd8g/n9qZ1qzt0\nxAAatX3XPh06ckySNJCvb5SUNYys0XfcG5L06t62QD/z/t0enytm/u+WYI2+1m/lH0AW9OPfLcEa\nAICUY80afSGoTvu0FWOBXZDc6CcOpEfcv9ssYWSNTNt8755qoPXTb+UfQBb0498twRp9z1vyFRXg\nAXRfL+5J3QyCNTIr7o48TmnXWGHE98q8n3b2AXqFU6q5Y8uFmR5ROyKDtTHmHGPM45Xbq40xPzfG\nPF7571LPc4eMMV83xjxpjHnaGHNxuw4cCNOqkixKuwCkQWiwNsbcIOkBSc7iwJmS7rTWvrfy3zc9\n3/IfJU1Ya98t6Tck3d3qAwbiaLS0wzuKjvs6m+/doytv/X5zBw0AAaKywfdLukTS1ytfnynprcaY\nD0l6UdKnrLWHXc9/WNJ/q9zOS5pt4bECbTV15Fi1k5mkumxTAOiW0JG1tfZR1QbcpyVtsta+R9JL\nkv7U8/wj1trDxpiCyoH7cy0+XiCWpDvybFu/RrNzxViv3S8JLQDSI2md9bestYcqt78t6S7vE4wx\nb5D0qKR7rLW74rzo+Hj2kwO6rd/O8e3XvVuX3/w9HTxUHikvXTyih276QPg35SSV6u/O56TKfgHK\n56TXLZmvs07/9ZrnDAyUt/PyO8/O9PiDWy9K9o+Ar377Xe4GznH6JA3WjxljNlhr90p6n6Rn3A8a\nY14v6fuS1ltrH4/7ohMTJOu00/h4oS/P8bUfPl237NxbvR11Dk5bXt9oIZ+TrvrgSu34yxckSYUF\n8zQ3V6p5re279ulfJl+TJG35ypN1G4HMzZU0OT2jyz//PZqqNKlff5c7qVfPca81Lkp6QRS3dMsZ\nb1wt6YuV7PDzJN0qScaYnZUR9WckLZZ0kytjfCTREQEtkrS0w9toIZ+TxgojOvdtyzRWGNFAPq9D\nlXVtp4QrKlvcSVgrlspr4klQ7w3E0w/llblSyWfer7NKvXgV10t69Uq5G376i2ndsnNvddpbklau\nGNOLPzuk45417bHCsCanj/q+zlhhWCcvXRDYEjHOxUOvjRQ6gd/l9uu1c+y9YJaS/Z11y/h4IZfk\n+TRFAVyWLytoIF/7Z/HCK5N1gVpSYKB29OPOQECn9cvfGcEafc1vqjluVrgkDQ7UXxxHbSjwq8NH\n9emv/JApbgCxEayBBo0VhvW53z8rcEMBv/IxSSqVpOlXk61fA/CXtEyzVxGs0beCklKC/vgLC4Zq\nvs7ncrr70ecDNxTwZoW7FUuqeV+SyYDG9MsOXARr9KWwLO6gP/6Nl55RDcpLRudVA/3DT+zXWGFE\nr1syv+4D4k0nLwo9Dud9k0y9u7mDPAEf/aofduAiWKMvRSWl+P3xO7tzDeTzevnAiWzZF16Z1OT0\njI7P1gfcrZedpXxEzufk9FGmxYEmRO2clwVJm6IAfcGp0fbjlxleLEn5gKh842VnV5uzFAMqJd3T\n4mHT5wD8Zb3EkZE1+lInk1LcV/0rA5LOHM60+K07n0nU5KEfmkIA/Yxgjb7UTFLK0ED9n81YYVhb\nrzgn8fv6mZw+qpcOTFW/jtpDe+rIMfbcBjKOYI2+1UhSyrb1a3T/5gt8A/2ppyxJ/L5xhTV5CGrY\nkrWmEEA/Y80afStsXTrKhrWrquvQfoHe2yrUvZ7mTItL0utPmt/wvtnO1DeA7GNkDTQgLPs0yfqx\n33S837r2WGFYxWKpehHg1w/Z+/ygiwjKu4DeQ7AGGrRt/Zq6DNSoXbj8vnds1BWsR4cD19OnXz2m\nyenySNqv9Mzt5KUL6i4iNt+7p/r9AHoLwRpooaSbCmzfta8mmeylA1O6/p6n9NEL3lyznr591z4V\nS+USrzjZ3iSZAdlCsAa6KCi47/jLFySV99N++In9daP1AZ8NRPxex32RMHXkWF2bUwC9gWANtFBY\n/fZVt+81mw8bAAATe0lEQVTWVbfvTvyafgF9dq6kXIJs8u279tVkjXdr5J10zZw1dqCMYA20UNL6\nbb/gns9JhQXzIt9rdP5Qdar8jSfXv747ycwvGS2t5V0EaKAewRpoMb/67aA1Z78mKcWS9OrMbPXr\nBSP+FZa/+763VDPSb7zs7MCLBKa8gd5HsAZazKnf3rHlwmqwDMsQ9yuxOj5XVLEkzc4VawK328NP\n/KTm66AmL0GZ47kM71AEZA3BGmhS1LRtVIZ4WIvTQ0eCd+Py7tSVdOehJaPDHd2hqJX9y+NMlUc9\nh+n23sHPimANtF3ARlux97AOyvyenStpcnom8nX81sULC4Y6Oqr2m1246vbddQlubEgC+CNYA01w\nB5cb72v8yt9vcxBHWOZ3eR38xOWA02zFPRLxrovnc9KXN7yro6Nqv9mFYkk1CW5+AZ3ADS7gygjW\nQIO8weVHL074lkMFVVgNugL0W96wOPS9RucPBT52+LXjkcfqrGdL8TLNWy1odsE9lR/Wlc1Z59//\n81/Fej8+4LMhbkfAfkCwBhoUt1tZ1N7Zcfp8b7z0jMDHvYHcL1A569lLF43oi9e9s+b5nVgPDJo5\nmJ0rxf7wnZw+qlu/9nToczbfu0ef3PYEH/AZkbQjYJYRrIE2i6q9DhtRup8btI/2xkvPqAbcdo5E\nmgnqixbOC9wS1Pnw9buo8RM1ambLUGQRwRpoUNSI2a2RvbOd73PMFWuDUGHBkO645nw9/MT+avAK\nan5yy869Xd/EI2r63a/m3G2sMKyli0cYNfeRJH9jWUewBhrkDS5LF48EdisLK6sK6mK2eOG8mueO\nFUa0eOGJgLfx0jMip9AdTkOW1580P9a/rR0GB/KBswPOh6/7oqawYKjmOXdcc75e/Fn9mrVzMeKI\neg8p/WvafrMY/Vi+tGnd6poZmaiOgFlGsAaa4A4uW684p6HX8JsmHyuM1CSgOcHl0JFjGhrIa+mi\nctCP2irTq50j0aBg4g6Mknw/fO9+9HltvndPzUXNxkvPaGg2YtHCeaHLDmlJWurH4NsIZ0Ym6e9B\n1hCsgSa4g8uppyxp+HXCpsm9waXc3azUcHBxr9+2e4Tpd+ylUrl7Wj4nFYulwIB196PPS1LNbMQ7\nTh2ve55fL/Ww89mppKWgYEyQTmZwIK98TrGb/WSVf9NhALFtW7+m6ec5LUr9hAWX01aM1U2DFxYM\n6fCrx5XLlae+gwSNMDesXVX3oegEdef2pnWrg1844thLKo8Sxgojvt/jnCe/gHbL1Wv0B//pMU1O\nH5UkDQ7kNDtX0qEjx3Tdl57UkUpr1oef2F99/eXLCtXXivuzQjq4f++mQrr59QNG1kAKOc1NovhN\noX95w7t00qLyaH9lSIJO3BFmnGnjpCP0xaPDev1J86vf4/0gdl6vWKr/kHaPlGfnTlyNHHH1UH/h\nlclqdzfvsQ3GWNNG9/nNyvRzMiHBGki5qIxYvylfJ9gn3bLTT1RQDwvmQcc+Njpc90F8cGpGt+zc\nG/kh/fAT+2Mdd7EkHZk5Xndsc8ViTaOa8rnLpWaKNerCxz2NnuUpdWqsaxGsgZSLCrhRG3gErd+2\nqiwm7EM16NhfPjDl+1ovH5gO3Xv7xvv2xMp+d7hH3o5i6URHtbh7h7eKOxD7BeSgC59bdz4TOAuB\n/kCwBnpAVJ122LS5d8tOR9xRd7NBvdEacz/P7p+I/dzBgA1QuiWozM49ExF04fOS6+Lm+FxRk9Mz\nmZ8Opsa6FsEa6AFBAbdZcQJpVFCP+lD1G/nH7Vbm93pxzc6VIgN2sVS/1Wgz3J3kvCPnsDK7pNO7\nxZL0ha8/U/cem+/do6tu393U1PiVt34/FVPrrVjCyRKCNdDH4u6BHRbUG/lQDetWNlYY9m2IsnxZ\nwbd0SwreLCVsx7J2mTpyzHcqOyQxvyrJRYx7it95j7jbrvaKRmdlsriWT7AGEKnRdfEw7p3AHE5g\ndm9c4n69W65eU3dhsHTRiE5aNBIYsEfnDwX2JW/HmnVQb/KwUb4zc7Bp3eqGp+8np49GzhL0WhCL\nezHZD6izBtA050PVue3lt57ufM/sXLEaZJzA7M74fviJ/TV13RvWrtLnH9pbvR3UmUw6EQTvfvT5\nam90p/bcabThiFuL7Q12cWu3nYsCp0Y8iF9SXFzOv82ZGo9bD4/0Y2QNoGu2rV+jL173zprRU1Rd\ntzeQO8L6SG9bv0ZjhREVFsxz9R6fp6kjx1qWZe1u4OHlXDQEzTq416wbnbX3jsizsslJ3J4DWUew\nBvpc3A/DTn1ohpWCeUu3vAEpqo/04EC+emHwb8YX1kxZu2u9kwrbUMV90RBnKjcoYS/MWGFYcz4j\n8nbVJad5Oj3tm7Q0imANoCU6Ecz9SrfcASlJH+mg7OyXD0wnHpGGZXp7LxqidgXzzhDkc9Id15wf\n+v4b1q6KlcCWdWnZpKUdCNYAus4d6Jupr3Wmu+O+V5BWjUjzufo1/Ps3XxCZPe9OenNu+wV5ZyvV\noK5uhQVDvuctbIe0f5l8rWdHpVnuekawBpAqYaVgfqVb3kDuF4z9pkaT1np7X8e9xh20J3lQprn7\neP2C6Reve2f19qEjx7R91z4tWjivbk3e2Uo1aGQ/OJCvbmQSNW2d5VFpFkQGa2PMOcaYxyu3Vxtj\nfm6Mebzy36UB3/NrxpifGWPe2uoDBpB9QaVgfqVbUTXdQUHooxecGlrrHWerUqeTmPcCw5mK99s0\nRCqPtp1/n9+xX/elJ2u+djYmWTA82LJucF5ZGJVmuetZaLA2xtwg6QFJzm/hmZLutNa+t/LfN32+\nZ0jS/ZKOtPpgAfSHsPrapDXdYUEorNbb+75+r1MsqRrM3Mfit1uYl/Pv89q+a1/NDmLu1xwaGqh+\nn7Pf97b1a2IFKXfmexYTsKRsdz2LGlnvl3SJTlQTnCnpt4wxPzDG7DDGjPp8zzZJX5V0oHWHCaDf\nBK0tt7JRhvNaixfOa3jE6tSIe9eNG+3hHZasFsSvI9zJSxdUS+EOTs34NmuJs0Nar41KW9mLPk1C\ng7W19lFJ7ku8pyVtsta+R9JLkv7U/XxjzOWSJqy136/cla5O+gAyIUnmeZwg5C7pCroACFrjnp0r\n6fp7nvIt3XKPvL2SZs/nIoLPktHa9fEXXpnUJ7Y9HrlLWdQOad7zkfbSqKx2PUvawexb1tpDldvf\nlnSX5/E/lFQyxrxf0hmSdhpjPmSt/eewFx0fz84JTSvOcWdwntsv6Tm+/bp36/Kbv6eDh8oNS5Yu\nHtFDN32g+vhApZnIg1svSvQ6bmFdyfL5XKJjfsdbxvWjF2tL1PI56Y5PvUennrJEA9/5h5rHnNd+\nxWcEH7cbmnOMN111rq7/0g8kSTdddW7dcfvVuW/+6h5tveIcnXrKkljv1QnOzzRLf49Jg/VjxpgN\n1tq9kt4n6Rn3g5URtySpkpT2yahALUkTE2QbttP4eIFz3AGc5/Zr9Bxf++HTq81Orv3w6TWv4TQT\nifO613749GqrU6/BgVxdcBwrDNe9n5tfi9MNa0/X9fc8VXMBcONlZ2vx8IAmJqY1N1fS1JFj1Wnt\nS274ru7ffIEaLbR2H+Pi4QG9bsl8zc2Vqu/n9uyL9XXuBw/N6OYd/8u3FjxuC9dWu+2T50lKd2xJ\neiERt3TL+TW4WtIXK4H4PEm3SpIxZqcx5g2J3hkAOqRVU6PLlxUCm5p87vfPallykzvxbfHCeTWv\n8fqT5td1Xrv+nqe04uT694lah+yVBKw0d0zrlMiRtbX2FUlrKreflfROn+dc5nPfe1twfADQVklH\nfYsWztPk9Ex10wwn4EnlIOuM4JtJbnIuLianZ+rKv4Ky26XydLlzXH4jfelEAI9aA/dz2oqxwM1S\nmtGtEXgvoSkKACTk3hDEHahandw0VhhJFMDcvdGD1qtzOemkRY0dYzdKo9Ke0NYpBGsASChO9niU\nqCAUlC0enJVe1Kc++g4trQTidpXitLo0Kuw80FXtBII1ACTQig1LmglCQXtUT796vKZMLKwFqvff\nkGRNOO7sQZwRcdR5yEJXtVYhWAPoC53a4jPO+3QiCPlNWYe1QG2luBcjBOP4CNYA0GNWxuw21siU\n9YNbL2r6oqZVQTgrXdVagWANAB3WbBDy7nntTvRyj+yTJLxNTs90vDwq6jxkudd3UgRrAOiwVgSh\nGy87O/GouVNLAXEvRsIuOhxZ7fWdVNIOZgAANV8T3GxNtjNqdm4HiTpOJxFMit4pLK5N61bXdGFz\n16J7FRbM06EjxwKDcdx/Z9YxsgaALkjDhhN+e3RfsuW7+ukvppvuGhZnRLx91z4dqlwgDOTzgeeh\nUzMCacbIGgB6VDsSwY7PFnXzzr1atGBeYOZ4nPeNGhH7XShcf89T2rB2VV+PoIMwsgaALknriLFU\nOrFPd7tQtpUMwRoA+lRQN7RWmZye0eR0/ZaiSI5gDQB9atO61RocqG9MOlYYrvYZbxdqqJMhWANA\nH/vc759V8/XSxSO645rzm+50tn3XPhVL5V3A/NqNUkOdDMEaAPrY8mUFLV54YreurVec0/Rr+rUb\nver23XXtRqmhjo9gDQB9bnAgr3yuvCXnqacsafr1/JLHiiXVJY+loXytVxCsAQBIOYI1AEBjhZGW\nlZEFbc/JVHfjcqVSqdvHUJqY6L+NxDtpfLwgznH7cZ7bj3PcfuPjBW35ypPVNeeVK8Zq9tB2uppF\nBXZ3u1GH97X62fh4oT4NPwQjawBA1Y337Ym1F3UUv1F0o68FgjUAwOXZ/RN19zXSWSwoYYwuZY0h\nWAMAYnF26Do4NeNbO432IVgDAKrecep43X1jhWEtGZ2XeHp8JV3KWoZgDQCouuXqNb6dxV45UB+U\no6a06VLWOgRrAECNVnYWo0tZaxCsAQA1/DqLNbrxBl3KWoNgDQCIxJR2dxGsAQCxMKXdPYPdPgAA\nQG9wprSd2+gcgjUAILZG+oe3qud4P2MaHACAlCNYAwCQcgRrAABSjmANAEDKkWAGAKhDUli6MLIG\nACDlCNYAAKQcwRoAgJQjWAMAkHIEawAAUo5gDQBAyhGsAQBIOYI1AAApF9kUxRhzjqTbrLXvNcas\nlvRdSS9WHv6qtfabnud/RtLFkoYk3W2t3dniYwYAoK+EBmtjzA2Sfk/S4cpdZ0q601p7Z8DzL5B0\nnrV2jTFmoaQbWnisAAD0paiR9X5Jl0j6euXrMyW91RjzIZVH15+y1h52Pf8iSc8bY74taZGkzS0+\nXgAA+k7omrW19lFJs667npa0yVr7HkkvSfpTz7eMqxzQPyLpakl/0bpDBQCgPyXdyONb1tpDldvf\nlnSX5/FfSvqxtXZW0j8aY2aMMa+z1v4y5DVz4+OFhIeBpDjHncF5bj/OcftxjtMnaTb4Y8aYsyu3\n3yfpGc/jP5T0G5JkjPl1SQslHWzqCAEA6HNxR9alyv+vlnSPMea4pAOSPiFJxpidkj5nrf0rY8y7\njTF/r/KFwHprbcn3FQEAQCy5UolYCgBAmtEUBQCAlCNYAwCQcgRrAABSjmANAEDKJa2zbpgxJi/p\nXkmrJB2VdJW19ieuxy+WdKPKTVi+Zq3d0aljy4oY5/h3Jf2xyuf4eZGtn1jUOXY9788kHbTWfqbD\nh9jzYvweny3pDkk5Sf8k6Q+stce6cay9LMZ5/rCkz6pcDfQ1a+19XTnQHufeX8Nzf6KY18mR9e9I\nmmetXSPpT1T+Y5MkGWOGJN0p6d9Leo+kTxhjfq2Dx5YVYed4vqRbJF1grX2npMWSPtiVo+xtgefY\nYYz5pKS360TJI5IJ+z3OSfozSZdba98l6W8lvbErR9n7on6Xnc/k8yVdb4xZ3OHj63mV/TUekDTs\nuT9xzOtksD5f0mOSZK19WtJZrsdOk7TfWnvIWntc5eYq7+7gsWVF2DmeUXmTlZnK14OSXuvs4WVC\n2DmWMWaNpH8n6X6VR35ILuwcv1XlRksbjTFPSFpirbUdP8JsCP1dlnRc0hJJ81X+XebiMzlnfw3v\nZ0HimNfJYL1I0pTr67nKNIzz2CHXY9Mqj/yQTOA5ttaWrLUTkmSMuU7SQmvt33ThGHtd4Dk2xpws\n6SZJ14pA3Yywz4rXSVoj6SuS3i/pfcaY9wqNCDvPUnmk/b8l/YOk71pr3c9FDD77azgSx7xOBusp\nSe6Gs3lrbbFy+5DnsYKkyU4dWIaEnWMZY/LGmO0qt4pd2+mDy4iwc/wRlYPJ/5C0RdLHjDF/0OHj\ny4Kwc3xQ5RGJrexB8JjqR4SIJ/A8G2P+rcoXncslrZD0emPMRzp+hNmVOOZ1Mlg/Jek3JckYc66k\n51yP/V9JbzHGjBlj5qk8HfB3HTy2rAg7x1J5anZY0odd0+FIJvAcW2u/Yq09q5JIcpuk/2qt/c/d\nOcyeFvZ7/JKkUWPMmytfv0vlkR+SCzvPI5LmJB2tBPB/UXlKHK2ROOZ1rN1oJTHEyTyUpD9UeTvN\nUWvtA8aYD6o8hZiX9KC19qsdObAMCTvHKm+68oykJ13f8mVr7bc7epA9Lur32PW8yyQZa+1nO3+U\nvS3GZ4VzMZST9JS19tPdOdLeFuM8f1rSx1TOd9kv6eOV2QwkYIxZofKF+5pKRU5DMY/e4AAApBxN\nUQAASDmCNQAAKUewBgAg5QjWAACkHMEaAICUI1gDAJByBGsAAFLu/wNsYAYd1k6NEgAAAABJRU5E\nrkJggg==\n", "text": [ "" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is your task: **Fit a Fourier model to this data** of the form\n", "\n", "$$\n", "y = a_0 + a_1 \\sin(\\omega t) + b_1 \\cos(\\omega t) + a_2 \\sin(2 \\omega t) + b_2 \\cos(2\\omega t)\n", "$$\n", "\n", "and note that we can go to as high an order as we like.\n", "First we'll do this using ``scipy.optimize.fmin`` (i.e. an iterative solution), and then we'll do it using the closed-form linear algebra solution.\n", "\n", "You should strive to make these functions as general as possible, i.e." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. Iterative Solution\n", "\n", "We'll start with an iterative solution to the problem.\n", "\n", "1. Create a function which evaluates the model given an array of times *t*, a base frequency $\\omega$, and an array of coefficients $\\theta$. **For this section, we will not treat $\\omega$ as a model parameter, but as a constant.** For your model, you can use $\\omega = 2\\pi$ (we'll relax this assumption in part II).\n", "\n", "2. Create a function which evaluates the log-likelihood as a function of $\\theta$, using the above function. **Keep in mind that ``theta`` must be a one-dimensional array (this is what ``optimize.fmin()`` requires)**\n", "\n", "3. Use ``scipy.optimize.fmin`` to maximize your log-likelihood (i.e. minimize the negative log-likelihood) to find the optimal model.\n", "\n", "4. Plot this model over the data to see how it looks.\n", "\n", "5. Use new data, from ``fig_code.sample_light_curve_2()``, and apply your code again (hopefully you've written the functions so that they can easily be reused, right?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. Direct Solution\n", "\n", "Now that you've done this iteratively, let's compute the direct solution using linear algebra.\n", "\n", "1. Create a function which will construct the design matrix ($X$) given *t*, $\\omega$, and a specified number of terms\n", "\n", "2. Create a function which, given this design matrix, will compute the maximum likelihood parameters using ``np.linalg.dot`` and ``np.linalg.solve``.\n", "\n", "3. Plot this answer over your previous answer. Hint: they should agree!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part II: Bonus \u2013 finding the optimal phase\n", "\n", "It's possible (though much more difficult) to use maximum-likelihood estimation to find the best phase. Loading the data this way gives the raw MJD values of the observations." ] }, { "cell_type": "code", "collapsed": false, "input": [ "t, y, dy = sample_light_curve(phased=False)\n", "plt.errorbar(t, y, dy, fmt='o');" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[10003298 10004892 10013411 ..., 9984569 9987252 999528]\n" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAFVCAYAAAD7Sga4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+QHOV95/HPrCQkJA1oI2+iUPZJFWN/gxwIOuITCP/A\nOIdzsX0+w5koTi4kEnaIsBQKBIoN2Bf/qEBJyAoYQ4w4ozgpb6EKuJK4CvsPgRWjs4PvFscx5AHZ\nwNmJSAQRaCVYoZXm/piZ3d7Z7pmeme7pfrrfryqK2d7Z0TPP9Dzf5/dTqdVqAgAA/hrKOgEAAKA/\nBHMAADxHMAcAwHMEcwAAPEcwBwDAcwRzAAA8N7fTE8xstaRbnHPvMrNVkv5G0tONX9/lnLu/5fkf\nl/R+SfMkfcE5tyvhNAMAgIC2wdzMbpD025KONC6dJ2m7c257xPMvknSBc26NmS2SdEOCaQUAACE6\ntcz3S7pU0lcaP58n6c1m9gHVW+fXOOeOBJ5/iaQfmNnXJJ0m6fqE0wsAAFq0HTN3zj0gaTJw6buS\nNjvn3inpx5I+1fInI6oH/P8u6SpJf5lcUgEAQJiOY+YtHnTOvdx4/DVJt7f8/gVJTzrnJiU9ZWYT\nZvY659wLUS9Yq9VqlUqly2QAAOC1RANft8H8ITPb5Jx7TNK7JX2v5ffflvSHkrab2RmSFkl6sd0L\nVioVHTw43mUy0I2RkSp5PADkc/rI4/SRx4MxMlJN9PXiBvPmaSxXSbrTzI5LOiDpo5JkZrsk3eic\n+7qZvcPM/l71LvwNzjlOcgEAIEWVHJyaVqMWmC5q2oNBPqePPE4feTwYIyPVRLvZ2TQGAADPEcwB\nAPAcwRwAAM8RzAEA8BzBHAAAzxHMAQDwHMEcAADPEcwBAPAcwRwAAM8RzIEMXf/Ffbr+i/uyTgYA\nzxHMAQDwHMEcAADPEcwBAPAcwXwAGBcFAKSJYB6C4AsA8AnBHAAAzxHMW2wbHdOLhyf04uEJbRsd\nyzo5AAB0RDAP2DY6pieePTT18xPPHtJ1dz6q554fzzBVAAC0RzAPeDIQyJsOjR/T7X/1Dz2/Ji19\nROHeAJAUgnmKaOkjCvcGgCQRzAPOWjE869pwdb42XXZOT6+XRksfxcC9ASBJBPOAzWtXabg6f+rn\n4ep83Xb1hVq+rJphqgAAaI9g3iLYCu+1Rd6UdEsfxcG9ASBJBPMWux/ZH/q4F7T0EYV7A0CSCOYB\naUxKSrKlj2LZdNk5GqpIQxXuDQD9IZgHpDEpafmy6lSB7UOri61sB2f5sqqGqws0XF3gxb0BIL8I\n5gAAeI5gHsCkJACAjwjmAWlMSto2OqaTNelkTezyBQBIBcG8xZLFp4Q+7oVvu3yxvSgA+IlgHrBt\ndEzPHJgOtM8cGO8r+Pq0y5dvFY+iODQ+oUPjE1knA4DnCOYBPgXfpJX5vQOA7wjmKWJCHVAsLN1E\nXhHMA5IOvj7t8kXFAwD8RTAPSCP4+rLLl08VjyJpbhrjI1qpQH4QzFskvf2qT7t8tb73LArrMgUI\nVg8ASArBvEWSB634JljZyHvFw3dFWD1Qtpn4VL6QZwTzgCIUsP0IFlAbd+yl4OrS+s9+M3avAqsH\n/FL2sgH5RzAPKHMB21pYHZ2YnHpMwYWyK3PZAD8QzAdg64Y12rphTdbJmKF1bDqssAqi4EpWEVYP\nNLcpBqKUaQ5M1gjmAUUoYHtFmTxYrB7wS5nLBviBYB5AARuNgit5Sa+cGKTgPArmVADZI5i38GVd\neNIqbX5HpSYdy5dVp+41n/K2jJPBGDNH3nUM5ma22swebjxeZWY/NbOHG/9dHvE3P2tmPzGzNyed\n4LT5tC68V2FLbMK6ESWpMsBKjc9Lf7aNjunfDr3qZdq7RWAD8qdtMDezGyTdI6nZ93yepO3OuXc1\n/rs/5G/mSfozSUeTTuyg5HHCWlKiWlUfuujMWUMMS09boJ8ZUKXG59Zer2nnrHt/MGbem7LtRZCl\nTi3z/ZIu1XQv7HmS3mtm3zKznWa2OORvtkq6S9KB5JKJpLRrVbWe5T7ISo3Prb1e0u5z5WXhgrmz\nrlUq0ocueuPUz8xiBgarbTB3zj0gaTJw6buSNjvn3inpx5I+FXy+mf2upIPOuW82LrUbikWOjL/y\nWqJnuaM9HyovUQH5lYnJWddqNWn3Iz8aRLIy8UTE57X9/sczSI0f6HkarNlV7PYedM693Hj8NUm3\nt/z+9yTVzOxXJZ0raZeZfcA596/tXnRkpJhj03nSzONfftOIHn/64IzfLT19gf795dldYYfGj+kL\nD/5A933yPV3/e+s/W6/P3XvTJbGeH5Wum9atzv390VPaKwpdDzg0VMnN+50zp14Xb01P1DLGk7Xa\n1HOj/jYJeckfSTry6vFcpScp/b6nm+/eN6vn6fq79ummdat15uuX9Js8hOg2mD9kZpucc49Jerek\n7wV/2WixS5Iak+Z+v1Mgl6SDB2n9pWlkpDqVx5suO1vX3fmoDo0fk1Qf99v6B2u0/pY9oX978mSt\np8/nhZdelRT/s41KVzevkZVe0n7W8uFZrb3h6nx97INn5+L9Nif0SdKWO/Zq89pVHf8meK+cOFEP\n+Um/l+C9nBd5S0+/ksjj77dUbiXpxZcn9Omd39FtV1/Y12sXRdKVwLhL05qV8askfb4RqC+Q9FlJ\nMrNdZvaGRFOG1IStb05ygk+v3Ws+Lwvsds14nvc06Hc83+dVCd1afOq8rJMASIoRzJ1zzzrn1jQe\nf98597bGTPYPO+eONK5f4Zz7Scvfvcs591Q6yUY/wk6GSyq49BMIfF4W2Mtpe3ndNKaf8XyfJ/a1\nszKisnvt5edmkJr8Y/b/4LFpTMm0K2yTCC5lnCjUawDz9bjdqFmtc+cMeTGxrxeb165SdeF0K7y6\ncF5uelIAiWBeOu0K2zSDy5FXjyf6enlStKVpnVpVZW11BVvhtMjbK2qlLs8I5inzZb3t+CuvpRpc\nGFucKc+FXachl3a/L3KgD7bCaZEjbwjmJRO14cfkidkLjnoJLv2OLfq4Y1ScTVR802kyYtSQTGug\nl1So7uhmnqC9Ilfq8opgnrK8BaejERt+oHe9bKKS98KuUyu03ZBMsyIgFW/XqOYETbSX59UaRUUw\nH4CTNXnR1d6ql+DSb/exj4VlVF1o8sTJyL/xobBrLi9s1Wm8v7kqQar3UKCc8rpao6gI5ilqrreW\npMNHX8s2MR0sWjA3keDSS2ArqhNhkTAgz+vqN+7YG/pYildh+/fD9d6ok7XZf++zIh/ClDTmGAwW\nwTwlra2X4ydO5ma2cpihoUqug4uPXj02u/s9KK/r6jfu2DtjOOboxKTW37pH3/nh87H/PliN6fbv\nAXSPYJ6SvM5WbrdGOK/BxVdxZvDnsaUXNa/inr99QlLn8f5Of49yCO7+V/SdAPOAYJ6SqA7Wl44c\ny7R1nveJVz5aFDKbvci7g/kw3o9s5XkfhaIimA9YraZMW+dpF8RRLf9OXc6+2jY6NqslWmkMUxQt\nuJ06f7rS0m5IJqxyU6lIH3nfytTTiHzIa89kkRHMe+DLRjBR2s0y7bfbN6zlL9XXsXeqmft4QEdY\noZV1ha1fURWy+fPmTD1uNyRzxzXvmDGLvVKR7t1ysc5/y7IUUgtAIpgPXCUHk8uWL6tq6WkLtPS0\n5MfGwzYNaTo0fkxbI4I03XL5kcRQTLAVTou8fBjOGzyCeUzB1nicjWCiWjdLFs8vXPdrq3Zf2LAN\nViR/u+WKWGglMRQTbIXTIgfSRzDvUtyzusMK+erCebkp5NOcRd2u0C/aJiJF3b6UDT/QD18r5z4j\nmHfh8NH4h5G0FvJDFelPN73d+0I+rm4nQfncwg2msSh1lX5P0GNZEjBYBPMYmhOzjofsYtauthks\n5KsLT0ktfXkUVmlZsni+fn7potDn+7zcKZjGe//o4gxTkoy48xeieneY/wCfK+e+Iph30FowhYna\nqnT5surU8p3Pb3xbGsnLrV662Xzdga5ordB+u0iT6mL1fdUIMEgE8w7CCqZu+HhwSFZ83IGOVuhs\n7M8PxswHj2COVCRxxrcPLbOwXptD48e0/f7HM0hNMtLqIu108AyA3hHMO4jaBCWuPO69PQi9nPFd\nJEdePZ51EnqW1vyFou4CWFbtKtuMmQ8ewbyDdpugNJ2kxYEWcQ5ZybM05i90kyc+7gZYNu322/B5\nQquvCOYxNAu2KGGnRJVdvzVzXwrzqCV4v/nuN2WQmuT0M39hZcRnH/fgGeYhFAN7FQwWwTyGZsEW\npWgboSShn5p5N+v5s1bk4YReh4g2r12l6sLpVnh14byuWmVMniqG4GoeWuTpI5jHtHXDGk6DihA1\ndtZLV+3WDWtCZz1TmPsl2Aov6lGw6IzVPINDMO8Cp0F1p9euWp+WNjHRJ1w/rbIkVkIAZUMw7xKn\nQc3ky9g2Bq/XVlmRhy6AtBDMu3T+W5ZNtTjK3iJPa6JS1BSEuXPyd7syvhutrMsyyyDugVMYnNn9\nWeho5xb/999OQrtAdtvVF0pST4X5WSuGZ23GQtd1efD551tUJX7TZefMGlKhMjc4+WvqoPR8WqPK\nmHnyfPr8y4jeqHwimKNnaQay4cXzQx/nDYEnHaxRzi+fJqiWCcEcPUsrkG0bHdOPDxye+vnHBw7n\ndp255O9pb3m2fFlVS09boKWn+XPoDpAlgjn6smTxKaGP++Hb4SXNJXivW3IqgQeF59ME1TIh99Gz\nbaNjeubAdGv5mQPjqbag83x4ydYNa3TvTZdknYxCabf3N7LDPJF8IpijZ4OeCOP74SXoDruHAfER\nzNGztCbC9HtQB4D0MJs9nwjmyJ1+D+qA/9hZEOgOwRyRog5QaUpzIgwHdZQXR6DmG2Pm+UQwR6g4\nLaM0v9TBVjgt8nKhGzff2FshnwjmmCUPLaNgBYJuViBf0liSiv4QzDFL3JZRWi2oPFQmkB26cfNt\n0EtSEQ/BHLlDNyuQX3w/84lgjlnitoxoQZVLpwmRSSFYAN3rGMzNbLWZPdx4vMrMfmpmDzf+u7zl\nufPM7CtmttfMvmtm708r4UhP3AkuaU2EoZIA5Bffz3xqG8zN7AZJ90hqltjnSdrunHtX47/7W/7k\ntyQddM69Q9KvSfpC0gnGYMQ9PCSNQ0aYLZs/g1z3TbDIN76f+dSpZb5f0qWaXlJ8nqT3mtm3zGyn\nmS1uef5uSZ8MvPZkYinFQDUPDxmutj+1Ku7zusVJZPkx6AmJm9eu0tw507sYzJ1TIVjkDN/P/Gkb\nzJ1zD2hmQP6upM3OuXdK+rGkT7U8/6hz7oiZVVUP7DcmnF4M0NYNa7R1w5q2z0mrxZZWJQHdG/QY\n9rbRMU2emN4sePJEjdnSOcP3M3/mdvn8B51zLzcef03S7a1PMLM3SHpA0p3OudE4Lzoyws2QtjTy\n+Oa7981qsV1/1z7dtG61znz9kr5f/75Pvafv1xi0Qt7LFYVuxD80VEnl/T75XHjl4QsP/kD3nX1G\nMfM4Z+Lk8ZxG7wmfRz50G8wfMrNNzrnHJL1b0veCvzSzn5P0TUkbnHMPx33RgwepcadpZKSaSh5/\n/+mDs669+PKEPr3zO7rt6gsT//fyLq18ztpZy4dnnTE/XJ2vj33w7HTeb8QJPidP1n9RxDzOk7j3\n8YkTMz+P5kqHTr15qEu6EhR3aVrz63WVpM83ZrdfIOmzkmRmuxot8o9LOl3SJwMz3jnDEPDYoCc8\nMQEO6F6lVos6yHJgatS005VWi7F1YpQ0XeiWcRytqC1zSXru+XF9ZtdjkqSbr3hr6p/vdXc+qkPj\nxyRNVx6kYudxXvSax7TMuzMyUo06q6onbBrTYlAbYxQBS1TKY9ATnpgtDXSn2zFzYIZNl50z1WKj\n0C22Qba4mpWH5mMA7dEyDxjkxhhFsXxZVTu3XKydWy6m0AVKirIzewTzBk7qAoDuUXbmA8G8gcMd\ngHyJs2kRskfZmQ8EcwAAPEcwb2BtKwB0j7IzHwjmDWkus2K5G4rgylv36Mpb92SdDOQMS1TzgWAe\nwNpWAOhesLyk7MwGwTwgjWVWLNlAEWzcsVcna9LJWv0xELT7kf2hjzE4BPMUsWQDRbBxx14dnZg+\nCfnoxKTW37pH3/nh8xmmCnlBOZcPBPMUsWQDRRAM5E21mnTP3z6RQWqQN5Rz+UAwBwDAcwTzFLFk\nA0WwaMHsIxwqFekj71uZQWqQN5Rz+UAwTxFLNlAEd1zzDlUChzVWKtK9Wy7W+W9Zll2icoAlp3WU\nc/lAME8Zy91QBMFWOC1ytGJpWvYI5ikb9DnQQBrOf8uyqUppnBY5rdZyWb6sOnV/UM5lg2A+ABwY\nAd9tGx2bWmfOfgnsH9GK+yN7BHMAbXW7jrjogY511TORH/lAMO9DkboSi/RekKxu1hGXoWBnXfVM\n5Ec+EMwBtFWLuD554uSsaxTsQDYI5j0qelcigHCsq56J/MgHgnkPitaVSMUE7VQirs+dM7v4KEPB\nzrrqmTavXTXjHqlIpc6PrBDMe1CkrsSiVUyQvG4CdFkCHftHTNu4Y++MoZiaxEE8GSCYl1y3FRMm\nypVPtwG6DIGO/SOmcRBPPhDMe1CGrsQwdMeXVzcBmkAHDF6lVouaqzowtYMH892lG+yKXrliWJvX\nrtJ1dz6qQ+PHJE23VJqaLde8bBQzMlJVVB63drNL0xWTYEEc93ll1i6fkQzyOH3d5nHreffS9EE8\nZd+/v52RkWrUdJSe0DLvIGpMeeH86ZOkliw+JYukJSJuF2qR5gkASbjy1j268tY9WScjcxzEkw8E\n8w6igtg/v3B06udnDox7PWmMQxIA9IODeLJHME+Iz63UOIcklHWeABBm4469U3uRb9yxN+vkZK7b\ng3iQPIJ5B2FBrIzKsuQI6KR1jPjoxCRLsZA5gnkHYUFsZQFbqc3Zx+2UYckR0AlLsZBHBPMYWseU\ni9ZKjbvkbPmyqnZuuVg7t1zs7XsFgCIimMew+5H9sx4XpZXKDnBAdxYtmDvrWnMpVpnF6d1Degjm\nHUQFu6DW9dg+bazCDnBAd1iKFW7rhjW52VujjAjmHbQLdq010aK3cn2rqABpOWPpotDHQFYI5j0a\nf+W1WYHNx41V4i45K3pFBYhr2+jYjH0m/vmFo3wXkDmCeQdhwW7unIomT0xvg9sMbJlvjNsDdoAD\nusN3AXlEMO9g89pVmjtneoBs7pyKTpyYHbYPjR+b8bwmH5asFWUyHwCUFcG8g22jYzNa4ZMnapEt\n8OrCU7xcshbnlKuFETN4P3TRG9NOHpAr7IaIPJpdQmOGsC61MMEv8x/f95ikYrVyX4nYKGP3Iz8q\n/SzeIsrbyX+DFHZKYlCnUxOBLNAy7yCqFR5cmhJsgcfZ5xzIszKvWog70ZOhKeQNwbxHp86fW6ov\nM12L5VD2VQtxJ7fFGZoCBolg3kHU6fHz580p1Je504YPRdvCFuGYqQ34qWMwN7PVZvZw4/EqM/up\nmT3c+O/ylucOmdndZrav8XvvZ0d1O/Fr2+jY1NGIReuipGsRRddNDxQ7niFP2k6AM7MbJP22pCON\nS+dJ2u6c2x7xJ/9N0inOuTVmtlrSbY1r3mo38WuoMrPdHtVFuemycwrRgm0etILiOmvF8Ix7WCrX\ncEo3k9vKPEkQ+dOpZb5f0qWa7m0+T9J7zexbZrbTzBa3PP9CSQ9JknPuu5J+JcnE5klRdoADghhO\nqfc6DVfnl6oSA/+1bZk75x4wsxWBS9+V9CXn3JiZfULSpyRdH/j9aZIOB34+YWZDzrmT7f6dkZH8\nFhS//KYRPf70wRnX5s0d0vHJ6bf0xLOHdP1d0YePDA1VMn+PSfz7N9+9byovzn3TiD5zFS2SVll/\nzkkYWXLqVMt0ZMmpuXtPaadnZKSqPz/7jLbPufnufXrx8IQk6fa/+kHhvgt5+8zRWbfrzB90zr3c\nePw1Sbe3/P6wpOBd0DGQS9LBg/mdKbvpsrNndbu91Hgc9OLLE7O2eW0+/2MfPDvT9zgyUu36329d\naytpRvfr408f1O/8z4cKM4SQhF7yOW+2jY7pqZ+8NPXzUz95KVefcx7yuHU4rWjfhX7zmOGHeJKu\nMHU7m/0hM3tr4/G7JX2v5fePSvp1STKz8yUVon85brebrzvASTOPNg0b+28dR5UYQigihoo6I4/a\nOzQ+oUPjE1kno3Titsybzc2rJN1pZsclHZD0UUkys12SbpT0oKT/bGbNA79/L8G0Zmb3I/unWua7\nH9nfcZLQZ3b5vQNc3F3vAAD50DGYO+eelbSm8fj7kt4W8pwrAj/+QVKJy4OwVupwdb6qC+dp/JXj\nkmbPeG2ece5Di7wfTBAqnrLPZo+DPEIesWlMB1FdapIKu+Y66thXFB+z2Tsjj6IVeZ+NvCOY92ju\nnOms8/1L3LoXd1hhFXXsK+OExcPmQJ2RR7OVfSvgrBHMOyj6nuRRX8APXfRG1tqWFPuOd0YezRY1\nSXb7/Y9nkJryIZh3ENWltvuR/ZHdSc1Wrg+ihhF2P/Ij3Xb1hVPdh0Wv1HQrOPsf5cR2rvEcefV4\n1kkoBYJ5DK1dap26k5rPLRLGCcuFQIWkLD51XtZJKAWCeQzNPcl3brlYy5dVC7XOtJsWN+OEdWU+\n7xuIsjKiLLn28nMzSE35EMwHIM9dst20uFsrNWXEJB8g3Oa1q1RdON0Kry6cR+/dABHMe9BNa9aH\nVlww3WVuccdRpF4ZIGnBVjgt8sEimPegXWs2uM5y4469XrTili+rTnWfU4sG0CvKkuwQzHsUNn7c\n2gV7NOQs9Dy24tjoIT5m9QPII4J5j8LWmYats8w7xoC7U5ZZ/VnM88jz3BLE1ywXMVjdHoGKgF6W\n7uStFdduDDi43zymbbrsHO8P0+mEU6/QK5Y0ZoOW+QAVtRVXNkXf/Wvjjr0z5n0AyD+CeYLC1llW\nF85TJcdrsxkD7p4PKxR6tXHH3hlzPY5OTGr9rXv0nR8+39Xr0GVeXnz22SCYJyhsneWfbnq7fibH\nrbhex4DL+oUt+hyDsEmbtZp0z98+Efs1ilzZAfKKYJ6w152+YNbjvG+N2e3ObmUurDlMor2iV3aA\nvCKYJ2jb6JieOTBdaD1zYNyLgqybMWAK63BFOUxi0YLZc2IrFekj71sZ6+973VSnzBXEojk0PsEE\nygwQzBNUht3ByvAee1GUwyTuuOYdqgQOCapUpHu3XKzz37IstX+TCmJxsGdFdgjmQBfKcJjEzy9d\nGPo4jl4mVFJBLAYqZdkimCeoDDPDy/Ae2yn6YRLbRsf0Ly+8MvXzv7zwSlcFclk21cFsVMqyRTBP\nkM8FWdxJej6/x6Rce/m5UxMGi9Qil5IpkLudUFn2CiKQBIJ5wspw5ncZ3mM7g9w0xoclgK1p7DZ/\nqCAWw8KIyZMfuuiNGaSmfAjmCSv67mBSOd5jJ3lfbtirblvJUbPQu82fslcQi+CViD0Kdj/yowxS\nUz4E8xY+tITyoKjBLE+yWK7VTSs5yQlPVBCB/hDMA1jrirzIcmZw3FYyE54QxNyHbBHMG1hWgTzJ\nMlAuX1bVzi0Xa+eWi2klIzbmPmSLYN5AKwPoDi0xtAp+9twHg8V55iko+lhysBdj5YphbV67KuMU\nFc9ZK4Zn7QOft0C5ee0qXXfnozo0fkzSdEsM5bV8WVVDlenHGBxa5g20MuJhOKIu7YmSvnRZJjkL\nnUmV/mM71+wQzBt8KTyzxnDE4PiwXItZ6Giiop8tgnnApsvO0XB1Pi1ytDWoVQ8+TERjBQiaqOhn\ni2AesHxZVbddfSEt8jbKPhxB62MaeQHkB8EcXSn7cAStj2nkBYLKXtHPGsEcXWM4AkCrslf0s0Yw\nR9fqy08qGqpUSvdFpfUxjbxAKx8mbRYVwTwm9myfVuZJT5vXrtLcOZWpn+fOqZS29UFLDK1Y3ZAd\ngnkMZQ5erco+6Wnb6JgmT9Smfp48USvV+29FSwzIB4J5B2UPXq3KPump7O+/lQ/L54AyYDvXDtoV\n3mxdCRTPttGxqe/9WWxX3DV28csGLXN0peyTngb9/pmrMVjNnriapJroiYM/COYdlD14tSr7pKdB\nvn8f5moUrbLBMAp8RTDvoOzBK8ySxaeEPi6LQUz68mGuhg+VDaAsCOYxMGN32rbRMT1zYDqgPHNg\nPHdBJm2DmPSV9xaiD5WNXtATB191DOZmttrMHm659mEzm9W3ZmZDZva/zOzbZrbXzCzJxGaFGbvT\n8h5kiqIWcX3yxMmBpiNKUe8DeuLgq7bB3MxukHSPpPmBa6skrYv4k0skLXLOvU3SpyV9LqF0AsBA\nsF0xfNSpZb5f0qWSKpJkZktVD9DXNK+1eFXS6WZWkXS6pNeSSyrygG7IwQj7cknS3Dn5GBkr8n3A\n6YnwUduSwTn3gKRJqd6FLuleSddKOhLxJ49KWiDpnyT9maQ7EkspBmrb6JjW3bJH627ZM2NyE92Q\ng5H3YMl9AORLpVaLGp2rM7MVkr4qaZOkL0s6qHrAXinpXufctYHnfkL1bvYbzez1kvZI+iXnXLsW\nevsEYOBuvnufHn/64IxrS09foJvWrdaZr1+i/T99Sdft+JYk6bZr3qkzX78ki2QW3u9++ht68eUJ\nSfX8v++T78k4RTNxHwB9ieqA6+3F4gZz59wFgWvLJY0GrzWuf07SYefcrWa2SNI/SlrpnHu1zT9R\nO3jQ7xmweTcyUlU3ebz+lj2hNaxm6wvhus3nTp57flyf2fWYJOnmK95Kq1fJ5zFmI48HY2Skmmgw\nj7uda2vZXgleM7Ndkm6UtFXSl83s7yTNk/TxDoEcQITmCVTNxwAQpWPLfABomaes25p26xpiaXq8\nlqASrZnPwfxbyd7eiaLVmD7yeDCSbpnnY2oscoXJTdO63a60qJupdKto27wCeUcwR6h2a22jZroX\nTS/blRZ1M5UoYUGbbV6BwSOYI1TUWtuytDzL8j77ERa0yTcgGwRzdKUsLc9e32fe14cnJSpot861\nkIp5fwB5QzAHElSW+QZRlR0A2SCYoytlaXn28z7LsLd31BqYOUOzJ+g284FJcUB6COboSllanv28\nzyT39vbTV/K3AAAMuklEQVQtAC5cMLcU9weQNwRzdK0MLU+Jc+zbaXcQDPkGDF7cHeCAKc2WZ9E1\nz7HPSnO2ePNxnjafOWvFcNuNhbLMt6SxCRB8QMscyKG8L/Eqy3BL3j8HoIlgDuSQD0sAu+1OPzQ+\noUPjEwNIWXJ8+BwAiW52AD1Kuju9OdFv64Y1ib0mUBa0zIEcKtoSwG2jYzpZk07W5NUWr0X7HFBc\nBHMgYdtGx7T+lj1a38fe9UUak/Z53LlInwOKjWAOJKgZuGqqb6zST+AqyhKvqC1et9//eAap6V5R\nPgcUG2PmQILaTZjqdjlf0ZZ4tTry6vGskxBL0T8HFAMtc6AN33Zg88niU+dlnQSgMAjmQIKSnjBV\nhMrEyog8ufbyc6d+5gx0oD8E84Bto2Nad8seretj4hLKLckJU0UJcJvXrlJ14XQrvLpw3ow88WGC\nXBEqVSg2gnmDDwUK/JDE3vVFux+DrfDgYyn/G7MUpVKFYiOYN+S9QMHg9VqIJ3FqWtHux+XLqlMz\nwn1a1lW0ShWKi2AeA93v5UMhnrzh6gINVxfMup7njVmKVqlCcRHMG6IKlOHF80tdqCexAYqPsi7E\n8xzgksbGLED/COYNUQXKMwcOz3puWWrmSW6AUiZJ9OSULcDldWOWMlWq4DeCeUASE5eKJOvWaZZ6\nLcST7J7Pa4BLw/Jl1alu+DxVWMpWqYK/2AEuoDlxKeisFcOztqMk2Bff5rWrdN2dj+rQ+DFJ04V4\nJ+wAF82X09CCFbKVK4a16bJz9Jldj0kqfqUK/qJlHhDWPVrmmnnZuxiD77Ms77nswnpWPveV702d\n+Lb7kf0Zpg6IRjBvaNc9Wtbu9zJXZKSZBXfcQrzsFSDfhfWsTJ6oTT1m3gjyqlKr1To/K121gwez\n/2Ksv2WPwnIibvdqno2MVNVtHgcrN3PnVFRdeIo2XXZOaQJ5a+VOmg7KUXnQzOdeuufLqrmrWtwu\n+F7u5W5ElQOtivy5pp3HqBsZqVaSfD1a5pilNZAFWyZl0c/kv7L25HQrjzurhfWsAD4gmDfQPTqt\nzLPYk5DEDnBFl9dNeVqHlubOmd14Kmu5gHwjmDeUfXwYM1G5S1eeK4zBnpUb/8evUC7ACwTzALpH\n6whkVO7KrLVnhXIBPmACXAn0MqGFSVzSc8+PT60vvvmKt3YM5Ewciq+XCYYSeTwI5PFgJD0BjmBe\nAr18OZ97fnyqy7NMs9j7QSHYnV4qjORx+sjjwUg6mLMDHEKF7YYHJImd1YDkEMwBZMKX7Wpbt3fd\nvHZVxikCZmMCHIBcuP6L+6Y2kcmLvC6hA1rRMgcibBsdm1pCdRYtstiK1JJN8uAcIE0E84AiFULo\nT1SLLM5kwDLfR/3kW14EK3GZTw8GYqKbvYHuNAT1uqlJ2e+jPG8GE0fz86spOpCz3hx5RDBv8L0Q\nQj5wH/kt7PMLYvMg5BXd7ECIs1YMR25qgmjd5tuVt+6RJP3i8mG9eHhCUr11nLehiUpFWrKYzx/5\n1bFlbmarzezhlmsfNrPQaadm9nEz22dmj5nZFUklNG1sYYqgXrdzLft91E2+bdyxVydr0smacjM0\nEfX5ffKKt9IiR661DeZmdoOkeyTND1xbJWldxPMvknSBc26NpIsk/UJSCU0be3GjVS97cnMf1fNt\nqCINVaI3g9m4Y6+OTkxGvkZWQxN8fvBVp272/ZIulfQVSTKzpZI+J+ka1YN8q0sk/cDMvibpNEnX\nJ5fU9G267JwZW5g2lXl2cpn1ugte1H1UFnE2g2kXyLNW9s8Pfuq4N7uZrZD0VUkXSnpA0h9JmpD0\nVefcBS3PvUfSGyS9T/VW+V87536xQxpyvTd7rwdC5Al7LQ8G+Rzfulv2tP191HeMPE4feTwYWe7N\nfp6kMyXdJWmBpJVmtt05d23gOS9IetI5NynpKTObMLPXOedeaPfCIyP5DYpPPhc+O/kLD/5A933y\nPRmkqDd5zuMiIZ/jWXzqPB159Xjo75aevqDtd4s8Th957J/Ywdw595ikX5IkM1suabQlkEvStyX9\noaTtZnaGpEWSXuz02rmuBUZ0XJw8Wct3ugOoaQ8G+Rzf7X/4dq2/dY+CHYNDjXbKxz54dmQ+ksfp\nI48HI+kKU9x15q0hrRK8Zma7zOz1zrmvSxozs7+X9NeSNjjnvN5Eqeyzk4G0fOR9K6cef/T9KzVc\nXaDh6gJvhq+APOE88xh6OXc5T6hpDwb53J/mIStbN6yJfA55nD7yeDCSHjNnB7gYelmiBADAoLAD\nXAy9LlECAGAQCOYAcqFd9zqA9uhmBwDAcwRzAAA8RzAHAMBzBHMAADxHMAcAwHMEcwAAPEcwBwDA\ncwRzAAA8RzAHAMBzBHMAADxHMAcAwHMEcwAAPEcwBwDAcwRzAAA8RzAHAMBzBHMAADxHMAcAwHME\ncwAAPEcwBwDAcwRzAAA8RzAHAMBzBHMAADxHMAcAwHMEcwAAPEcwBwDAcwRzAAA8RzAHAMBzBHMA\nADxHMAcAwHMEcwAAPEcwBwDAcwRzAAA8RzAHAMBzBHMAADxHMAcAwHMEcwAAPEcwBwDAcwRzAAA8\nRzAHAMBzBHMAADzXMZib2Woze7jl2ofNbF+bv/lZM/uJmb05iUQCAIBoc9v90sxukPTbko4Erq2S\ntK7N38yT9GeSjiaURgAA0Eanlvl+SZdKqkiSmS2V9DlJ1zSvhdgq6S5JBxJKIwAAaKNtMHfOPSBp\nUpLMbEjSvZKuVaClHmRmvyvpoHPum41LUQEfAAAkpFKr1do+wcxWSPqqpE2SvizpoKQFklZKutc5\nd23gud+SVGv8d64kJ+kDzrl/TSPxAACgi2DunLsgcG25pNHgtZC/e1jS7zvnnkoorQAAIETcpWmt\nEb8SvGZmu8zsDYmlCgAAxNaxZQ4AAPKNTWMAAPAcwRwAAM8RzAEA8BzBHAAAz7XdzrUXZvZ/Jb3c\n+PEZSbdLukP1zWeOSfod59y/mdmfSrpQ0njjuf+18Zy/kDTSuH6Fc+4FMztf0o7G77/pnPt00un2\nSUse/1jSbZK+1Pj5aUlXOudOmNlHJH1U9Xz7rHPu62Z2qsjjjrrIY+7jHrXmsXNufeP6hyV9zDm3\npvEz93Efushn7uUehZQXd0j6uqTm0uwvOud2p3kvJ9oyN7MFkuSce1fjv3WNxFztnHuXpAckbWk8\n/T9KuiTw3HFJfyDp+865d0j6c0k3NZ57t6TfdM69TdJqMzs3yXT7JCSP16u+xe4fNfJHkt5vZssk\nbZS0RtJ7JP2JmZ0i8rijuHnc+D/3cQ8i8njW2Q/cx/2Jm88N3Ms9iMjjX5F0W+Da7rTv5aRb5r8s\naaGZfaPx2p+Q9BuBHeDmSXrVzCqS3iTpHjP7OdV3kvuy6rXCWxvPfUjSzWZWlXSKc+6ZxvVvSPpV\nSY8nnHZfhOXxpc65WuPGWCbpJUn/SdKjzrnjko6b2X5J54g8jiNWHje2OOY+7k1YHu/X9NkP9zSe\nx33cn1j5zL3cl9Y8vlH1ipGZ2QdU78m7Rinfy0mPmR+VtNU59x5JV0n6S9W3f5WZrZF0taTPS1qk\nevf7b0n6NUkbzOxsSadpuqtiXNLpjWuHA/9G83pZheVxxcz+g6QfSloq6R8kVTWdl1J4fpLH4eLm\n8UJxH/eqNY9HJd2n2Wc/BPNS4j7uVtx85l7uXWse/4Wk/yNps3Punap3u39KKZfJSQfzp1Qv+OSc\ne1rSi5LOMLPfUP0ktV93zr0o6RVJtzvnJpxzRyTtUb12c7jxJqT6G3+pca0a+DdOa1wvq9A8ds79\nP+fcm1Q/fna7ZudbWH6Sx+Hi5jH3ce9a83i5pF9UvZz4qqSVZrZd9cKP+7h3cfOZe7l3YeXFN5xz\nY43fPyhplVIuk5MO5r+n+kQhmdkZjQRcpHqL/CLn3LON55mkb5vZUOP887epXpN5VNKvN57zXyTt\nbYzbvGZmv9Donr9E0t6E0+2TsDz+kpmd2fj9EUknJP29pLeb2XwzO13SWZL+UeRxHHHz+M3iPu5V\nax7/kyRrzK1ZK+mJxiFOj4n7uB9x85kyuXdh5cWDZvbWxu9/VdL3lHKZnPSY+b2Svmxme1Xfu329\npL+R9JykB8xMkh5xzv2xmf25pP8t6bik+5xzT5rZs5J2mdnfqT7z/cON1212dc5RvcbzWMLp9kkw\nj6X6jVSRdJ+ZvaZ6l8+Vzrl/NbPbJf2d6pW2TzjnjpnZXSKPO+kmj7mPe9Oax+uccycbj6fOfnDO\nPc993Je4+fwk93LPwsqLVyXdaWbHJR2Q9FHn3JE072X2ZgcAwHNsGgMAgOcI5gAAeI5gDgCA5wjm\nAAB4jmAOAIDnCOYAAHiOYA4AgOf+P/9HfQ8DtZHSAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model we want to fit is:\n", "\n", "$$\n", "y = a_0 + a_1\\sin(\\omega t) + b_1\\cos(\\omega t)\n", "$$\n", "\n", "except this time our parameter vector is $\\theta = [a_0, a_1, a_2, \\omega]$. Since we're fitting for $\\omega$ itself, this is **not a linear model**, so the closed-form solution will not work. Furthermore, this is a **non-convex** problem, so standard optimization will not work either!\n", "\n", "Instead, we'll do a bit of a hack: for each value of $\\omega$, we find the best $(a_0, a_1, b_1)$, and then report the value of the log-likelihood as a function of $\\omega$.\n", "\n", "Your task is to plot $\\omega$ vs. $\\log L_{max}(\\omega)$ and find the $\\omega$ which best-fits the data. Once you've done this, divide *t* by the phase, and re-produce the plots you did above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll see later in the week that there is a more efficient way of doing this procedure, called the *Lomb-Scargle Periodogram*. At its core, however, it's essentially doing exactly what we do here, just much more efficiently." ] } ], "metadata": {} } ] }