{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyzing 9-year NANOGrav data\n",
    "\n",
    "In this tutorial we will use `enterprise` to analyze the [NANOGrav 9-year data release](https://data.nanograv.org) for a stochastic GW background. We will reproduce the power-law GWB limit from [this paper.](http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:1508.03024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "% matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "%load_ext line_profiler\n",
    "\n",
    "from __future__ import division\n",
    "\n",
    "import numpy as np\n",
    "import glob\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.linalg as sl\n",
    "\n",
    "import enterprise\n",
    "from enterprise.pulsar import Pulsar\n",
    "import enterprise.signals.parameter as parameter\n",
    "from enterprise.signals import utils\n",
    "from enterprise.signals import signal_base\n",
    "from enterprise.signals import selections\n",
    "from enterprise.signals.selections import Selection\n",
    "from enterprise.signals import white_signals\n",
    "from enterprise.signals import gp_signals\n",
    "from enterprise.signals import deterministic_signals\n",
    "\n",
    "import corner\n",
    "from PTMCMCSampler.PTMCMCSampler import PTSampler as ptmcmc\n",
    "\n",
    "datadir = enterprise.__path__[0] + '/datafiles/ng9/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to convert PAL2 noise parameters to enterprise parameter dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_noise_from_pal2(noisefile):\n",
    "    psrname = noisefile.split('/')[-1].split('_noise.txt')[0]\n",
    "    fin = open(noisefile, 'r')\n",
    "    lines = fin.readlines()\n",
    "    params = {}\n",
    "    for line in lines:\n",
    "        ln = line.split()\n",
    "        if 'efac' in line:\n",
    "            par = 'efac'\n",
    "            flag = ln[0].split('efac-')[-1]\n",
    "        elif 'equad' in line:\n",
    "            par = 'log10_equad'\n",
    "            flag = ln[0].split('equad-')[-1]\n",
    "        elif 'jitter_q' in line:\n",
    "            par = 'log10_ecorr'\n",
    "            flag = ln[0].split('jitter_q-')[-1]\n",
    "        elif 'RN-Amplitude' in line:\n",
    "            par = 'log10_A'\n",
    "            flag = ''\n",
    "        elif 'RN-spectral-index' in line:\n",
    "            par = 'gamma'\n",
    "            flag = ''\n",
    "        else:\n",
    "            break\n",
    "        if flag:\n",
    "            name = [psrname, flag, par]\n",
    "        else:\n",
    "            name = [psrname, par]\n",
    "        pname = '_'.join(name)\n",
    "        params.update({pname: float(ln[1])})\n",
    "    return params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get par, tim, and noise files\n",
    "Here we collect the tim and par files as well as noise files made from the `PAL2` code. These are the same par, tim, and noise files used in the 9-year analysis papers. We use the convienience function above to convert from `PAL2` noise files to `enterprise` parameter dictionaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "parfiles = sorted(glob.glob(datadir + '/*.par'))\n",
    "timfiles = sorted(glob.glob(datadir + '/*.tim'))\n",
    "noisefiles = sorted(glob.glob(datadir + '/*noise.txt'))\n",
    "\n",
    "# 18 pulsars used in 9 year analysis\n",
    "p9 = np.loadtxt(datadir+'/9yr_pulsars.txt', dtype='S42')\n",
    "\n",
    "# filter\n",
    "parfiles = [x for x in parfiles if x.split('/')[-1].split('_')[0] in p9]\n",
    "timfiles = [x for x in timfiles if x.split('/')[-1].split('_')[0] in p9]\n",
    "noisefiles = [x for x in noisefiles if x.split('/')[-1].split('_')[0] in p9]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load into Pulsar class list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psrs = []\n",
    "for p, t in zip(parfiles, timfiles):\n",
    "    psr = Pulsar(p, t, ephem='DE421')\n",
    "    psrs.append(psr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get parameter dict from noisefiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {}\n",
    "for nfile in noisefiles:\n",
    "    params.update(get_noise_from_pal2(nfile))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up model\n",
    "\n",
    "When setting up the model for our upper limit run we fix all of the white noise (EFAC, EQUAD, and ECORR) parameters to the values obtained from the noise files. This is done by using `Constant` parameters. In this case we do not specify a default value for all instances of that parameter but instead will set them, based on their initialized pulsar and backend specific name, later via the `set_default_params` method of `PTA`. \n",
    "\n",
    "Speaking of white noise parameters here, we also use the `Selection` object.\n",
    "\n",
    "Another feature to notice is that we do not use a uniform prior on the log of the red noise or GWB amplitude. Instead we use a `LinearExp` prior (short for linear-exponent prior), that is a prior of the form $p(x)\\propto 10^x$. This is how we can still use the log of the parameter to sample but place a uniform prior on the parameter itself. We do this for both the red noise and GWB amplitude parameters.\n",
    "\n",
    "Next, in order to save on computing time we do not include spatial correlations here. Instead we model the GWB as a common red process across all pulsars. In `enterprise` we can do this with a simple trick. We pre-initialize the parameters before passing them to the `Signal` model. In this way the *same* parameter instance is used for all pulsars. Lastly, we fixt the spectral index of the GWB to be 13/3 (4.33) using the `Constant` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the maximum time span to set GW frequency sampling\n",
    "tmin = [p.toas.min() for p in psrs]\n",
    "tmax = [p.toas.max() for p in psrs]\n",
    "Tspan = np.max(tmax) - np.min(tmin)\n",
    "\n",
    "# selection class to break white noise by backend\n",
    "selection = selections.Selection(selections.by_backend)\n",
    "\n",
    "##### parameters and priors #####\n",
    "\n",
    "# white noise parameters\n",
    "# since we are fixing these to values from the noise file we set\n",
    "# them as constant parameters\n",
    "efac = parameter.Constant()\n",
    "equad = parameter.Constant()\n",
    "ecorr = parameter.Constant()\n",
    "\n",
    "# red noise parameters \n",
    "log10_A = parameter.LinearExp(-20,-12)\n",
    "gamma = parameter.Uniform(0,7)\n",
    "\n",
    "# GW parameters (initialize with names here to use parameters in common across pulsars)\n",
    "log10_A_gw = parameter.LinearExp(-18,-12)('log10_A_gw')\n",
    "gamma_gw = parameter.Constant(4.33)('gamma_gw')\n",
    "\n",
    "##### Set up signals #####\n",
    "\n",
    "# white noise\n",
    "ef = white_signals.MeasurementNoise(efac=efac, selection=selection)\n",
    "eq = white_signals.EquadNoise(log10_equad=equad, selection=selection)\n",
    "ec = white_signals.EcorrKernelNoise(log10_ecorr=ecorr, selection=selection)\n",
    "\n",
    "# red noise (powerlaw with 30 frequencies)\n",
    "pl = utils.powerlaw(log10_A=log10_A, gamma=gamma)\n",
    "rn = gp_signals.FourierBasisGP(spectrum=pl, components=30, Tspan=Tspan)\n",
    "\n",
    "# gwb (no spatial correlations)\n",
    "cpl = utils.powerlaw(log10_A=log10_A_gw, gamma=gamma_gw)\n",
    "gw = gp_signals.FourierBasisGP(spectrum=cpl, components=30, Tspan=Tspan)\n",
    "\n",
    "# for spatial correltions you can do...\n",
    "#orf = utils.hd_orf()\n",
    "#crn = gp_signals.FourierBasisCommonGP(cpl, orf, components=30, name='gw', Tspan=Tspan)\n",
    "\n",
    "# timing model\n",
    "tm = gp_signals.TimingModel()\n",
    "\n",
    "# to add solar system ephemeris modeling...\n",
    "#eph = deterministic_signals.PhysicalEphemerisSignal(use_epoch_toas=True)\n",
    "\n",
    "# full model is sum of components\n",
    "model = ef + eq + ec + rn + tm + gw\n",
    "\n",
    "# intialize PTA\n",
    "pta = signal_base.PTA([model(psr) for psr in psrs])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set white noise parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pta.set_default_params(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set initial parameters drawn from prior and evaluate likelihood to fill caches\n",
    "\n",
    "Evaluating the likelihood is not necessary, the caches will be filled the first time it is called within the sampler if not called here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1396202.32558\n",
      "-32.2501201076\n"
     ]
    }
   ],
   "source": [
    "xs = {par.name: par.sample() for par in pta.params}\n",
    "print pta.get_lnlikelihood(xs);\n",
    "print pta.get_lnprior(xs);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dimension of parameter space\n",
    "ndim = len(xs)\n",
    "\n",
    "# initial jump covariance matrix\n",
    "cov = np.diag(np.ones(ndim) * 0.01**2)\n",
    "\n",
    "# set up jump groups by red noise groups \n",
    "ndim = len(xs)\n",
    "groups  = [range(0, ndim)]\n",
    "groups.extend(map(list, zip(range(0,ndim,2), range(1,ndim,2))))\n",
    "groups.extend([[36]])\n",
    "\n",
    "sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, groups=groups, \n",
    "                 outDir='chains/nano_9_gwb/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sampler for N steps\n",
    "N = 1000000\n",
    "x0 = np.hstack(p.sample() for p in pta.params)\n",
    "sampler.sample(x0, N, SCAMweight=30, AMweight=15, DEweight=50, )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "chain = np.loadtxt('chains/nano_9_gwb/chain_1.txt)\n",
    "pars = sorted(xs.keys())\n",
    "burn = int(0.25 * chain.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAusAAAIQCAYAAAA4rCDnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAWJQAAFiUBSVIk8AAAIABJREFUeJzt3Xu8b1VdL/zPl1BDhI1YlLfaSlw8x8yHTUpRKnIiCs/J\n7G7e6IpaqOV5Dh5MxdJ4dQVNX9lFt9lJy07q00MWhWgqUUfQeqqNILrzgpciEAQUjfH8MedPfizX\n2ntd5l5rrLXe79drveZec8w5fmMN5mJ9fuM35pjVWgsAANCfgza6AQAAwOKEdQAA6JSwDgAAnRLW\nAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgH\nAIBOCesAANApYR0AADp18EY3YD1V1YeSHJ5k7wY3BQCArW1nkptaaw9aSyXbKqwnOfyQQw458iEP\neciRG90QAAC2rj179uS2225bcz3bLazvfchDHnLkFVdcsdHtAABgC9u1a1euvPLKvWutx5x1AADo\nlLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBT\nwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBApw7e6AYAAMzbec5Fa65j7/lnTNAS2HhG1gEAoFPC\nOgAAdEpYBwCATpmzDgB0bTnzz6eY5w49MrIOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAA\nOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADo\nlLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBT\nwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTk4T1qvreqnp5Vb2z\nqm6qqlZVv7/Kuh5QVa+uquuq6nNVtbeqLqiqe0/RVgAA2CwOnqie5yf5hiSfSfLRJMevppKqOjrJ\nZUmOSvKWJFcleUSSZyU5vapObq1dP0mLAQCgc1NNg3lOkmOTHJ7k6Wuo55UZgvrZrbXHt9bOaa09\nNsmvJzkuyUvW3FIAANgkJhlZb61dOvt3Va2qjnFU/bQke5O8YkHxC5P8RJInV9XPttZuWV1LAYDt\nYOc5F63o+L3nn3GAWgJr09MNpqeM24tba3fMF7TWbk7y7iT3THLSejcMAAA2wlRz1qdw3Li9eony\nazKMvB+b5JJ9VVRVVyxRtKq59AAAsBF6GlnfMW4/vUT5bP8R69AWAADYcD2NrE+mtbZrsf3jiPsJ\n69wcAOAAW+mc85XOaYeN0tPI+mzkfMcS5bP9N65DWwAAYMP1FNbfP26PXaL8mHG71Jx2AADYUnoK\n67PlH0+rqru0q6oOS3JykluTXL7eDQMAgI2w7mG9qu5WVceP66p/UWvt2iQXJ9mZ5JkLTjsvyaFJ\nXmeNdQAAtotJbjCtqscnefz47VeP22+qqt3jv/+ttfbc8d/3T7Inyb9kCObznpHksiQvq6pTx+Me\nmWEN9quTnDtFewEAYDOYajWYhyd56oJ9Dx6/kiGYPzf70Vq7tqpOTPLiJKcn+c4kH09yYZLzWms3\nTNReAADo3iRhvbX2oiQvWuaxe5PUPso/kuTMKdoFAACbWU83mAIAAHOEdQAA6JSwDgAAnRLWAQCg\nU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6NTBG90AAGBz2XnORRvdBNg2jKwDAECn\nhHUAAOiUsA4AAJ0yZx0AWJO955+x0U2ALcvIOgAAdEpYBwCATgnrAADQKWEdAAA6JawDAECnhHUA\nAOiUsA4AAJ2yzjoAsO3tPOeiFR1vbXnWi5F1AADolLAOAACdEtYBAKBT5qwDANvOSuecr3ROO0zF\nyDoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J\n6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHTq4I1uAADAZrPznItWdPze8884\nQC1hqzOyDgAAnRLWAQCgU8I6AAB0ypx1AID9WOmc85XOaYelCOsAsM0JltAv02AAAKBTwjoAAHRK\nWAcAgE6Zsw4A3IUH+EA/jKwDAECnhHUAAOiUsA4AAJ0S1gEAoFPCOgAAdGqysF5VD6iqV1fVdVX1\nuaraW1UXVNW9V1jPt1TVW8bzP1tVH66qP6uq06dqKwAAbAaThPWqOjrJFUnOTPJ3SX49yQeTPCvJ\n31TVfZZZz9OTvDPJqeP215O8I8mjk7y1qs6dor0AALAZTLXO+iuTHJXk7Nbay2c7q+rXkjwnyUuS\nnLWvCqrqbkl+Mclnk+xqrb1/ruylSd6b5Nyq+pXW2ucmajcAAHRrzSPr46j6aUn2JnnFguIXJrkl\nyZOr6tD9VHVkkh1Jrp4P6knSWtuT5OokhyS511rbDAAAm8EU02BOGbcXt9bumC9ord2c5N1J7pnk\npP3U86kk/5rk2Ko6Zr6gqo5NckyS97XWrp+gzQAA0L0ppsEcN26vXqL8mgwj78cmuWSpSlprraqe\nmeT3k1xRVW9Kcl2S+yf57iT/lOQHl9OgqrpiiaLjl3M+AAD0YIqwvmPcfnqJ8tn+I/ZXUWvtjVV1\nXZLXJ3nKXNEnk7wmw02rAACwLXS1znpVPSnJX2VYCeYhGabPPCTDiPxvJHnDcuppre1a7CvJVQeo\n6QAAMLkpwvps5HzHEuWz/Tfuq5JxXvqrM0x3eXJr7arW2m2ttauSPDnD0pDfV1WPWXuTAQCgf1OE\n9dnKLccuUT67WXSpOe0zpyW5W5J3LHKj6h1J/nr8dtdqGgkAAJvNFGH90nF7WlXdpb6qOizJyUlu\nTXL5fuq5x7j9yiXKZ/tvX00jAQBgs1lzWG+tXZvk4iQ7kzxzQfF5SQ5N8rrW2i2znVV1fFUtXJnl\nneP2e6vqYfMFVfXwJN+bpCV521rbDAAAm8FUTzB9RpLLkrysqk5NsifJIzOswX51knMXHL9n3NZs\nR2vt76rqNUnOTPJ/xqUb/yXDm4DHJ7l7kgtaa/80UZsBAKBrk4T11tq1VXVikhcnOT3Jdyb5eJIL\nk5zXWrthmVX9aIa56U9L8u1JDktyU5J3Jfnt1tqyVoMBAICtYKqR9bTWPpJhVHw5x9YS+1uS3eMX\nAABsa12tsw4AANxJWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTkz0UCQDo\nw85zLtroJgATEdYBAA6wlb6B2nv+GQeoJWw2psEAAECnhHUAAOiUaTAAsMWZUgGbl7AOADCxlb5B\nclMwSzENBgAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWs\nAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAO\nAACdEtYBAKBTB290AwCAfdt5zkUb3QRggxhZBwCATgnrAADQKWEdAAA6Zc46AGwye88/Y6ObAKwT\nI+sAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADol\nrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKcmC+tV9YCq\nenVVXVdVn6uqvVV1QVXdexV1nVBVf1BVHx3r+mRVvaOqnjJVewEAoHcHT1FJVR2d5LIkRyV5S5Kr\nkjwiybOSnF5VJ7fWrl9mXT+V5MIkNyS5KMnHkhyZ5KFJvjPJ703RZgAA6N0kYT3JKzME9bNbay+f\n7ayqX0vynCQvSXLW/iqpqtOSvCzJXyb53tbazQvK7zZRewEAoHtrngYzjqqflmRvklcsKH5hkluS\nPLmqDl1Gdb+c5LYkT1wY1JOktfb5tbUWAAA2jylG1k8Ztxe31u6YL2it3VxV784Q5k9KcslSlVTV\nQ5M8LMmbk/x7VZ2SZFeSluR9SS5dWD8AbDY7z7loo5sAbCJThPXjxu3VS5RfkyGsH5t9hPUk3zhu\nP5Xk7UketaD8/6uqJ7TWPrDKdgIAwKYyRVjfMW4/vUT5bP8R+6nnqHH7oxluKj0jybuSfFWSFyR5\nUpKLqurrW2u376uiqrpiiaLj99MGAADoRk/rrM/a8mVJfrC19mettZtaa9ckeUqS92QYnf+ejWog\nAACspylG1mcj5zuWKJ/tv3E/9czKP9Fa+5v5gtZaq6q3JDkxw5KQr99XRa21XYvtH0fcT9hPOwBg\n3ew9/4yNbgLQsSlG1t8/bo9dovyYcbvUnPaF9SwV6m8Yt4css10AALCpTRHWLx23p1XVXeqrqsOS\nnJzk1iSX76eeyzMs87hziWUeHzpuP7SGtgIAwKax5rDeWrs2ycVJdiZ55oLi85IcmuR1rbVbZjur\n6viqusvNnq21W5P8bpIvT/ILVVVzx399kqcl+UKSP15rmwEAYDOY6gmmz0hyWZKXVdWpSfYkeWSG\nNdivTnLuguP3jNtasP/nMizZ+Owk3zSu0f5VSZ6QIcQ/e3xzAAAAW94kq8GMAfrEJLszhPSfTXJ0\nkguTnNRau36Z9dyU5FuTvDTJkUl+KsnjMizh+O2ttQunaC8AAGwGU42sp7X2kSRnLvPYhSPq82Wf\nyTASv3A0HgAAtpWe1lkHAADmCOsAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0\nSlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnTp4oxsAAJvZznMu2ugmAFuYkXUAAOiU\nsA4AAJ0S1gEAoFPmrAPAhPaef8ZGNwHYQoysAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcA\ngE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTwjoAAHRKWAcAgE4J6wAA0KmDN7oB\nAADc1c5zLlrxOXvPP+MAtISNZmQdAAA6JawDAECnhHUAAOiUOesAABtsNfPNVzOvnc3HyDoAAHRK\nWAcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACdEtYBAKBTHooEAHM8aAboiZF1AADolLAO\nAACdEtYBAKBT5qwDwD7sPf+MjW4CsI0ZWQcAgE4J6wAA0ClhHQAAOiWsAwBAp4R1AADolLAOAACd\nEtYBAKBTwjoAAHRqsrBeVQ+oqldX1XVV9bmq2ltVF1TVvddQ56Oq6j+qqlXVL0zVVgAA2AwmeYJp\nVR2d5LIkRyV5S5KrkjwiybOSnF5VJ7fWrl9hnYcleW2SW5Pca4p2AgDAZjLVyPorMwT1s1trj2+t\nndNae2ySX09yXJKXrKLOC5PsSPKLE7URAAA2lTWH9XFU/bQke5O8YkHxC5PckuTJVXXoCur8riRn\nJjk7yXVrbSMAAGxGU4ysnzJuL26t3TFf0Fq7Ocm7k9wzyUnLqayqjkry20ne3Fr7/QnaBwAAm9IU\nc9aPG7dXL1F+TYaR92OTXLKM+n47w5uIs1bboKq6Yomi41dbJwAArLcpwvqOcfvpJcpn+4/YX0VV\n9SNJ/luSH2itfXKCtgEAwKY1yWowU6iqnUkuSPLG1tofraWu1tquJV7jiiQnrKVuAABYL1PMWZ+N\nnO9Yony2/8b91PPqJLclecYEbQIAgE1virD+/nF77BLlx4zbpea0z5yQYfnHfx0fgtSqqiV5zVh+\n7rjvzWtrLgAAbA5TTIO5dNyeVlUHza8IMz7Y6OQMDza6fD/1/F6GVWMWOibJo5K8L8kVSd675hYD\nAMAmsOaw3lq7tqouzrDiyzOTvHyu+LwkhyZ5VWvtltnOqjp+PPequXrOXqz+qnpahrB+UWvt+Wtt\nLwAAbBZT3WD6jCSXJXlZVZ2aZE+SR2ZYg/3qJOcuOH7PuK2JXh8AALacKeasp7V2bZITk+zOENJ/\nNsnRSS5MclJr7fopXgcAALaTyZZubK19JMmZyzx22SPqrbXdGd4EAADAtjLJyDoAADA9YR0AADol\nrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSw\nDgAAnTp4oxsAAAfSznMu2ugmAKyakXUAAOiUsA4AAJ0S1gEAoFPmrAOwrew9/4yNbgLAsgnrAABb\nwEpvpvbGdXMwDQYAADolrAMAQKeEdQAA6JQ56wAAm9BK55x7QNjmZGQdAAA6ZWQdgE3F6CCwnRhZ\nBwCATgnrAADQKWEdAAA6Zc46AJuapzACW5mRdQAA6JSwDgAAnRLWAQCgU8I6AAB0SlgHAIBOCesA\nANApYR0AADolrAMAQKc8FAmADbXznIs2ugkA3TKyDgAAnRLWAQCgU8I6AAB0ypx1ALqy9/wzNroJ\nAN0wsg4AAJ0S1gEAoFPCOgAAdEpYBwCATgnrAADQKWEdAAA6JawDAECnhHUAAOiUsA4AAJ0S1gEA\noFPCOgAAdEpYBwCATgnrAADQKWEdAAA6JawDAECnhHUAAOiUsA4AAJ06eKqKquoBSV6c5PQk90ny\n8SRvTnJea+2GZZx/aJLHJzkjyQlJHpjkjiTvT/L6JC9vrd0+VXsBmN7Ocy7a6CYAbCmThPWqOjrJ\nZUmOSvKWJFcleUSSZyU5vapObq1dv59qvjXJ7yf59ySXZgj6907y35L8SpInVNWprbXPTtFmAADo\n3VQj66/MENTPbq29fLazqn4tyXOSvCTJWfup4xNJnpTkjfMj6FX13CRvT/LNSZ6Z5FcnajMAAHRt\nzXPWx1H105LsTfKKBcUvTHJLkieP01yW1Fp7X2vtfy2c6tJauzl3BvTHrLW9AACwWUwxsn7KuL24\ntXbHfEFr7eaqeneGMH9SkktW+RqfH7dfWOX5AGyAveefsdFNAJaw0ntM/D5vjCnC+nHj9uolyq/J\nENaPzerD+o+M2z9fzsFVdcUSRcev8vUBAGDdTbF0445x++klymf7j1hN5VX1UxlWmHlfklevpg4A\nANiMJlu68UCoqickuSDDzaff01r7/H5OSZK01nYtUd8VGZaFBACA7k0R1mcj5zuWKJ/tv3EllVbV\n45O8IcmnkpzSWvvg6poHAMBK55x7bkIfppgG8/5xe+wS5ceM26XmtH+Jqvq+JG9M8skkj26tvX8/\npwAAwJYzRVi/dNyeVlV3qa+qDktycpJbk1y+nMqq6oczPLH0ugxB/ZoJ2ggAAJvOmsN6a+3aJBcn\n2ZnhoUXzzktyaJLXtdZume2squOr6ktWZqmqpyb5vSQfTvIoU18AANjOprrB9BlJLkvysqo6Ncme\nJI/MsAb71UnOXXD8nnFbsx1VdUqG1V4OyjBaf2ZVLTgtN7bWLpiozQAA0LVJwnpr7dqqOjHJizMs\ns/idST6e5MIk57XWblhGNV+bO0f6f2SJY/4lw+owAACw5U22dGNr7SNJzlzmsV8yZN5a251k91Tt\nAWDtrAYBsLGmuMEUAAA4AIR1AADoVNdPMAVgWqa1AGwuwjoAy7bSJyACsDamwQAAQKeMrANsYqa1\nAGxtwjrANmZaC0DfTIMBAIBOCesAANAp02AAthDTWgC2FiPrAADQKWEdAAA6JawDAECnhHUAAOiU\nsA4AAJ0S1gEAoFPCOgAAdEpYBwCATgnrAADQKWEdAAA6JawDAECnhHUAAOiUsA4AAJ06eKMbAMCd\ndp5z0UY3AYCOGFkHAIBOCesAANApYR0AADplzjpAx/aef8ZGNwGADWRkHQAAOiWsAwBAp4R1AADo\nlLAOAACdcoMpwAHkIUcArIWwDrBMgjcA6800GAAA6JSRdQAA9mulny56TsQ0hHWAdeSPFwArIawD\nrJLgDcCBJqwD25YbRgHonbAOAMCXWOmnhwZADgxhHdgy/KEAYKuxdCMAAHTKyDrAyA2jAPTGyDoA\nAHTKyDqwZRkpB2CzM7IOAACdMrIOdMvqLgBsd8I6sG6Eb4DtY6X/zzd1cXHCOrAqgjcAHHjmrAMA\nQKeMrAObho9IAdhuhHVgwwjfAFvHSv+fbjrl8gjrwCQEbwCYnjnrAADQKSPrsEX5eBEANj8j6wAA\n0ClhHQAAOmUaDGwSprUAwPYjrMNE1hqmD/RqKlZrAaBnK/07ul3+rk0W1qvqAUlenOT0JPdJ8vEk\nb05yXmvthhXUc2SSFyR5fJL7Jrk+yZ8neUFr7aNTtZfNbb1HmdfjfwhGzgGAhSYJ61V1dJLLkhyV\n5C1JrkryiCTPSnJ6VZ3cWrt+GfXcZ6zn2CRvS/KGJMcnOTPJGVX1Ta21D07RZg6s3keZV0qQBgA2\nwlQj66/MENTPbq29fLazqn4tyXOSvCTJWcuo56UZgvqvtdZ+dq6es5NcOL7O6RO1ed30Ngo8RXt6\nC9Pbkf8GALD1VWttbRUMo+ofSLI3ydGttTvmyg7LMB2mkhzVWrtlH/XcK8mnktyR5L6ttZvnyg5K\n8sEkXzu+xqpG16vqihNOOOGEK664YjWnr9pWDOvs34H+7yCsA7CdbLbBxl27duXKK6+8srW2ay31\nTDGyfsq4vXg+qCdJa+3mqnp3ktOSnJTkkn3Uc1KSQ8Z6bp4vaK3dUVV/keQnxtczFWYftmMYn/qX\nbz36UNgGAPZnirB+3Li9eonyazKE9WOz77C+nHoy1sMms9lG+wVpAKAHU4T1HeP200uUz/YfsU71\npKqWmufyDXv27MmuXWv6NGLFPv6xu/5ID73/jiWOXJ1//NhSXbY8y2nPWl9j11++YJ/lC/to6voB\ngM1trVkhWd+8sGfPniTZudZ6tts66/9x2223ffrKK6/cu86ve/y4vSpJrvzkOr/6fqxHeyZ+jbv0\n5wGof7v5kv5kzfTptPTntPTn9PTptA5Yf65zXtiZ5Ka1VjJFWJ+9zVlqeHa2/8Z1qidrncg/tdlI\nf2/t2qz057T05/T06bT057T05/T06bT0510dNEEd7x+3S80lP2bcLjUXfep6AABgS5girF86bk8b\nl1j8onHpxpOT3Jrk8v3Uc3mS25KcPJ43X89BGW5SnX89AADY0tYc1ltr1ya5OMO8nGcuKD4vyaFJ\nXje/xnpVHV9Vx88f2Fr7TJLXjce/aEE9PzXW/xeeYAoAwHYx1Q2mz0hyWZKXVdWpSfYkeWSGNdGv\nTnLuguP3jNtasP9/JnlMkp+pqocn+bskD0nyXRkemLTwzQAAAGxZU0yDmY2un5hkd4aQ/rNJjk5y\nYZKTWmvXL7Oe65N8U5KXJfm6sZ5HJnlNkl3j6wAAwLZQrbWNbgMAALCISUbWAQCA6QnrAADQKWEd\nAAA6JawDAECnhHUAAOiUsA4AAJ0S1gEAoFPC+gpV1d2q6llV9Zqqel9V3V5Vrap+bD/nHVVVv1RV\n/1hVN1fV9VV1RVX996o6bBXteGpV/V1VfaaqPl1Vb6+qx63+J9s4q+nTqto7HrOvr59b5us/bT/1\nnDXdT3vgbXR/ztW5Ja7R1f7Oj+fuqKoXV9U/jP1w0/j/gFdV1d2W+frb/vqcO3fN/TlX15a4PpNV\n/85Pdl25Rg9MH2yVa3Qtv/Nzddxj/F1vVfXRFb7+i/bz3+b0lf9U6+vgjW7AJnRokgvGf38yySeS\nPHBfJ1TVziR/m+SoJG9P8tYkX57ktCS/lORJVXVSa+225TSgqn4lw9NdP5rkt5PcPckPJvnTqvrp\n1tpvrOgn2ngr7tPx+CMW2V9J/meGa/utK2zHW5K8b5H971lhPRttw/tzi12jq+nPVNXxSS5Ocv8k\nf5Wh/+6WZGeS78/QP59fQTu28/U5aX9useszWWWfjqa8rrb1NTqapA+22DW6lv6ceWmSr11jO16b\nZO8i+z+wxnoPvNaarxV8ZfiF+Y4k9x2/f1GSluTH9nHOK8ZjXrhg/5cluWQse8oyX/+bx+M/kOTe\nc/t3Jrk+yWeT7NzofjrQfbqPur59PPfKFZzztPGcp210X2yR/txS1+gqf+fvmeTqJDckOWmR8oMz\nPkF6Ga+/7a/PiftzS12fa+jTya4r1+jk/bmlrtG1/k1K8pgkdyQ5azzvoyt8/dnrPWaj+2K1X6bB\nrFBr7fbW2ltbax9fwWkPHrf/z4K6/iPJReO3X7nMumYfpb2ktXbDXF17M7wpuEeSM1fQtg23yj5d\nyk+M21dNUNem1EF/bqlrdJX9eVaSY5I8r7V2+SJ1fqGNf0W2mw76c0tdn8nkv/PbXgf9uaWu0bX0\nZ1UdnmR3kktaa785eeM2CWF9ffzTuD1jfmdVHZTh3eYdSd62zLoeO27/fJGyty44Zlupqq9K8l+T\nfCbJH6yiiodX1bOr6pyqenJVPWDaFm4ua+hP12jyxAwjOW+oqp1V9fSqel5V/XBV3WeVdW7n63PK\n/nR93tWU19V2vkZnpugD1+idXpbk3kl+dIK6vqWqnltV/6OqfqCqvmKCOteFOevr45eSPC7Jz1fV\nKUmuzPCx0GlJvjrDR0Hv3V8lVXVohvman1niHeo14/bYSVq9+fxIhjmsu1trN6/i/Gct+P4/qup3\nkjy7tfbZNbdu81lxf7pGh5upknxDkn9N8uMZ5lrO/7/2lqo6u7X26hVWvS2vzyn70/W5qCmvq215\njS6wpj5wjd6pqr47yVMzZKQPT1Dlzy/4/nNV9ctJXtD7J51G1tdBa+1TSU5K8qYM74afm+TsJMcl\n+aMMN0stx45x++klymf7F7tRcEurqkoyu7P8t1Z4+oeS/HSG/x6HJrlfhhvW9ib5ySQrDVWb3hr6\n0zWaHJkhTN4nyS9m+APxwCRfkaFPW5Lfqarljoxt9+tzyv50fd5pyutqu1+jyXR94BrNFz/Z/a0k\nb22t/e4aq/v7DINPD05ySIYbVX88yY1Jnp/kJWus/8Db6EnzG/GV4ZenreDr9/dR14uy/xtPdib5\nhwzviL8d/jzGAAAMpElEQVQjyeEZRtR/MsMUg08kedAy2n2/7OPmigyjoC3J57Z6ny5yzreN51wx\n4c/0wCT/Ptb7DfpzWed1eY2uZ38mue9cPb+5SPlPj2V/4fpc3/7s9fpc7z5dr+tqu1yjU/ZBr9fo\nevdnhnv8bkhyvwX7l+ybVfxMJyS5ffz6ivXsz5V+bddpMNdmuJt6ua5b4+vtTvL1GX5R/2Hcd1OS\nV1XVl2dY0uiFGe4m35fZO+odS5TP9t+46pau3nr36UKzGyFXOqq+pNbaR6rqz5L8cJJHZXh3vl42\na3/2eo2uZ3/Oj4i9aZHyN2WYh/mINbzGdro+p+zPXq/PZON/55NMe11to2t0Savog16v0XXrz6p6\nSob7pZ7aWjsg/12SpLV2ZVX9XZKTk3xTkj89UK+1VtsyrLfWTl2v16rhgUePTvLvc0F93qXjdtf+\n6mqt3VJVH0ty/6q6b/vS+WzHjNurV93gVVrPPl2oqo5K8l1Z/Y2l+/Kv4/bQievdp83an71eo+vZ\nn621W6vqIxlG1Bb7gzpb3eGQCV5uy1+fU/Znr9dnsrG/84uY8rra8tfoMiy7D3q9Rte5P08Yt6+t\nqtcuUn7/qmrjv+/dWlvLG5cNuT5Xypz1A+/u4/bwqrr7IuWzJRtvX2Z9s1VjFnvi1ncsOGa7ODPD\nR4Ovb6u7sXRfHjluPzhxvT1ba3+6Ru+8D+Whi5TN9n1ogtfZLtfnlP3p+ty/Ka+r7XKN7stK+2C7\nX6N/k+R3l/hKklvnvv/cal9kvHl99sag6+tTWD/AWmvXJ9mT4VOMuzyufZwC8/zx20sWlN23qo6v\nqoUfhc3WGT23qu49d/zOJM/McOG+Zqr2927BjZD7XAt8qT6tqhMXOfagqnpeho/G/i2LL6G15UzR\nn3GNJsNayHckOaeqvvgMhfF3fnYz0+vnT3B97tNk/RnXZ5LVXVeu0aVN2Z/Z5tdoa+0PW2s/ttjX\neMgNc/u++OT3qvqasT/vObfvsKo6buFrjIOnFyT5miRXpfOn7NY4yZ4VqKpzkhw/fvvwDMuKXZY7\nl1R6V2vtd+aO/y8ZHn509yR/Ox57SIZ3yF+b4SllJ43BfnbO7gxLFp3ZWtu94PV/NcnPZHgM8R+P\n9f5AhtUSNttjiJOsvE/nzjs1w6jbla21fU4lWqpPx4/T/jHDXMKPZZgTeHKGEbtbk3x3a+3i1f5s\nG2Ej+3Ms21LX6Gr6s6pekOS8JJ/KcLPUZzM8EfaY8dxT29xSbq7P9enPsWxLXZ/Jqv4urfi6co2u\nT3+OZVvqGl3t36RF6mlJPtZa+5L166vq7RmmHZ/SWnv7uG9nhlHz92QYOP14hhkNpyR5UIY3Ud/W\nWnvf6n6ydbLRd7huxq8kb8++74Levcg5D0vyuiQfzjDl5bYMD0t6aZIjFjl+91jX05Zow9OS/J8k\ntyS5Ock7kjxuo/tmPft0PO8Px/KfXMZrLNqnSX557L/rMgSAWzO80/6NJA/e6L7ZbP25Fa/RNfTn\nE5L8dYYbyj87/s6fm+Qers+N6c+teH2upk9Xc125RtenP7fiNbra3/lF6llyNZi513jM3L7DM9yA\nfnmGlfduz3A/1t8nOT/JURvdN8v5MrIOAACdMmcdAAA6JawDAECnhHUAAOiUsA4AAJ0S1gEAoFPC\nOgAAdEpYBwCATgnrAADQKWEdAAA6JawDAECnhHUAAOiUsA6wRlW1s6paVe3e6LYAsLUI6wCbVFU9\noKrOrao3VtUHquqO8U3D1+3nvEOq6ryqen9VfbaqPlVVf1RVD5mwbSePbWlV9RNT1Quw3QjrAJvX\niUl+Icn3JKkkn97fCVV1jyR/meQFSW5KcmGSv0ry3UneU1WPnKhts4De5v4NwAoJ6wCb13uSPCrJ\nEa21o5P8/TLO+ZkkJyf54ySPbK39j9baE5N8b5J7Jnl1Va3pb0NVHZHk+5Jck+RPkuyqqv9rLXUC\nbFfCOsABUlX3rapXVNXeqrq9qv61qv6kqnYtcfyOqrqgqj46Tk+5qqp+pqoevNic+NbaR1tr72yt\n3bTM9lSSs8Zv/+/W2h1zdb0lyTuT/Kckj17NzzvnSUkOSbJ7/EomGF2vwbOq6p/H/vlYVf3G2G97\nq2rv3LHfPvbZSxbUccrc9JwHLij7w3H/g9faVoCpCOsAB0BVPSjDyPczklyb5FeT/EWSM5JcVlWP\nW3D8lyd5W5JnJflUhukpb09y7njuFI5O8jVJrm6tfWiR8reO28eu8XV+PMkdSX4vyZ8n+USSJ1bV\noWus9xVJLkiyI8lvJXl9ktMyTOu524Jj35nk9iSnLth/6mL/Ht/InJJkb2vtg2tsJ8BkhHWAA+M3\nk9wvyfNba6e21p7XWntShkD4ZUleW1X3mjv+vyc5Ickbkuwap6ecleThSb5lojYdN26vXqL8mnF7\n7GpfoKpOSvKwJH81jvx/Icn/SnJ4kh9cQ73fmuTpGdr+n1trZ7fWnpvkoUk+m6Gvv6i1dmuSv01y\nYlXtmCs6Ncl7k1yfuwb3hyX5ygxvmAC6IawDTKyqHpBhxPfDSX5pvqy1dlmGEeEjkzxhruipGUaj\nn9daa3PHfyTDaPIUZqF1qRtRZ/uPWMNr/Pi43T23b/bvtUyFeeq4fUlr7cbZztba7Umet8Q5l2R4\nY/ToJKmqwzLclPuXSS7NXT9BOHXuHIBuCOsA05vdTPnO1trnFyl/2/xxVXV4hikqH2ut7V3k+HdN\n3sIDYPw5fiDJjUneNNvfWvvHJFckeURVPWyV1c/6dLG+uDzJFxbZP+vnWRB/dJKDMwTytyW539xy\nlY9dcA5AF4R1gOnNRrA/vkT5bP9sBPvwcfvJJY5fav9KzUbOdyxRPtt/4xLl+/PDSQ5N8oettc8u\nKNs9blc7uj5r25f0RWvtPzJMa1no8iS35M6wfmqGeezvyp0j6KdW1cEZVtX559baJ1bZPoADQlgH\nmN4sFH/1EuX3XXDcbDWXr1ri+KX2r9T7x+1Sc9KPGbdLzWnfn9kUmJ+cW3GlVVVL8vKx7Ier6pBV\n1L1kH1XVlyW5z8L946ca70ryn6vqqzOE9b9prd3aWrs6yUeT/Jckj0hyWIyqAx06eKMbALAFvXfc\nfktVHTzeZDnvlHF7ZZK01m6qqg8m2VlVOxeZCjPVDabXZphHf2xVPWiRFWG+Y9yuOLRW1YkZpqpc\nlztXlVnoGzPcyPn9SV67wpd471j/tyRZuFrLSVn679klSb49yQ9luBn1hXNlb0vyXbnzv5f56kB3\njKwDTKy19tEMNzHuTPLs+bLxCaFPTHJD5uZ1Z1jm8KAkvzguIzg7/oEL61hDu1qGVWqS5JfmH35U\nVd+V5FuT/HOSd6yi+tn0lgtbaz+22FeGBzLNH7sSvzduz51f3aWq7p7kpfs4b/bG45wMT3m9ZEHZ\njgzLa96RYalMgK7U3KIDAKxCVe1M8qEkr22tPW3c9+Ak784wFeYvM6y5/sAMT/Y8KMn3jQ8imtVx\nSJLLMizV+N4kF2cIkt+f5K+TPD7Ja1prP7LgtXfPfXt6hmkif5Lk5nHf77TW3jV3/D0yhNRvHtt0\nSYa1178vw3zux7bW/naFP/+9MszDv0eSB7bWFp1jP74J+UCSByd5aGvtn1b4Oq/KEPQ/luR/J/l8\nkv+aYTrR/ZN8rrX24AXnHJTk35LcO0OfHDn7pGNctecj46Hvaa1940raA7AejKwDHADjg3VOzDCS\nfVyS52aYZvLnSU6eD+rj8bdlmB7z8gwB/znj9y9N8ovjYYs9qfSpc1+z+dxPmNv3dQte53NJvi3J\nz2e4wfU54/dvTvKNKw3qox9Kcq8k/+9SQX187Zbkd8dvVzO6/vQMo/OfyfAk1icm+asM7T88i/TP\n+JTWS8dv3zk/JWn8BGQ2P998daBLRtYBOldVP57hiZ1ntdZetdHt6U1VHZMhdL+htfZDG90egCkZ\nWQfoRFXdb5F9X5Pk5zKsI/6n696ojlTVV8/Psx/33TN3PjTqTV96FsDmZjUYgH7876q6W4YHCN2Y\n4QbVxyW5Z4Ynm163gW3rwbOT/FBVvT3DHPnZcowPyLACzRs3rmkAB4ZpMACdqKpnJHlyhvXOd2SY\nm/3eJL/RWvuTDWjPi5Z56Jtba+9bYd07kzxtmYdf0Fq7sapOzTD3/+FJjszwacPVSf5gPGaxp8UC\nbGrCOgCLGh9mtBxnttZ2r7Dux+TOGz/350GLrD0PsC0I6wAA0Ck3mAIAQKeEdQAA6JSwDgAAnRLW\nAQCgU8I6AAB0SlgHAIBOCesAANApYR0AADolrAMAQKeEdQAA6JSwDgAAnRLWAQCgU8I6AAB06v8H\nXvsA7eQkeZQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1086fed10>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 264,
       "width": 373
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(chain[burn:,-5], 50, normed=True, histtype='step', lw=2);\n",
    "plt.xlabel(pars[-1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Upper limit value\n",
    "\n",
    "We see that the upper limit agrees perfectly with the published value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.49899289556e-15\n"
     ]
    }
   ],
   "source": [
    "upper = 10**np.percentile(chain[burn:, -5], q=0.95)\n",
    "print(upper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}