{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Refitting PyMC models with ArviZ\n", "\n", "ArviZ is backend agnostic and therefore does not sample directly. In order to take advantage of algorithms that require refitting models several times, ArviZ uses {class}`~arviz.SamplingWrapper` to convert the API of the sampling backend to a common set of functions. Hence, functions like Leave Future Out Cross Validation can be used in ArviZ independently of the sampling backend used." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below there is an example of `SamplingWrapper` usage for [PyMC](https://www.pymc.io/projects/docs/en/stable/learn.html)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import arviz as az\n", "import pymc as pm\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the example we will use a linear regression model." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(4)\n", "\n", "xdata = np.linspace(0, 50, 100)\n", "b0, b1, sigma = -2, 1, 3\n", "ydata = rng.normal(loc=b1 * xdata + b0, scale=sigma)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7kUlEQVR4nO3deXxj1Xn4/8+RZHmT93332B5mnwFmGHYCJJBAIJCEkCahof0mgTbpkqZNmub7ev3SJfl+06ZtkrZJvyU0LaFpAoHskAAhULYZmBlmgdnwPuNVtrxKtiRbOr8/7r2y5VW2pbFlP+/Xi9dIV1e6R2AeP/Pc55yjtNYIIYRIPrbVHoAQQojlkQAuhBBJSgK4EEIkKQngQgiRpCSACyFEknJcyIsVFhbq2traC3lJIYRIekeOHOnXWhfNPH5BA3htbS2HDx++kJcUQoikp5Rqn+u4lFCEECJJSQAXQogkJQFcCCGSlARwIYRIUhLAhRAiSUkAF0KIJCUBXAghklRMAVwp1aaUekMpdUwpddg8lq+UekYp1Wj+mZfYoQohxNrVNTTOT452XtBrLiUDv0FrfbHWep/5/PPAs1rrzcCz5nMhhNiQ/vX5Zj79yDGeOtlzwa65khLKHcBD5uOHgDtXPBohhEhSB1s8AHzxpycZ9U9ckGvGGsA18LRS6ohS6j7zWInWutt83AOUzPVGpdR9SqnDSqnDfX19KxyuEEKsPf3eAI1uL7fuKqV31M/f/ersBblurGuhXKO17lRKFQPPKKXOTH9Ra62VUnPuzaa1fgB4AGDfvn2yf5sQIun0jQYIhTWlOWlzvm5l3/ddV09xVhoPHWjjzkvK2VuTn9BxxZSBa607zT/dwI+B/UCvUqoMwPzTnahBCiHEavr84yf4nf94bd7XD7Z4cKU62FmezZ+9cwtl2Wl8/vE38E+EEjquRQO4UipTKZVlPQZuBt4Efgbca552L/DTRA1SCCFW01vuUc70jNLa75vz9YMtA1xWm4fDbsOV6uD/vG8XjW4v//fJ0wkdVywZeAnwklLqOPAa8ITW+lfAV4CblFKNwDvM50IIsa5MhMJ0DfkBeObU7A4T96ifJreXK+oKIseu31LMx6/ZxEMH2vnVm92z3hMvi9bAtdYtwJ45jnuAtydiUEIIsVZ0DY0TChu375451ct919VHvf5qywAAV9YXRB3/3Lu2cqhtgM8+doId5TlU5WfEfWwyE1MIIRbQ7hkD4Mq6Ao60D+LxBqJeP9jiISvVwfay7KjjToeNf/7QpaDhD79/lIlQOO5jkwAuhBALODdgBPCPXbOJsIZnz0T3axxo8XDZpnwc9tnhtLogg6+8fzeNvaOc6R6N+9gkgAshxALODYzhdNi4YWsx5TlpPHOqN/Kae8RPS5+PK+sK5n3/u3eX8cLnbmBXZU7cxyYBXAixZj18oI13fu2FVR3DOc8Y1fkZ2G2Kd2wv4cXGPsaDRnvgwVaj/n3FAgEcoMCVmpCxSQAXQqxZLzd5ONs7ylhwctXG0D4wRo15A/Lm7aX4J8K81NTPT4528lc/O0mhy8n28uxFPiUxLuiu9EIIsRTNfV4APN4gGfkXPlxprTnn8XH5JmNG5eV1+WSlOfjMI8cYDUxycVUuX3n/Luw2dcHHBpKBCyHWqMlQmDaPMXGmf0bnx4Uy4AviC4aoKTAy8BS7jXfvKkMDf33HDh7//avYWro62TdIBi6EWKPaB8aYCBn91x5vcNXGAFA9rYf7r+/YyRdv30G6074qY5pOArgQYk1qcnsjjz2+1cnAz5sB3MrAwejvXivWzkiEEGKa6QG8f7UycHMST2Ve/GdRxoMEcCHEmtTs9lKanYYr1bFqNfBzA2OUZqeRlrL65ZK5SAAXQqxJzX1eGopdFLqcq1YDP+cZo7pgbWbfIAFcCLEGaa1p7vNRX5RJgSv1gtXA3SN+vIGpnvP2AV/UDcy1RgK4EGLN6TEDaUOxi4JMJ/2jFyYD/+ADB/m9h4+gtcY/EaJ3JBCZxLMWSQAXQqw51g3M+mIXhVmzM3BvYJIn34jvOtvuET+t/T5eaurnhcb+SAeKlFCEEGIJrADeUOyiMNPJgC8YWZMb4Mevd/DJ773OW73xW+HveMcwAOkpdr7yyzO0eWb3gK81EsCFEGtOk9tLdpqDIlcqBa5UwhoGx6bKKFZwPWEG3Xg4fn4Iu03xl+/ZzunuEb71fBMANQWZcbtGvEkAF0KsOVYHilKKApcTiJ6NaZU33uyMYwDvGGJLSRYf2FvFjvJsjp4bwpXqIC8jJW7XiDcJ4EKINafJ7aOh2AVAobkU6/SdcM4PjgPwRpwCeDisOX5+iD1Vudhsis/fshUwyidKrc5CVbGQqfRCiDVleGyCfm9gWgA3MvA+M4BrrekwM/BTXSOEwnrFqwG2eXyM+Ce5uMrYdOHazUW895IKqvLSV/S5iSYBXAixpjT1GTcmrQBekGll4EYJZXh8gtHAJNvLsjnVPUJLn5fNJVkruubxjiEA9lTlRo597YMXr+gzLwQpoQgh1pRIB0qREZRz0lNw2FSkldDao/LWXaVAfMoox88Pk+G0s7l4Zb8ILjQJ4EKINaXJ7cXpsFFhli9sNkV+5tR0+vMDRv37+i3FpKXY4hLAj50fYmdFzqptzLBcEsCFEGvKa22D7CzPjgqmBa7UyIJW5wenlnjdXpYd1Yky6p/gq0+dYdQ/EfP1gpNhTnWNcPG08kmykAAuhFgzhsaCnOgY4rqLiqKOF7qckSVlzw+MkZuRQlZaCrsqcjjZNULYnOTzHy+38c3nmvnlmz0LXqexdzQS5M/0jBAMhdlTmRv/L5RgEsCFEGvGS039aG10gUxXOG1Bq3MDY5HZkTsrchgLhmjp9+GfCPGfr7QB8GrLwLzXcI/6ufWfXuSOb75Mu8fH8fNDAFxcnRv375NoMXehKKXswGGgU2t9m1JqE/ADoAA4Avy21np11nwUQqwLL77VT1aagz2VOVHHpy9o1TE4zvYyYx/KXeZ5b3YOc6C5nwFfkOr8DF5t9cx7jadO9jIR0rhHArz3W69QU5BBoSuV8py0BH2rxFlKBv7HwOlpz/8W+JrWugEYBD4Wz4EJIdaGnmE/9zz4atREmkTQWvNiYx9X1xfisEeHpgJXKuMTIbyBSToHx6nMN25wNhS5SHXYOHZ+iG+/2Mol1bn87tW1dAyO02HWymd68kQ39UWZ/PwPryE7zcHRc0NcXJWzpifszCemAK6UqgTeDTxoPlfAjcBj5ikPAXcmYHxCiFV2uH2Al5r6Odw+mNDrNPd56Rr2z6p/w9RknlNdRr26ytzizGG3sa0sm0cOnefcwBj3X1fP5ZsKgLnLKP3eAK+2erh1VxmbCjP58Sev5n2XVPDhy6sT+M0SJ9YM/OvA54Cw+bwAGNJaWyufdwAVc71RKXWfUuqwUupwX1/fSsYqhFgF/aNG5t3W70vodV54qx+AazcXznrNmk5/9JzxS6Rq2gqBuypyGJ8IUVeYyU3bS9hamkVOesqcZZSnTvYQ1nDrrjIA8jKd/OMHL+bGrSVx/z4XwqIBXCl1G+DWWh9ZzgW01g9orfdprfcVFc3+zSqEWD3NfV6+e6BtwXOs7g9rBcBEebGxj02FmVHB2WItaHX03BAQvcTrrgqjDn7fdXXYbQqbTbF/Uz4H58jAf/lGD5sKM9lamlwTduYTSwZ+NfAepVQbxk3LG4FvALlKKesmaCXQmZARCiES5tHD5/n/fnqS4fH5+6at/utEZuCByRAHWwa4bo7sG4waOBgTbpSC8typG4637Snjr+/Ywfv3VkaOXb4pn3MDY3QPj0eODfiCHGjxcOuu0qSsd89l0QCutf4LrXWl1roW+C3gN1rrjwDPAXeZp90L/DRhoxRCJMSAmV2fWyC7tgJ4uyc6gA+PTXDD3z/Pobb5W/ZidaRtkPGJ0Kz2QUtBppGB94z4Kc1OI9UxtUt8htPBR6+sJWXajc8r6mbXwZ8+2UMorLllZ9mKx7tWrKQP/M+BzyilmjBq4v8enyEJIS6UAZ9VHpk/u+4zg3zXsB//RChy/Oj5QVr7fbzY2L/icbzQ2I/DpriivmDO19NS7GSlGn/ht25gLmRbWTZZaY6oOviTb/ZQU5DBjvLsFY93rVhSANdaP6+1vs183KK13q+1btBaf0BrfWG2jRZCxI3HCuALlEf6RwOkpRihwlpICuB0t7FqYGMctjV7sbGPS2vycKXOPzXFqoNbLYQLsdsU+2uNOvhEKMyDL7bwclM/t+wsWzflE5CZmEJsaFMZ+NwlFK01/d5AZJp567RAf7p7BGDF+1J6vAFOdo1wbcPc9W+L1YkS6x6Vl9fl09rv451ff4EvPXGaazcXct91dSsa61oj64ELsYFZAXxmfdviDUwSmAyzrzaPV1sHos6zAnibZ4zgZBinY3n54CvNRpnjmnluYFqsDDyWEgrA1eYvhIlQmAc/uo+3byteV9k3SAAXYsPymzMbYf4auNVCWF/kIi8jJZKp+yeM9Ueq8zM4NzBGa7+PLctszXup0Zg+v3uRxaSsTpS52gznsqM8h1/84TU0FLtIS7Ev/oYkJCUUITYoK/uuzs+g3xuccwlWqwOl0JVKTUFmpFbe2OslFNa8Z0+58dy9vDKK1pqXmvq5qr5g0bW4C81OlKoYauCWnRU56zZ4gwRwITYsK4DvrckDoH2OOrg1C7PQlcqmwszIOVb55N27y7ApeKvXu6wxtPb76Bwa55p52genu3V3GfdfV0dpdvItOpUoEsCF2KCsDpRLzWVU5wzgVgae5aSmIIOu4XH8EyFOdY+QnmLnopIsqvMzlt2J8lKTOX1+kRuYAFtLs/mLW7etuzr2SkgAF2KDGjDX176k2sjA56qD93mDKAX5GU42FWaitbGhwunuEbaUZmG3KTaXZNHoXl4G/lJjP5V56dQUxFbXFtEkgAuxQVl7TFblZVCclTpnL7jHGyAvw4nDbqOmIBMwyh6nu0fYbk6IuajERVu/j+BkOPK+nxztnLezxTIZCnOg2cO1mwslq14mCeBCbFADviAOmyI73UFtQea8JRRrKddNZgB/pdnDiH+SbeamCpuLs5gM60gG39Ln5dOPHONbzzUveP3jHcOMBia5pkEWuVsuCeBCrFP+iRDjwdC8rw/4guRlOlFKUVuYMWcJpd8bjEygyclIITcjhadOGvtNbi8z2gY3l7iAqQk9PzzSAcDr5xZeP/ylxn6UgqvmmT4vFicBXIh16s9+eJzf+6/5V4H2+IKRRaJqCjJxjwbwBSajzjEy8NTI89qCTLqH/QBsKTUy8PoiV6QTZTIU5vEjHdgUNLq9C65y+HJTPzvLc8gzxyCWTgK4EOvUqe4Rjp0fQms95+sDviD5ZvCsNcsjM8so/aMzA7hxs7GmICOybklaip3q/Aya3KO80NiHezTAR6+sBaY2YJhpLDjJ0fODkdmSYnkkgAuxDmmt6RoaZ3h8ItIuONP0AG51gUy/8TgeDOELhijMmsqQawuNQL+tNHpFv80lWbzV6+XRQx0Uupz8yTsuwqbg9Xm2YXutdYCJkJbyyQpJABdiHfL4gvgnjK6Q5nla/DzeQKSEYgXm1mkBfPosTIuVqVs3MC2bi1209vv49ele3ntJBTkZKWwry+bIPBn4gWYPTruNy2rzl/P1hEkCuBDrUOfg1E40TX2zA/hEKMyIf5L8TCM4u1IdFLpSae+fKqH0mQG8aFoAtwL3pTW5UZ93UUkWobBmMqz5wL4q45zqPI6dGyIUnl3Cebm5n0uqc0l3rt9p7heCBHAh1qHOoakA3uye3V0yaJZV8l3TyiMF0Z0o06fRW7aUZvHCZ2+YtXNOQ7HRiXJxVS4XlRjdKXtr8vAFQ5ztiZ6lOTQW5GTXCFfVS/17pSSAC7EOWRl4VX46zXNk4FZdvGBaB0jNjF5wayXC6TVwgOo5Zk02FLuoL8qMWm/bWmNlZhnlQLMHreHqBql/r5QEcCHWoc6hcVypDi6pyqNpjhq4NQszPzM6A+8Z8UeWmLVq4AWZqbPeP1Naip1n//R6bt01td9kZV46ha7UWTcyX2n2kOG0s6cqd8nfS0STAC7EOtQ5NE5FbjoNxS46h8ZnTejx+KzyyFQAv2yTcUPxf872AUYAz0lPWfZGDUop9tbkzprQ83JzP/s35UdtQiyWR/4NCrEOdQ6OU56bRn2RUZueWUaxlpLNn5ZdX1abT6HLyZNvdgPR0+iXa29NHu2eMfrMenrPsJ+WPh9XS/07LiSAC7FGaa15/EhH5IbjUnQOjVORlx65uThXALcpyE1PiRyz2xQ37yjluTNu/BMh+keDUTcwl+NSc6VDKwt/2Vw+9iqpf8eFBHAh1qinT/Xypz88zr8817Sk93kDkwyPT1CRm0FNQQY2Bc190Z0oHl+QvAwnthm74Ny6s4yxYIj/eavPyMCzVhbAd1bkkGJX/M0vTvGnjx7n4YPt5GWkzJoIJJZHArgQa1A4rPn6rxsBePz1DvwT8y9KNZPVgVKRl05aip2q/IxZk3kGvMGoG5iWy+vyyctI4ZdvdNPnDUS2MVuutBQ7X7x9B5sKM3mhsY9j54e4YWvxrF8cYnlkU2Mh1qCnTvZwunuEu/ZW8tiRDn71Zg93XlIx67xwWPPx7x7myroCPmG28HUOGa2AFbnG3pENRa45SyhzBfAUu42bt5fyixNdxjT6FZZQAO65ooZ7rqiJXNdaQ0WsnGTgQqwxVvZdV5TJ/33fLmoKMvjv187Nee5jr3fwmzNuHn+9I3LMysAr84wAXl/soqXfFzUj0uMLUDDPDcpbdpXiM7tWVlpCmSk/07nsrhYxm/ybFGKNefLNbs72jvLHb99Mit3Gh/ZX81rrwKx+bl9gkq8+dRabgjM9o3jMvu2OoXGcdltkCnx9USbByTAdg1OTdObLwAGuqi8kO83IkuORgYvEkQAuRAIdaR/gM48cIzzHeiAWb2CSL/3iFF9+4hT//Gwj//D0WzQUu7htdzkAd+2tJMWu+P6MLPxfn2+mbzTAF27dBhgr/AF0Dfkpy02L1JlndqKEwpqh8YmoFsLpnA4bN20vBVhxG6FIrEUDuFIqTSn1mlLquFLqpFLqr8zjm5RSryqlmpRSjyil5L+0EDM8c8rNj452zrnbjeUnRzt58KVWHj7Yzj888xat/T7+7OaLsJsBuNCVys3bS6NuZnYMjvHtF1u44+Jy7r2qlvQUOwdbPAB0Do5RnpMe+fy6QjOAm2uiDI4F0Tp6Gv1MH768moZiF3VmH7lYm2K5mxAAbtRae5VSKcBLSqlfAp8Bvqa1/oFS6v8BHwP+NYFjFSLpuEeM3WtOd4/OGwx/caKLhmIXz/zJdUyENIHJEFlpKVHnfGh/NU+80c1vPXCQ6vwM2j0+lII/f9dWUuw29tXmcbDFyMA7h8ajFpvKy3RSkOmMlGCmJvHMH8D31uTx68+8bflfXFwQi2bg2mAV31LMfzRwI/CYefwh4M5EDFCIZOY2ZyCe7h6Z+/URP6+2DvDuXWUopXA6bLOCNxj7Rn7k8mpsCo53DNHa7+PT77iIcrPT5Iq6As72jtIz7Mc9Goh0oFjqi1yc6RlBax1ZB2WhDFwkh5j6eZRSduAI0AB8E2gGhrTW1gZ6HcDsHifjvfcB9wFUV1evdLxCJJXeSAY+dwB/8o1utIbb95TN+brFZlN8+b275n39SnNnm58c60Rrowd8uhu2FvO3vzrDt55vjmzKUCA3KJNeTAFcax0CLlZK5QI/BrbGegGt9QPAAwD79u2b/06OEOuQFcBPzRPAn3ijm62lWTQUZ63oOrsqcshw2nnc3BG+ckYGfv91dZztGeGrT51lv7kLzkIlFJEcltSForUeAp4DrgRylVLWL4BKoDO+QxMiufknQoz4J8nNSKF72M/QWPSaJt3D4xxqG+S23Qtn37Ew6uD5NJp17pkZuM2m+Lu79nDdRUW81mbUyvMyZpdqRHKJpQulyMy8UUqlAzcBpzEC+V3mafcCP03QGIVISu4Ro/5t3VCcmYU/ccJY9e/dZrvgSl1RZ2TWSkFZTvqs150OG//6kUvZU5lDWU4aDlnONenF8l+wDHhOKXUCOAQ8o7X+BfDnwGeUUk1AAfDviRumEMmnd9Qon9ywxQjgp7ujtxb7xYludpRns8ncUHilrqwz6uDFWanzznbMTHXwyP1X8vjvXxWXa4rVtWgNXGt9ArhkjuMtwP5EDEqIRBn1T3DPg6/y5ffuYmdFTkKvZdW/t5dnU+hK5VTXVAZ+fmCMY+eH+PN3xXw7aVE7K3LIdNpndaDMlJZij3SviOQmf4cSG8rZnlGOdwxzoNkT18/9+fEu/u5XZ6KOWSWUkqw0tpVlRXWiPHakA6WIS/3bkmK3cf/b6nnfpZVx+0yxtkkAFxuKtVv79F3b4+GbzzXx7y+1Rk2Z7x3147TbyM1IYXt5Nk1uLxOhMBOhMN9/7Rxvu6iIqvzZGwSvxB+9fXNk5T+x/kkAFxtKx2D8A/j5gTHO9IwSmAzTNTz1ue6RAMXZqSil2F6WTTAUprnPy9Mne3GPBvjolRJoxcpIABcbSiQDH4xfAH/6VG/kcWv/1Jon7lE/xeZyrNvKjB1oTneP8PDBNirz0nnbRcVxG4PYmCSAiw2lMwEZ+DOneiKBenoA7x0JUJKdBkBdYSZOh42fHuviYMsAH7m8JrJYlRDLJQFcbChW4B4en8AbmFzk7MUNjQU51DbI3fuqyHTaaembHsD9kQDusNvYUpLF82f7cDpsfPCyqhVfWwgJ4GLD0FrTOThOkZktd8UhC//NGTehsOam7SVsKsqkxczAx4MhRv2TkWsBbCszpsvftqtMprGLuJAALjaMAV+Q8YkQ+zcZMxbjUUZ55lQvJdmp7KrIoa7QRWu/MZXdbU7isTJwMNYrAbhHbl6KOJEALjYMK2BfbgXwFd7I9E+E+J+3+njHthJsNsWmwkw6BscJTIbotXrAs6cy8A/sq+LR+6/k0uq8FV1XCIsEcLFhWAH74qpcUuxqxRn4K839jAVD3LS9BIC6oky0hnOesUgGXpw1lYGnpdgj2b8Q8SABXGwYVsCuysugNCdtwRr48NgEgcnQgp/3qzd7cKU6ImtxW2uatPT75szAhYg3CeBiw+gYHCfTaSc3I4WK3PR5SyiToTDv+sYLfPVXZxf4rDF+crSL2/eUk+qwA1BrBfA+H+4RP06HjZx0WbJVJI4EcLFhdA6NU5GXjlKK8tz0eTPww+2DdA/7eampf97P+qdnGwH4wxsbIsey01IodKXS2u/FPRqgxJyFKUSiSAAXG0bn4Hhkpb7K3HR6RvxMhMKzznv6pDGz8mzvKMPjE7Neb+7z8vjrnXzkiupZq/rVFWXS2u+jd8QfVf8WIhEkgIsNw8rAwdixJqyhZ9gfdY7WmqdP9VDocqI1HD03OOtzvvbMW6Q6bHzy+oZZr9UVTgVwqX+LRJMALjYEb2CS4fEJKnKN1f+szHlmGeV09ygdg+N86oYG7DbF4bboAH6qa4RfnOjmd6+ujZqkY9lUmEm/N8j5gXHJwEXCSQAXG4J1wzKSgZsBfGYr4dOnerApuH1POdvLsjncPhD1+jeefYvsNAf3XVs/53WsTpRgKEyxZOAiwSSAiw2hc2gMmArcVgY+sxPlqZO97KvJp9CVyt6aPI6dH4rUyd0jfn592s2HL68hZ54NgeuKprZHK5EMXCSYBHCxIVjrgFeZGXhaip1ClzNq/e7zA2Oc7h7h5h3GxJx9tXn4J8KRrdB+dLSTUFhz9775d7ypys/AWmRw+jR6IRJBArjYEDoHx3HabRS6psoa5bnpkcAOU+t6WzMr99UYsyYPtw+iteaHh8+zryaPuiLXvNdJddgju+xICUUkmgRwsSF0DI1TnpuGbdoa3BUzesGfPtnD1tIsagqMMkhpThoVuekcaR/g9XNDNPf5+MAC2bfFqoNLCUUkmgRwsSF0Dk61EFoqctPpHBpHa80Lb/XxausAt+yM3mR4X20eh9sG+eHh86Sn2Hn37vJFr7W1NJuc9BSy0x1x/Q5CzCQBXGwInUNTk3gs5bnp+CfCnOkZ5U8eOcaWkizuu64u6px9NXm4RwP86Ggnt+4qw5W6eFD+gxsb+MmnrpZZmCLhJEUQ655/IkTfaCDSA26xMvJPfPcwY8EQ//LhS0h32qPO2WvWwYOT4QVvXk7nSnXEFOiFWCn5KRNrlnvUjz8YprogY/GTF9DuMVsI5yihgNGh8tW7drO5JGvWe7eUZpGV6iDf5ZSlYMWaIwFcrElvdg7z0e+8Rlaag//57A3L/pzu4XF+/7+OkJ5iZ19N9EYKVfkZpNgVt+8u5669c2fXdpviL9+zg6IsWZhKrD0SwMWac6R9gN/5j0N4A5MM+IL0ewNR7X+xOucZ48MPHmR4bIKHP7Y/styrJSc9hac+fR3V+RkLBuf3zxPchVhti97EVEpVKaWeU0qdUkqdVEr9sXk8Xyn1jFKq0fxT9okSK3ag2cM9D75GoSuVf/jAHgDe6Bhe8ud0DY3zgX97BW9gku994nL21c5d/qgrcuGwy718kZxi+cmdBP5Ua70duAL4lFJqO/B54Fmt9WbgWfO5ECvypSdOUZKdyiP3X8HNO0pRCo53DC35c77yyzMMjU3w/U9cwe7K3LiPU4i1YNEArrXu1lq/bj4eBU4DFcAdwEPmaQ8BdyZojGKDGPQFOdU9wvsuraQ4Kw1XqoOGIhcnlpiBHzs/xM+Od/GJa+vYVpadoNEKsfqW9HdHpVQtcAnwKlCite42X+oBSuZ5z31KqcNKqcN9fX0rGatY515t9aA1XGXuMQmwqzKHEx3DaK1j+gytNV9+4hSFLie/d/3cKwYKsV7EHMCVUi7gceDTWuuR6a9p4/+uOf8P01o/oLXep7XeV1RUtKLBivXtQLOH9BR7VMljT2Uu/d4A3TM2XpjPUyd7OdQ2yJ/cdJH0Yot1L6YArpRKwQje39Na/8g83KuUKjNfLwPciRmiSFZDY0EePXyefm8gpvNfafZw2aZ8nI6pH8vdlTkAnIihDh6cDPOVX55mc7GLD+6rWtaYhUgmsXShKODfgdNa63+c9tLPgHvNx/cCP43/8EQyOt09wuceO87l/+dZPvfYCR56pW3R9/SNBmh0e7myriDq+LaybBw2FVMd/JvPNdHmGeMLt26TzhKxIcTyd8yrgd8G3lBKHTOPfQH4CvCoUupjQDtwd0JGKJJK19A4d3zzZexK8f69lTx/xs1bvaOLvu9AiweIrn+DsW73ltKsRQP4K039/NNvGnnfpRXcsLV4+V9AiCSyaADXWr8EzDfL4e3xHY5Idr854yY4GeapT1/HltIs7n/4ME1u76LvO9DcT1aqgx3ls7tGdlfm8sSJLrTWc064cY/6+aMfHKO+yMWX7twZl+8hRDKQv2eKuHr+rJuq/HQuKjE2PWgodtHmGSM4GV7wfQeaPVxelz9n6WNPZQ4j/snImiaj/gmeO+PmcNsAzX1ePv2DY3gDE3zrI5eS4ZQbl2LjkJ92sSxaaw61DXJZbV4kKw5Mhni5ycNdeysjxzYXZxEKa9o9vjkXiwKj7NLmGeOeK2rmfN3qSjneMUS60849D75K44ys/qt37eaieT5fiPVKArhYltdaB/jgAwf5x7v38L5LKyPHxidC3LB1ql20odjIxBvd3nkD+IFmq/5dOOfrm0tcpDpsPHWyh394+i083gD/8uFLyE5LYXAsSHZ6Cjdskbq32HgkgItlOTdglDO+/WIr772kAqUUz53pw+mwcWXdVCCuL3KhFAvWwV9p9pCXkcLW0rkDfIrdxo7ybJ58o4ec9BS+94kruLgqN67fR4hkJDVwsSw95sSa090jvGJm0M+/5eaKuoKoTRHSnXYqctNnlTwsE6Ewz591c1VDYdR+lTO9fVsJ5TlpPHK/BG8hLJKBi2XpHvGTk55Cil3x7RdbqMxLp6XPx2/PUcduKHbNm4E/f7YPjy/I+y6pWPB6n7y+nk9eXy9rcgsxjWTgYll6hv1U5qXz0Stref5sH99+sQVgzlr05mIXzX1eQuHZqy08duQ8ha5Urrto4WUWlFISvIWYQQK4WJbuYT9lOWncc0UNaSk2/uvgOWoLMmZtmgBGBh6cDNMxOBZ13OMN8OxpN++9pJwUmTkpxJLJ/zViWXpH/JRkp5Gf6eT9ZhfK9fN0gjQUGzcnZ5ZRfnqsi8mw5q69sm6JEMshAVwsmX8ixIAvSFlOGgCfuLaOspw07ri4fM7zp7cSTvfYkQ52V+awZZ7uEyHEwuQmpliy3hGjA6U0x9jVvbYwkwN/Mf+qCjnpKRRnpUZl4Ce7hjnVPcJf37EjsYMVYh2TDFwsmbU2t5WBx6Kh2BWVgT9+pBOn3cbtu+fO2oUQi5MAvo4caR+I9GcnknWN0iUE8M3FLprdXrTWNPaO8sihc9y0vYS8TGeihinEuicBfJ0IhzX3fucQ33yuKeb3/P1TZ7n7/x1Y8BxfYJLfnOmNOmZl4KXZS8vAvYFJzvaO8vHvHibd6eB/v3tbzO8XQswmAXyd6BwaxxuYpGck9gz81VYPxzuGFtxv8l+ea+J//edh2vp9kWO9I36y0hxkLmHLMqsT5X/9xyG6h/z8229fSnlueszvF0LMJgF8nbBuEMa6fRlAc5+PwGSY4fGJOV8PTob54eHzAJzsmtoGtXt4fEn1b5jqROka9vOlO3eytyZ/Se8XQswmAXydaHQbu970jcYWwAd8QQZ8QYB5s/ZnT/fS7zXOOdU9tSNOz7A/0oESq0KXkx3l2fze2+q5+zLp+xYiHqSNcJ2YnoHPt3PNdC19Ux0hPcN+tpbO3gnnv187R3lOGpmpDk5FZeBzn78QpRRP/NG1S3qPEGJhkoGvE1aLnn8ijC8YWvT85mkBvHeODPz8wBgvNfVz92VV7KrI4VS3EcAnQmH6vIEldaAIIRJDAvg6oLWmqddLdprxF6pYyijNfT6cDuM/f8/w7PMfOXQeBdy9r4rt5dn0jgTo9wZwjwbQemk94EKIxJAAnkBdQ+McaR9M+HXcowFGA5Ps32Ts6B7Ljcxmt5e6wkwKMp2zauCToTCPHj7P9VuKKc9NZ3uZUS453T2yrB5wIURiSABPoH95ron7Hz6S8Os09hrlkCvrjQAeWwbupb7YRUl22qwSynNn+3CPBvjQ/moAtps7xZ/skgAuxFoiATyB+kcDDPgChOdYBzuerA6Uq+pjy8ADkyHODYxRX+SiNCdt1uzNV1s8pDps3LDFWKM7N8NJRW46p7pG6B4eB6AsW3q4hVhtEsATaGhsgrCGUf9kQq/T6PaSk57CRSVZ2JTxi2Mh7Z4xwhrqizLnzMBb+31sKszEMW2N7m1l2ZwySyjpKXay06WBSYjVJgE8gQbHjB7qofFgQq/T5PayudiF3abIz3TSNyMDD4U1E6Fw5Hmz2bFSX+SiNDsNjy9IYHKqc8UK4NNtL8+mpc9La7+Pspw02R1HiDVAAngCDZkzHAfH5p7pGC9Nbm9kpmOhK5W+0ehfGPc/fITf/6/XI8+tFsK6okxKc1KBqbr5ZCjMuYGx2QG8LJuwhgMtHql/C7FGLBrAlVLfUUq5lVJvTjuWr5R6RinVaP6Zl9hhJh+tNUNWBj6WuAzc4w0w4AtGAnhRVuqsGvix84P8+nQv7R5jPZPmPh8VuelkOB0UmwtSWWWUjsFxJsN61tZoO8wbmWPBkARwIdaIWDLw/wTeNePY54FntdabgWfN52IaXzDERMi4eTnfWiPxYM3A3FxiLBZlZOBTAXzUPxGZDv+DQ8a6Js19XuqKjABtrSho9YK3motW1c0I4JV56WSZi1ctZRVCIUTiLBrAtdYvAAMzDt8BPGQ+fgi4M77DSn6DvuCcj1fKF5jkmr/9DQ+au8BbMzBnZuDWCoPtHmMj4axUBz883EFwMkyz20t9kXF+JICbGXiLGcBnllCUUmwzs3CZxCPE2rDcGniJ1rrbfNwDlMRpPOvG9Kx7aEYG3tg7ypefOLWs9sLXzw3SMTjOl544zY+PdtDk9pLptFNuBtVCl5PAZJjRgNH5YmXUv3d9Pf3eAP/9aju+YIh6MwPPzUjB6bBFSiht/T6y0xzkz7HRglVGWepCVkKIxFjxTUxtpHrzRiKl1H1KqcNKqcN9fX0rvVzSGJxW9x6acRPzl2/28O0XWznbO7rkzz3UNohNwf7afD77wxM8fbKHhmJXpCuk0GXclLRaCa11vO+9qpaynDS+8WwjQCQDV0pRmj3VC97a72NTkWvOLpOd5TkAVMg63kKsCcsN4L1KqTIA80/3fCdqrR/QWu/TWu8rKipa5uWSz/TOk5k3Ma2bjMuZZn+4bYDt5dk8+Dv7aCh20TXsp94sn4BRQjGuYVyzzTNGaXYarlQHd++rioxr+ntKs9MiJZTWfh+bCjLmvPbte8r5t9/ey7Yy2UVeiLVguQH8Z8C95uN7gZ/GZzjrhxW0y3PSZpVQrAD++hID+EQozNFzQ+yrySc7LYX//N397CjP5votxZFzrAzcupHZ5vFRYwbkuy+rwqbAleqg2Az0ACU5xmQe/0SIzqFxNhW6mIvTYeOdO0qlB1yINWLR6XRKqe8D1wOFSqkO4IvAV4BHlVIfA9qBuxM5yGRklU1qCjJnlVD6zT7tI+eWFsBPdo0wPhHislpjN5vSnLRZa2xPZeBTJZSbthu3KCpy07llVxkj4xNRQbg0O5Wnh/20mW2Gm4qib2AKIdamRQO41vpD87z09jiPZV0ZHAviSnVQlJXKiY6hqNesmZLtnjH6RgORoLuYw21GM9C+2vnb7vMynMZ0em+AEf8EHl8wqqf7Gx+8eNZ7SrLTCEyGOXbOGOfMFkIhxNokMzGX4BPfPRxp31vM0NgEuRkp5GakzC6hjAbYXWncEFxKHfxQ2wDV+RmULNCHbUynN3rB2/uNFsLagqmA7LDbotY4gamVBQ+0eIzzJYALkRQkgMfIG5jk16d7+eHhjpjOHxwLkpfhJDfDyfD4RKRl0D8RYjQwyfUXFeG023g9xjKK1prDbYOR8slCrF7wVs/cPd0zWb3gB1s8FGWl4lrCbvNCiNUjATxGZ7pH0BrO9o7iHp17E+DpIhl4egpaw4jfyMKt2nR5bjq7KnNizsBb+314fEEuW6B8Yil0OenzBiMthNX5c3eVWEoi0+kDiwZ7IcTaIQE8Rienbep7oNmz6PlDY0FyM5zkZqSYz60AbtzALHSlsrcmjzc6hqNWApzP4TYj0O+LJQN3pdI/GqDNXDkw3Wlf8Pzi7KkavNS/hUgeEsBjdLJrmLyMFHLSU3i5qX/R8wfHJsjLSCEvw5jRaNXBrQk2hVmpXFqdRzAU5s3OkXk/x3KobYC8jJTIDMqFFGWl0meWUKbXv+eT6rBHZl5KBi5E8pAAHqNT3SPsKM/hyroCXm7yRNYamUsorBnxT5Cb4STHzMCtmZlWCaXQ5WRvjVEOiaUf/FDbAPtq82PqwS50pRKcDHO6eyTmG5JWGUUCuBDJQwJ4DCZCYd7q8bKjPJurGwroHBqPLBI1l5HxCbSG3HSjBg4wPBZdAy90pVKUlUpNQcaCdfDAZIjvvNRKm2cspvo3QGGWkU37J8JsKly4/m0pNcsoEsCFSB7SbhCDxl4vwVCY7eXZ7Kww2v9ebu6PZLf+CaOGnZZi1JqtbDsvc1oJJZKBB8lKc0TO3Vudx4tN/Wito7JrrTU/OdbJ3z/1Fp1D41xVX8Bde6tiGm+Ra6rNMJYSChithEpB9TzT6IUQa49k4DE41W3UqHeUZ1NXmElpdhqvNBk3Mn2BSe785st88ntTO95Y643kZjjJTk+JOtbnDVDkmrppuK82n77RAC82RtfVv3ugnT955Dh5mSk8/LH9fO/jl8+5QuBcrAwcYu/p/sjlNfzl7TtIdSx8w1MIsXZIAI/Bya5h0lPsbCo0Vum7uqGQV5r7CYU1n3vsBGd6RjnRMRw5f9jcAzMvw4ndpshOc0SWl+0fDUTWKwG485JyNhe7+Myjx3CbC0qd6hrhy0+e5satxfzsU9dw7eaiJa0/Yv2CUGrxFkLLzooc7r2qNuZrCCFWnwTwGJzsGmFrWRZ2mxFEr24oYHBsgs/+8DhPvNFNQ7GLfm8gUuce9JkZuJl952U6p5VQAlEZcobTwbc+cim+QIg/+sFRRv0T/MH3Xyc3PYWv3rUbm23pC0dZvzjKc9IjpRohxPojAXwRWmtOd42wvSw7cuzqhkIAfnS0k3fvKuPz79oKQHO/sTtOpAZu1r9z01MiJZR+bzAqAwdjO7S/uXMnB1sGuOUbL9La7+PrH7yYAldsa6TMZDN3p6+N8QamECI5yU3MRZwfGGc0MMkOczMDMFrutpVlEwqH+bu7dkeWbm12e7m0Oo+hsQlsCrLSjH+9ORlOhsYnCEyGGB6fmBXAAe7aW8nBFg+PHengUzfUc5X5S2K5fvfq2pjLJ0KI5CQBfBEnu4zatrWdmOW/PrYfp8NGZqqDVIcNp91GU5+RgQ+NG7MwrfJHXkYK7R4fnmmzMOfypTt3cvP2Em7cWjzn60vxyesbVvwZQoi1TQL4Ik51j2C3KbaURu9CM7284bDbqC3MoNltrD0yaK6DYslNT2FobCJqEs9c0lLs3LyjNN5fQQixTiV1DfzouUG+/UJsy7su18muEeqLMhe9GVhf5KLFysDHgpEbmGCUUEb8E7hHpqbRCyHESiV1AH/4YDtffvI05wfmnxW5Eie7hjl6bjCq/j2f+iIX7QNjBCfDDPomIjcwwSihaA0t5k3OomXenBRCiOmSOoBb09l/fqIrrp/rHvHzuceOc9s/vwTAh/ZXL/qe+uJMQmHNuQEfw+PGOigWq5zS2GsE8Plq4EIIsRRJXQO31rv++fHuuN20C0yGuO2fX2JwLMjHr9nEH9y4mZxp5ZD51BcZGwE3uX3mZg7Ta+BGMG/q85LptC+6vKsQQsQiaQO4td9jRW46p7tHaHJ7aSieezf1pTjnGcM9GuDv7trN3ftiW3sEoM4M4Ke7RxgLhqJvYpqPm9xeqX8LIeImaUso1n6P97+tDqXgF3Eqo7SZZZmLSrIWOTOaK9VBaXZaZGXB6BKK8XjUPynlEyFE3CRtALf2e9y/KZ/9tfn8/HjXgmt0x6rd/NzaZazKV1+cyVFzj8vpNzGnd6TM10IohBBLlbQBvN2sf9fkZ3L7nnKa+3yc7h6N+f1aa558o5vgZDjqeGu/z9xNfumBtr7IhS9oLC07vQaenZ6CtRaVZOBCiHhJ2gDe6vFRmm3s93jLzlLsNrWkMsoLjf188nuv88s3u6OOt3vGqIlxDe2ZrBuZQGQnHsBckdB4LgFcCBEvSRvA2z1jkcWaClypXN1QyM9PxF5GeeZUD2As3Tpdm8e3rPIJRAfwvBkZvHUjU25iCiHiJWkDeFu/L2r7r+s2F3J+YJwBX3DR92qt+fUpNwCne6bKLoHJEF1D4zHvYjNTffHU+2YHcON5kdTAhRBxkpQB3GohnF7qsILuuRhmZZ7sGqFnxE9WqoMz3VMZ+PmBccKaZS/DWpqdRobTTqrDNqvX27qRKSUUIUS8rCiAK6XepZQ6q5RqUkp9Pl6DWozVQjg9U64xyx4zA7jWGo+5iJTlmVO92BTcc2UN7tFA5HWrA2W5NXClFPVFrlnZN0wroUgAF0LEybIDuFLKDnwTuAXYDnxIKbU9XgNbiNVCOD1TrswzA/iM3eKfP9vHZV/+NS9N23Py16d72VuTx9X1xprbZ8wySqvZ2bJpmQEc4LqLCrm0JnfWcSuoSw1cCBEvK8nA9wNNWusWrXUQ+AFwR3yGtbDpLYSWdKed4qzUWRn48Y4hwhq+8OM3GA8aNe6TXSO8Y1sJW8uMyTqnzTJKu2eM7DRH1CzKpfrsO7fyrY/snXV8/6Z8rt1cSKZMoxdCxMlKptJXAOenPe8ALp95klLqPuA+gOrqxReFisX0FsLpagoyaJ8RwJv7fGQ47ZwbGOPrz75FRW46AO/YXkKhK5WirNRIBt7mMW6MLmUD4VjduquMW3eVxf1zhRAbV8LXQtFaPwA8ALBv376VT5UkuoVwuqr8DA40e6KONbm9XL4pn+KsNB58sZWaggzqCjMjLX9bS7MiGXibx8clVXnxGKIQQiTcSkooncD01Z4qzWMJN7OF0FKTn0nPiB//hDEbMhTWtPQZi1x94dZt5GU4aenz8Y7tJZH3bCvLprHXy3gwROfg+LJ7wIUQ4kJbSQA/BGxWSm1SSjmB3wJ+Fp9hzW+uFkJLdUE6WkPH4DgAXUPjBCbDNBS7yMlI4W/u2IFNEVXK2FqaRTAU5oXGPrOFcPk3MIUQ4kJadglFaz2plPoD4CnADnxHa30ybiObx1wthJbqfKsX3EdDsYsmt7GBglUuuWVXGce+eHNkWjvA1lJjs+JfvWnMzFxuC6EQQlxoK6qBa62fBJ6M01hiMlcLoaU6P7qVcGYAB6KCN0BDsQuHTfHr073G50oJRQiRJJJuJuZcLYSWQpfT7DgxSijNfV4KMp3kZc4/fd3psNFQ7GLUP0lWmoP8Bc4VQoi1JOkC+Budw5TlzG4hBGMmZHV+BucGjCDf5PZSH8MuPVtLjX7w2oLEtBAKIUQiJFUA7xgc49kzbm7fUz7vOVX5GZwbGENrTVNfbNusbS0z6uByA1MIkUySKoA/9EobAPdeVTvvOTVmAO/3Bhkam4iqf89nKgOX+rcQInkkTQD3Bib5wWvnuXVXWWQ25VyqCzLwT4Q52GJM6IklA99dmUtWqoNLq2USjxAieSTNrvSPHjrPaGCSj12zacHzrE6U584Y633HEsDzM50c/+LN2GxS/xZCJI+kyMBDYc13Xm5lX00eF1flLniuFcCff6uP9BQ7ZdlpMV1DgrcQItkkRQB/+mQPHYPjfPzahbNvMJaVVQoGfEHqizMlMAsh1q2kCOD/8UobVfnp3LS9dNFznQ4b5TlGjbwhhhuYQgiRrJKiBv6Pd++ha8iPPcZsujo/g86h8Zg6UIQQIlklRQZemZfB/k35MZ9v1cFjuYEphBDJKikC+FJVm/3csczCFEKIZJUUJZSles+ecsaDIamBCyHWtXUZwKvyM/izd25Z7WEIIURCrcsSihBCbAQSwIUQIklJABdCiCQlAVwIIZKUBHAhhEhSEsCFECJJSQAXQogkJQFcCCGSlNJaX7iLKdUHtC/z7YVAfxyHkwzkO28M8p3Xv5V+3xqtddHMgxc0gK+EUuqw1nrfao/jQpLvvDHId17/EvV9pYQihBBJSgK4EEIkqWQK4A+s9gBWgXznjUG+8/qXkO+bNDVwIYQQ0ZIpAxdCCDGNBHAhhEhSSRHAlVLvUkqdVUo1KaU+v9rjSQSl1HeUUm6l1JvTjuUrpZ5RSjWaf+at5hjjSSlVpZR6Til1Sil1Uin1x+bx9fyd05RSrymljpvf+a/M45uUUq+aP9+PKKWcqz3WeFNK2ZVSR5VSvzCfr+vvrJRqU0q9oZQ6ppQ6bB6L+8/2mg/gSik78E3gFmA78CGl1PbVHVVC/CfwrhnHPg88q7XeDDxrPl8vJoE/1VpvB64APmX+d13P3zkA3Ki13gNcDLxLKXUF8LfA17TWDcAg8LHVG2LC/DFwetrzjfCdb9BaXzyt/zvuP9trPoAD+4EmrXWL1joI/AC4Y5XHFHda6xeAgRmH7wAeMh8/BNx5IceUSFrrbq316+bjUYz/uStY399Za6295tMU8x8N3Ag8Zh5fV98ZQClVCbwbeNB8rljn33kecf/ZToYAXgGcn/a8wzy2EZRorbvNxz1AyWoOJlGUUrXAJcCrrPPvbJYSjgFu4BmgGRjSWk+ap6zHn++vA58DwubzAtb/d9bA00qpI0qp+8xjcf/ZXpebGq9HWmutlFp3PZ9KKRfwOPBprfWIkZwZ1uN31lqHgIuVUrnAj4GtqzuixFJK3Qa4tdZHlFLXr/JwLqRrtNadSqli4Bml1JnpL8brZzsZMvBOoGra80rz2EbQq5QqAzD/dK/yeOJKKZWCEby/p7X+kXl4XX9ni9Z6CHgOuBLIVUpZydR6+/m+GniPUqoNo/x5I/AN1vd3Rmvdaf7pxvhFvZ8E/GwnQwA/BGw271o7gd8CfrbKY7pQfgbcaz6+F/jpKo4lrsw66L8Dp7XW/zjtpfX8nYvMzBulVDpwE0bt/zngLvO0dfWdtdZ/obWu1FrXYvy/+xut9UdYx99ZKZWplMqyHgM3A2+SgJ/tpJiJqZS6FaOOZge+o7X+8uqOKP6UUt8HrsdYdrIX+CLwE+BRoBpjGd67tdYzb3QmJaXUNcCLwBtM1Ua/gFEHX6/feTfGzSs7RvL0qNb6r5VSdRjZaT5wFLhHax1YvZEmhllC+TOt9W3r+Tub3+3H5lMH8N9a6y8rpQqI8892UgRwIYQQsyVDCUUIIcQcJIALIUSSkgAuhBBJSgK4EEIkKQngQgiRpCSACyFEkpIALoQQSer/BzERFafMf8sJAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.plot(xdata, ydata);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will write the PyMC3 model, keeping in mind that the data must be modifiable (both `x` and `y`)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "with pm.Model() as linreg_model:\n", " # optional: add coords to \"time\" dimension\n", " linreg_model.add_coord(\"time\", np.arange(len(xdata)), mutable=True)\n", "\n", " x = pm.MutableData(\"x\", xdata, dims=\"time\")\n", " y_obs = pm.MutableData(\"y_obs\", ydata, dims=\"time\")\n", "\n", " b0 = pm.Normal(\"b0\", 0, 10)\n", " b1 = pm.Normal(\"b1\", 0, 10)\n", " sigma_e = pm.HalfNormal(\"sigma_e\", 10)\n", "\n", " pm.Normal(\"y\", b0 + b1 * x, sigma_e, observed=y_obs, dims=\"time\")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [b0, b1, sigma_e]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [6000/6000 00:04<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 5 seconds.\n" ] } ], "source": [ "sample_kwargs = {\"chains\": 4, \"draws\": 500}\n", "with linreg_model:\n", " idata = pm.sample(**sample_kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have defined a dictionary `sample_kwargs` that will be passed to the `SamplingWrapper` in order to make sure that all refits use the same sampler parameters. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will create a subclass of `az.SamplingWrapper`. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from scipy import stats\n", "from xarray_einstats.stats import XrContinuousRV\n", "\n", "\n", "class PyMCLinRegWrapper(az.PyMCSamplingWrapper):\n", " def sample(self, modified_observed_data):\n", " with self.model:\n", " # if the model had coords the dim needs to be updated before\n", " # modifying the data in the model with set_data\n", " # otherwise, we don't need to overwrite the sample method\n", " n__i = len(modified_observed_data[\"x\"])\n", " self.model.set_dim(\"time\", n__i, coord_values=np.arange(n__i))\n", "\n", " pm.set_data(modified_observed_data)\n", " idata = pm.sample(\n", " **self.sample_kwargs,\n", " )\n", " return idata\n", "\n", " def log_likelihood__i(self, excluded_observed_data, idata__i):\n", " post = idata__i.posterior\n", " dist = XrContinuousRV(\n", " stats.norm,\n", " post[\"b0\"] + post[\"b1\"] * excluded_observed_data[\"x\"],\n", " post[\"sigma_e\"],\n", " )\n", " return dist.logpdf(excluded_observed_data[\"y_obs\"])\n", "\n", " def sel_observations(self, idx):\n", " xdata = self.idata_orig[\"constant_data\"][\"x\"]\n", " ydata = self.idata_orig[\"observed_data\"][\"y\"]\n", " mask = np.isin(np.arange(len(xdata)), idx)\n", " data_dict = {\"x\": xdata, \"y_obs\": ydata}\n", " data__i = {key: value.values[~mask] for key, value in data_dict.items()}\n", " data_ex = {key: value.isel(time=idx) for key, value in data_dict.items()}\n", " return data__i, data_ex" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Computed from 2000 posterior samples and 100 observations log-likelihood matrix.\n", "\n", " Estimate SE\n", "elpd_loo -255.84 6.31\n", "p_loo 2.70 -\n", "------\n", "\n", "Pareto k diagnostic values:\n", " Count Pct.\n", "(-Inf, 0.5] (good) 100 100.0%\n", " (0.5, 0.7] (ok) 0 0.0%\n", " (0.7, 1] (bad) 0 0.0%\n", " (1, Inf) (very bad) 0 0.0%" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "loo_orig = az.loo(idata, pointwise=True)\n", "loo_orig" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, the Leave-One-Out Cross Validation (LOO-CV) approximation using [Pareto Smoothed Importance Sampling](https://arxiv.org/abs/1507.02646) (PSIS) works for all observations, so we will use modify `loo_orig` in order to make `az.reloo` believe that PSIS failed for some observations. This will also serve as a validation of our wrapper, as the PSIS LOO-CV already returned the correct value." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "loo_orig.pareto_k[[13, 42, 56, 73]] = np.array([0.8, 1.2, 2.6, 0.9])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We initialize our sampling wrapper" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "pymc_wrapper = PyMCLinRegWrapper(model=linreg_model, idata_orig=idata, sample_kwargs=sample_kwargs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And eventually, we can use this wrapper to call `az.reloo`, and compare the results with the PSIS LOO-CV results." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/oriol/Public/arviz/arviz/stats/stats_refitting.py:99: UserWarning: reloo is an experimental and untested feature\n", " warnings.warn(\"reloo is an experimental and untested feature\", UserWarning)\n", "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [b0, b1, sigma_e]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [6000/6000 00:04<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 5 seconds.\n", "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [b0, b1, sigma_e]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [6000/6000 00:04<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 4 seconds.\n", "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [b0, b1, sigma_e]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [6000/6000 00:03<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 4 seconds.\n", "Auto-assigning NUTS sampler...\n", "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [b0, b1, sigma_e]\n" ] }, { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", " \n", " 100.00% [6000/6000 00:04<00:00 Sampling 4 chains, 0 divergences]\n", "
\n", " " ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 1_000 tune and 500 draw iterations (4_000 + 2_000 draws total) took 4 seconds.\n" ] } ], "source": [ "loo_relooed = az.reloo(pymc_wrapper, loo_orig=loo_orig)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Computed from 2000 posterior samples and 100 observations log-likelihood matrix.\n", "\n", " Estimate SE\n", "elpd_loo -255.82 6.30\n", "p_loo 2.69 -\n", "------\n", "\n", "Pareto k diagnostic values:\n", " Count Pct.\n", "(-Inf, 0.5] (good) 100 100.0%\n", " (0.5, 0.7] (ok) 0 0.0%\n", " (0.7, 1] (bad) 0 0.0%\n", " (1, Inf) (very bad) 0 0.0%" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "loo_relooed" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Computed from 2000 posterior samples and 100 observations log-likelihood matrix.\n", "\n", " Estimate SE\n", "elpd_loo -255.84 6.31\n", "p_loo 2.70 -\n", "------\n", "\n", "Pareto k diagnostic values:\n", " Count Pct.\n", "(-Inf, 0.5] (good) 96 96.0%\n", " (0.5, 0.7] (ok) 0 0.0%\n", " (0.7, 1] (bad) 2 2.0%\n", " (1, Inf) (very bad) 2 2.0%" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "loo_orig" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.10" } }, "nbformat": 4, "nbformat_minor": 2 }