{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inference plots - Autocorrelation plot\n",
    "\n",
    "This example builds on the [adaptive covariance MCMC example](https://github.com/pints-team/pints/blob/master/examples/sampling-adaptive-covariance-mcmc.ipynb), and shows you a different way to plot the results.\n",
    "\n",
    "Inference plots:\n",
    "* [Predicted time series](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-predicted-time-series.ipynb)\n",
    "* [Trace plots](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-trace-plots.ipynb)\n",
    "* __Autocorrelation__\n",
    "* [Pairwise scatterplots](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-pairwise-scatterplots.ipynb)\n",
    "* [Pairwise scatterplots with KDE](https://github.com/pints-team/pints/blob/master/examples/plot-mcmc-pairwise-kde-plots.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up an MCMC routine\n",
    "\n",
    "See the adaptive covariance MCMC example for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import pints\n",
    "import pints.toy as toy\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load a forward model\n",
    "model = toy.LogisticModel()\n",
    "\n",
    "# Create some toy data\n",
    "real_parameters = [0.015, 500]\n",
    "times = np.linspace(0, 1000, 100)\n",
    "org_values = model.simulate(real_parameters, times)\n",
    "\n",
    "# Add noise\n",
    "noise = 50\n",
    "values = org_values + np.random.normal(0, noise, org_values.shape)\n",
    "real_parameters = np.array(real_parameters + [noise])\n",
    "\n",
    "# Get properties of the noise sample\n",
    "noise_sample_mean = np.mean(values - org_values)\n",
    "noise_sample_std = np.std(values - org_values)\n",
    "\n",
    "# Create an object with links to the model and time series\n",
    "problem = pints.SingleOutputProblem(model, times, values)\n",
    "\n",
    "# Create a log-likelihood function (adds an extra parameter!)\n",
    "log_likelihood = pints.GaussianLogLikelihood(problem)\n",
    "\n",
    "# Create a uniform log-prior over both the parameters and the new noise variable\n",
    "log_prior = pints.UniformLogPrior(\n",
    "    [0.01, 400, noise*0.1],\n",
    "    [0.02, 600, noise*100]\n",
    "    )\n",
    "\n",
    "# Create a posterior log-likelihood (log(likelihood * prior))\n",
    "log_posterior = pints.LogPosterior(log_likelihood, log_prior)\n",
    "\n",
    "# Perform sampling using MCMC, with a single chain\n",
    "x0 = real_parameters * 1.1\n",
    "mcmc = pints.MCMCController(log_posterior, 1, [x0])\n",
    "mcmc.set_max_iterations(6000)\n",
    "mcmc.set_log_to_screen(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Autocorrelation\n",
    "\n",
    "The [autocorrelation](https://en.wikipedia.org/wiki/Autocorrelation) in a Markov chain indicates how much each sample in the chain differs from the next. Checking for (lack of) autocorrelation is an easy way to check if your MCMC routine is converging. It can easily be plotted using Matplotlib's [acorr method](https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.acorr)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running...\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "print('Running...')\n",
    "chains = mcmc.run()\n",
    "print('Done!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select chain 0 and discard warm-up\n",
    "chain = chains[0]\n",
    "chain = chain[3000:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAE49JREFUeJzt3X+QH3d93/HnyzIOxTY2GQmGSLKltDLBpY1NDseBTHrENpFdsJIpAaslcShFJMHGMbSNIa1pnc6QHxRCU4eiJC6OQ2wUh6RqKlDAWDgwtUcydgDJUSNkgw+bWoAxBDcIwbt/fFfLl/PpbvVjb9Hd8zHznfvufj/f3feONPe6z352P5uqQpIkgBOGLkCS9N3DUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLLUJAktQwFSVLrxKELOFxLly6tVatWDV2GJB1X7r777i9U1bK52h13obBq1Sp27NgxdBmSdFxJ8pku7Tx9JElqGQqSpJahIElqGQqSpJahIElq9RYKSW5I8kiSTx3i8yT5L0n2JPlEkuf2VYskqZs+ewrvBtbO8vnFwJrmtQF4Z4+1SJI66C0UquoO4EuzNFkH/EGN3AmcnuSZfdUjSZrbkGMKy4EHx5anmnVPkGRDkh1Jduzbt++IdjY5Ocnk5OQRfVeSFoshQyEzrKuZGlbVxqqaqKqJZcvmvEtbknSEhgyFKWDl2PIK4KGBapEkMWwobAZ+trkK6Xzgsap6eMB6JGnR621CvCQ3A5PA0iRTwJuBJwFU1X8DtgCXAHuAx4FX9lWLJKmb3kKhqtbP8XkBr+1r/5Kkw+cdzZKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWoZCpKklqEgSWr1GgpJ1ibZnWRPkmtm+PyMJLcnuSfJJ5Jc0mc9kqTZ9RYKSZYA1wMXA2cD65OcPa3ZvwM2VdW5wGXA7/RVjyRpbn32FM4D9lTV3qraD9wCrJvWpoCnNu9PAx7qsR5J0hxO7HHby4EHx5angB+e1uY/AH+R5ErgZODCHuuRJM2hz55CZlhX05bXA++uqhXAJcBNSZ5QU5INSXYk2bFv374eSu3X5OQkk5OTQ5chSXPqMxSmgJVjyyt44umhVwGbAKrqfwNPBpZO31BVbayqiaqaWLZsWU/lSpL6DIXtwJokq5OcxGggefO0Np8FLgBI8mxGoXD8dQUkaYHoLRSq6gBwBbAVuI/RVUY7k1yX5NKm2RuAVyf5K+Bm4OeqavopJknSPOlzoJmq2gJsmbbu2rH3u4AX9FmDJKk772iWJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSy1CQJLUMBUlSq9PjOJMsB84cb19Vd/RVlCRpGHOGQpJfB14O7AK+2awuwFCQpAWmS0/hJ4FnVdXX+y5GkjSsLmMKe4En9V2IJGl4XXoKjwP3JrkNaHsLVfW63qqSJA2iSyhsbl6SpAVuzlCoqhuTnASc1azaXVXf6LcsSdIQulx9NAncCDwABFiZ5HIvSZWkhafL6aP/DLyoqnYDJDkLuBn4oT4LkyTNvy5XHz3pYCAAVNX/wauRJGlB6tJT2JHk94GbmuV/AdzdX0k6ViYnJwHYtm3boHVIOn50CYVfAF4LvI7RmMIdwO/0WZQkaRhdrj76OvC25iVJWsAOOaaQZFPz85NJPjH91WXjSdYm2Z1kT5JrDtHmZUl2JdmZ5I+O7DAkScfCbD2Fq5qfLz6SDSdZAlwPXARMAduTbK6qXWNt1gBvBF5QVY8mefqR7EuSdGwcsqdQVQ83b3+xqj4z/gJ+scO2zwP2VNXeqtoP3AKsm9bm1cD1VfVos89HDv8QJEnHSpdLUi+aYd3FHb63HHhwbHmqWTfuLOCsJB9LcmeStTNtKMmGJDuS7Ni3b1+HXUuSjsQhTx8l+QVGPYLvnzaGcCrwsQ7bzgzraob9rwEmgRXAXyZ5TlV9+Tu+VLUR2AgwMTExfRuSpGNktjGFPwLeD7wFGB8k/mpVfanDtqeAlWPLK4CHZmhzZzOX0v1JdjMKie0dti9JOsZmG1N4rKoeqKr1zTjC/2P0l/4pSc7osO3twJokq5sJ9S7jibOt/hnwQoAkSxmdTtp7BMchSToG5hxTSPKSJH8D3A98hNHEeO+f63tVdQC4AtgK3AdsqqqdSa5LcmnTbCvwxSS7gNuBf1NVXzyiI5EkHbUudzT/J+B84ENVdW6SFwLru2y8qrYAW6atu3bsfQGvb16SpIF1ufroG81f7yckOaGqbgfO6bkuSdIAuvQUvpzkFEZzHr0nySPAgX7LkiQNoUtPYR2jQeargQ8AnwZe0mdRkqRhdJkQ72tjizf2WIskaWCz3bz2Vb7zZrM0y2E0RvzUnmuTJM2zQ4ZCVZ06n4VIkobXZUyBJD+a5JXN+6VJVvdbliRpCF1uXnsz8MuMprgGOAn4wz6LkiQNo0tP4aeAS4GvAVTVQ4wmxZMkLTBdQmF/c+dxASQ5ud+SJElD6RIKm5K8Czg9yauBDwG/229ZkqQhdLlP4a1JLgK+AjwLuLaqPth7ZZKkeTdrKDTPWd5aVRcCBoE6m5ycBGDbtm2D1iHp8Mx6+qiqvgk8nuS0eapHkjSgLhPi/R3wySQfpLkCCaCqXtdbVZKkQXQJhf/VvCRJC1yXMYWLquoV81SPJGlAXcYUljXPWJYkLXBdTh89AHwsyWa+c0zhbX0VJUkaRpdQeKh5nYDTW0jSgtbl5rX/CJDk1NFi/W3vVUmSBtFlltTnJLkH+BSwM8ndSf5h/6VJkuZbl7mPNgKvr6ozq+pM4A0495EkLUhdQuHkqrr94EJVbQOcKVWSFqAuA817k/x74KZm+RXA/f2VJEkaSpeewr8ElgHva15LgVf2WZQkaRhdrj56FHCeI0laBLpcffTBJKePLT8tydZ+y5IkDaHL6aOlVfXlgwtNz+Hp/ZUkSRpKl1D4VpIzDi4kOZPmec2SpIWlSyj8CvDRJDcluQm4A3hjl40nWZtkd5I9Sa6Zpd1Lk1SSiW5lS5L60GWg+QNJnguc36y6uqq+MNf3mmm3rwcuAqaA7Uk2V9Wuae1OZTSQfdfhFi9N52NApaPTpacA8HxgsnmdP2vLbzsP2FNVe6tqP3ALsG6Gdr8K/AajJ7xJkgbU5eqjXwOuAnY1r6uSvKXDtpcDD44tTzXrxrd9LrCyqv68c8WSpN50uaP5EuCcqvoWQJIbgXuYe1whM6xrB6iTnAC8Hfi5uQpIsgHYAHDGGWfM0VqSdKS6nj46fez9aR2/MwWsHFtewei5DAedCjwH2JbkAUanpTbPNNhcVRuraqKqJpYtW9Zx95Kkw9Wlp/AW4J4ktzP66//HgDd1+N52YE2S1cDngMuAf37ww6p6jNGUGQAk2Qb866ra0bl6SdIx1eXqo5ubX9jPYxQKv1xVn+/wvQNJrgC2AkuAG6pqZ5LrgB1VtfnoSpckHWtzhkKS26rqAmDzDOtmVVVbgC3T1l17iLaTc1YrSerVIUMhyZOBpwBLkzyNbw8cPxX4vnmoTZI0z2brKbwG+CVGAfDxsfVfYXRTmiRpgTlkKFTVO4B3JLmyqn57HmuSJA2ky9VHjyX52ekrq+oPeqhHkjSgLqHwvLH3TwYuYHQ6yVCQpAWmyyWpV44vJzmNbz+vWZK0gHS9o3nc48BZx7oQSdLwutyn8D/59pxFS4BnA5v6LEqSNIwuYwpvHXt/gNH9Cuv7KUc6PvkcBy0UXcYUPpLkHEbzFr0MuB/4k74LkyTNv9nuaD6L0SR264EvAu8FUlUvnKfaJEnzbLaewl8Dfwm8pKr2ACS5el6qkiQNYrarj/4Z8Hng9iS/m+QCZn5wjiRpgThkKFTVn1bVy4EfALYBVwPPSPLOJC+ap/okSfNozvsUquprVfWeqnoxo6en3Qtc03tlkqR5d1g3r1XVl6rqXVX1430VJEkazpHc0SxJWqAMBUlSy1CQJLUMBUlSy1CQjnOTk5Pt3EvS0TIUJEktQ0GS1DIUJEktQ0GS1DIUJEktQ0GS1DIUJEktQ0GS1DIUJB0Rb5pbmAwFSVKr11BIsjbJ7iR7kjzhwTxJXp9kV5JPJLktyZl91iNJml1voZBkCXA9cDFwNrA+ydnTmt0DTFTVPwZuBX6jr3okSXPrs6dwHrCnqvZW1X7gFmDdeIOqur2qHm8W72T0uE9J0kD6DIXlwINjy1PNukN5FfD+mT5IsiHJjiQ79u3bdwxLlCSN6zMUMsO6mrFh8gpgAvjNmT6vqo1VNVFVE8uWLTuGJUqSxp3Y47angJVjyyuAh6Y3SnIh8CvAP6mqr/dYjyRpDn32FLYDa5KsTnIScBmwebxBknOBdwGXVtUjPdYiSeqgt1CoqgPAFcBW4D5gU1XtTHJdkkubZr8JnAL8cZJ7k2w+xOYkqeWNc/3p8/QRVbUF2DJt3bVj7y/sc/+SpMPjHc2SpJahIElqGQqS1NGQYxnztW9DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLUMhQkSS1DQZLU6jUUkqxNsjvJniTXzPD59yR5b/P5XUlW9VmPJGl2vYVCkiXA9cDFwNnA+iRnT2v2KuDRqvoHwNuBX++rHknS3PrsKZwH7KmqvVW1H7gFWDetzTrgxub9rcAFSdJjTZKkWaSq+tlw8lJgbVX9q2b5Z4Afrqorxtp8qmkz1Sx/umnzhUNt93vPfHZd9KYbDruee//qXgDO+cFzDvu7R2uofXvMC3+/Q+7bY55fR7vvTT///LuramKudn2Gwk8DPzEtFM6rqivH2uxs2oyHwnlV9cVp29oAbAA45Zl//4cuefNNvdQsSQtV11A4sccapoCVY8srgIcO0WYqyYnAacCXpm+oqjYCGwEmJibqva/5kV4KlqSFatPPd2vX55jCdmBNktVJTgIuAzZPa7MZuLx5/1Lgw9VX10WSNKfeegpVdSDJFcBWYAlwQ1XtTHIdsKOqNgO/D9yUZA+jHsJlfdUjSZpbn6ePqKotwJZp664de/93wE/3WYMkqTvvaJYktQwFSVLLUJAktQwFSVLLUJAktXq7o7kvSfYBnznCry8FDjmFxgLlMS8OHvPicDTHfGZVLZur0XEXCkcjyY4ut3kvJB7z4uAxLw7zccyePpIktQwFSVJrsYXCxqELGIDHvDh4zItD78e8qMYUJEmzW2w9BUnSLBZNKCRZm2R3kj1Jrhm6nr4lWZnk9iT3JdmZ5Kqha5oPSZYkuSfJnw9dy3xIcnqSW5P8dfNvveAfNpLk6ub/9KeS3JzkyUPXdKwluSHJI83TKQ+u+94kH0zyN83Pp/Wx70URCkmWANcDFwNnA+uTnD1sVb07ALyhqp4NnA+8dhEcM8BVwH1DFzGP3gF8oKp+APhBFvixJ1kOvA6YqKrnMJqWfyFOuf9uYO20ddcAt1XVGuC2ZvmYWxShAJwH7KmqvVW1H7gFWDdwTb2qqoer6uPN+68y+mWxfNiq+pVkBfBPgd8bupb5kOSpwI8xei4JVbW/qr48bFXz4kTg7zVPa3wKT3yi43Gvqu7giU+hXAfc2Ly/EfjJPva9WEJhOfDg2PIUC/wX5Lgkq4BzgbuGraR3vwX8W+BbQxcyT74f2Af89+aU2e8lOXnoovpUVZ8D3gp8FngYeKyq/mLYqubNM6rqYRj90Qc8vY+dLJZQyAzrFsVlV0lOAf4E+KWq+srQ9fQlyYuBR6rq7qFrmUcnAs8F3llV5wJfo6dTCt8tmvPo64DVwPcBJyd5xbBVLSyLJRSmgJVjyytYgF3O6ZI8iVEgvKeq3jd0PT17AXBpkgcYnR788SR/OGxJvZsCpqrqYA/wVkYhsZBdCNxfVfuq6hvA+4DnD1zTfPm/SZ4J0Px8pI+dLJZQ2A6sSbI6yUmMBqY2D1xTr5KE0bnm+6rqbUPX07eqemNVraiqVYz+fT9cVQv6L8iq+jzwYJJnNasuAHYNWNJ8+CxwfpKnNP/HL2CBD66P2Qxc3ry/HPgffeyk12c0f7eoqgNJrgC2Mrpa4Yaq2jlwWX17AfAzwCeT3Nuse1Pz3GwtHFcC72n+2NkLvHLgenpVVXcluRX4OKMr7O5hAd7ZnORmYBJYmmQKeDPwa8CmJK9iFI69PN/eO5olSa3FcvpIktSBoSBJahkKkqSWoSBJahkKkqSWoSB1lORvh65B6puhIElqGQrSUUjykiR3NRPSfSjJM5r1y5o57z+e5F1JPpNk6dD1SnMxFKSj81Hg/GZCulsYzdIKoztQP1xVzwX+FDhjoPqkw7IoprmQerQCeG8zQdlJwP3N+h8Ffgqgqj6Q5NGB6pMOiz0F6ej8NvBfq+ofAa8BDj4acqbp2qXveoaCdHROAz7XvL98bP1HgZcBJHkR0MvzdKVjzQnxpI6SfIvvfA7H24BPA29nFAx3As+rqskkTwduZhQGHwFeDqyuqq/Pb9XS4TEUpB4k+R7gm8207T/C6Olo5wxdlzQXB5qlfpzBaO77E4D9wKsHrkfqxJ6CJKnlQLMkqWUoSJJahoIkqWUoSJJahoIkqWUoSJJa/x/qcFhLHLYw4wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f03012464e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.acorr(chain[:,0] - np.mean(chain[:,0]))\n",
    "plt.xlim(-0.5, 10.5)\n",
    "plt.xlabel('Lag')\n",
    "plt.ylabel('Autocorrelation')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}