{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 6. Model comparison (part II)\n", "\n", "Slides introducing this section are available [here](model_comparison_part2_slides.pdf)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model\n", "\n", "We again consider the sine model with gaussian measurement errors.\n", "\n", "$$y = A_1 \\sin\\left(2 \\pi \\left(\\frac{t}{P_1} + t_1\\right)\\right) + B + \\epsilon$$\n", "\n", "where $\\epsilon \\sim \\mathrm{Normal}(0, \\sigma)$\n", "\n", "We want to test if this is preferred over pure noise.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from numpy import pi, sin\n", "\n", "def sine_model1(t, B, A1, P1, t1):\n", " return A1 * sin((t / P1 + t1) * 2 * pi) + B\n", "\n", "def sine_model0(t, B):\n", " return B + t*0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "The model has four unknown parameters per component:\n", "\n", "* the signal offset $B$\n", "* the amplitude $A$\n", "* the period $P$\n", "* the time offset $t_0$\n", "\n", "## Generating data\n", "\n", "Lets generate some data following this model:\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "np.random.seed(42)\n", "\n", "n_data = 50\n", "\n", "# time of observations\n", "t = np.random.uniform(0, 5, size=n_data)\n", "# measurement values\n", "yerr = 1.0\n", "y = np.random.normal(sine_model1(t, B=1.0, A1=0.9, P1=3, t1=0), yerr)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualise the data\n", "\n", "Lets plot the data first to see what is going on:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEECAYAAAArlo9mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaw0lEQVR4nO3db4wcZ30H8O/3ro5zl1Aa7MRFquw7IgVOSm0iBpRCQowoEiUvEkCVItxAQuoF1BLxBqE2KQxSAwikNoBUpOFFIIqjFKkt4k8CKEi2oMRBa0SmEUcQ9BwIqmLiiJBgxzH2ry9mN3fe27nbPzPPM8883490ut292Z1ndud++8zz/J7noZlBRETab8Z3AURExA0FfBGRSCjgi4hEQgFfRCQSCvgiIpEINuCT7Pgug086fh2/7zL4pOOf7PiDDfgAov7AoePX8cdNxz+BP6q6FFXYvn27LSwsbLjN/Pw8kiSJdhCBjl/Hr+PX8Q8+fuTIkafM7OKy5zUy4C8sLKDb7W64TZIkm27TZjp+Hb+OX8c/iOTjGz0v2CadTifuKzodv44/Zjr+yY6fTZxaIUkSi/nbW0RkEiSPmFlS9vdga/giIjIeBXwRkUgo4IuIREIBX0QkEgr4IiKRUMAXEYmEAn5T5ClwL1d/8tRveUSkdRo50jZKu1Pg2MHi9l8e9FgQEWkr1fBFRCKhgC8iEgkFfBGRSCjgi4hEQgFfRJopT5W5VjFl6YhIMylzrXJeavgkzyOZk1zwsX8RkRj5atL5MIB5T/sWEYmS84BP8pUAFgB83/W+RURi5jTgk5wB8CkA/1jy9w7J7vLyMpIkQZZlLosnIhKkLMuQJAkALJHskhy6BqLTJQ5Jvg/AC2Z2F8kvAUjN7OjgdtEucfjg3uK3OqhECvqfGMtmSxy6ztJ5E4CXk3wPgFcBeBXJt5nZ047LISISHadNOmZ2g5ldY2Z7AXwLwA0K9iXyVDnIIlIpb3n4ZnaTr30HQTnIIlKxdo+0zVPVkkVEeto90jbkWnKeAscOFbfvJXD5x4rjERGZULsDfsh2p8Bjdxa3//q3HgsiIm3R7iYdERF5kQK+iEgkFPBFRCKhgC8iYcpTZeGNSZ22TbFyAHjqMHD2FPDVBWDPHb5LJNJsIWfheaIafhOsHAB+2CmCPQCceLy4f+YFv+USkVZRwG+CR24Dzpw497EzJ4Czz/spj4i0kgJ+E5z4Zckf3M1kKiLtp4DfBPM7S/5Ap8UQiVKeRtP5q4DfBHvuAGYHVnycnQcwA5x+ZuMTMU+jOVlFarE7BS65pvh5l7V6ChNl6TTB4r7i98O3FB2387uKL4Hu3wF/eA64+KryLARlKkhbDctc6/+vyERUw2+KxX3A9iuLWsb1R3ViS9zKMtdWDvgtV+AU8EWkecoy1x65zU95WkIBX0SapyxzrTSjTUahgC8izVOWuVaa0SajcBrwSW4jeT/Jh0h2SV7hcv8iEoiyzLXQpxzJU69Zda6zdG4F8CMzu53k9QBSANc5LoOMK0+BRz++el+rb0ndyjLXQk9m8JxV5zrgPwDgVyQJYCeAJxzvXyah1E/xYXEf8IsvFrd13lXCaZOOmR02s18D+AaATwO4a+3fSXZIdpeXl5EkCbIsc1m8Zlk5AJz+HWBnilxkpaOJSIksy5AkCQAs9ZrLO8O2c1rDJ3kxgN+a2bUkXwvgKyQvNTMDADPLAGRJkli323VZtGbp5yD359I5e6p3H+Ff0opI5TqdDjqdDkgum1lStp3rLJ3PALild/s4gK217q0/Uu/YoWKkXii1ZOUgi0gNXLfhpwDuJrm/t+/9/dp95cpG6gHNryUrB1lEauA04JvZUQBvdLKzjWrJTQ/4W14GnD4+/HERkQm1d+BVyLXkslmRmzhbcp5qts62ylN9ti3T3tky53cWzTjDHm+6F54e73GflLLZXtN8tnmqsRsN1N6Av+eO3rqwa5p1fIzUy9PxT/yNvqyGvZ5I06gi0EjtDfhNGak3yYm/0ZfV4r71r/fg3ipKKhIWzZc/tvYGfMD9SL08reYytn/SPnQjAANmtgKvy3Qyi/SFnIXnUXs7bX2ocqm0xX3Alj8GOFssjKKTWGSVxqpMRAFfRMITchaeR+0J+Hk6eQrZNM8VEfc0X/5E2tOGP01WgDIKRMJSZRZeRJ2/7anhu5Cnbq8EzrxQzJa50VxAoc4XJDKNxX1FIsNMbzqu+V2TJTZEtlh6e2r4LtR5JZCnRdAGii+TP3sHcPbk6t/XZiH0KVMhXnmqgU1VZOGFPAXLBBTwm2J3eu4/7FcX1m/TPxEv7P0tspNV1lAzZDUi6/xVwG+qjU7EfsCP7GQVrK/Zz+/yVZJ2CHkKlgmoDb+pRslCUKZCM+Rp/X07/X30g/38rmK8R//LXybT1sXSSyjgN9WwE27wRIzsZG2sKgfcjbIPBfrqVNX5GwgF/KZa3AfMzK3eH3YiRnayitRicV8xmv2Sa4Drj9b7/+M5q04B35c83bwZYPa8YmqFjU5EFyerUj/Fhzwtzrljh9oxILIBKaDqtPUllCwLpX6KL4OZa6FrQFad0xo+yTmS95D8Acmc5NUu9y8TcD1JVZ5qmgtppwZk1bmu4e8H8KyZvZ7kGwB8DsBrHJehkKfr09vUEbae65M0lCsfkXE1IAXUdRv+YwA+37v9BIBZx/tfNZhZoWA/XOnJaKu18OeOuiyRVClPdUXlSgOy6pwGfDP7tpn9hOQrANwD4Pa1fyfZIdldXl5GkiTIssxl8WSYspP0Ja/Sl2UbuEgplUKNWXVZliFJEgBYItkl2Rm2nfNOW5LXAfgIgA+a2ZG1fzOzDECWJIl1u13XRZNhypaK7M9hIiKjq2kVvk6ng06nA5LLZpaUbee60/ZaADcC2DsY7KXBNkr9DDVlM0/VlLGZUD/bvjzVZzzAdQ3/vQAuA/AdkgCwYmY3Oy5DGPIUOP1McfvYoeJ+0y63Tz4ZbsqmOoc3Nm06bhPmmNdnvI7rNvx3mtmfm9ne3k/7g/2ktaR+2+qWlxa/+8G+SbWu369oXdG2miYdtwEDjGS4do+0zVO/I/WqPvHLXu/kk9WUd1z9cgzSbJ3hmyYdVwuMN1a7R9oOjtTL0yLw91Wdez94GXv6uclH1uXp6oIo/Sad//3S8Nf7/Qowt6OKIxjPzNbhQV+zdYZvmpzxBgwwkuHaXcMfVGfu/bDa9+njw7cd5cTfnRZl7Jd3d1r+vLKadt0uWPSeVyw1mSZnXNN2b+y5o6utDo5bHuIK+HUadhlbZtITv+x5/bxe1+Z2aLZO1/rBov9T16C3aXLGGzDAqNEuXFjtm3M89qE9AX+azswqOkJHvVyd5sQv+0e6YHGy16uCy6llm8xVra0fLPoBA6gv9XDSz1bTdjdWO9rwp0khqyq1sKzN87xtwB+eO3fQ0qQnvgZBNdeFC8Dp3wIXvdptCuCFC6tNk01KPaxpgJFMpx01/GmyAsZJLdzoSqCs9v2az1ZbA66yRp2nGpgyqjzVe9U0edqu+fIdaEcNf5qsgFFTCze7iiirfa+t6TSNBqaMrmnv1ZkXVjPCZrb6bdbzpW3z5TvQjhr+NFkBZR2eg88d5SpC7dniwpkXgLMnVysfZ08Bz/5MA5tkU+0I+NNkBYyaWqjc4rDkaTubYE4+WQT7dc5qYJNsqh0Bf5qsgFFTC6vOLc5TtT/WyeW0v/2+ndPPFL/rqGmvHAB+89/Asz8t30aVD9lEOwI+MF1zyijPrTq3eHe6ml6necjDNdi3c/ZU9fPG9Pdhf9h4Ow1smkyeTlb5ytPgriLb0WnrwkadsnXJ09XpFe4lcPnH1n8xjLKN1MfFwtQjDeqb0cCmSU3a+du0jvwRKOCPo8rc4pUDwOnfAbDy6WNHORGVqeCXi76dUV7rJZcpSUA21Z4mnZD0L9FhxX1NHxsuF/PGbPhaM8Vykz4mz5PgKOD7oOljyzVpvv9RuJg3Ztg+ih0VNXsFexlRvE06g1MZz865+8dRiudw066y5MNg387M1urnjRncBwjMnA9sf111+5AoxFnDHzZ/zrM/c7eQSOjTxw7Wwqt630K98ulneW15afG7ji+ntZlkW/4YmD2vmtcN7YpKpuI84JOcJXkTyftc7/tFw+bPwdnicRdCnj525QBw+OaBL8vHqgn6uvJxS0sRupenxZfr6WdWFzZyyGnAJ7kNwKMAPuVyv+uUzZ/jaiGR/kAx9Fbfavr0sXm6mqf80N8AdnpgA6vmyzL0Kx+gmCY5lNzsaa+o8jScY22KYQsbOeR6EfPjZrYE4AaX+12nbP4clwuJLO4rLs23vLT58+70T9J3GV78khpUxZelyyufupoyLlxwN8J3WtNeUbkczTyJ0Jqr8rT2L9BGteGT7JDsLi8vI0kSZFlWz46GzZ+DmThnHBxXnatuuVo4w1dTRp42q0bchiuqMiE2V03xBZplGZIkAYAlkl2SnWHbNSpLx8wyAFmSJNbtduvb0dwO4PLbzx016zJLJzR5Cjz68dX73DLQrFPhl6WLhTNcjI4dpmkjM/fcUQTBte9FKH1Jm/H1GXvS6XTQ6XRActnMkrLtGlXDd2pw/hwF+3KDNY8r78I5/Q+h5YKrc7jQ5qUI9RkPFVfAryudMDaD/Q8hBXug3U0Z42rrGg76jIfyEvDN7KCZue24Hdam5zL3Xpoj5LRYGU2TP2OPncnx1PCHzjjoMPe+rfI0vHn9Q27KOPlkMelef+79OisseRreZ9vX1M/Yc2dyozpta1XWducq976tQp2ts+rO4X4gPnaovjVmTz5ZXJX2J907e6pYEOW5nxdz5ZfNujqpaT7bPPU/bbeLBIBxee5MjqeGX2c64UbytNpaUtWvJ9NbOTAkENewxuzvVwCcXf94f2GUJqUerh270cQcfV88dybHE/CHzjjoIPe+6hO/ytcLbWCKTxu9V4/chvWBuIY1Zke5Gg1h7qGYee5M3jTgk/wWybeTDLv5Z1ibXtXphHkaTu07xIEpvmz2XrmqtY16NRp56mGjee5MHqWG/yEAVwA4SPITJMMdjlp37n1Il7Ghzkzpw2bvlata2wWLGOlfNvLUw0bz3Jm86dljZj81s48CeHNv+5+Q/CbJPbWXTurT9oEpVTZXbfZe7bkD6/+Valhjdm5HcVX64nxGs1g3t1FTUg+lnMexD6M06fwFyS8AOAzgWQCLAN4H4F9rLtt48nTy5pRpnjup/pq2p5/x037e5oEpVTdXbfZeLe47NxDPbK1vjdm5HauD3i65CnjJK5uXeiiNNUq7/EcB/BuAvzezM/0HSTYr4E+TQuY6tbBsTVvA3T+r5lEZ3Sjv1dwO4Pn/Ay569URFntjcjtWmyaakHsZicNW8KlNiazJKk85fmdnX1wb73uNfr69YDZWn1VwJNKH9vKkDU6pQdXNVm98rmUygSQ9hZ964VtWVQFPaz5s4MKUK8zuLf8Bhj0+qae/VmedXBzYBxZeQuBPobJzx5OE3SZvbz5ugyfOoVGX2/NVssEuuKRZeEXeaUmkbkwK+DzEEJJ/a3ATz3NGio//0M80f69FmgVba1KTjQz/wPHQjACsCUgAdPkFx2QSTp6vNK8cOFZ9n1TXutfsAin1cf7S4/eDeavclmws06UE1fF9CWtNWNja4MHUdzStrB/WpCce/QK8iVcOX8fXHEMCqSUfL03OXUKyjhixStaZ15I9ANXxf8jTMttiyMQTTpKMNLqEYcrA/+WQ8E9Jp8r3gqIbvS9MWtB5VoOloTrw4X31v5kwfA+pcKctDB9p3rK44GMjltIbPwp0kHyZ5mORbXe5fKlB3OlqehjPj6KBh89WvHVDXphpxEwYPtomjgVyua/hvAbAbwJUALgXwXZILZmaOy9GMFXlCsbbmwVng3EHXharS0UJdQQson6/+xC/bVyMONA+9sRxdObtuw78CwPes8HMUs01td1yGQkhTGfs0GKiGBXtuaXw6mhNl89XP72xujThPJ7uiCjQPvbEcfYG6DvjnA3h6zf3jAOb6d0h2SHaXl5eRJAmyLHNcPFln6OLva8zvAq68K8xaatWGzVffz81uao140oqPBg9Wa8ov0CzLkCQJACyR7JLsDNvOdcB/HsC2Nfe3AzjZv2NmmZklS0tL6Ha76HSGlllc2iggeZjPu9H689UPy81uW4040Dz0xpryC7TT6aDb7QLAspklZja0tuw64P8YwFW9zttLUeT2PeW4DDIOX4u/1y1P6+kcntsxfHGLNtaIPS7kMZI8DScBwNEXqOtO2+8AuBbAw7377/PSYSujKxtCHmrNtM9153D/H/fhW4r+EE2nUb/QEgAcDORyGvB7wf1Wl/uUKZUFqv6JKaMLcGSmtItG2srmmn7pLiIjiSvg52k4bXoiIhWLa2qF0Nr0RKR98tTboM+4Ar5I6PJ0fbCQsHiseCrgi1QhT92sMTsYLNbuV1OEyCYU8EWqsDYQu1yBSs2UMoa4Om3lXHmqTuxR5aneKwmeavgxU+1wdHqvpAVUwxcRiYQCvrRTnqoJRmSAmnR8yVNlV9RJTTAi6yjg+6KAJBK2PA2u0qaAL/45WLxZpHIBVtrUhi9+OVq8WUQU8MW3pq71KtJCCvjiV1PXehVpIQV88atta72KNJgCvvjVxrVeRRrKacAnOUvyJpL3udyvVCxPqxvU5GLx5jwtytn/ydPqXlskIM7SMkluA/B9ABcBOOhqv1KDqtPR6l7rdXcKHDtY3+uvlafB5WZPLE/jOdaWcBbwzew4gCWSewG839V+RZwKMDd7YjEdqwt5WvsXaKPa8El2SHaXl5eRJAmyLPNdJJHq5Knm95Fyu1PgXbb6M0awz7IMSZIARaW6S7IzbDuaWSVlPedFi529e+Dhu80s69fwzeyGsucnSWLdbrfycsmE8hR49OOr9+u4dO8vGlJXk0vdry/SACSPmFlS9vdamnTMLAOg6nlb6NJdpBUa1aQjIiL1cT55mpkdhLJ0REScUw1fRMSnPHU2TkTTI4uI+ORwnIhq+CIikVDAFxGJhAK+iEgkFPBFRCKhgC8iEgkFfBGRSCjgi4hEQgFfRCQSCvgiIpFQwBcRiYQCvohIJBTwRUQioYAvIhIJBXzxL0/rXet15QDw1OHi9b+6UNwXiZCmRxb/6lxCceUA8MMOcPZUcf/E48V9AFjcV88+RRpKNXxpt0duA86cOPexMyeKx0Ui4yzgk5wjeQ/JH5DMSV7tat8SsRO/HO9xkRZzWcPfD+BZM3s9gA8AuNPhviVW8zvHe1ykxVwG/McAfL53+wkAsw73LbHacwcwO3/uY7PzxeMikXEW8M3s22b2E5KvAHAPgNsHtyHZIdldXl5GkiTIssxV8aStFvcBr8uAma3F/fldxX112EqLZFmGJEkAYIlkl2Rn2HY0s8p33tvZuwcevhvAkwA+AuCDZnak7PlJkli32628XBKxB/cWv2teJFpkIhWdnySPmFlS9vda0jLNLANwTvWc5LUAbgaw18xeqGO/IiLB6Y8TOXuqGCey547arkBd5uG/F8BlAL5DEgBWzOxmh/sXEWkWx+NEnAV8M3unq32JiARho3EiNQR8DbwSEfHF8TgRBXwREV8cjxNRwBcR8cXxOBEFfBERXxyPE9FsmSIiPi3uA37xxeJ2zeNEVMMXEYmEAr6ISCQU8EVEIqGALyISCQV8EZFIKOCLiERCAV9EJBIK+CIikVDAFxGJhAK+iEgkFPBFRCKhgC8iEgkFfBGRSDgL+CS3kbyf5EMkuySvcLVviVyeAscOFT/3srgvEiGX0yPfCuBHZnY7yesBpACuc7h/idXutPgRiZzLgP8AgF+RJICdAJ5wuG8Rkeg5a9Ixs8Nm9msA3wDwaQB3DW5DskOyu7y8jCRJkGWZq+KJiAQryzIkSQIAS70m886w7Whmle+8t7N3Dzz8dQD/YmanSb4WwL8DuNSGFCBJEut2u5WXS0SkkR7cW/yecsUrkkfMLCn7ey01fDPLzOyqtT8AlgDc0tvkOICtdexbRESGc9mGnwK4m+T+3n73D6vdi4hIPZwFfDM7CuCNrvYnIiLn0sArERGf8tTZOBGXTToiIjLI4TgR1fBFRCKhgC8iEgkFfBGRSCjgi4hEQgFfRCQSCvgiIpEINuDHPrGajl/HHzMd/2THr4AfKB2/jj9mOv7Jjr+W2TKnRfI3AB7fZLMlAMsOitNUOn4dv44/XmXHv8vMLi57UiMD/ihIdjeaBrTtdPw6fh2/jn/c5wXbpAMg7ms6Hb+OP246/gkEW8MXEZHxhFzDFxGRMSjgi4hEIqiAz8KdJB8meZjkW32XyQeSsyRvInmf77K4RHKO5D0kf0AyJ3m17zK5RnIbyftJPtRbrPoK32XygeR5vXNgwXdZXOrFwKMkD/Z+bhvn+aHNh/8WALsBXAngUgDfJbkQ01KJJLcB+D6AiwAc9Fsa5/YDeNbMXk/yDQA+B+A1nsvk2q0AfmRmt5O8HsXSodd5LZEfHwYw77sQHuwCcMjM3jPJk4Oq4QO4AsD3rPBzAASw3XOZnDKz42a2BOAG32Xx4DEAn+/dfgLArMey+PIAgC+QJICdKN6HqJB8JYAFFBWf2LwawKUkD5H8Zu+9GFloNfzzATy95v5xAHOeyiKOmdm3AYDkKwB8GcDtfkvknpkdBgCS3wTwZgBX+S2RWyRnAHwKQAfAZzwXx4enAXwSwP0A3oHi/+DKUZ8cWg3/eQDb1tzfDuCkp7KIBySvA3APgA+Z2Td8l8c1kheT3GJm1wK4GsBXerX9WOwH8DUz+43vgnhyBMC3es3Y/4Witj9yHA8t4P8YwFW9jotLARiAp/wWSVwheS2AGwHsNbMjvsvjyWcA3NK7fRzAVo9l8eFNAG4ieRDAWwHcR/Jlfovk1CcBfLB3+3IAT5jZ2VGfHNTAq15N5rNYvYT5mJk94LFI3pDcC+D9ZhZNWz7J/wBwGYpABwArZnazxyI518tKuRvABSiaZP/BzO73WihPSH4JQGpmRz0XxRmSLwdwL4A/QdHi8QEz+/HIzw8p4IuIyORCa9IREZEJKeCLiERCAV9EJBIK+CIikVDAFxGJhAK+iEgkFPBFRCKhgC8yApJ/S/I/e6O8P0Hyn3yXSWRcGnglMoLeKO8vAzgF4E8BXDfOkHaRJghttkwRL8zMSH4awP8AeLuCvYRITToiIyB5HoAvoFhs5BMkL/FcJJGxqYYvMpp/BnC/mX2N5IUA7ib5NtX0JSRqwxcRiYSadEREIqGALyISCQV8EZFIKOCLiERCAV9EJBIK+CIikVDAFxGJxP8Dygm5HJN49poAAAAASUVORK5CYII=\n", "text/plain": [ "