TestOscillationFitv1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Two level oscillation tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy.linalg as npl\n",
    "from IPython.core.pylabtools import figsize\n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "create an artificial data set with several points per time value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x106fc7908>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAERCAYAAABowZDXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4XGW9L/DvNEmTtEnTlDPKw/EgF8uPculuKsjVCuwt\neAMvHCh6pBxaYoU2HIRKURS2KFI2BcWI9OwIsnm2h9I+yinKoV5O2lMobkFNQaW8pbjLHTf2Ti9p\n0875Y2aSycyaZK1ZvzXrXe98Pz48JjPNyps1a6339nt/byqTyYCIiGg0Y+IuABERJQMrDCIi8oUV\nBhER+cIKg4iIfGGFQUREvrDCICIiX+rj+sUicgqARcaYs4te/xKAOQDezr001xizodrlIyKi4WKp\nMETkegCfB/COx9vTAVxqjOmrbqmIiGgkcQ1JbQTwGQApj/feD+CrIvKEiNxQ3WIREVE5sVQYxpif\nAhgo8/ZDAOYCOAfAmSLy8aoVjIiIyrJx0vtuY8wWY8x+AI8B6Ii7QEREFOOktxcRaQPwnIgcB2A3\nsr2M+0b6mUwmk0mlvEa2iIhoBIEfnHFXGBkAEJHPAmgxxvTk5i1WAegH8GtjzMqRDpBKpfD22zuj\nL2mNSKdbeT4V8Xzq4bnUlU63Bv6ZlAPZajO8iPTwptTF86mH51JXOt0auIdh4xwGERFZKO4hKSIi\nqqLuvh6YrRuRQebgspn3Buo0sMIgIqoR3X09eGHri/lvOSRFRETezNaNoX6eFQYREfnCISkHLV7a\nh/WbtgIAphzRjgWXcO0jEQHS/r7CIanA2MNwzOKlfXh+01ZkkF3k8vymrbjunrV4+S2GIxLVuq6O\nTkxsbKv451lhOCbfsyi0dWc/vveT52IoDZXT3deD+b0LMb93Ibr7euIuDtWQuVMvy1carwf9WVYY\nRFWWj1TJ5P73wtYXcePaW/HKztfiLhrVgMNb34Nbz7gRy2be+56gP8sKwzFTjmgvea29tRFXXzg1\nhtKQF68x5G3923HPuhHTphHFjhWGYxZc0oH21sbB79tbG3HnvDPw3kOD542h6tq1f3fcRSAaESsM\nB1194VS0tzayZ5Ew4xvGx10EohExrNZB7z20FXfOOyPuYlAZx7ZPLhmWmtjYhrlTL4upRHZhWLi9\n2MMgqrKujk60NLQMft/S0IJbz7gRh7cGnoN0DsPC7cYKgygG86bNxsTGNkxsbMO8abPjLo41GBZu\nNw5JEcUgH9pIlCSsMBzEMWBKqilHtOP5ol4GgzfswR33HJMfAy6Uv+H8hNZyVzNd5c5nfk8CIJvf\np6ujs9pFs9Z196zF1p39AIbCwgFem9q44x5xDDgBuNJ7ZAwLtxeHpMgTh7Wi47Unwbb+7fifz/0L\n5zXAsHCbsYfhGI3UIAxtJCIvrDAco5EahMNa0ZL295W8xoV7lASsMBzEMWC7dXV0oj41NBpcn6rn\nwj1KBM5hOCjsGDBDG6PV3deDgczA4PcDmQHcuPZWzJ16GSsNshp7GFSCGW+jNdKkN5HNar6HwXh4\nb1dfOHVwzoI9i+F4zXhjZJ37anrhXj4evlB+8rFWhwa4OGpkQa8Zr/Pp4nUXdsGoH7w2dXHhXkAc\nGqCgNK6Zro7O/J7KALKVRdInvRlZVxtqusIgisvcqZcNZqtlOC0lRU3PYUj7+7iRDQWidc24lq2W\nkXW1oabnMADgxrW3Ylv/dgBDQwO1jOPEowtyzdTS+SyXNFBLkHPJCfjR1eQcxsyHr0J3X0/FP8+h\nAQqK14w3WxaMMrVNdBLfw7j44SszQPKjTGxRSy3iauD51OP3XM5Z1Auvp1oUvZ4kq6SH4cwchgvZ\nPtmNJiKbJX5IqtDAwYHR/5Gl2I0m0qGRsZm8OVVhAIF7WNZgHDuRDqa2iY4zQ1IAUD+mLu4iUASY\nioOCYmqbaDhTYSQ9YoVx7N6K02jktzOtJMDBporHprK4iLv2RcOJISkXUiuwG+1NK32LTfto21QW\noiASX2FMap6Y6J5FIVvi2F1kU94wm8pio8VL+zBnUS/mLOrF4qV9cReHCiR+SGrJBbc5E+duUzfa\nlhBfpm9JDo1rpjjrbT5a8OoLpyKdru3etg0S38MgfTaF+GpldrVpH22byqJF65phtKDdWGE4KGyX\nvnjyHcjetHctW6dRvMA0UnHYlFLcprJo4YO+NrDCcEyUvYN39uwPfYxK5DO7hn2o2pQDyqay2ISL\n7uyW+DkMGm6kll7Y+ZHmxuCXiy1zIYBdKcVtKouGcU312LV3eKaFVAq46KyjAx1nwSUdkWe9pcqx\nh0G+1Y0JtpLeprkQitbuvaVpeTIZYPnqlwIfi9GC9mIPwzEaCwBTgGe2z/q6YO2LKHs75C6bogW1\nuLJQM/E9jLD7YXT39WB+70LM710Y6ji20FgAyHFkCorXTHkuLdSMrcIQkVNEZJXH6+eLyNMi8pSI\nXDHaccJ8AC59kIXCdum1Vp3b9hBxrXFgE2YqKM+lhZqxVBgicj2AHgCNRa83ALgLwIcBfAjAF0Tk\nXX6OWckH4NIHWSjfpQ9zw2qMI9v0EHG1cWATzj0kw5fX3Ix5vdcPbj4XRFw9jI0APoPSfORTAGw0\nxmw3xuwH8CSAGdUuHOlUOoA9DxFXGwc20bpmXGPTQs0vr7kZuwf2VPzzsUx6G2N+KiJHeLw1AcD2\ngu93Amjz+HclKvkAmHYiepzATAabwp812PT3dHV04sa1t2Jbf/bRll+oGYcwlQVgX5TUdgCFTZNW\nAKWhNkUmNU/EkgtuC/zLbjn3Wnzx0a9gy55toY7jGhdz9pzwbsEf//rCsNcmNU/E9WdeifQk/3/v\nN1ffXZJu/eu/+TauP/NKHDXpcM+fsf18fn3JUyX5m75871P42uxT8L73TIyxZKX8nEsb/54bZlyF\nf3ryXgAIfM3ZxLYK4wUAk0WkHcAuZIej7hjpByY1T0TnCbMqTkDYecKswWGJMMdxRTrd6uQ5+OLx\ns3HjtuGtvG+e9lXgAAL9vX/6qyl5bcuebVi05geercYknM9nX3y75LXN2/filh/+m1W9Q7/n0sa/\npxXt2esNCHzNaRpX3xyqlxF3WG0GAETksyLSmZu3uBbALwA8BeA+Y8ybIx1gyQW3hUoXoZV2guzH\ndBxU6+6Y8Q2kQmxlHVsPwxizCcDpua8fKnj95wB+HlOxSJlNY8ka6ThcnPdybbdH1/4ebbOOm4kH\nn38YGWQOBP3ZVCYTOLLKNhnbu/xJojmEUry3ATB04yY5iibIBGYShqQAqORvirpxEORcMh/V6NLp\n1sBdjbiHpMhhrqa8dnFoK2z4s2beMI0d92wJ53YNexg0jGaLeM6iXs+cVLXU4ktKDyMsrc96pF7p\nSSceVhPnslrYwyCr2JYahOznaq/UFbaF1VYdF2FFx9W9DVy8ZsLiRHNtqOkhqXx+oUL5cemkhtiG\nnWjWHkJ5+a2dg63DpE92A8GvmVoZkgJ0Jpq7vrvGcyOmzk8ch/PPmlwz59KLVkMlf5wMMpllM+8N\nNMpU0xXG/N6FyHiMvMa5dD+ssGPJLj/gNG64oNdM1OfTpt6ORuNgpOv3wX/8SNWvTVt661qN2+Lj\nLJt5b6B5DM5hUE1wMVutbX+TTckHNSKtbNoxUit5ptdxgqjpCsOmLJJaONHsTeuGa65vKnkthRQ+\nefRHKy5bpVzMwKtx/Wo96DkBX6qmK4yujk5MbBxKhpsfVkjq/AVg1x4ULvLKw5NBBj95kckJNGhc\nvy4+6LUat17HCaKmKwyAi7CSQGOnvKh7k7v271I5ThAu9pCB8NdvuVnZgQMHAx3Hpt66LY3bmp70\nplK2TXprRrJp7Ekwr/d6z9dbGsbj9g/eXPJ61OfTln0WqsHvuZy9qNfz9dZxDbj76g8G+p02hYW/\nsvO1weHGSiM5i4M2OOlNVgk7+ag5Tq/Rmzy2fXLJaxMb2zBv2pyKjheWiz1km9jUW1+x8XFs79+B\n7f07sGLj47GUgQv3LApLdE3xmpD85GNc6zE0stV2dXRi4RPfwDu5IaiWhvGxtuo1/iYqz5YdI4t7\n2vmIuKA9Da9sy0Ekvocx8+GrKh7X1gpL1Ajhc5HG5KON4/Tzps0ZbNXH1bPI05jfcU25MZb6uuQ+\n7rR62sVzIUEl9wzmhHnQa3wItmXpdI0tk32FbNl0y7Z1GLawabLaRvlhTACvB/3ZxFcYedv6t+Oe\ndfdX/fdqhfDZtEhIi9aNy3F6by6uw9CgGVpuSyNOs6edb/Asm3lv4NaOMxUGEDy0UWsR1lh5Bk0n\nr0TTySsxVp4J9LN5tsWOa9woWjfussfexltPnIa3njgNyx4r3a+5VnmlKAGAgYOBN1JzjsZktU2N\nuK6OTtSnhqac61P1sfRwnaowxjeMC/Tv9wzsLXktgwxWvOQ/AmHitD7UtW1GKpVNklbXthnNHatx\n8cfSgcpiE80bxaaNeWrFgczA6P/IcRppSmxqxHX39WCg4HMdyAzEMvzoTIUR1wRk/9j/KH2xYS8e\nfWN5oONojruG7R1o3ihhb1ybbtqk2DvQH3cRSJktw49OVBj50Mag3TObInAWXNKB+rqh+I76ulRF\nD1m2yCloT9s2tswbcPK8VOIrjEnNEyvuWWhE4GhVOouX9mHgwNCY9MCBTGwJ02zq7fCmLc+2RYQa\nbGrw2JSXzZbGbeIrjCUX3BZq4idsBI5W2KfW0ItGHh2tG0Xj5rfpprVNV0cnWhrGD35faU/bJrYN\nQdqy0tuW8PKaX+mtsVJ27tTLhuV4qYRWwjQt7S2NQzl0WhpH+dfeinf+A7I3/13L1gXK6XP1hVOH\nbczjAq0MA4c0TxpcdX5I8yS18lGW1kpvjc9b4zkTFpMPWqJcwrRxTfX4/jUzfB9HY8e9hd1rQm3z\nmlfub0qlgPsWnuP7OFriSAPjlTAvqt3TKj2OTUbaYvikEw+zKjGmX7Z+Tul0a6DEg4ADQ1Ku29Mf\nLERSY8w/6mGBluYGleMEYdOq6Ch3T0v6wj0XhyBd+pxYYSiIMqoj6MPVphvuuDKV17UXT6t6WVy6\naV1ny7wBlWKFEZJWVIfmwzXsDacVmbTgkg60jhuq8FrHNSS+taghyt3TXEifYtPe4Bpc+px8zWGI\nyPEAJqEgEaQxZk2E5Qoi1jmMsHMGhf7H957Azt37AVS22YuG/Ji71sYxL7+1c9iEdVwPgLjGkctt\n+qO18RE3UIrW4qV9g0O0U45ox4JLOio6jo2fUyVzGKNGSYnIPQDOB/AXDA/mOTvoL7ORTfthXHvx\nNGuigbQik2zZT8A2WhEvNkTOuEpzP5dxb56KrS2rs19vOVWzmFU1ag9DRF4EMNUYs6c6RQqs4h6G\nRqtzpKiOoBeVDZWXbVu0AjqtvOKtKfOibu3ZeD6TqtrnUmv0QPMZoSH/nMkgk1k2895A0xJ+/vFf\nfP67xNGYCNWaZLYpiscmNq38JaqETYsRC58zKL/XVFl+KoKtAJ4XkYdE5Ee5/6q/8YTF3n3Sc2g+\neSWaT16Jd59U2UXAKB5vWjebSxOPVB3jmkpH7FMp4KKzjo6hNDq8njNB+FnpvTL3X753lkL5hcmJ\n4rW/bdCHSHdfD17e9e+DdfXLu/69or12KVpdHZ1WTjzaQmty1yW79w5grDyDMRM2AwAO7jgE+8zJ\nWL76JZx6/KG+jzPliPayQ1JB2DBkPWoPwxjzAIDfA5iAbKTUOmOME81ejfwsXhuqZ3f/uy9QWdgC\n9qaZfJA793njsJ+3BnmmZK+bpmmrcLBpW6DjaAxbaw1Zez1nghi1whCRSwH8bwBHAngvgEdEJLnp\nMItE9RDZtX93oH/fb05CZt/QRZXZ14i9685CZlflG7a7QHMhoi17cdvGpjF2m9TlehaFUmP7Meao\npwMfK+zaKK0h6+JGclB+hqQWAPiAMWYzAIjItwD8PwDBmtCW0kg+6KWpPljCvvWbtgLjpqPxmD8A\nAPo3TMfe3dmbttbDUl1MPkj2y8B7VnjvweABozaFl+dDsbf1b3896M/6qTDG5CsLADDG/E1EuGnw\nKOpSwRMBZ3a3Ye86J5a3qLLpZgPsGEvWpDXGbpNI52QGxgb+kbDXjMZ8a16+kZxOB+9m+4mSek5E\nvisiJ4rIVBG5G8CzgUtZc4LFBXCjoGRwMfzZpvxjGrTmZJr3v7v0xf1NmHXMrEDHcema8VNhdALY\nB+B+AD/KfX1VlIVyQ7AQZ9duWle5Gv5cuOdJpfuf2EJrTubOj1wH7C/oTewfi3vOuwWnHHlMoONo\nXDO2XHejjpsYY3YDuL4KZanIzIevinVYIIWU5wri+jF1gY/FsfpouTaUpGXx0j785c0dg9//5c0d\nFafAcM2sY2bhwQ3/mvv68zGXJn5lexgi0pf7/4Me/1kzhxF3F08zHNa1LJ02iTIsMelhuq5FSWkO\n755y5DG457xbKupZ5GlcM7Zcd2UrDGNMR+7/xxT/B2Bc1UroU1zDArbstUsjiyoskZ+3fWwb3nXp\nmvGzDuM3Rd/XAfhdZCVKIC4Ii1Z3Xw/m9y7E/N6F6O7ribs4zn3eLgZc2LYJU9hrxvo5DBFZBeBD\nua8PFrx1AMCKiMsVWJw3b1RrOSrl0lh9cUbh/FBS0NQrUYQlumLBJR1q+5/YwrZQ7BUbH8f2/h2D\nXyf1nhxpSOrs3PDT94uGpBqMMf+1imUcVZK7eNpcCuEDOJRULRoJNMmbxj1p/RxGgetF5NMicqmI\nzBKRy0XklshL5tOk5olODAtosaXraiPXhpK0DEugmRpKoBm0kRHl3vZJpnFP2tLg8bMc+acAmgFM\nBrAGwAxYNCS15ILbuEGNxWxc4UrDjfRA83u+NHenI2827K7op4chAM4B8AiAOwB8AMDhlf5CERkj\nIktE5CkRWSUiRxe9/yUR+VPuvVUiUlksW42ypesK6HTFuzo6UXegefD7ugPNHEqykGuhuZq07kkb\nkmf6qTD+aozJAHgB2a1a3wDgPxl8qU8BGGuMOR3ADQDuLHp/OoBLc3MoZxtjNoT4XTXHlq4roNMV\nX7y0D7vWT0NmXyMy+xqxa/00pt5WZlMjw0U23ZNh+akw/iwi3QBWA7hGRL4CIEzugDOQ3ZAJxpjf\nAjip6P33A/iqiDwhIjeE+D01y6Wx+vWbtg4mZdy77mxkdrex5apM44GmGZrr4lyIK/eknwrjiwCW\nGWP+DOBmZHsXnwvxOycA2FHw/QERKSzHQwDmIjsMdqaIfDzE76pJNnRdAbZckyTsA01rsZyrmznZ\nck+GNdI6jA+hYFtWEZkBYDuyk+CTQvzOHQAKr6IxxpjCdR53G2N25MrwGIAOAI+NdMB0mpNqmrTO\n5y3nXosvPvoVbNmT3aFsUvNELLngtkDH+LvJaax78e1hrx3S1oSvzT4lMZ97EsqZTk/BPx+1KNwx\nJjYPruVIT2yu6O9e/7L3XMj3H/kjHjjxsEScS5eNFCX1DYyco7vSjRvWAjgfwHIRORXA4NiCiLQh\nm079OAC7ke1ljLpRE6Ok9KTTrarns/OEWYNzFp0nzAp87KsvPBFfuGMVBg5kL8X6uhTuuPJ0AMn4\n3LXPp60WL+3DhleHti7d8Oo2zPrHlcGjpMo8cQ4ezL5R7XPp0iLYYpVUvqlMJti+DWGJSArADwDk\nBzcvR3beosUY0yMinwXwJQD9AH5tjPnGKIfM1MINWS22PeCKwzWBobHxJIRr2nY+ywn7YJyzqNfz\nWR901fhIn/dJJx5W1XNZnGUAGBpSDTqsFOmGThVKp1uD7cEAHxVGLkVIsYwx5pygvywirDAU2faA\nm72o1/P11nENuPvqD1a5NMHZdj69aDwYtSoMAGXTlFT7XM7vXei5dUE+KMAvWxs9lVQYfia9v1Hw\n37cB/BnZYSXKcTGqw3bv7NkfdxGcoRH+rBklZVviwLBcWqPiZwOl1UUv/UpEngbw9UhKlDBc4Rq9\nsfIMxkzIbit/cMch2GdORktzQ8ylokKaCQxtSRyomWXAFaNWGCJSuKo7BeAEhIuScspIrYegF73L\nE2yVap/Wh71jNw9+X9e2Gc0dq/E57n6mRuvB6NqOkV0dnbhx7a3Y1r8dQPChqLwpR7SXHZJKGj+5\npNZgKHYhA+BvALoiK1GN0krjbZuwk339Y/+j9MWGvXj0jeU45UjmhbKJLT0DTRr5m1xKH+9nSOqI\nKpQjsbRaDxoJ4GzD4broafRKi3sXQPbau2fdfbj9gzeHLmOSaSWsdKX35WfHvfeKyAoR2SkiW0Xk\nxyKSrkbhksC27SBtojHZ17jvXaUv7m/CBYddFKZoToh675Nd+3erHIeAR99civ4pj6J/yqN49M2l\nFR3DhuAaP1FSPwbwKwD/GcCRyG7Pys0VCmhEdTCNhrdt6zqQ2TdUIWf2NWJP31lY9n/eHuGnvNlw\nw2mKeu+T8Q3jVY5T6zQqdltSpviZw2g1xny/4PvviMh/j6g8VaexoEZj7FZrgs0mWsN1/Rumo/GY\nPwx+XQkOj5V3bPtkRgOVoTLkt+XFbLhQgW3923FP349w+wx/waaawTVh+OlhrBORS/LfiMh5AP4Y\nXZGqR6vW1mq5upLRMk9juG7KEe0l2WorqXRcioXP0+qVdnV0oqWhZfD7loaWipLkdff1YH7vQszv\nXYjuvp5AP2sjrSG/ckujk7iWyE+F8fcA/peIbBORvwF4HMBFIrJHRBI9yKnxENHsKrqS0bJQ2OE6\nG+eIbBna6jcnlQzX7V13FjK72kb4KW/zps0ebKzMmzY78M+7tpc8oDfkd3DHISWvZfY1Yuxrp/g+\nhubCyDBGrTCMMe8xxowxxkw0xvyn3NeNxphmY8y4ahRyJBcsWBHrTetiy1VTfrguzENeY45I64az\nqVe6ftNW9G+YPri5VP+G6RVfeys2Po7t/TuwvX8HVmx8PPDPcy/58t7Xf25Jxd780kdwzcdn+D6G\nLQ0nP1FS40Xkn0Tk9yLyrIh8R0SsmQ3LZCq/aW2ptfNc69Jr0ah0tG4423qlxcN1lXCxd6BBMxCl\nuGKvhA0pU/wMSX0fwDhks8peBmAsgCVRFqoSlbSstMbYi1XygfKmLU+rIrXhhgP0eqVa194LW7zX\nYXT//n7fx3Axyk9ra1WtXSM1Gk5h+akw3m+MmW+Mec4Ys84YMw+l26omli1j7Jpdepd6KpoVqcYN\nZ1OvVOvaKzcpu2uv/0nZro5O1KeGgi7rU/VOzMVNHDvB8+ta5afCSInI4F2S+9q66f1Kb1pbxti1\n2NZTCVt52TY2rvGQHtdUGs2eSgEXnXV04PJoXHvlJmX3veh/6KS7rwcDmYHB7wcyA4nvIXf39WDT\nzlcHv9+089WK/iabGhlh+akw7gLwtIjcKSJ3AXgGwN3RFiuYuCNnNCodrS69TQ9Y2yovLWEf0rv3\nDpS8lskAy1e/FPhYGtdew8unlUzK9j97Nq445zTfxxgpvUhSad1LtkxYa/CzcO/HAA5HNp15Ctnd\n8H4UZaGCOKStCfM/fWLcxQjNxYV7GvmxmuubsHtgz7DXUkjhk0d/VKWMlXAtyV73NTNwRfdmjJ2c\nXRy578XpuG+hzv5ocaUX+fKamwevm3H1zbhjxmgbd0arZnJJAfghgOMBfBrAhQA+BOA7URYqiAdu\nOi+RNbUXjYV7rk0+7hnYW/JaBhmseCl46KctbByiuOKc09D/7NmBexajaapvUjuWX4WVBQDsHtiD\n+b0L8fRbfwh0HM17yYYJaw1+ehgfADDFGJMBABF5FNld95xg0167GpkxbeqpcAMabzamuz71+ENx\n6vGHqh+3LuWnTaqruEcKZBsZDz7/MD5wqP95ma6OTlz9q5txoC57vLoDzbH2+m14Vvn5NF8DcFTB\n9+8C8EY0xQkuzMI9WxJ6abMlskMjLNG1HlOeTYESGlLFyZJy6sf4aZPaafHSPuxaP21w/cSu9dNi\nez7Y8qxKZTLlguqyROTXAE4F8GsAAwDOBvB67r+MMeZjURdyJOdftyIDVLapuubG9Ro0Ep0Vb8QE\nDD1g/Tyo0+lWvP223kX4ys7Xhm1AU0mYpVaPKY4Wmvb5tFXY684Pv+eya9UNOJg5WPL6J448Fx89\n8h98/z6bng9RlCWdbvWu5Ufgp4fxLQCfAPBdZBfxXQjgagC35/6zQtLTcWhFFNkUJQXo5MfSmNux\npYVG0SvXCH7yjd9WuSTu8bPj3uoqlCMWNu21y13PytOY27ElPbQ2G8a1Ad0dI23Z296m54MtZan+\njFREKjl5SYiPDhqW6OqYP5XS7DXZkh1Ao6etdQ/Y9HywpSxOVBhhTp7tk49Bdz3Tyn/jGhtDWcPS\nykll00NaY0hV8x6w6flgQ1mSG8KQE3bhni2LsDR3PZs79bJhE81kZyhrWOXCVfr3Hwh0HI3hJJvC\nuYFsdOBgWUJECtryfADsKEviK4wHbjovVBSKLWPAXR2dWPjELXhn/zsAhnY9q4TGmL+LXFltO5o9\n/aWpR6pBo6GisXanXA4ozYitONjwrBo1rDYBMpVWGMX7PAOVhedq0QhBBcJNGtZKGGi1RHk+Zy/q\n9Xw9lUKg1B7VCIkNolxPxe+5nN+7EBmP/lfcvZ4wonhWRRVW6yzbdsvTCEF1NeEflRpfJutt5yeO\nC3Qc2+a9XNvbXoMtz6rED0nRcJrhjWS37mtmYM7tvcgPEgTtWRSyad4r7JCqjSlpwg4nlRsHGjhQ\nukAxSjXdw9Dc5zns/sxElej8xHFIpSrrWRTS6N3awrYek0uLRmu6h6EROVM8tpi/GOKaB7GxdWUL\nGyYNtUWVNLAStiy4A+zqMWksGk3Bu5dRX1fdNn9N9zCA8LHNtowt5tnWurKFS608G2nOnWksInSp\nxwTYs44o8RVGmGy1gDt56gvZNGloywri4ggTIFux37VsXQylcY9WDjMXgzamHNGOsfIMmk5eiaaT\nV2KsPBP4Yc+V3koymXhbi7bU/IVsaV0l4eZ/Z49129PXNNuSZ2polN+hrm3z4FxTXdtmNE1bjdT4\n7YGOY8NK78RXGHlxDQPZUvPbKAk3f0tzQ9xFcAJzmJWndR/YMBriTIURp4ktYz2/JnscV6YneO3F\n0wIfi1GZAbC1AAANIElEQVRxpbTmzljx2M2ZCiOubtripX14fULv4Pjk6xN6OZmaY9PNv+CSDrSO\nG+pNtI5rqKilxsnz8jTmzlwM2rDpPgjLibDa+rpUxUm5woZabmz8JeraNg9+X9e2GXuaV+K7j23D\nd+Z8vKIyucK2hHTXXjwtdC4prX01ur67Brv2ZnM+jW+qR/c1Myoqj020cpjZFBKrod+chMzEV5Aa\nmw3fz+xrxN71ZyFzdBuQsJFrJ3oYAwcyFbXyNFqLYyZsLnktNbYf+97D3b0AuyK2bBgDBoZXFgCw\na+8A5tzei3/781uxlckmtgRtaFm/aSv6N0wf3Bu8f8P0xO4Q6kQPA6islae1oMYLJ1OzXMucq7Hz\nWWFlkZfJAD0/f96aRXiUpbUYMbO7DXvXna1ZtFg40cOI07GTJpe8NrGxDfM6Lo+hNBS1BZd0oL5u\nqJmQHw5lVJx7tMLCbQy9r5QzPYxKPgCN1mJXRycWrvkm3hnIDmO11LdW3KLWas3YlKLBNYuX9mHg\nwFCShvxwaJBUMOOb6kt6GWFzQYVlU9oUW65frUSeLm3e5UQPo9K1D1prKOZ1XD44Tl9pz0KrNZOE\nxXJJppEKpvuaGUgVjGXms8zGNRxlU+SXq9evDYvuNCS+wjikrSnUB6DxQWpM0mkt7knCYjnSyzKr\nwaZ8aMWJM4Hs9XvPuvuqXhbNcFiNgAut9T/545x/3YrAudETPyQVdovW5as3Yluuq7h89UYnMphq\nsGVYwDbjygwnXXTW0YGOY1OW2STYtX931X+nTWHhWlmxi47DHfeCsKkrrtWa0TiOq8MCGnaXiXBa\nvvqlGEqjIwmTsuMbxsfye20JC9fqBXodJ4iarjBs6oprrXDt6uhEfWqo41ifqg98HA5r1Rab8qEd\n214m6nDa7EDH0cqS7NqakLCqXmGIyBgRWSIiT4nIKhE5uuj980Xk6dz7V4x2vLDpzW2i0Zrp7uvB\nQGaoFTyQGWDvQFESWuOVsGVStqujEy0NLYPftzS0BH5Yu9hD1rruvI4TRBw9jE8BGGuMOR3ADQDu\nzL8hIg0A7gLwYQAfAvAFEXnXSAcLk97ctpvflslzl3LfaLOpNa7JllXwADBv2uyhqMOAPQvAzR6y\n1nVXfJyg4qgwzgCwEgCMMb8FcFLBe1MAbDTGbDfG7AfwJABfSXaSutTeRi4mgNNkS2vcVRwG8qZ1\n3eWPA+D1oD8bR5TUBAA7Cr4/ICJjjDEHc+8V7iqyE0AbIqKVSM4mWnt6tzVOGIwOaWucoFY+F+Rb\n42QnV/e117ru8sdJp4PXxnFUGDswPEdjvrIAspVF4XutAHxN6x/S1oSvzT4F6XSALlqZndXHjEkF\nO45Fbjn3Wnzx0a9gy55tAIBJzROx5ILbAh1jyZ/vx8s7Xh38/uUdr+Lrv/k2rj/zShw16XDV8taC\npF5LNvJzLjXuAfIWR4WxFsD5AJaLyKkACseRXgAwWUTaAexCdjjqjtEO2N7aiDuuPB0AAq3JmPJe\n79Qg8z99Yqi1HXHrPGHW4Hht5wmzAv0t6XQr/vRXU/L6lj3bsGjND5xKJFgN6XRroq8lmwQ5lxPq\nW7EF2wa/5mdQqpKGTBwVxiMAPiwia3PfXy4inwXQYozpEZFrAfwC2fmV+4wxb450sEPamjD/0ydW\nVBCXcrwUci1DLFEQ3X092LRzqIe8aeeruHHtrZg79TLOiYSUymQ8xmSSJROm9fDyWzuHbaoTd4RI\n3NLpVtz0y7vKjgEHveFqfcU4exh6/J7L+b0LkfEYa457Ay/bpNOtXOkdlE3hhLbQipJyMR6eqJbV\nfIVB3gojoyqNkrIpkRygl7yN7MZ1RNFJfPJB0tfd11MSJaU5BhxHIjmt5G35Y9myd4RNZbGFTUkD\nXcMeBpWIeqVsHInktPKG2ZSw0qay2MaWpIGuYQ+DSnhNGALAwMHSTK0jObZ9snMLqLQWe2r0DFxc\neKqFkYLRYA+DAggWVKGRSE6LTXnD2DOgpGKFQSVSZSqG+jF1gY8VNpGcbTQqHq3hMZsqQaoNrDCo\nhGaUiS2J5LQe0jZlq7WpLJQcYbZoZYVBJZitdmRhs4Zq9gxczJyrtfkRlSocDkUFW7TW/EpvGi6/\nmvaVna8NRkW5kFKhOKwWGHpIR9kiL7c62cWUNBryiz0L5Xu37z9qClfNhzR7Ue+w73925ycDVRqM\nkiJPWlEmtqQGsS1v2NUXTh2WkoayRgrp/uejFsVQIirEISmKjG2pQWwavmFKGkoiVhgUGdu2yuRD\n2n5M62E3DkkRBcR0HNFhWo9oldkzzjf2MCgyLrYWueguekzrER2vCL0gGCVFw+SjerQmq11rLc5Z\n1OvZQis3iV4uSoq9lOC4t4iOwuCPoFFS7GFQCc3JarYWS7GXQnHKB38AeD3oz7KHQcOk062Y+fBV\n3LGsjKDrObxaxUF7KZRV6z0M7V4pd9wjihjTcVAcNHulTA1CqlycrNZkU2oQ7iJYG6LYzwUVpAZh\nhUElmEtqZGHXc2j1UjgXQkF5VTxBsMIgT5ysjpbGqnOtVifZz5ZU9ly4R564Y1m08r0UIj+0cqFN\nOaK9JGgjCPYwiBLKllYnVYdGr7R4ODSoxFcYFyxYwck+qkmM2KotWrnQanodxvnXrcgA1dnboBbU\neqy7tqjP58tv7RyWJt2F679clgFem7oqWYfhTIUBcOGTBt6Uung+g+EGStXDhXtElGi2pcSn4Zyp\nMDjZR0QULScqDE72EbmBWQbslvgK45C2JvYsiBzBLAN2S3yF8cBN57FnQeQQZhmwF1d6E5FVmGXA\nXonvYRARUXWwh0FEVEPyGzFlgIM/u/OTgToNrDCIiGpE0Y6RgRfuscIgT+XSMxBRcnE/DFKXT8+Q\nyf3vha0v4sa1t+KVna/FXTQiihErDCrB9AxEbvJKiR8EKwwiohpR8/thkD6mZyByV03vhwEgw5TH\nevLpuG9ceyu29W8HMJSegYJjenM9PJe6mN6c1DA9AxEVY1gteWJ6BiIqxh4GERH5wgqDiIh8YYVB\nRES+VHUOQ0SaAfwrgDSAnQAuM8b8rejf3A3gjNz7GQCfMsbsqGY5iYioVLUnva8E8Kwx5hYRmQng\nawCuKfo30wGca4zZUuWyERHRCKo9JHUGgJW5r1cC+IfCN0VkDIDJAHpE5EkRubzK5SMiojIi62GI\nyByU9h7+CiA/vLQTQFvR++MAfA/AXbmyrRKR3xlj/hhVOYmIyJ/IKgxjzH0A7it8TUR+AiC/AXcr\ngG1FP7YbwPeMMXtz/74XwN8BYIVBRBSzas9hrAXwMQDPAPgogDVF7wuAh0RkOoA6AGcCeGCUY6bS\n6dZR/gkFwfOpi+dTD89lvKpdYdwL4F9E5AkA/QA+BwAi8iUAG40xPxORBwH8BsB+AA8YY9ZXuYxE\nROTBheSDRERUBVy4R0REvrDCICIiX1hhEBGRL6wwiIjIl8Tsh5FbBf4DAFORjbC6whjzUsH75wP4\nOoABAPcbY34YS0ETwMe5/BKAOQDezr001xizoeoFTRAROQXAImPM2UWv87qswAjnk9dmACLSAOB+\nAO8F0AjgW8aYnxW8H+j6TEyFAeBTAMYaY07PXUx35l7Ln5S7AJyE7OK/tSLyqDHmP2Irrd3Knsuc\n6QAuNcb0xVK6hBGR6wF8HsA7Ra/zuqxAufOZw2szmP8G4G1jzKUi0g5gHYCfAZVdn0kakhrMQ2WM\n+S2yf2TeFGTXcWw3xuwH8CSAGdUvYmKMdC4B4P0AvioiT4jIDdUuXAJtBPAZAMV7JPO6rEy58wnw\n2gxqOYCbcl+PQbYnkRf4+kxShTEBQ3moAOBAbmgl/972gve88lTRkJHOJQA8BGAugHMAnCkiH69m\n4ZLGGPNTDL8R83hdVmCE8wnw2gzEGLPLGPOOiLQiW3kU7rsc+PpMUoWxA0N5qABgjDHmYO7r7UXv\ntQLYWq2CJdBI5xIA7jbGbMm1Oh4D0FHV0rmD16U+XpsBich/AdAL4EFjzNKCtwJfn0maw1gL4HwA\ny0XkVADPFbz3AoDJuTG6Xch2q+6ofhETo+y5FJE2AM+JyHHIjmueg6IkkuQbr0tFvDaDE5F3A/gl\ngKuMMauK3g58fSapwngEwIdFZG3u+8tF5LMAWowxPSJyLYBfINtrus8Y82ZcBU2A0c7lDQBWIRtB\n9WtjzMpyB6JhMgDA61KN1/nktRnMV5EdZrpJRPJzGT0AxldyfTKXFBER+ZKkOQwiIooRKwwiIvKF\nFQYREfnCCoOIiHxhhUFERL6wwiAiIl9YYRBFQETaROSRuMtBpIkVBlE02gFMi7sQRJpYYRBF43sA\nDhORn8RdECItrDCIotEF4A1jzIVxF4RICysMomh47eVAlGisMIiIyBdWGETRGECyskETjYoVBlE0\n3gLwioj837gLQqSF6c2JiMgX9jCIiMgXVhhEROQLKwwiIvKFFQYREfnCCoOIiHxhhUFERL6wwiAi\nIl9YYRARkS//H3YL34HVqhaaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106f9f390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "sigma =0.2;\n",
    "Omega =0.5;\n",
    "Nt = 20;\n",
    "tmax =2;\n",
    "Nrep = 5;\n",
    "tlin = np.linspace(0,tmax,Nt);\n",
    "t_1 = tlin[:];\n",
    "t_2 = tlin[:];\n",
    "n1_simu = np.sin(2*np.pi*Omega*tlin)**2;\n",
    "n2_simu = 1 - n1_simu;\n",
    "\n",
    "n1_noise = 0.2*np.random.randn(Nt);\n",
    "n2_noise = 0.2*np.random.randn(Nt);\n",
    "\n",
    "n1_exp = n1_simu+n1_noise;\n",
    "n2_exp = n2_simu+n2_noise;\n",
    "\n",
    "for jj in np.arange(Nrep):\n",
    "    n1_noise = 0.2*np.random.randn(Nt); \n",
    "    n2_noise = 0.2*np.random.randn(Nt);\n",
    "    n2_exp = np.r_[n2_exp, n2_simu+n2_noise]\n",
    "    n1_exp = np.r_[n1_exp, n1_simu+n1_noise]\n",
    "    t_1 = np.r_[t_1, tlin]\n",
    "    t_2 = np.r_[t_2, tlin]\n",
    "\n",
    "nt_exp = np.r_[n1_exp, n2_exp];\n",
    "t_all = np.r_[t_1, t_2];\n",
    "plt.figure(1)\n",
    "plt.clf;\n",
    "plt.plot(t_1,n1_exp, 'o');\n",
    "plt.plot(t_1,n2_exp, 'o');\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('population')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the simulated datas let us simulate them with pymc. \n",
    "\n",
    "The key is to put the mean value function into the Deterministic symbol, then pymc unstands that it is supposed to be a variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape.0\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "zero-dimensional arrays cannot be concatenated",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-d09398a614f1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     17\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mn1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m         \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m         \u001b[0mrate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mswitch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mNswitch\u001b[0m\u001b[0;34m>=\u001b[0m \u001b[0midx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr_\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mr_\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m     \u001b[0;31m# Likelihood (sampling distribution) of observations\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/fnj/myPy3Env/lib/python3.4/site-packages/numpy/lib/index_tricks.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m    336\u001b[0m                 \u001b[0mobjs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mobjs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfinal_dtype\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    337\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 338\u001b[0;31m         \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_nx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobjs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    339\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_retval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mres\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    340\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: zero-dimensional arrays cannot be concatenated"
     ]
    }
   ],
   "source": [
    "basic_model = pm.Model()\n",
    "\n",
    "\n",
    "with basic_model:\n",
    "    # Priors for unknown model parameters\n",
    "    sigma = pm.HalfNormal('sigma', sd=1)\n",
    "    Omega = pm.Normal('omega', mu=0.55, sd=0.1)\n",
    "    amp = pm.Normal('Amplitude', mu=0.55, sd=0.1)\n",
    "    \n",
    "    \n",
    "    # Expected value of outcome\n",
    "    n1 = amp*pm.sin(2*np.pi*Omega*t_1)**2\n",
    "    n2 = 1-n1\n",
    "    Nval = len(nt_exp)\n",
    "    Nswitch = len(n1_exp)\n",
    "    idx = np.arange(Nval)\n",
    "    if n1.shape:\n",
    "        print(n1.shape)\n",
    "        rate = pm.switch(Nswitch>= idx, np.r_[n1, n1], np.r_[n2, n2])\n",
    "    \n",
    "    # Likelihood (sampling distribution) of observations\n",
    "    Y_obs = pm.Normal('Y_obs', mu=rate, sd=sigma, observed=nt_exp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now sample it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Nsamples  =5000\n",
    "with basic_model:\n",
    "    # obtain starting values via MAP\n",
    "    start = pm.find_MAP()\n",
    "\n",
    "    # instantiate sampler\n",
    "    step = pm.NUTS(scaling=start)\n",
    "\n",
    "    # draw 500 posterior samples\n",
    "    trace = pm.sample(Nsamples, step, start=start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and here is the final result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.stats.mstats import mquantiles\n",
    "Nplot =200\n",
    "omega_sample = trace['omega']\n",
    "amp_sample = trace['Amplitude']\n",
    "ts = np.linspace(0, tmax,Nplot)\n",
    "n_om = np.zeros((Nsamples, Nplot))\n",
    "for jj in np.arange(Nsamples):\n",
    "    amp = amp_sample[jj]\n",
    "    omega = omega_sample[jj]\n",
    "    n_om[jj,:] = amp*np.sin(2*np.pi*omega*ts)**2\n",
    "\n",
    "n_av = n_om.mean(axis = 0)\n",
    "\n",
    "qs = mquantiles(n_om, [0.025, 0.975], axis=0)\n",
    "print(n_av.shape)\n",
    "plt.figure(1)\n",
    "plt.clf;\n",
    "plt.scatter(t_all,n1_exp)\n",
    "plt.fill_between(ts, *qs, alpha = .4)\n",
    "plt.plot(ts,n_av)\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('y')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pm.traceplot(trace);\n",
    "pm.summary(trace);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(n1_exp.shape)\n",
    "print(n2_exp.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}