{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING (theano.configdefaults): install mkl with `conda install mkl-service`: No module named 'mkl'\n"
     ]
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from pymc3.gp.util import plot_gp_dist\n",
    "\n",
    "#https://docs.pymc.io/notebooks/GP-Marginal.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fcf4339ec18>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 18876 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(5)\n",
    "x = np.linspace(0, 10)[:, None]\n",
    "y = np.sin(x)\n",
    "\n",
    "X = np.random.choice(x.flatten(), 7, replace=False)[:, None]\n",
    "Y = np.sin(X)\n",
    "\n",
    "plt.plot(x, y)\n",
    "plt.plot(X, Y, 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/theano/tensor/basic.py:6611: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  result[diagonal_slice] = x\n",
      "Only 100 samples in chain.\n",
      "Sequential sampling (1 chains in 1 job)\n",
      "NUTS: [f*, ls]\n",
      "100%|██████████| 150/150 [03:43<00:00,  1.05it/s]\n",
      "The acceptance probability does not match the target. It is 0.9064175024104284, but should be close to 0.8. Try to increase the number of tuning steps.\n",
      "The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n",
      "Only one chain was sampled, this makes it impossible to run some convergence checks\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as model:\n",
    "    ls = pm.HalfCauchy('ls', beta=5)\n",
    "    \n",
    "    # zero mean\n",
    "    mean = pm.gp.mean.Constant(0)\n",
    "    # input_dim = the total columns of X\n",
    "    cov = pm.gp.cov.ExpQuad(input_dim=1, ls=ls)\n",
    "    \n",
    "    # The GP\n",
    "    gp = pm.gp.Marginal(mean , cov)\n",
    "        \n",
    "    # Observed data\n",
    "    f = gp.marginal_likelihood('y_obs', X, Y.flatten(), noise=0)\n",
    "    \n",
    "    # Conditional GP\n",
    "    f_star = gp.conditional('f*', x) \n",
    "    step = pm.NUTS()\n",
    "\n",
    "    trace = pm.sample(chains=1, step=step, tune=50, draws=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [00:01<00:00, 891.40it/s]\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    pred_samples = pm.sample_ppc(trace, vars=[f_star, f], samples=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fcf281c8748>]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16904 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(pred_samples['f*'].mean(0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fcf28183d68>]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 74164 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the results\n",
    "fig = plt.figure(figsize=(12,5))\n",
    "ax = fig.gca()\n",
    "\n",
    "plot_gp_dist(ax, pred_samples[\"f*\"], x)\n",
    "plt.plot(X, Y, 'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/theano/tensor/basic.py:6611: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  result[diagonal_slice] = x\n"
     ]
    }
   ],
   "source": [
    "with model:\n",
    "    prior_pc = pm.sample_prior_predictive(samples=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
