Revision 1db48f3a735eb0fba06a7d503f080a7ead512604 authored by Artem Artemev on 11 July 2018, 12:50:44 UTC, committed by GitHub on 11 July 2018, 12:50:44 UTC
1 parent 707b195
mcmc.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fully Bayesian inference for generalized GP models with HMC\n",
"--\n",
"\n",
"*James Hensman, 2015, 2016, 2017*\n",
"\n",
"It's possible to construct a very flexible models with Gaussian processes by combining them with different likelihoods (sometimes called 'families' in the GLM literature). This makes inference of the GP intractable since the likelihoods is not generally conjugate to the Gaussian process. The general form of the model is \n",
"$$\\theta \\sim p(\\theta)\\\\f \\sim \\mathcal {GP}(m(x; \\theta),\\, k(x, x'; \\theta))\\\\y_i \\sim p(y | g(f(x_i))\\,.$$\n",
"\n",
"\n",
"To perform inference in this model, we'll run MCMC using Hamiltonian Monte Carlo (HMC) over the function-values and the parameters $\\theta$ jointly. Key to an effective scheme is rotation of the field using the Cholesky decomposition. We write\n",
"\n",
"$$\\theta \\sim p(\\theta)\\\\v \\sim \\mathcal {N}(0,\\, I)\\\\LL^\\top = K\\\\f = m + Lv\\\\y_i \\sim p(y | g(f(x_i))\\,.$$\n",
"\n",
"Joint HMC over v and the function values is not widely adopted in the literature becate of the difficulty in differentiating $LL^\\top=K$. We've made this derivative available in tensorflow, and so application of HMC is relatively straightforward. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exponential Regression example\n",
"The first illustration in this notebook is 'Exponential Regression'. The model is \n",
"$$\\theta \\sim p(\\theta)\\\\f \\sim \\mathcal {GP}(0, k(x, x'; \\theta))\\\\f_i = f(x_i)\\\\y_i \\sim \\mathcal {Exp} (e^{f_i})$$\n",
"\n",
"We'll use MCMC to deal with both the kernel parameters $\\theta$ and the latent function values $f$. first, generate a data set."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:08:57.256029Z",
"start_time": "2018-06-13T10:08:56.212731Z"
}
},
"outputs": [],
"source": [
"import gpflow\n",
"from gpflow.test_util import notebook_niter\n",
"import numpy as np\n",
"import matplotlib\n",
"%matplotlib inline\n",
"matplotlib.rcParams['figure.figsize'] = (12, 6)\n",
"plt = matplotlib.pyplot\n",
"\n",
"X = np.linspace(-3,3,20)\n",
"Y = np.random.exponential(np.sin(X)**2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"GPflow's model for fully-Bayesian MCMC is called GPMC. It's constructed like any other model, but contains a parameter `V` which represents the centered values of the function. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:08:57.312496Z",
"start_time": "2018-06-13T10:08:57.257138Z"
},
"scrolled": false
},
"outputs": [],
"source": [
"with gpflow.defer_build():\n",
" k = gpflow.kernels.Matern32(1, ARD=False) + gpflow.kernels.Bias(1)\n",
" l = gpflow.likelihoods.Exponential()\n",
" m = gpflow.models.GPMC(X[:,None], Y[:,None], k, l)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `V` parameter already has a prior applied. We'll add priors to the parameters also (these are rather arbitrary, for illustration). "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:08:57.322680Z",
"start_time": "2018-06-13T10:08:57.314029Z"
}
},
"outputs": [],
"source": [
"m.kern.kernels[0].lengthscales.prior = gpflow.priors.Gamma(1., 1.)\n",
"m.kern.kernels[0].variance.prior = gpflow.priors.Gamma(1., 1.)\n",
"m.kern.kernels[1].variance.prior = gpflow.priors.Gamma(1., 1.)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Running HMC is pretty similar to optimizing a model. GPflow only has HMC sampling for the moment, and it's a relatively vanilla implementation (no NUTS, for example). There are two things to tune, the step size (epsilon) and the number of steps [Lmin, Lmax]. Each proposal will take a random number of steps between Lmin and Lmax, each of length epsilon. \n",
"\n",
"We'll use the `verbose` setting so that we can see the acceptance rate. <- this is broken :("
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:08:57.883746Z",
"start_time": "2018-06-13T10:08:57.324047Z"
}
},
"outputs": [],
"source": [
"m.compile()\n",
"o = gpflow.train.AdamOptimizer(0.01)\n",
"o.minimize(m, maxiter=notebook_niter(15)) # start near MAP"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:09:12.888491Z",
"start_time": "2018-06-13T10:08:57.885123Z"
},
"scrolled": false
},
"outputs": [],
"source": [
"s = gpflow.train.HMC()\n",
"samples = s.sample(m, notebook_niter(500), epsilon=0.12, lmax=20, lmin=5, thin=5, logprobs=False)#, verbose=True)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:09:15.681006Z",
"start_time": "2018-06-13T10:09:12.889934Z"
}
},
"outputs": [],
"source": [
"xtest = np.linspace(-4,4,100)[:,None]\n",
"f_samples = []\n",
"for i, s in samples.iterrows():\n",
" f = m.predict_f_samples(xtest, 5, initialize=False, feed_dict=m.sample_feed_dict(s))\n",
" f_samples.append(f)\n",
"f_samples = np.vstack(f_samples)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:09:15.795838Z",
"start_time": "2018-06-13T10:09:15.682398Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(-0.1, 3.2763890749508007)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XmcXFWZ+P/Pc29tvaX3zh5CCCGEPYQl7iJREBBEFNRRmRllnMH56oy7Ds7o/HAZ/TrOGEdHkEEdfgMyAgKiGCQqypawBLLvSafTa1V1177ce8/3j6rudDq9VHfXnvN+vfKil9vVp+mqp899znOeI0opNE3TtOpilHoAmqZpWv7p4K5pmlaFdHDXNE2rQjq4a5qmVSEd3DVN06qQDu6apmlVSAd3TdO0KqSDu6ZpWhXSwV3TNK0KuUr1jdva2tTSpUtL9e01raxZjsJlSKmHoZWhF154YUAp1T7VdSUL7kuXLmXz5s2l+vaaVrZCiTQ7u8NctLQZER3gteOJyKFcrtNpGU0rM5G4RSieJmk5pR6KVsF0cNe0MhOIpUjbjg7u2qzo4K5pZUQpRTCWwus2SKbtUg9Hq2A6uGtaGYmnbRxH4TVNQol0qYejVTAd3DWtjESTNgrwuAwiSavUw9EqmA7umlZGgrEUXtPEbRqEEzq4azOng7umlRF/JInXbWAagu0okpbOu2szo4O7ppWJRNomaTm4zWMvS10xo82UDu6aViZiqRNn6cm0Du7azOjgrmllYjCWwmUce0m6DIOwrpjRZkgHd00rE4FoCp/72EvSoxdVtVnQwV3TykDadggnLbwuc+RjHpcO7trMTRncRcQnIs+LyBYR2SYiXx7nGq+I3Ccie0XkORFZWojBalq1SloOY1uEmYZgOQ4pvaiqzUAuM/ckcJlS6jzgfOAKEbl0zDV/CQSVUsuBfwW+kd9halp1S48T3IfpckhtJqYM7iojkn3Xnf2nxlx2LfDj7Nv/C7xFdK9STctZ2nFOeFEN0+WQ2kzklHMXEVNEXgb6gA1KqefGXLIQ6ARQSlnAENCaz4FqWjVLpGzMceZDLsMgovPu2gzkFNyVUrZS6nxgEXCxiJw9k28mIreIyGYR2dzf3z+Th9C0qpSwHFzmiS9HtymEdY8ZbQamVS2jlBoENgJXjPlUF7AYQERcQCPgH+frf6iUWqOUWtPePuUpUZp20oilrHGP1fOYBhFd667NQC7VMu0i0pR9uwZYB+wcc9nDwIeyb98APKmUmiiFqGnaGIm0M25wd5kGScvBdvTLSZueXM5QnQ/8WERMMn8MfqaUelREvgJsVko9DPwI+KmI7AUCwE0FG7GmVRmlMg3Cat3muJ8XhLTtYBrjf17TxjNlcFdKvQJcMM7HvzTq7QTw7vwOTdNODmlb4ThMeBi2QpGyHXwTBH9NG4/eoappJZa2HcbJyBx/jS6H1KZJB3dNKzHLnjyfbiAk9Hmq2jTp4K5pJZayJ97ABOAyZdx2wJo2GR3cNa3EkpaNMcmGbpdh6OCuTZsO7ppWYvGUPW4Z5DC3KcR1WkabJh3cNa3E4mn7uEM6xnKZBom0jaNr3bVp0MFd00oskbJxmZOXywiZ5mKalisd3DWtxBKWjTlFLaQC3dddmxYd3DWthNK2g6OYdEH12LU6LaPlTgd3TSshy1YTHtIxmgApfWiHNg06uGtaCaXs3FItuhxSmy4d3DWthKwpNjANc5lCVAd3bRp0cNe0EkpZDpJDYsZtGsT1oR3aNOjgrmklFEuPf0jHWC5DSFg2+pgELVc6uGtaCSVSzpQ17pBpB+woXTGj5U4Hd00roal2p46VznEBVtN0cNe0Ekqkp97ANCxTDqmDu5YbHdw1rURsR2E5KufgrtAzdy13OrhrWomk7dwqZYbpWndtOnRw17QSSdsOKqcq9wyXIcTSuhxSy40O7ppWIlMdrzdWptZdp2W03OjgrmklMt0Wvi5Tz9y13Ongrmklkkw7mDl0gxxmiGA7Si+qajnRwV3TSmQ6ZZDDdDmklqspg7uILBaRjSKyXUS2icjHx7nmTSIyJCIvZ/99qTDD1bTqkbCcaQd30OWQWm5cOVxjAZ9USr0oIg3ACyKyQSm1fcx1Tymlrs7/EDWtOiWnuTsV9IlMWu6mfGYppbqVUi9m3w4DO4CFhR6YplW7maRlXIZBJA/dIVOWQ1If/lHVpjVtEJGlwAXAc+N8eq2IbBGRX4nIWRN8/S0isllENvf39097sJpWLZxp7k4d5s5TX/f9AxGe3uvn0EBU3wlUqZyDu4jUAz8HPqGUCo359IvAKUqp84DvAg+N9xhKqR8qpdYopda0t7fPdMyaVvHSjjONvanHuE2D2Cxn7vGUTc9QgqYaNwf8UZ7b7ycQSc7qMbXyk1NwFxE3mcB+j1LqgbGfV0qFlFKR7NuPAW4RacvrSDWtilj2dPamHpOPvu5dg3FMQ3CZBq11Xjwug0OB2IwfTytPuVTLCPAjYIdS6tsTXDMvex0icnH2cf35HKimVZNcD8YeS0RAQXKGqZSkZdMVjNHgdY98rMZtMhhP6xx8lcmlWua1wAeAV0Xk5ezHvgAsAVBK/QC4AfhrEbGAOHCT0kfGaNqE0k5uZ6eOZ7g7pM9tTvtre4cSKDgu1z/8B2MolqZjzvQfUytPUwZ3pdQfYfJJhlJqPbA+X4PStGqXzvHs1InMZBHUsh0OB2I0+twnfK7WY9ITStAxxzfjMWnlRe9Q1bQSiM+gDHKYgZBITz+F0h9OYjkKl3niy77GbRKIpnTlTBXRwV3TSiBh2TkdjD0et2nMqBzySDBGvXf8m/XskhmhRHpGY9LKjw7umlYCifTMWg9ApjtkdJrlkJGkRSRp43VNnFP3ukz6QokZjUkrPzq4a1oJpNKzm7lP90Sm/nAClzn596v1mPRHkli6d01V0MFd04pMKUXSnvnM3TSEtO3kHIQdR9EVTFDvmbx+whBBKQgldM/4aqCDu6YVmeUolDqW554JAdI5nuQ0FE9jOc64C6ljuQ2DQFTvVq0GOrhrWpHNdAPTWLlWtvSEEnjN3OrX3aZBPK3TMtVAB3dNK7LZbGAapoCkPXXePWU59IYS1HpzC+4uU2bdu0YrDzq4a1qRTfdg7PG4DINYcurgHowmUSqTT8/tcWffu0YrDzq4a1qR5aMaxW0KsVRmhp20bLZ2DdIZiI08tuMoOgMxtneHafDl0mUkQ7KLqrnm87XylftvXdO0vEik7WkdjD0el5Eph4ynbF45MkjKdhiIpDgciHFqax3doTihuEVzrWfGR/l5XHruV8l0cNe0IktazrSP1xvLbQqD8TQvHg5iIDTVeIBMUN7dF8brMmmr987osQV9Tms10MFd04psNn1lhokIaVvhdQm1o+rX3dke7bOhz2mtDjq4a1qRJdPOlLtFc9HRUJgOjsMpH62y6aSaphVZ0pr9zL2QXIYQS+tyyEqng7umFZGdPRg719LEUnCZQjyp0zKVTgd3TSuitD2zg7GLyWUYxPXMveLp4K5pRWQ5MzsYu5hMQ7AcpbtDVjgd3DWtiILRVKmHkJPpNCbTypMO7ppWAEopHEcd9/7+/gh7esMjNenlLqVn7hVNl0JqWgH0h5Ps6g3T0eClrd5LXzhJ91Cc1npvWS+mDlPojUyVTgd3TSuAWMpGOeCPpOgeSmCI0FbnnVUP92IyEJIzOIRbKx86uGtaASQtB7fLmPBA6nLnMkVvZKpwU+bcRWSxiGwUke0isk1EPj7ONSIi/y4ie0XkFRFZXZjhalplSMzijNRyoHepVr5cphUW8Eml1Isi0gC8ICIblFLbR11zJXB69t8lwPez/9W0k1Iybed0rF25cptCXKdlKtqUzz6lVLdS6sXs22FgB7BwzGXXAj9RGc8CTSIyP++j1bQKkSjzFgNTMQ0hqQ/tqGjTmlqIyFLgAuC5MZ9aCHSOev8IJ/4B0LSTgu0onGmcflSOhg/t0OWQlSvn4C4i9cDPgU8opUIz+WYicouIbBaRzf39/TN5CE0re5kWA5Ub2EfTG5kqV07BXUTcZAL7PUqpB8a5pAtYPOr9RdmPHUcp9UOl1Bql1Jr29vaZjFfTyl6mxUB1BMW07utesXKplhHgR8AOpdS3J7jsYeCD2aqZS4EhpVR3HsepaRXDqoDmYLnQJzJVtlyqZV4LfAB4VURezn7sC8ASAKXUD4DHgLcDe4EY8Of5H6qmVYa0XR3zdtMwiCZ1xUylmjK4K6X+CJNPRFRmSf3WfA1K0ypZ2nIwqmDu7jENQokUUFfqoWgzULmFuJpWpuKWVdFlkMN8boNgLK3PU61QOrhrWp4lrfyckVpqw31wwol0iUeizYQO7pqWZ4m0U9E17qN5TZO+cLLUw9BmQAd3TcuzZIX3lRmtxmMyEEliO9WwRHxy0cFd0/JIKUXKdqoi5w6ZNgS2o4gk9JmqlUYHd03LI8tRKEXF9G3PhcswGIjo1Eyl0cFd0/LIslUVFEEer9Zj0hdO6CZiFUYHd03LI8txqmID02hu0yBpOUR1f/eKooO7puWRVaWNtgwRBmOpUg9DmwYd3DUtj9JOdW74qfWY9AwlSj0MbRp0cNe0PEpWUY37aF6XSSRhEdepmYqhg7um5VGln506KUGnZiqIDu6alkdJq3pq3Meq87g4qlMzFUMHd03Lo8zMvTpfVj63SSieJqEPzq4I1fksrEJp22Fb15A+PKHMVfPMHUAEhnRqpiLo4F4hYkmbzmCMvX2RUg9Fm4DjKNJV1HpgPDVuk56QTs1UAh3cK0QkmabW46J7ME6ffnGVJeskaK5V4zYJxNIkLZ2aKXc6uFeIYCyN12XQVOthR3dIl6SVIcupjrNTJyMioCAU143Eyp0O7hVAKcVgLIXXZeI2DTymyc6ekO71UWaq5ezUqdR6TI4OxfXzr8xVZHBP2w6BaIrD/iivdg1VfZoiaTnYDiO53Hqfi6FYmpievZcVa4aL3ffefQf+gf6R9/0D/dx79x35Glbe1bhN/OEkfSHdKbKcTXlAdjnqGUqwsyeELzuT3RFJ0uBzU+MxSz20goilbMbOCUUgGEtR563IX2FVspzpd4S89+47+Optn+a+n/6IO+97BIAP33gN+3bvBOCmmz+S51HOnohk0oM9Iep9Lv0cLFMVOXNXSlHjdtFU66HO68JlGOzuDVftbWIonj6hdrrOm1lc1cpHOHHi72kq666+jtNWrGTf7p1cf/larr98Lft27+S0FStZd/V1BRrp7LlNA6/LZEd3SJ/SVKYqMriPNafGjT+apHuwOtMzgVgKr+v4X5XXZRJN2nphtUwopegNJamd5t1ja1s7d973CM2tbQT9AwT9AzS3tnHnfY/Q2tZeoNHmR73XRSRpcWAgWuqhaOOYMriLyF0i0iciWyf4/JtEZEhEXs7++1L+hzm15hoPu/vCxFLVtYpvO4pw3DohuAO610cZiaZs0raDy6yK+VLOmms9dAZiRJPV9bqrBrk8E+8GrpjimqeUUudn/31l9sOaPpdp4DYMDg7Exv28U6G3jvF0Jt8+3rFttW6X3lBSJoZiqRmVQfoH+vnwjdeMzNiHZ/AfvvGa4xZZy5Uhgts0OBwY/3Wnlc6UwV0p9QcgUISxzFqDz0VvKHFC74uBcGYBthJNdidS4zEZjOleH+WgJ5Sk1jP9hcUNjz40kmN/4IlneOCJZ0Zy8BsefagAI82/OT4XPUNxInr2Xlbytcy9VkS2AEeBTymltuXpcadFRDAE+kIJlrTWAZkZ+77+KPG0zYoKvG0ejKXxmhPncUUyC64+d3VWClWCpGUTTqRprfNO+2uHq2HWXX3dSI79zvseYcOjD5Vlpcx4RASvy+SwP8qqBY2lHo6WlY9I9yJwilLqPOC7wITTDRG5RUQ2i8jm/v7C3HLO8bk5HIiNrOD7o8mRevBKnFkEoim87ol/TbrXR+mFE7N7Xt1080eOWzxtbWuvmMA+rN7rojeUJJxIl3ooWtasg7tSKqSUimTffgxwi0jbBNf+UCm1Rim1pr29MJUALtMg7TgEokkcR7G/P0p9tlxyIFJZmy7CiTTxlI17kruNGreJP5KquoXkStIfTuJ1ndx3TiKCz2VyyK8rZ8rFrIO7iMyT7GqfiFycfUz/bB93Nurcbg4FYgxEMrN2n9ukzmPSH05WRC28UoquYJzNB4NTltZlbokN9vRGKuJnqzaOoxgIJ6nRaTHqfS76wym9BlQmpsy5i8j/AG8C2kTkCPCPgBtAKfUD4Abgr0XEAuLATarEUabGYzIQSbAvHaHBl/kRXaZBMp7Zsl/OO+rStsOO7hD+SIrmWk9O7WMbfG4Gokl6hxLMa6opwii1YeGkha1UVbf5nQ6RTPpTrwGV3pRRTin13ik+vx5Yn7cR5YnP5SJp29T73Mc+KDAUS5dtcFdKsac3TDCaoq1+eotzjT43e/oiNNV59AuriPyRJGYVHog9Ux7TwB9JTvv5q+VfZZWOTEO9z3VC9UKt20VvuHwXH3uGEvSGkrTMoOrCbRqIwP5+fZhHscRTNocDMRpGTyBOcrWeTGpGpwhLr2qD+3h8boOheJqUVX5H1UWSFjt7wjTVzDxQzPG56QklGCjjP2DV5KA/isswdEpmFNMQLMchqttilNxJFdyHd3mWW0mklT0ftdZjzqoOX0SY43Ozqzeiz1otsKF4mu6hOHN85ZniKyUBQjFdEgmZ9iFJyyaWsor+mjzpnpke06B7MI7bFOo8LowymHUdDsRIpO0ZpWPGyjQUszg4EOX0uQ15GJ02llKKvX1h6jyucdtCnOxq3C76IwkWNJ/ci/vdg3F290Yy7boVzKlxsfqUlqJ9/5MuuNd5XfijKfrDSVymMK/Rx6lt9SW7tY4kLQ75Y7TUefL2mM21HjqDcToafDTW6nxwvvWHk4QSFm15+GNcjXxug2A8hVWBO8LzZSieZmdP+LiKt4FIgnAiXbQ1mpPu/7whQnOth9Z6L/VeN0eCcbZ2DZUkD5+ZAUbwuU2MPM4ARYR6j4udPaGyXF+odN1DCWp1RdKERASlyi/9WSyJtM2rXYPUe13HTRrdpsnRIp7BcNIF99FMQ2it8xJKpHmlc7Domy/6w0mC0RT1BSjNrPGYpG2HzYcCui1wnkWS1km/I3UqLsMgED35nne2o9jRHcJATihJbvC56BlKkLSKE2dO6uA+rKnGQ8p22NI5WLRFj5TlsLsvTOMsqmOm0ljjwW0YvHh4kP39EX1iTh5YtkPKdnSFzBRqPSZ94cpq95EPXcEYg7HxUy/Dd+f9Rfr/ooN7VoPPTcKyOewvTl/q3lACy1aT9o3JB5/bpLXOw+FAjM0HAwzFdRXDbCQtR79ocuA2DRJp+6Q6xCNp2Rzyx2iunXj9rMHnpjMQK8r5Evp5OkpTjYfDgWjBO9vZjuKQP0pjsRZWJJN+EoQXDgXZ2xfB0qWSM5K0HPT9T25MEfwV1qxvNo4OxnGmaEXhNg0SlsNgESZZOriPYohQ43axqzdc0L+s/kgSy1FFrySo8WRm8Z2BKAf1uZczEk9ZGDM6c+nkU+d1cWQwXrGnoE1HIp2ZtTfWTF315nOZ9AwVfmFVB/cx6rwuwok0PUOF2eWplOKwP0bdDE7tyYfhWfyhQEwvtM5ANGXjHu88W+0EbtMgZdmz7ndfCY4EY5iG5LQW4zKEtK3TMiXR5Msctj1UgF12oYRFuMRd84Z3su7oDh23gKyU0j1BphCOW7hNPXPPldsw6QkVr/yvmJRSpG2HoXiazkCcOWXWY+ik28SUC5dpUO918eLhIOcuaqQ1jx3uuoJxvGUw8/O5TfxRi0MDUZa21dEfTnLQH2VxSy2LmmtLPbyypJQimkrndOutZdT7Mic0LWt3Cl48UCy2o9jXH6ZnKIGjMu0WavK8VyUfquP/dgF4XSZzfG5eOTJEf54acSXSNr2hREHq2meipdbDoUCMZ/f72dkTBoTuIm6yqDRJy0Epyu5FPJV7774D/8CxYy39A/3ce/cdRfnehgiOUlVTpZWyHLZ2DXF0MEFTjYfWOi8tdd6ybCNefiMqIx6XQWNNJsAva7dY0lI3q/rmg/4oLlPKph+JiNCSLdsaXtz1R5LEUzY1U5wAdTJK2ZVXKXPv3Xfw1ds+zX0//RF33vcIAB++8Rr27d4JUJSzWmvcJl3BWMX3eI+nMjtPk5Yzo8PQi00H9ym4TYO2ei8HB2L4IynOnD9nRn+lA9EUXcE47WX2BB9bsSMCQ/EUNZ6Tu+nTeJLpyisfXXf1ddz30x+xb/dOrr98LQBB/wCnrVjJuquvK8oYaj0uBiJJEmm7Yg+ScRzFtq4hbDtTMl0JdFomB4YIbfVe0rZi08EAR4LT24SQsjJH5zXWuMtm1j6RGreLnpDuBz+eaNLCZVTWS6a1rZ0773uE5tY2gv4Bgv4BmlvbuPO+R2htK8wh9eMRoaI3NB0dihNJWdRXUIvnynqmlli910VTjYc9fWFe7gzm/GTdPxDBcVRF9CPxuQ2CsfI80KTUwkldKTNTbsMgWKGlt4m0zb6+SMXM2Ifp4D5NpiG01flIWorNBwNTBvhgNh1TyB4y+SQimcMWCrxLtxJFExaeCqv48A/08+EbrxmZsQ/P4D984zXHLbIWWubQ+soM7vv6IpgVeOJWZT1Ty0i9N3PQx8AUTYD2DURo8JZ/OmY0n8ukT6dmjmM7ioRlV1x/8g2PPsS+3Ts5bcVKHnjiGR544hlOW7GSfbt3suHRh4o2juFeM8XuvDoTadshlrKIJC16huL0hRMVMzkbrXISSGWozuPi6FCcJa214wbvUCJNOJ6mrd5XgtHN3PAs62Q+bGGspGVXZNOB4WqYdVdfN5Jjv/O+R9jw6ENFqZQZK5Yq70XVlOXw4qEgScseqYya46usdMwwHdxnwW0amQCetMbdnXZ0MI7HLN8n8kSGa5PDCYvmPJ4QVckqsVJm2Ngg3trWXpLAnunxnszrqWP5NHx8Ytpx8nLkZalNOS0TkbtEpE9Etk7weRGRfxeRvSLyioiszv8wy5cp46dmEmmb7qFERa2uj+Z1mXQGi9P+uBIU64CFalbrMfGXcd69dyhBbyhZcQunE8nlnvtu4IpJPn8lcHr23y3A92c/rMpR73VxdDBxQmlkXyiBQeXtZhxW783UJgdPwtN0xhNOWFWzfX4qhdrR6jYN4mm7LP9QRpMWu/oiNFVgbn0iUz5blVJ/AAKTXHIt8BOV8SzQJCLz8zXAcucyDSzn+M53tqM4HIiVXSOh6WrwutndG9YnOJE58NhTBj2BCm14R+twNc1wtc1Xb/t03loWRJPlFdwtO7MPxWMaVbXGlI+fZCHQOer9I9mPnTRchkl/JFNdopSiL5QoSb/2fPO5TeJpm94CtT+uFPGUTSRxcpybuu7q60aqaa6/fC3XX752pNomHztaXYZRVneDSil294aJJK2y6fmUL0WNPiJyi4hsFpHN/f3Fq7EttHpv5uDbvlCCFw4G2dUbqvhZ+7BGn5t9/ZGyvJUulsFYiooslZmBQu9orfWYDJTR6UyH/TF6Q8mK6BUzXfkI7l3A4lHvL8p+7ARKqR8qpdYopda0txdv63OhmYZgOYptR0MooLXOVzX5WZdpoIDOIp0tW46ODiVKdrhKtXGbBrFU6fPujpO5w943EJn0zNNKlo8I9DDwwWzVzKXAkFKqOw+PW1Fa67y01XvLuoZ3phpr3BwOxE7KXauJtE04ka7K3+t4irGjVYSSnM7kOIo9vWGe3+/nqb39bD06RFONp+J2nuYql1LI/wGeAc4QkSMi8pci8lER+Wj2kseA/cBe4A7gbwo2Wq0kDBHqvC529xT2bNlyNFSkfijBWKosdm8WY0er12XSP8XO7kIYiCTpDMZxmQbNNR7a66vnDns8U95rKqXeO8XnFXBr3kaklaVajwt/NEnPUIIFzSdPO+CjQwlq3YVNyRwdjPN/7n0Jr8vgQ69ZyuVnzi1ZCW0xdrQO591tRxVt1qyU4sBAlAavq6oD+mgnx0+p5UWjz83e/nBZzDCLIWnZDMXT+NyFfZls3NVH0nIIJSy+++RePnX/Fvb1Rwr6PSdz080fOW7xNN87Wg0RHEcRKWJqJhBNEavgfvIzoYO7ljOXaWAaBju6QxXdmztXwwekF7Lpm1KKP+zO5LKvO38hLXUe9vRF+MKDr9IVrN4jD03DKFrVjFKKA/1R6gp8B1ZudHDXpmWOz00i5bDpYID9/RHSduX2XJlKTyhBTYFnevv6oxwdStBU6+bm1yzlB++/kEtObSGWsrn9se3EUtX5R7TWY9IXTpDJ6hbWYCxNJJU+6Y6O1MFdm7Z6n4vmWg+HAzFePjxYlQE+kbbxR1IFD+6/z87aX7e8DdMQajwmf79uBYtbaukMxvnOE3tmHQAzDbEiPLGjF3+Z1Ji7TYOk5RBNFTbFp5TigD9KzUk2awfdFVKbIUOE1jovgWiSPb1hzpw/p6J61k9lIJJEpLApGUcp/rg3E9zfcPqxHHetx8U/vP1M/v5nL/PMfj/3v3CE96xZPNHDTCgQTfHLV7t5ak8/3dldxobAJae2csXZ87hgcVNJf2dCZoNYoXaGxlIWu3rChOLpqtykNBU9cy9zhWrilC/NtR56Q0mODlZPiwKlFJ1F6A20/WiIgUiKjgYvK+c1HPe5BU01fPKtZyDAfz97iC2dg9N67IMDUf7+Zy/zs82ddGfTPsPB/Jn9fv7x4W18e8Pukt511XoyO7vzbfj39/yBAImUc1IGdtAz97I23MTpvp/+iDvvewSAD994Dft27wRO7NNdCiJCc62H3b1h6n2uijyxZqxQ3CJp2dR7C/uz/GFP5o/2609vH3cGfdHSFm68aDH3burkW7/ZxXduPJ/W+qkD1atdQ9z+y+1EUzZnzmvg/ZeewtkLGjENIRBNsWF7D//74hF+t7ufgUiSL759VUlaU/vcmZLIeMrOWz7ccRR7+sJ0DcZpqfVW7QalXOiZexkrdBOnfDENod7r4uXDQY4G4xW/0enoUOEPWbFshz/uHQDgjSvaJrzupouWcP7iJgbjaf7l8V1Yk8y0bUfxq63dfOkXW4mmbNYua+X/u+4czlsrTNUbAAAgAElEQVTUNBLkWuo83HjREr5+/bm01HrYejTEp3++hb7wzGfQs7m79LoMXuoMMhSf/e5nx1Hs6g3TPZSgre7kDuygg3tZK3QTp3zyuU3m+Nzs7guz6WAg02yrAqUsh95QgroCdwh87kCAcMJicXMNS1vrJrzONIRPrltBS52H7d0h/uvpgycssCqleGbfAH/7Py/yH7/bh+UorjpnPp+9YuWEbYpPa6/nW+8+j1NaajkSjHP7L3fMqN/LbFsEN/jcuA2DFw4F6My2uOgLJ9jbF55WuwvbUezsCdEzlKC1zltV6z8zpYO7ljcu06C1zoshwsudgxXZi2a4mqSQO0Qdpbhvc6ZL9lXnzJ8yEDXVevjM287AEHh4y1Fuf2zHyEx3R3eIT//vK3z1VzvpDMaZN8fHp956Bn/1hmVTzlzbG7x8/fpzmd/oY/9AlDv+sH/aP0s+7i59bpOWWi/7+iO8eCjI9qMhjg4m2No1lNOagOModvWE6AsnacshbXWy0MG9jBWjiVMh+NwmtW4X24+GKqpM0skeslLovt7PHwhwYCBKS52Hdavm5fQ1Zy1o5HNXrKTOa/LcgQD/596XuP2x7Xzm56+wqzdMU42bj75hGf/x/tW8ccX4Ofzx1PtcfO6KlbhN4fHtvfx2R++0fpZ83V2aRqb6avhfc60Hy1ZT7tRVSrG7L1y1bXtnQwf3MlaMJk6FUuMxSabtkm6jn66eoQSxdGEP5VBKce+mwwC8a/WiaZ3utPa0Nv79xgs4c/4cAtEUz+4P4HEZ3HjRYv7zAxdy1bkLZtQ3ZVl7PR9942kA/Mfv97G/TH5nTTVuuoLxCWvzlVLs64vQPRSntUwP3S4lXS1TxorRxKmQmms9HA3Gaan10DHHV+rhTCqRttnbH6bJV9ggselgkH39UZpr3bztrLnT/vqOOT6+9s5z+MXLXQxEkly/elFeUhHrzpzL9qMhfruzjy88+CqfvWIlFyxpnvLrxt5dAiN3l7NdGxIRmmo87OgJcdHSlhP+6B72xzgcjNOmc+zj0jP3MlfoJk6FJCI01WYWAsfm3x1H0RWMlU0Tsv39EUyjsGdojp61X7960YzvEExDuH71Im55w2l5yzGLCH/9ptNYu6yVaMrmnx7Zxi9fOTrl1xX67tLjMhCErV1Dxy349ocT7B2I0FLr0YF9AnrmrhWU2zSo87jYcniQC05ppt7rws4ugHUNxmmscXPe4qaSnk8aiKboCWXK5wrFUYr//7nD7OmL0FTj5oqzcsu1F5PXZfK5K1fy388e4v4XjvCDP+zn1aMhPnjpKSxoGr/NczHuLuf43IQSaV4+PMg5ixqxHMXWrhDNVXzQRj5IMRr3jGfNmjVq8+bNM/raw/4ohwPxqtgwc7KIpSzSjsO5i5o45I/ij6RorfMSSqTxmgbnLG4saoAPJdLEkhahhEX3UJw6j6tg3z+StPj2hl1sOhhEgE9cfjqXrZx+SqaYNu7qY/2Te0nZDoZk0jbvWbO4pOm1SNLCdhxEBLdhVGwjsETaxuMyOG9x04y+XkReUEqtmeo6PXPXiqLW4yKatNh8MIgpjFQ2DM/KXj0yxLmLmqa1wDgTSikO+aPsH4hiiOA1TRq87oId4NAVjPOVR7dxdChBvdfFp996BqtPmTqXXWpvPqODsxbM4d5Nnfx2Ry+Pb+9lw45eLjm1lavOnc+5CxuLng6p97pIpG0cpSo2sBeTnrmXOdtRPLPfz+PbekjbDsva6ljWXs+5ixrpaCjvRcrxTHT6zmAsxZwaN+csbMQo0K2242RK6zoDcVrrPQU/7WhXT5gvP7qNcMLi1LY6vnDlmcxrrLzf2ZFgjHs3dfLHvQPY2d3Hyzvq+avXL2Pl/DklHl3lKdbMvaKDu89t0B9OYohk2qW6TeZUQcBXKlNv/eyBAL/e2jPuoQYuQ3j7OfO5cc3iqviZIdOJcWGzj9M7GvI+K4ynbA4MRLL10IVfhNt8KMDXf7WTpOWw5pRmPnvFyoo/BSgQTfH4th5+tbWbYPYgkzef0c7NrzmVFl2KmDMd3Cdwz3OHePjlo+zvjzIQSTJ69AJce/5C/uK1SytmBV0pxe7eCF2DMfojKXqHEmw5MkjfqAOEFzbVcM15C5g3x8f+/gi7esM8fyCAAuq8Ju+5cDFXnTu/pIuS+aCUYiCaZMXcBhY118768RxH0R9OcGQwQTiexjQy1TuF9nLnIP/48FYcBW9Z2cHH3ry8oFU4xZZI29z/whEefOkIaVvRXOvmG+86l/mNJ8/ZurOhc+4T2NMb4bkDASBTEtZe70UELEcRiKZ46OUuTEP40NpTyj7AH/JH+eFT+3nlyNAJn2uqcbP6lGZef3obq5c0j6QQLszmaw8MRLjrTwd5uXOQ/3r6IL94+SjvuWgxb101t2IPABYRWmq97OoJ4zYM5o5JYSilSNsKJzshcZvGpNUSnYEYe/sjNHjdOXVTzJd7Nx3GUXDNufP5yOuXlf3zcLp8bpMPXHoK686cy3d+u5ttR0P8w0Nb+Zd3nVvU/8/a5Cpu5r61a4gtnYO4TYPTO+qPmxE9d8DP1361E9tRvO/iJbz34iX5HHLeJNI2P37mII+92o2jMgtFq5c00d7gpa3ey4q5DSzvqJ8yJ6yU4sXDg/z02YPs648C0FLr4aJTW7hoaTPnLWqqyFRA2nYYjKdZNa+BedkSvFjKYmd3iFDCIvN/RXC7hOXt9bTVe0/I04cSaV44GKS5trjlcgcHovztvS9R4za5+88votZTcfOnaYmlLG77xVZ290ZY3FzD164/96RYC5sNnZaZxGQLqn/cO8A3H9+Jo2DV/DmsmNvAirn1rF7SXPBOf7mwbIfbH9vB5kNBDIErz57P+y5eMqu8uVKKp/f5uef5w3QGYiMfd5vCeYuauOTUVi4+taWi8qKW7RCIpThz3hxMQ9jZE8Jjmsf1HU9ZDqFkilq3i+Ud9bRkc+mW7fDi4SBKUfTg+r2Ne/n1th6uOmf+yJb+ahdOpPn8A69yKBBjYVMN779kCa85rU3XoE+grIK7iFwB/BtgAncqpb4+5vM3A98EurIfWq+UunOyxyxktcyTO/v47pN7sEb1Fa/1mFxz7gLecd6Cki1AKqX4zhN7eHJXHw0+F/987dmc1l6f18ffPxBl08EAmw4G2NMbGVmTMA3hposW8+4LF1fMi852FP5oEkOgscYzYbopkbaJJNM01no4ra2e/kiSI8FY0RtJRZIWN//X8yQth/9432oWt8x+3aBSBKIpvvDgq3QNxoHMOtH7Ll7CG1aUV2vqclA2wV1ETGA3sA44AmwC3quU2j7qmpuBNUqpj+U6wEKXQg7F0+zpDbO7N8yWI0Ns7w4B4HMbnLWgkZZaD021bs5e0MgFS4pzluTdTx/g5y924XUZfPWd57BibsPUXzQLwWiKTYcCPLvfz6aDQQDOmNvA369bMeGOw3KjlMr5dxNLWURTFkqRSdUUOdf98JYu7njqAOcuauT2684p6vcuB2nb4YkdvfzvC0dGCgJuWL2ID1bA+lcxlVNwXwv8k1Lqbdn3Pw+glPraqGtupsyC+1g7ukPcu6mTFw8HT/jc0tZa3rV6Ea8/vb0gs9pwIs1dfzrAEzv6MA3htqtWjSyMFsvLnYP82293MxBJ4TaFCxY3s2ZpMxee0lyR9fKTcZQqemB3lOJv7nmRrsE4n79yJa85beLTlaqdZTs8vq2HHz61H0fB5Wd28LE3n14xd4yFVk7VMguBzlHvHwEuGee6d4nIG8jM8v9OKdU5zjUlc+b8OXz5HWdxJBjj6GCcQDRz4stvd/Rx0B/j/27YzR1P7Wf1Kc1cuKSZsxc20lI3u40uSime2jPAHU/tZzCexmUIn7h8RdEDO8D5i5v47k2r+c8/7ON3u/t5/mCA5w9mqo7WrZrLzWuXVk29fLEDO8CLh4N0DcZpq/dwyamtRf/+5cRlGlx17gLmzvHx9V/v5IkdfQzF03z+yjMrtpKrEuUyc78BuEIp9eHs+x8ALhk9SxeRViCilEqKyF8BNyqlLhvnsW4BbgFYsmTJhYcOHZrRoPO5QzVtO2zc1ceDL3VxJBg/7nNuU+ho8LGouYZV8+dw5vw5LGiqIZq0iCQtfG6TJRPkVdO2w/qNe3lyZx8AZy2Yw8fevDwv9duz5Y8keeFwkM0Hg2w6GMByFA0+Fze/ZimXnNrKHJ8LEWEwluKlzkF29oS57IwOzphX2DRSpXr+QIBv/WYX8bTNBy49hfesWVzqIZWNnT0hvvLIdsJJi3VnzuVvL1tetSmatO3w6pEhXGbm4JGWOs+4bRIqKi0z5noTCCilGid73HJrP6CU4shgnBcOBtl8KMAhf4zBHA7tXXNKMx9au5SlbcfOwRyKp7n9sR3s6A7hdRl8+HXLeOtZc0syo5zKkWCM7/9+33G19rUek6YaN0eHjh2aXOc1+dYN55XFH6dyoZTiwZe6uPvpgyjgdcvb+LvLVxS8P06l2dMb5nMPvkrKcvjz1yzl+tWLSj2kvOsMxvi/v9k1UpI8bOW8Bi5b2cHrlrfR4MvEq3IK7i4yqZa3kKmG2QS8Tym1bdQ185VS3dm33wl8Vil16WSPW27BfTzxlE1vKMH+gSjbu0Ns7w4RiCSp97lo8LrpGowTT9sIsPa0VlpqPRiG8Ox+f/Y8Rw//cNWqvFbEFIJSit/v7ufhLUc5Esz8TJC5czlnYSNJy2Hb0RDzG31864bzeOy+u49r8eof6K+YA0TyZWd3iHueP8zLnYMAvP+SJdy4ZnHVzkpn6097B/j6r3ciwOfffiZrl1VH6sqyHR7f3stdfzpAynJoq/cyd44XfyTFQCQ5UrHnMoTFLbV0NGRm9Ke21fGZK1bO6HvmuxTy7cB3yJRC3qWUul1EvgJsVko9LCJfA94BWEAA+Gul1M7JHrMSgvtUhuJp7t10mF9v7Tmu7BJgxdx6vvj2VRVVWw6ZQB9KWJk+L001+Nwm8ZTN5x54hf0DUer2PcH2//0Op61YyZ33PQLAh2+8hn27d/KFf/5m1Qf4I8EYd/7xAC8cyizM13pM/s9lp/Pa5SfvAmqu7t/cyU+ePYTHNHjX6oVcd8HCitnkFUtZmbv5WIrBeJreUJKdPSH29EVIWZlzgi9b2cFfvWHZyM+USNs8s9/Pkzv72NI5eFyrlNM76tnw92+c0VhO2k1MpdAzlOClzuDI1vg6j8kbVrSP2+vl3rvvqMhZ70AkySfv30J/Xx9DP/8Hwt0HjztWbTjYz+ZYtXLXF07wqfu3EIyl8bkNrjl3Ae+8YOHI7bY2OaUU3//9Pn61tQeAOT4XN160mKvOWVCWlTTbjg7x+9397OgOccgfY6JIuai5hvdfcgqvm+QPfDiRpnsoQX84ydGhTOz65FvPmNG4dHAvQ/fefQdfve3TFTvr3dsX4R8eepVQ0E/PXbdixTJ5+ubWNh544pmqDuzRpMVnfv4KhwMxzlnYyGevWFlxz79yse3oED9+5hA7sntPzlnYyKfeekbZ3OUe8kf5yTOHRqrJIJNWOaW1ltY6L021blrqPJze0cAZ8xqm/Twop1JILU/WXX0d9/30R+zbvZPrL18LHJv1rrv6urx9n0LdHSzvqOff33sB33jweUafrpmyHHoGE7S05r7hqJKkbYev/WoHhwMxFjfX8IW3n0l9GbSyqFRnLWjkG9efw/MHA6zfuJdXu4b4+L0v8cm3nsH5Mwx4sxVNWmw6GOBP+wZ4/kAARzFyd3bhKc0s76ivuK6r+hlaRK1t7dx53yNcf/lagv4BIDPrzWc6Y/ju4L6f/uiEuwNg1gHeTIbZd/dncGJDGLWZgqjoUIAP3nAVZ33kW1y4cinnLWpi9ZLmqqibt2yH7zyxhy1HhmiqcfOP15ylA3seiAiXnNrKio4GvrVhF68cGeK2X2xlfqOPsxc2ct6iJl63vDD9aVKWw/MHA+ztC9MTStIXSnBgIDqybmYawlVnz+PGixbTXIQW0YWin6VVptB3BxsefYj9ezKn3X/rv37OK0eG+PbfvZ+howc4/MKTBNXVPLGjD5/b4J3nV9ai2ViJtM03fr2TzYeC+NwGt129irklPEO0GjXXefjKO87mZ5s7efClLrqHEnQPJdiwvZen9vTz2StW5mXjk1KKPX0RntjRyx/29BNN2sd93pBMemjtslZec1prVbQu1jn3IvIP9I/Mogu5GOkf6D/h7iCfOfHx0j6PP/ogl1z1PrZ0DvLC4eBI3XxjjZs/u+QU3nbW3IpK2YQTab7y6HZ29oRp8Ln4p2vOKngvoJOdnT0GcWvXEPe/cIRI0uKipc2z2tkaTqT53a5+frO9h4P+Yx1Tl7fXc/GpLcxv9NExx8eippqi3WnqnHsV2vDoQ+zbvXPcBdVSVsxMN0c/9uOtbe287+ZbADitvZ7rVy9i29Eh7n76IDt7wnzvd3vZ1Rvib960vCK2n/sjSb708DYOB2K0N3j58jvOYrHevFVwpiHZFt0NnLe4idt+sZVNB4Pc/tgOPn/lynFz3pbt8Pvd/Ty+LVOBM7fRR0eDj4FIkj29YY4E4yNVLnN8Li5b2cFbVs49btNhtdIz9yIrdCnkdO8OClnBo5Tid7v7Wb9xLynL4awFc/j8lWdS6zFJWQ5u0yi73Zw9oQS3PbSVnlCCxc01fOXas2mrglv0SnRgIJqpzkpYzJvj4yOvP5WLs317BmMp/rTPzwMvHjnuSMqxXIZw9sJG3nbWPC45taUsJhdls0O1UE7W4F5o0w3WxUgV7e2L8M+/3E4gmjru4y5DWL2kmdcub+OSU1tKfpjK4UCM236xlUA0xfKOev7pmrP0c6zEDvmjfPPxXRzKHkJz7qJG4imbPX2RkWsWNtVww4WL6Gjw0htK0BdO0ljjZsXcBk5tqyuLgD6aDu6T0MF9ctO9Oyh0jh4yqY5v/mYX246GcBmC2zRIpO2RW2af2+Ajr1/GujOLn5tPWjYPvdTF/S8cIWk5nL1gDrddvapiF4KrjWU7PLa1m3ueO0wsdXxrjMvPnFtxpz7p4D4JHdzzqxjBfZjtqJEXYiCa4pl9Azy1d4BtRzMbWl5/ehu3vml5UWbxtqN4et8Adz99cOTW/o0r2vnby5ZXXE3zySAYTfHHvQPMa/RxzsLGijwfGHRwn5QO7vlTrAqeqTy5s5fv/34fibRDR4OXj7x+GZec2nLcLD4YSxFL2sTTNo5SLGurO+6A9Ik4SrG/P0o8beN1GbhNg+cPBvj11m4GIplU0dLWWj7y+mWcu6g0m2i0k4eultGKolwqeC5bOZcz5s7hm7/Zyb7+KLc/toOzFszh2vMXsqc3zNP7/CPncw5b2FTDn116Cq89rfWEVI7tKLZ3h3h63wDP7PPjH5PvH/0Y156/gLeumldRt/aVrFL7K1UaHdxPcsMvqNEvtjvve6QkL7aFzTV884bz+PXWHv5n02G2HQ2NpGsAatwmzbVufB6TcMKiazDON369k+Xt9Zy/uInmOje1bhdbjw7x/MEA4YQ18rVt9R7mzvGRTDskLJtFzTW8/ez5nLe4qSz77FerQu+g1o7RaRmt4GYyU4smLe5/4QhbOgc5Y14Da09r5ewFjSOza8t22LCjl3uf7yQQG39WPr/Rx2tOa+U1p7Vxekd9RW2iqlblkgYsJZ2W0arCTGdqdd7MsX8TcZkGV549nzef0cHT+/z0hxMEY2lCiTRLW+u4dFkri5trqjKgV3Jaoxj9lbQMHdy1gip0rxuf2+SylR2zfpxKodMaWq7Kq7pfqzrDM7Xm1jaC/gGC/gE9U5uFdVdfx2krVo78sbz+8rUjC+L5bBtdKMNpmeHnwfDz4sM3XoN/oL/Uw6sqOrhrWp7de/cdxwUq/0A/9959R14eu9L/WI6uznrgiWd44IlnRv5YbXj0oVIPr6rotMwsVXL+sxjGztSAkZlapQSk6dBpk8mVU3VWtdMz91kYfiEP31IOB7Kv3vbpvM3UhqVth0jSIpG2SdsOSmXOa7UdRdp2SFkOScsmkbaJpSxC8TTBWIpANIk/+28wliKWsnCKWCE1k5laIWe+hVbotEk1pDVuuvkjx/1Rb21r14G9APTMfRZyWSxMWjaWfXwwnaiAw1GAyuyotJXCcRRkr/Vla7yTlkM8bRNJ2hgimIZgimAYgiGCYYLHZeIxDWrcJm7TQEQwJFOC5Y+mGIqnsbOnzohkvucwBbgMg1qPmZeGS9OdqVX6zLfQ1SD52HRmOwrLcRj9N95RCpV97jkKbMdBMfL0G3mKjH7fYxr43Pl5nkxG3x3PjA7uszDuC7mljX/98QNQM4eBSJIGn4um2mP1+AqFApTDcS8eEY4F62xjLY/LwGUIdV5X3vpoLGyuRSlFynZI2wrLzrzITTPzRyJtO/gjKfrCCUKJ9AkvbgBThHqvK6et/zB+//eJXpjFOme2Uk3nj2XSsoml7JE7NUFQKFyGUON24XYdm2UYhuA2BJcpeEwDr9vEZQgu08BtZp6ThgiWnXnuJC2bwVgafyTJUDw9MmExRfC5TbwuIy9lqJX+x76U9CamSTgqk/JI25n0h+1kQrMIqOykOuDv5y+ueSODAT8ALa1t/OoPz7NsyQIafK6KbkCVshws59jPPzy7iyYtjgTjWI6Ny8gEgeE/SPnYwl/MRma5mM7McTabdDLPNQdr+P83w0EZTMOgzmMe9wdVZWfZo9NzaXt4+gB1HpO5c3zU+9y4zczvJ1+/o9GSlp1NCzrEkhaBWIqhmIXKvkgMyfzB8GT7+kxHNW560puYJiEiJC2bQNTJfoDjp5ZjHJshy8gTP/Ne9guz0XrsbaghUOdx01zrosaTmY24zMxs2jCEwEA/f339exgM+GlvzzzJ+vv7+fN3X83GjRvx1ld2/bXHZeAZZ1mmvcHLkpZawgkr08wrbZFIZ9YELNsByW9qZyrTvW2fzvXTnTlOJ21i2Q4JyyGRzrSxrfOY1Hpd1LgN6jyZOyPTyDxLA9EUveEEqezd1PCz2G0amWZoLqGp1kudx4XPY1LrMYvWNdHrMvG6TBoA6r0saa3DcRSxtE0ynbl7CCfThOIW4WR6ZPCZ15hk7hAMA5eZeXv0jH8maS6dxsnIKbiLyBXAvwEmcKdS6utjPu8FfgJcCPiBG5VSB/M71GPmNfpoqnWP5KgVx247jx/48flt21ZYjkIplQnS2VQEmYdBKTVy+2mI4DZl0lvLn/7iQXZs386qVavYuHEjAG9+85vZvn07999/P7feems+f+yyYhhCY62bxtEpJ6Uys7eUzWAsRX84STiRHpWvzeT+PS4Dr8scdwY53eqa6Qbf6V6fS5poeAatlOJdH/hL0rbDZVdei6+hGUcp/vXHD/Dbx37Bund/kEA001pYAV6XQWOtm+X19TT4Jk+9Ndd5WNZeRyLtIMLI3VK57sA1jEzqrt7ronXUxy07M8MfvstIWjbxVOYPQDxtE05kJ2zZV7XLMIinrHG+w/h0GueYKdMyImICu4F1wBFgE/BepdT2Udf8DXCuUuqjInIT8E6l1I2TPe5s0jLl5Hvf+x7vfve76ejIzNL7+vqqPrBPRyJtE01a2Eph25nUQThhMRhPk7adkVno8H8fvOcu/v2fP8/S5Wfwbz95ABHh4x94Jwf27uIzX/4X3nvzRzCEkaA23dv2mdzmj00TNbW0cdcjv6O5NXPtcNrEZWTaIhhGdjaazV97sukQt8vAnZ2hDt8FasdznGxOP+2QtG0Odnbzvne+nX27d9LUkj1iL+Bn6fIz+Ncf/5ymUb8DyKRJ/+5D7+LQ3l00tbRlrx/g1OVn8N17HqKltX1kQqjUiTf8SqlxkwDCsefccMzMTAgzjzX6eTz8ufEeA8BWivZ6L+fMsL103vq5i8ha4J+UUm/Lvv95AKXU10Zd83j2mmdExAX0AO1qkgevluCuzVzKcjJrGdnAP1wldOd/fp+rr30njS3tpG2H3t4+Hv3FA1z//r8kZduksovAw2m1oL+fPx+17tHU0sbdj/x+JPiOvqETIOjv50NXv5HBwHCwbuWuUdePfoGOt67S2tbGU8+/yOIF8/HpIF1Q3/ve9/jYxz7GqlWrePLJJ7Ecxbq3vIUdO7bz7e/8G7d89G9GKnyc7LpQT18vb7xkzUhpaEtrG7/83XO0tLfjZA+LGS5gEGTk7n74RlJEjrvjV+pYQB9tuPjBEBmpSDOyXzt8Jzf8HDJGPaYgeFwGNZ6Zpc3ymXNfCHSOev8IcMlE1yilLBEZAlqBgTGDugW4BWDJkiU5fGutmk10OPbnPvmJ495f1l7P2rM/NfK+Upn0WqakT9HjTeMyjj2WacDStlra2o+dcD88ywKoV1FGx2PTMFg5fw4d7Y0jKbnMQiAEBgbGXVe54eor2LhxI/Udlb2uUu6G74BH3x3/7ncbJ707tmM+Rmf8TENY1lFPR0dDwcdbVpRSk/4DbiCTZx9+/wPA+jHXbAUWjXp/H9A22eNeeOGFStNmq7e3V61atUoBqr29XbW3tytArVq1SvX29s76+vXr1x/3+dFfv379+mL8iNo0TPf3W4mAzWqKuK2UymmHaheweNT7i7IfG/eabFqmkczCqqYV1P3338/27KL21q1b2bp1K6tWrRpZ1J7t9bfeeivr169n48aNdHR00NHRwcaNG1m/fr1eVylD0/39VrWpoj+Z1M1+4FTAA2wBzhpzza3AD7Jv3wT8bKrH1TN3LV/Wr19/3Kyst7d30ln1dK/XKku1/37Jceae0yYmEXk78B0ypZB3KaVuF5GvZL/JwyLiA34KXAAEgJuUUvsne0y9oKppmjZ9ed3EpJR6DHhszMe+NOrtBPDu6Q5S0zRNKwxdw6VpmlaFStZbRkT6gUMz/PI2xpRZlolyHReU79j0uKZHj2t6qnFcpyilpmyqU7LgPhsisjmXnFOxleu4oHzHpsc1PR2cFrIAAAPYSURBVHpc03Myj0unZTRN06qQDu6apmlVqFKD+w9LPYAJlOu4oHzHpsc1PXpc03PSjqsic+6apmna5Cp15q5pmqZNouKDu4h8UkSUiLSVeiwAIvLPIvKKiLwsIr8RkQWlHhOAiHxTRHZmx/agiMysmXSeici7RWSbiDgiUvKqBhG5QkR2icheEflcqcczTETuEpE+Edla6rEME5HFIrJRRLZnf4cfL/WYAETEJyLPi8iW7Li+XOoxjSYipoi8JCKPFvL7VHRwF5HFwFuBw6UeyyjfVEqdq5Q6H3gU+NJUX1AkG4CzlVLnkjl85fMlHs+wrcD1wB9KPZDswTTfA64EVgHvFZFVpR3ViLuBK0o9iDEs4JNKqVXApcCtZfL/KwlcppQ6DzgfuEJELi3xmEb7OLCj0N+kooM78K/AZ5j0BNXiUkqFRr1bR5mMTSn1G6XU8Hllz5Lp7llySqkdSqldpR5H1sXAXqXUfqVUCrgXuLbEYwJAKfUHMn2byoZSqlsp9WL27TCZgLWwtKOCbH+tSPZdd/ZfWbwORWQRcBVwZ6G/V8UGdxG5FuhSSm0p9VjGEpHbRaQTeD/lM3Mf7S+AX5V6EGVovINpSh6sKoGILCXTOPC50o4kI5v6eBnoAzYopcpiXGQaMH4GcKa6cLZyahxWKiLyBDBvnE99EfgCmZRM0U02LqXUL5RSXwS+mD2S8GPAP5bDuLLXfJHM7fQ9xRhTruPSKpeI1AM/Bz4x5s61ZJRSNnB+dm3pQRE5WylV0vUKEbka6FNKvSAibyr09yvr4K6Uuny8j4vIOWT6y2/JHlq7CHhRRC5WSvWUalzjuIdMN82iBPepxiUiNwNXA29RRayBncb/r1LL5WAabRQRcZMJ7PcopR4o9XjGUkoNishGMusVpV6Mfi3wjmwLdR8wR0T+Wyn1Z4X4ZhWZllFKvaqU6lBKLVVKLSVz+7y6GIF9KiJy+qh3rwV2lmoso4nIFWRuB9+hlIqVejxlahNwuoicKiIeMgfPPFziMZUtycysfgTsUEp9u9TjGSYi7cPVYCJSA6yjDF6HSqnPK6UWZWPWTcCThQrsUKHBvcx9XUS2isgrZNJGZVEeBqwHGoAN2TLNH5R6QAAi8k4ROQKsBX4pIo+XaizZBeePAY+TWRz8mVJqW6nGM5qI/A/wDHCGiBwRkb8s9ZjIzEQ/AFyWfU69nJ2Vltp8YGP2NbiJTM69oGWH5UjvUNU0TatCeuauaZpWhXRw1zRNq0I6uGuaplUhHdw1TdOqkA7umqZpVUgHd03TtCqkg7umaVoV0sFd0zStCv0/xQD7BhkC3sMAAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f526bc794e0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rate_samples = np.exp(f_samples[:, :, 0])\n",
"\n",
"line, = plt.plot(xtest, np.mean(rate_samples, 0), lw=2)\n",
"plt.fill_between(xtest[:,0],\n",
" np.percentile(rate_samples, 5, axis=0),\n",
" np.percentile(rate_samples, 95, axis=0),\n",
" color=line.get_color(), alpha = 0.2)\n",
"\n",
"plt.plot(X, Y, 'kx', mew=2)\n",
"plt.ylim(-0.1, np.max(np.percentile(rate_samples, 95, axis=0)))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:09:15.817633Z",
"start_time": "2018-06-13T10:09:15.796952Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>GPMC/V</th>\n",
" <th>GPMC/kern/kernels/0/lengthscales</th>\n",
" <th>GPMC/kern/kernels/0/variance</th>\n",
" <th>GPMC/kern/kernels/1/variance</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>[[-1.21828836385], [-0.319760275612], [0.07058...</td>\n",
" <td>0.561140</td>\n",
" <td>0.056633</td>\n",
" <td>2.669075</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>[[-1.21828836385], [-0.319760275612], [0.07058...</td>\n",
" <td>0.561140</td>\n",
" <td>0.056633</td>\n",
" <td>2.669075</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>[[-2.222449023], [2.85239361506], [1.037498786...</td>\n",
" <td>0.048587</td>\n",
" <td>0.034132</td>\n",
" <td>1.617707</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>[[-0.586121746297], [0.587378397531], [-1.3408...</td>\n",
" <td>0.450445</td>\n",
" <td>0.033564</td>\n",
" <td>0.523221</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>[[-1.33457942445], [1.0519000803], [0.08972336...</td>\n",
" <td>0.248194</td>\n",
" <td>0.293047</td>\n",
" <td>1.035331</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" GPMC/V \\\n",
"0 [[-1.21828836385], [-0.319760275612], [0.07058... \n",
"1 [[-1.21828836385], [-0.319760275612], [0.07058... \n",
"2 [[-2.222449023], [2.85239361506], [1.037498786... \n",
"3 [[-0.586121746297], [0.587378397531], [-1.3408... \n",
"4 [[-1.33457942445], [1.0519000803], [0.08972336... \n",
"\n",
" GPMC/kern/kernels/0/lengthscales GPMC/kern/kernels/0/variance \\\n",
"0 0.561140 0.056633 \n",
"1 0.561140 0.056633 \n",
"2 0.048587 0.034132 \n",
"3 0.450445 0.033564 \n",
"4 0.248194 0.293047 \n",
"\n",
" GPMC/kern/kernels/1/variance \n",
"0 2.669075 \n",
"1 2.669075 \n",
"2 1.617707 \n",
"3 0.523221 \n",
"4 1.035331 "
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"samples.head()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"ExecuteTime": {
"end_time": "2018-06-13T10:09:15.921849Z",
"start_time": "2018-06-13T10:09:15.818902Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7f526b837a90>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAD8BJREFUeJzt3X+MZWV9x/H3p6wosimLYiZkl3RJSmwoW22dEBqSZhbaFMUIfxiDpXaxJJsmVLGS6NL+QfsHCabBHyWtzUSo23TjlqrNEtFWgtyQJgVllbrAat3gIkuQ1SDoKNGufvvHHJ0Z3GV2z5nZO/PM+5VM9p5zz7nPd74Mn/vMc8+9k6pCktSuXxl3AZKk5WXQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhq3btwFAJx11lm1efPmXuf+8Ic/5PTTT1/aglYpezHHXsyxFwu11I+9e/d+t6pes9hxKyLoN2/ezEMPPdTr3NFoxNTU1NIWtErZizn2Yo69WKilfiR54niOc+lGkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIatyLeGTvEvqee55odd/c+/+Atly9hNZK08jijl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNWzTok9yR5HCSR+bt+9skX0vy1ST/nmTDvPtuTHIgydeT/OFyFS5JOj7HM6P/OHDZi/bdA1xQVb8F/C9wI0CS84GrgN/szvmHJKcsWbWSpBO2aNBX1f3Asy/a9/mqOtJtPgBs6m5fAeyuqh9X1TeBA8CFS1ivJOkELcUa/Z8Cn+tubwSenHffoW6fJGlMBv0pwSR/BRwBdvU4dzuwHWBiYoLRaNSrhonT4IYtRxY/8Bj6jrsSzczMNPX9DGEv5tiLhdZiP3oHfZJrgDcDl1ZVdbufAs6Zd9imbt8vqappYBpgcnKypqametVx26493Lqv//PVwav7jbsSjUYj+vaxNfZijr1YaC32o9fSTZLLgPcBb6mqH8276y7gqiQvT3IucB7wxeFlSpL6WnQqnOQTwBRwVpJDwE3MXmXzcuCeJAAPVNWfVdWjSe4EHmN2See6qvrpchUvSVrcokFfVW8/yu7bX+L4m4GbhxQlSVo6vjNWkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1btGgT3JHksNJHpm371VJ7knyje7fM7v9SfJ3SQ4k+WqS31nO4iVJizueGf3HgctetG8HcG9VnQfc220DvBE4r/vaDnx0acqUJPW1aNBX1f3Asy/afQWws7u9E7hy3v5/rlkPABuSnL1UxUqSTlzfNfqJqnq6u/1tYKK7vRF4ct5xh7p9kqQxWTf0AaqqktSJnpdkO7PLO0xMTDAajXqNP3Ea3LDlSK9zgd7jrkQzMzNNfT9D2Is59mKhtdiPvkH/TJKzq+rpbmnmcLf/KeCcecdt6vb9kqqaBqYBJicna2pqqlcht+3aw637+j9fHby637gr0Wg0om8fW2Mv5tiLhdZiP/ou3dwFbOtubwP2zNv/J93VNxcBz89b4pEkjcGiU+EknwCmgLOSHAJuAm4B7kxyLfAE8Lbu8M8CbwIOAD8C3rkMNUuSTsCiQV9Vbz/GXZce5dgCrhtalCRp6fjOWElq3OCrbla7zTvu7n3uwVsuX8JKJGl5OKOXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJatygoE/yF0keTfJIkk8keUWSc5M8mORAkn9NcupSFStJOnG9gz7JRuDdwGRVXQCcAlwFfAD4UFX9OvA94NqlKFSS1M/QpZt1wGlJ1gGvBJ4GLgE+2d2/E7hy4BiSpAFSVf1PTq4HbgZeAD4PXA880M3mSXIO8Lluxv/ic7cD2wEmJibesHv37l41HH72eZ55oV/9Q23ZeMZ4Bj6GmZkZ1q9fP+4yVgR7McdeLNRSP7Zu3bq3qiYXO25d3wGSnAlcAZwLPAf8G3DZ8Z5fVdPANMDk5GRNTU31quO2XXu4dV/vb2OQg1dPjWXcYxmNRvTtY2vsxRx7sdBa7MeQpZvfB75ZVd+pqv8DPg1cDGzolnIANgFPDaxRkjTAkKD/FnBRklcmCXAp8BhwH/DW7phtwJ5hJUqShugd9FX1ILMvun4Z2Nc91jTwfuC9SQ4ArwZuX4I6JUk9DVrcrqqbgJtetPtx4MIhjytJWjq+M1aSGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDVuUNAn2ZDkk0m+lmR/kt9N8qok9yT5RvfvmUtVrCTpxA2d0X8E+I+q+g3gdcB+YAdwb1WdB9zbbUuSxqR30Cc5A/g94HaAqvpJVT0HXAHs7A7bCVw5tEhJUn9DZvTnAt8B/inJV5J8LMnpwERVPd0d821gYmiRkqT+UlX9TkwmgQeAi6vqwSQfAb4PvKuqNsw77ntV9Uvr9Em2A9sBJiYm3rB79+5edRx+9nmeeaHXqYNt2XjGeAY+hpmZGdavXz/uMlYEezHHXizUUj+2bt26t6omFztu3YAxDgGHqurBbvuTzK7HP5Pk7Kp6OsnZwOGjnVxV08A0wOTkZE1NTfUq4rZde7h135Bvo7+DV0+NZdxjGY1G9O1ja+zFHHux0FrsR++lm6r6NvBkktd2uy4FHgPuArZ1+7YBewZVKEkaZOhU+F3AriSnAo8D72T2yePOJNcCTwBvGzjGirV5x929zz14y+VLWIkkHdugoK+qh4GjrQ9dOuRxJUlLx3fGSlLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcYP+OLhWp8077u597sFbLl/CSiSdDM7oJalxg4M+ySlJvpLkM932uUkeTHIgyb8mOXV4mZKkvpZi6eZ6YD/wq932B4APVdXuJP8IXAt8dAnGaYrLJ5JOlkFBn2QTcDlwM/DeJAEuAf6oO2Qn8NcY9EvqWE8SN2w5wjUDnkAktWno0s2HgfcBP+u2Xw08V1VHuu1DwMaBY0iSBug9o0/yZuBwVe1NMtXj/O3AdoCJiQlGo1GvOiZOm53J6uT0ou9/p5NtZmZm1dS63OzFQmuxH0OWbi4G3pLkTcArmF2j/wiwIcm6bla/CXjqaCdX1TQwDTA5OVlTU1O9irht1x5u3edVojAb8svdi4NXTy3r4y+V0WhE35+p1tiLhdZiP3ov3VTVjVW1qao2A1cBX6iqq4H7gLd2h20D9gyuUpLU23JcR/9+Zl+YPcDsmv3tyzCGJOk4Lcnv+VU1Akbd7ceBC5ficSVJw/nOWElqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa1zvok5yT5L4kjyV5NMn13f5XJbknyTe6f89cunIlSSdqyIz+CHBDVZ0PXARcl+R8YAdwb1WdB9zbbUuSxqR30FfV01X15e72D4D9wEbgCmBnd9hO4MqhRUqS+ktVDX+QZDNwP3AB8K2q2tDtD/C9n2+/6JztwHaAiYmJN+zevbvX2IeffZ5nXuhXd2smTmNF92LLxjNO2lgzMzOsX7/+pI23ktmLhVrqx9atW/dW1eRix60bOlCS9cCngPdU1fdns31WVVWSoz6TVNU0MA0wOTlZU1NTvca/bdcebt03+Ntowg1bjqzoXhy8euqkjTUajej7M9Uae7HQWuzHoKtukryM2ZDfVVWf7nY/k+Ts7v6zgcPDSpQkDTHkqpsAtwP7q+qD8+66C9jW3d4G7OlfniRpqCG/518MvAPYl+Thbt9fArcAdya5FngCeNuwEiVJQ/QO+qr6LyDHuPvSvo8rSVpavjNWkhpn0EtS41butXhqzuYdd/c+9+Atly9hJdLa4oxekhpn0EtS4wx6SWqcQS9JjfPFWK0KJ/pC7g1bjnBNd44v5Gqtc0YvSY1zRq/mDbmsE/yNQKufM3pJapxBL0mNM+glqXEGvSQ1zqCXpMZ51Y2kJeWH1608zuglqXEGvSQ1zqUbSQsMfYOZVh5n9JLUOGf00jIa1wuT88ed/wFvWpuc0UtS45Yt6JNcluTrSQ4k2bFc40iSXtqyLN0kOQX4e+APgEPAl5LcVVWPLcd40nIa14uTa/FF0dV6Df5Kr3u51ugvBA5U1eMASXYDVwAGvaQVp/Un1eVautkIPDlv+1C3T5J0ko3tqpsk24Ht3eZMkq/3fKizgO8uTVWr27vtxS/YizlrpRf5wHEfuqL6cQJ1H82vHc9ByxX0TwHnzNve1O37haqaBqaHDpTkoaqaHPo4LbAXc+zFHHux0Frsx3It3XwJOC/JuUlOBa4C7lqmsSRJL2FZZvRVdSTJnwP/CZwC3FFVjy7HWJKkl7Zsa/RV9Vngs8v1+PMMXv5piL2YYy/m2IuF1lw/UlXjrkGStIz8CARJatyqDXo/YmFOkjuSHE7yyLhrGbck5yS5L8ljSR5Ncv24axqXJK9I8sUk/9P14m/GXdO4JTklyVeSfGbctZxMqzLo533EwhuB84G3Jzl/vFWN1ceBy8ZdxApxBLihqs4HLgKuW8M/Gz8GLqmq1wGvBy5LctGYaxq364H94y7iZFuVQc+8j1ioqp8AP/+IhTWpqu4Hnh13HStBVT1dVV/ubv+A2f+p1+S7smvWTLf5su5rzb4ol2QTcDnwsXHXcrKt1qD3Ixa0qCSbgd8GHhxvJePTLVU8DBwG7qmqNdsL4MPA+4CfjbuQk221Br30kpKsBz4FvKeqvj/uesalqn5aVa9n9t3pFya5YNw1jUOSNwOHq2rvuGsZh9Ua9It+xILWriQvYzbkd1XVp8ddz0pQVc8B97F2X8u5GHhLkoPMLvVekuRfxlvSybNag96PWNBRJQlwO7C/qj447nrGKclrkmzobp/G7N+H+Np4qxqPqrqxqjZV1WZm8+ILVfXHYy7rpFmVQV9VR4Cff8TCfuDOtfwRC0k+Afw38Nokh5JcO+6axuhi4B3Mztge7r7eNO6ixuRs4L4kX2V2cnRPVa2pywo1y3fGSlLjVuWMXpJ0/Ax6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIa9/8IwDPEdvNpagAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f526bc55d68>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"samples['GPMC/kern/kernels/0/variance'].hist(bins=20)"
]
}
],
"metadata": {
"anaconda-cloud": {},
"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.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Computing file changes ...
