07 Chapter 7/01 Wines and Judges.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import pandas as pd\n",
    "import pymc3 as pm\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import numpy as np\n",
    "from pandas.plotting import scatter_matrix\n",
    "\n",
    "plt.style.use('seaborn-whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize(alpha=0.05):\n",
    "\n",
    "    def _summarize(ser):\n",
    "        hpd_ = pm.hpd(ser, alpha=alpha)\n",
    "        return pd.Series({\n",
    "            'mean_': ser.mean(),\n",
    "            'lower': hpd_[0],\n",
    "            'upper': hpd_[1],\n",
    "        })\n",
    "\n",
    "    return _summarize\n",
    "\n",
    "\n",
    "def agg_prediction(trace_var, source_data, group_columns, aggregate_fc):\n",
    "    df = pd.DataFrame(trace_var).stack().reset_index()\n",
    "    merged = pd.merge(\n",
    "        df, source_data,\n",
    "        left_on='level_1', right_index=True\n",
    "    )\n",
    "\n",
    "    return (\n",
    "        merged.groupby(group_columns)[0].apply(aggregate_fc)\n",
    "        .unstack().sort_index().reset_index()\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "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>judge</th>\n",
       "      <th>flight</th>\n",
       "      <th>wine</th>\n",
       "      <th>score</th>\n",
       "      <th>wine_amer</th>\n",
       "      <th>judge_amer</th>\n",
       "      <th>flight_red</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>76</th>\n",
       "      <td>Jamal Rayyis</td>\n",
       "      <td>white</td>\n",
       "      <td>G1</td>\n",
       "      <td>12.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>155</th>\n",
       "      <td>Daniele Meulder</td>\n",
       "      <td>red</td>\n",
       "      <td>F2</td>\n",
       "      <td>15.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>158</th>\n",
       "      <td>Daniele Meulder</td>\n",
       "      <td>red</td>\n",
       "      <td>I2</td>\n",
       "      <td>10.0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>Robert Hodgson</td>\n",
       "      <td>white</td>\n",
       "      <td>B1</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>141</th>\n",
       "      <td>Linda Murphy</td>\n",
       "      <td>red</td>\n",
       "      <td>B2</td>\n",
       "      <td>14.0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               judge flight wine  score  wine_amer  judge_amer  flight_red\n",
       "76      Jamal Rayyis  white   G1   12.0          1           0           0\n",
       "155  Daniele Meulder    red   F2   15.0          1           0           1\n",
       "158  Daniele Meulder    red   I2   10.0          1           0           1\n",
       "41    Robert Hodgson  white   B1   11.0          1           1           0\n",
       "141     Linda Murphy    red   B2   14.0          0           1           1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wines = pd.read_csv('../data/Wines2012.csv', sep=';',\n",
    "                    names='judge,flight,wine,score,wine_amer,judge_amer'.split(','),\n",
    "                    skiprows=1)\n",
    "wines = wines.assign(flight_red = (wines.flight == 'red').astype('int'))\n",
    "\n",
    "wines.sample(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: 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",
      "  rval = inputs[0].__getitem__(inputs[1:])\n",
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: 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",
      "  rval = inputs[0].__getitem__(inputs[1:])\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, beta]\n",
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: 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",
      "  rval = inputs[0].__getitem__(inputs[1:])\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:28<00:00, 418.42draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_01:\n",
    "    beta = pm.Normal('beta', mu=0, sd=1, shape=(2, 2, 2))\n",
    "\n",
    "    score_i = pm.Deterministic(\n",
    "        'score',\n",
    "        wines.score.mean()\n",
    "        + beta[tuple([wines.wine_amer.values,\n",
    "                      wines.judge_amer.values,\n",
    "                      wines.flight_red.values])]\n",
    "    )\n",
    "\n",
    "    sigma = pm.HalfCauchy('sigma', beta=5)\n",
    "    score_obs = pm.Normal('score_obs', mu=score_i, sd=sigma,\n",
    "                          observed=wines.score.values)\n",
    "    \n",
    "    m_01.name = 'full interaction'\n",
    "    m_01.trace = pm.sample(5000, tune=1000, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEcCAYAAADUX4MJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X14lPWd7/F3QKpgEC1NZXcL5kD1a/a4ra3G0l0r0AI+VEtL3erhUgHFau3WIq2WWqvV4+lliwX0WEXlIVTKYltrkSoWtDGLoAfWh9buxi9IDCAoBMpTmISi5Pwxd8YxzEyGkMz9G/i8rsuLyfzuue/P3Cbzne/vNw8lLS0tiIiItKdb3AFERKQ4qGCIiEheVDBERCQvKhgiIpIXFQwREcmLCoaIiOTlqLgDiBwKMzsamAMMBHYB33T3NWb2FeBuYEO06W3AS8ATQE/gGnf/s5mdDfyLu/8ky/6/DlwG7Ad6AD9w9+c6kPNaoB8wA7jV3a8zs3rgVHdvTttuXHTd5Dz2WQUscPens4wPAD7p7osONm87xx0KXOvul3bmfiV86jCk2F0NNLr7YOBbwH3R9WcAN7n70Oi/GmAkyYJxHXCVmZUA3wbuybRjM7sUGAF8wd2Hkiwcj5jZRzoa1t3fcffrOnr7g/R54F8KdCw5AqjDkGL3j8BiAHd3M6uIrj8D+JSZTQRWAt8DGkl2Fz2BPcAY4PH0Z/htXANMcvd90f7fNLPT3X1b9Oy+b/TfF4GbgM8B3YGp7v7rqHu5B9gOvAu8aGblJLuCwdExHoyu2wyMTT+4mX0rytgS3ebeTCGjZ/zfA/5GstNaANwFTAZ6mdkK4E3gXqAE2AZcCXwK+El0u2eAoe4+LNrn74EfAoOAb5LsrlqAr7Q59hzg49E5vcfdH8lyLuUwoA5Dit2rwIVmVmJmg4F/MLPuwFKSHcc5QClwLckHxRNJFoKHSD74/cnMHjSzmzLs+++BuvQr3H1b2o9/dPd/BgYD/8PdzwaGAT8ws+OBB4D/5e7DST5gZ/KAuw8B6kl2SwCY2T8ClwBnkyxEXzYzy3EeTgK+GmW5yd3fI1k05rv7E8DDJKfrhgJPkSxwAMe4++fc/XbgGDM7ycz+DviIu78CnAJ8Mbpv/w2cm5axN8nzOxo4D3gvRz45DKjDkGI3G6gAlgHLgZfc/T0zm+3uOwDMbCHwVXffD0yMrvs+yWf/t5AsLD8ys1PcfXXavtcB/YGdrVeY2bnAn6MfPfr3n4AzzOy56OceQDlwYtr+lpN8Jp7ub+7+YnR5Bcnpr1XRz6eRLALPRj+fAJycdsy2XnP3d4F3zawpw3gFcH9Uc3oAa9rcB4BZwBXAXpLrQgBbgLlm1gicCrzQurG77446uIeA44B5WbLJYUIdhhS7SuDZ6Bnwr4G6aG3iz2b2sWibL5Bc8AbAzD4KmLsvA3qRfGbcAhzbZt+zgR+a2VHR7U4BZvL+M+n90b+vA9XRs/fPA78C1gIb06bIKjNk/5CZnR5d/hzwl7QxB/4LGBbtt4r3C1UmmT4Ubj/v/407cEW0r5uA37e5D5CcyrqQZOc138z6ALcDlwITgCaSU1oARJ3IGe7+FZLTcj9tPVdyeFLBkGK3BphoZi8A/5vkmkMLyQe435pZDcmi8HDabW4B7owu3w/8Afg74E/pO3b3BcCLwPNm9h8kn3Vf5u5b2mRYBDSa2TKShanF3XeTnPr6hZk9S7JbaGsv8K1o3x8lWYxaj/0nkt3F82b2nyS7i435nxYAXgNGRYv334iyPE9yquqA4uPujdE5qI3y7yLZGb1AsoNrIjlN1+odoF+0RrIUuDvqcuQwVaJPqxURkXyowxARkbyoYIiISF5UMEREJC8qGCIikhcVDBERycsR/Zrpl156SS8RExFp44wzzijJdP2R/rLa4O98bW0tFRUV7W8YAGXtXA8//DBvv/02t956a9xR8lIM57SVsuaUsVjAEd5hiITs6quvpra2Nu4YIilawxARkbyowxAJ1NChQ0kkEqxcuTLuKCKACoZIsMaNG8emTZvijiGSooIhEqhx48ZpDUOCojUMkUDt27ePffv2xR1DJEUdhkigRowYoTUMCYoKhkigJkyYwMaNB/sVGCJdRwVDgrV+W4Kr5q6irmEPA8uOZdbYSgb07RV3rIK57LLLtIYhQVHBkLyMn7OSam/IY8u6Ljn+mi2NnDOlupP3ml/WYVbGnPFndfKx25dIJGhqyvT13CLxyLtgmFkVsMDdn85j2wHAJ919UR7bjgPuAKYDvwDmAz2BTcB4d0+kbduN5FdqfpLk11tOcPc32uzvapJfjfkuya/hfI7k9xef6u792stzMEZOq2H15sbO3GUWXfMgLPmp9gbKJz8ZX4BfxXjsyCknlrLkhiFxx5CYdVWH8XngVJLfdZyP+e4+1czujS5Xmdlkkg/809K2+zJwjLt/1swGAz8DRrUOmlk/4HrgTOAY4HngTHcfambvHPK9aqMQf0BH8mfejJhaw9qGRva3QEkJHNWthP376ZTpqWI4r48++igbN25k0qRJcUcRAQ6+YFxnZjdGt7vK3d8ws28BY0h+kN8C4OfAZKBX9OXwO4HbSL6EtxQY4+6rs+z/bODH0eXF0eVpbcafBnD3F83szDa3PwtY7u57gb1m9gbwCWDVQd5PySC/aamu6YZaWmDfe8nPiuy86amOZS3UFNUll1yiNQwJysEWjBXufpeZXQD81MxuAS4h+UAOsBT4A3AXySmgJ8zsOuAyd99kZjcD/wr8nyz7P45kgQHYDfTJMQ7wnpkd5e7vZhnPtI+U+vr6Ds8RX7twA+t2FOo18pqSCknhp6ji//9/0vE9mDGqf85tmpubi6bAKWt2uTrvgy0Y/xH9uwKYApwGnAQ8G11/AnBym9tsBO41s0bgH4DlOfa/C+gNNEX/7sgy3qpbWrHINJ5pHynl5eU5ouRWU6DpjGKYOmnVlVnTp6e6lcCgslKWTur4lGAxnNdi+yypYjinrZS1Yw72nd6tffjngL8ADvwXMMzdhwJVwJ+B/Wn7fpjk4vU4kgvZWT9rnWQxuSC6fD6wLNt4tIbxWpvxlcDnzOwYM+sDVEQ5pcjNGlvJoLJSupeUMKislFljK+OO1OWuv/56LrvssrhjiKQcbIcx2Mz+SHK94kp3X2dmzwLPm9nRJB+wN5J8IP+Bmb0MzAOWmdkeYDPw9zn2fycwN3ql01aSayOY2QJgIvA4MCJaGykBxkfjk4FX3f3paOF8GcmC9QN3bz7I+ygBGtC31yF1FMVo9OjRRTNtIkeGvAtG1CFkun4KyempdK8AFl1ecBDH2Aycl2GoDtjt7vuBazOM/zfwt2gfD5PsakSK2tatW9m+fXvcMURSQnnj3hgz2+LuU7OMz3D3PTlu/6q7r880YGalJN+HIVJULr744qJaw5DDX+wFw92rSK595NomYzHIZ9zdG4GhHYgmEqvvfOc7bNiwIe4YIimxFwwRyeyiiy7SGoYERd+HIRKod955h4aGfD6/S6Qw1GGIBOrSSy/VGoYERQVDJFCTJ09m/fqcy3ciBaWCIRKo8847T2sYEhStYYgEasOGDbz99ttxxxBJUYchEqjLL79caxgSFBUMkUDdcsstrFu3Lu4YIikqGCKBGj58uNYwJChawxAJVF1dnd7pLUFRhyESqCuvvJJEIsHIkSPjjiICqGCIBOv222/XGoYERQVDJFBDhgzRGoYERWsYIoFyd9588824Y4ikqMMQCdQ111xDIpHgggsuaH9jkQJQwRAJ1I9//GPq6+vjjiGSooIhEqh//ud/5oQTTog7hkiK1jBEAvWXv/yFNWvWxB1DJEUdhkig/u3f/o1EIsGXvvSluKMcYP22BFfNXUVdwx4Glh3LrLGVcUeSAlDBEAnUlClT2n2V1Pg5K6n2eL+Vb82WRs6ZUh39VFfQYw+zMuaMP6ugxzySlbS0tMSdIU7B3/na2loqKirijpGXkLKOnFbD6s2NcccQ6RSnnFjKkhuGFOpwJdkGgugwzGwccAdwHzAc6AlsAsa7eyJtu27A/cAngb3ABHd/I8P+yoDlwCfcvdnMfh/t93h3b+7iuyMBaO+PK6Tils2rr75KXV0do0ePjjvKAUZMrWFtQyP7055yDejTg3lfP5sBfXvFFyxPxfD/v1VIWYMoGJH5wMeA+e5eZWaTgWuAaWnbfBk4xt0/a2aDgZ8Bo9J3YmbnAncB/Vqvc/cLzay+a+OLdK6JEyeSSCTaLRghTEsBrN+5L21qqjA0JVVYIRUMgLOBH0eXF0eXp7UZfxrA3V80szMz7GM/yW7ipfYOVl9fT1NT0yEF7mrNzc1F8/EQIWS9duEG1u3Yl+fWhZ1vP2iDbwSgfPKTMQcJV7U3HML5Cfz/f5qTjt/AjFH9C3KsXN1MaAXjOGBndHk30CfHOMB7ZnaUu7/beoW7LwUws3YPVl5efihZCyKkdrQ9IWStyfP4IWTNR+g506emSoCPf7SUpZMKNtfeYaGf13QhZQ3tfRi7gN7R5d7AjhzjAN3Si4XI4WTVqlW89tprccfIadbYSgaVldK9pIT+fXro5bWHudA6jOXABUAVcD6wLMP4RcCvojWMsP+aRA7BjTfeSCKR4Gtf+1rcUbIa0LdXqqOora0tigVv6bjQCsadwFwzuxrYCowBMLMFwETgcWCEma0g2QGPj8YnA6+6+9OxpBbpAvfddx91dcUzzy6Hv6AKhrtvBs7LMFQH7Hb3/cC1Gcb/G/hbm32Vd3pAkQI67bTT6N69e9wxRFJCWsMYY2aTsozNcPc9OW77qrv/Mdtg9D6MftnGRUK0YsUKXnnllbhjiKQE0WG4exXJdYts4+vbuX174xd2KJhIjG6++WYSiQRjxoyJO4oIEEjBEJEDPfjgg6xduzbuGCIpKhgigTIz9u/fH3cMkZSQ1jBEJE1NTQ2rVq2KO4ZIijoMkUDddtttJBIJrrjiirijiAAqGCLBmj17Nm+8ccCHMYvERgVDJFADBw5k7969cccQSdEahkignnnmGVasWBF3DJEUFQyRQN155508+OCDcccQSdGUlEigHnnkEdasWRN3DJEUFQyRQPXv35/GRn0vuYRDU1IigXr66adZtqztJ/yLxEcFQyRQd911FzNnzow7hkiKpqREArVgwQJWr14ddwyRFBUMkUD169eP7du3xx1DJEVTUiKBWrRoEdXV1XHHEElRwRAJ1M9+9jOqqqrijiGSoikpkUD95je/0RqGBEUdhkigPvKRj3DCCSfEHUMkRR2GSKB++9vf8tZbb1FRURHL8ddvS3DV3FXUNexhYNmxzBpbyYC+vWLJImFQwRAJ1L333ksikeD6669PXTd+zkqqvaHgWdZsaeScKfkswNd12jGHWRlzxp/VafuTQ1fS0tISdwbMbBxwB3AfMBzoCWwCxrt7IsP2nwF+4u5DM4xdBNwKvAvMdveHzez30X6Pd/fmtM3jv/PtqK2tje0Z5sEKMevIaTWs3qyP15CwnHJiKUtuGJLXtjH8XZVkGwipw5gPfAyY7+5VZjYZuAaYlr6Rmd0EXA7sabsDM+sRbV8ZjS83syfc/UIzq+/a+BKibH+UIRa3TOLMOWJqDWsbGtmf9rTq5I+WZp2aKpZzCsWVNSQhFQyAs4EfR5cXR5entdlmLTAaeCTD7SuAN9x9O4CZPQ+cA/y6S9KKdKFHH32UjRs3ZnxgC3dqSlNSh7PQCsZxwM7o8m6gT9sN3P0xMyvP4/ZZ99Gqvr6epqamjiUtkObmZmpra+OOkZeQsl67cAPrduxrZ6vOe3DrGqWAce/kJ+MOEotqb6C8S+97vP//Tzq+BzNG9W93u0L/XeXqvEIrGLuA3kBT9O+ODt6+Vc59lJeXH+TuC6+YWueQsta0kyOkrNkkEglef/11Pv3pT8eaI31qqlsJDCorZemkA6f6iuGctlLWjgntfRjLgQuiy+cDB/vZzrXAyWb2YTP7EMnpqBc6MZ9IwfTq1YuePXvGHYNZYysZVFZK95ISBpUl1zDkyBRah3EnMNfMrga2AmMAzGwBMNHd38l0o+hVVkSL5ZOAP5AshrPdfWMhgot0tnnz5mVdwyikAX17Zewo5MgTVMFw983AeRmG6kiuR7RuVw8MTht/GTgzGlsELOq6lCKFMXPmTBKJBN/73vfijiIChFUwxpjZFnefmmFshrsf8DLaNH8F5mQbjN6H0e9QA4oU0tKlS4N5EYEIBFIw3L0KqMoxvr6d27/VzviFHQomEqMePXrQo0ePuGOIpARRMETkQFVVVWzatCn2NQyRVqG9SkpEIlVVVfzud7+LO4ZIijoMkUA999xzWsOQoKjDEBGRvKhgiATq4Ycf5te/1segSThUMEQC9eijj7J48eK4Y4ikaA1DJFDPPPOM1jAkKOowREQkLyoYIoG6//77+fd///e4Y4ikaEpKJFCLFi2isVFfLyvhUMEQCdTixYu1hiFB0ZSUiIjkRQVDJFD33HMPjzyS6avrReKhKSmRQD377LNaw5CgqGCIBOqJJ57QGoYERVNSIiKSFxUMkUDdfffdzJ49O+4YIimakhIJ1AsvvMDu3bvb31CkQFQwRAL12GOPaQ1DgqKCISIFtX5bgqvmrqKuYQ8Dy45l1thKBvTtFXcsyYMKhkig7rrrLrZs2cLUqVO79Djj56yk2hs6aW91B7X1mi2NnDOlupOOndswK2PO+LMKcqzDVUlLS0vcGTCzccAdwH3AcKAnsAkY7+6JDNt/BviJuw/Nsr9ewFLgKnd/3cxmAhcDg9399bRN47/z7aitraWioiLuGHk53LOOnFbD6s16X4QUh1NOLGXJDUM6ctOSbAMhdRjzgY8B8929yswmA9cA09I3MrObgMuBPZl2YmZnAjOifQHg7hPM7ONdFVyODB384zskh2MRHjqlmvpt7z8PLO/bi+duHNaV0Q5wOJ7XQgipYACcDfw4urw4ujytzTZrgdFAts9MOBr4So5xEUkT55QUQP22BOWTn+yk4+emaalDE1rBOA7YGV3eDfRpu4G7P2Zm5dl24O7LAcys3YPV19fT1NTUoaCF0tzcXDSvlDncsl67cAPrduwrUKJcDv5BWDKr9oa04lRM5zX/rCcd34MZo/p3+Ei5upnQCsYuoDfQFP27oysPVl5e3pW77xQhtaPtOdyy1sR8Xy677DJ27tzJokWLYs2Rr3z//4+YWsPahkb2t0C3EhhUVsrSSYWd7jvcflcLJbR3ei8HLogunw8sizGLSKzmzZvHT3/607hjdLpZYysZVFZK95ISBpWVMmtsZdyRJE+hdRh3AnPN7GpgKzAGwMwWABPd/Z1MN4peZYW7VxUmpoh01IC+vQreUUjnCKpguPtm4LwMQ3Uk1zRat6sHBqeNvwyc2WZfQzs/oUjh3HrrrWzdupX7778/7igiQFhTUmPMbFKWsRnunvFltJG/AnOyDUbvwzj9UMKJFNqGDRt4++23444hkhJEhxFNJVXlGF/fzu3famd8AjChI9lE4jJnzpyiedWZHBlC6jBERCRgKhgigfr+97/f5Z8jJXIwVDBEArVt2zZ27tzZ/oYiBRLEGoaIHOihhx7SGoYERR2GiIjkRR2GSKC++93vsm3bNubMyfqKcZGCUochEqimpib27t0bdwyRFHUYIoH6+c9/rjUMCYo6DBERyYs6DJFATZw4kb/+9a/84he/iDuKCKAOQ0RE8qQOQyRQ06dP1xqGBEUdhoiI5EUdhkigvvnNb7J9+3bmz58fdxQRQB2GSLB69uzJ0UcfHXcMkRR1GCKBuvvuu7WGIUFRhyEiInlRhyESqK9//evs2LGDX/3qV3FHEQFUMESC1bdvX1paWg64fv22BFfNXUVdwx4Glh3LrLGVDOjbK4aEcqRRwRDpROPnrKTaGzppb2dDXyif/GTWLdZsaeScKdWddLzchlkZc8afVZBjSZiCKBhmNg64A5ju7lPNbCLQz90nZ9j2NuCLwLvARHdf2Wb8IuDWaHy2uz9sZr8HhgPHu3tz194bORgjp9WwenNjTEevi+m4xanaG3IWr6RiOqfvZz3lxFKW3DAkxizFIYiCEZkPPGBmvwTOAh5ru4GZfRoYAnwG6B9tU5k23gOYFl23B1huZk+4+4VmVt/Vd0AOXlx/pLW1tVRUVMRy7HyNHz+eHTt28Pjjj3/g+qFTqqnflkj9XN63F8/dOKzQ8Q5QDOe0VTFlDUlIBQPgGGAusBQ4NcP42cASd28B1pvZUWZW5p6aA6gA3nD37QBm9jxwDvDrro8uR7LOnYqKnHgxnJh7Sgqgflsij2f+h0bTUQKBFYzogX5JNEWVyXHAtrSfdwN9gIa08Z0ZxjOqr6+nqampw3kLobm5uWhei3+oWa9duIF1O/Z1YqL2FNP0Sbzym46CYjmn/Y87iofiDpGnQj8G5Oq8gioYedgF9E77uTew4yDGP6C8vLwzs3WJYmqdDzVrTQHvZ7Gc10w5R0ytYW1DI/tboFsJDCorZemk+Offi+WcgrJ2VLG9cW85cK6ZdTOzAUA3d9+aNl4LnGxmHzazD5GcjnohjqAih+qyyy7jpptuOuD6WWMrGVRWSveSEgaVlTJrbGWGW4t0vqLoMMxsOlDl7i+Z2TKSRaAb8M1ofByAu1eZ2STgD9H4bHffGE9qkUNjZjQ0HLguMqBvryA6CjnyBFkw3L2qzVVrgcZo7EfAj9qMvwycGY0vAhZ1aUCRAvjhD39YNOtXcmQIaUpqTNQdZLLQ3d/Icdu/AnOyDUbvw+h3KOFERI50QXQYUUdRlWN8fTu3f6ud8Qs7FEwkRpdeeim7du3iqaeeijuKCBBIwRCRA51++uls2bIl7hgiKSoYIoGaPHmy1jAkKCGtYYiISMDUYYgE6qtf/Sq7d+9myZIlcUcRAVQwRIL12c9+ls2bN8cdQyRFBUMkUN/97ne1hiFB0RqGiIjkRR2GSKC+9KUv0djYyB//+Me4o4gAKhgiwfrCF76gNQwJigqGSKC+/e1vaw1DgqI1DBERyYs6DJFAnX/++TQ2NrJs2bK4o4gAKhgiwbrooot455134o4hkqKCIRKo6667TmsYEhStYYiISF7UYYgEavjw4ezZs4cXXtDX0ksYVDBEAnXJJZfw9ttvxx1DJEUFQyRQV199tdYwJChawxARkbyow5CCWr8twVVzV1HXsIeBZccya2wlA/r2ijtWkIYOHUoikWDlypVxRxEBVDAEGD9nJdXe0El7q8t7yzVbGjlnSnUnHTe3YVbGnPFnFeRYnWXcuHFs2rQp7hgiKSUtLS1xZ8DMxgF3ANPdfaqZTQT6ufvkDNveBnwReBeY6O4HPP0ys17AUuAqd3/dzGYCFwOD3f31tE0P6c6PnFbD6s2Nh7ILkQ455cRSltwwJO4YH1BbW0tFRUXcMfKirDmVZBsIqcOYDzxgZr8EzgIea7uBmX0aGAJ8BugfbVPZZpszgRnAx1qvc/cJZvbxzg5ciD/Yw+0Xe+iUauq3JVI/l/ftxXM3DuvqaAcohvO6b98+amtr+cQnPhF3FBEgrIIBcAwwl2R3cGqG8bOBJe7eAqw3s6PMrMz9A/MpRwNfAR7p8rSHkc6blsp/SgqgfluC8slPdsJxcyvGKakRI0ZoDUOCElTBcPftwJJoiiqT44BtaT/vBvoAqUc6d18OYGbtHq++vp6mpqaOxgXg2oUbWLdj3yHto30H9yAsB6r2hgyFKfDzOvhGgJwF9aTjezBjVP9CJcqpubm5aF4GrKzZ5eq8gyoYedgF9E77uTewo6M7Ky8vP9Q81HTxtEYxTJ20yifriKk1rG1oZH8LdCuBQWWlLJ1U+Ln4YjmvxZITlLWrhJS12N6HsRw418y6mdkAoJu7b407lORv1thKBpWV0r2khEFlpcwaW9n+jY5QiUTikDtgkc5UFB2GmU0Hqtz9JTNbBrxAsth9MxofB+DuVXFllPwM6Nsrlo6iGF1wwQVaw5CgBFkwMjzwrwUao7EfAT9qM/4ycGabfQztknAiBfKNb3yDjRs3xh1DJCWkKakxZjYpy9hCd38jx23/CszJNhi9D+P0QwknUmiXXHIJ559/ftwxRFKC6DCijqIqx/j6dm7/VjvjE4AJHckmEpedO3eye/fuuGOIpARRMETkQKNGjdIahgRFBUMkUNdffz1vvZWzeRYpKBUMkUCNHj26aN5cJkeGkBa9RSTN1q1b2b59e9wxRFLUYYgE6uKLL9YahgRFBUMkUN/5znfYsGFD3DFEUlQwRAJ10UUXaQ1DgqI1DJFAvfPOOzQ0dNY3IYocOnUYIoG69NJLtYYhQVHBEAnU5MmTWb8+54cciBSUCoZIoM477zytYUhQtIYhEqgNGzbw9ttvxx1DJEUdhkigLr/8cq1hSFBUMEQCdcstt7Bu3bq4Y4ikqGCIBGr48OFaw5CgaA1DJFB1dXV6p7cERR2GSKCuvPJKEokEI0eOjDuKCKCCIRKs22+/XWsYEhQVDJFADRkyJOcaxvptCa6au4q6hj0MLDuWWWMrGdC3VwETypFGBUOki4yfs5Jq74zPgqprd4s1Wxo5Z0p1Jxwrs2FWxpzxZ3XZ/qU4qGBI0EZOq2H15sYu2nv7D8SSVO0NlE9+Mo8twzunp5xYypIbhsQd47AQRMEws3HAHcB0d59qZhOBfu4+Ocv2Hwced/d/yjA2GLgHeBdY4u63m9lM4GJgsLu/3lX3QzpfV/2h19bWUlFR0SX77iwrVqygvr6eMWPGHDC2fluCL0x9jn3vtQBQUgIfLytl6aT4HhiL4ZzKoQmiYETmAw+Y2S+Bs4DHMm1kZpcD3wbKsuxnBvBVkk91njSzT7n7hKjIiBRE501H9eHmP7f/zL6lJTktlV8XcHA0HSWtQioYAMcAc4GlwKlZttkODAHWth0ws+OAo919bfTzH4DhwCuZdlRfX09TU1MnxO46zc3NRfPmra7Keu3CDazbsa/T9xvi9EmI8p+OgpDP6UnH92DGqP6A/q5yydUlBlUw3H07sCSaosq2ze8BzCzT8HHArrSfdwMDs+2rvLy8IzELqpja/K7KWtMF+yyG8zp06NCsnyU1YmoNaxv//HHEAAAI70lEQVQa2d8C3UpgUMzTUVAc57SVsnZMUAWjE+wCeqf93BvYEVMWkUMyZcoU3nzzzYxjs8ZWHvCSWpGudlgVDHffZWZ/M7NBJHvjc4HbY44l0iGVlZWUlpZmHBvQt1fsHYUceYris6TMbLqZnZ5j/Dwza31F1bXAL4GVwCvu/v8KkVGks7366qtFM88uR4YgOwx3r2pz1Vqgsc02/dJ+fBn4dHT9i8DgrswnUggTJ04kkUgwevTouKOIAGEVjDFmtsXdp2YYW+juub7cuAS4O9tg9D6MrB2KSIimT59OXV24rzqSI08QBSPqKKpyjOcqFrj75nbGJwATOpJNJC6nn346Rx99dNwxRFKKYg1D5Ei0atUqXnvttbhjiKQE0WGIyIFuvPFGEokEX/va1+KOIgKoYIgE67777tMahgRFBUMkUKeddhrdu3ePO4ZIitYwRAK1YsUKXnkl48egicRCHYZIoG6++WYSiUTGjzcXiYMKhkigHnzwQdauPeBDmUVio4IhEigzY//+/XHHEEnRGoZIoGpqali1alXcMURS1GGIBOq2224jkUhwxRVXxB1FBFDBEAnW7NmzeeONN+KOIZKigiESqIEDB7J37964Y4ikaA1DJFDPPPMMK1asiDuGSIoKhkig7rzzTh588MG4Y4ikaEpKJFCPPPIIa9asiTuGSIoKhkig+vfvT2NjY/sbihSIpqREAvX000+zbNmyuGOIpKhgiATqrrvuYubMmXHHEEnRlJRIANZvS3DV3FXUNexhYNmxzBpbyYIFC1i9enXc0URSVDBEDsL4OSup9oYuPcaaLY2cM6X6/SueerJT9z/Mypgz/qxO3accGWIvGGY2DrgDmO7uU81sItDP3Sdn2f7jwOPu/k9ZxrsDjwIz3f1pM7sTmACMc/enu+ROSIeNnFbD6s1xLewemd9mV+0NlE/u3CL0vq49p6ecWMqSG4Z06TEku9gLRmQ+8ICZ/RI4C3gs00ZmdjnwbaAsy/gg4BfAx4CZAO5+i5l9rCtCy6GL64+/traWioqKWI6dydAp1dRvS6R+Lu/bC568nUQiwcqVK2NMlr/Qzql0vlAKBsAxwFxgKXBqlm22A0OAbF8SUEqym/hep6eTI14hpqNa1W9LwOAbATq9G9CUlHRUMAXD3bcDS6Ipqmzb/B6S3xOQZfxPucbbqq+vp6mp6WCjFlRzczO1tbVxx8hLrqzXLtzAuh37CpyoPZqS6nxdf05POr4HM0b1P6R9HC5/V10hV5cYTMGIQ3l5edwR2lVMbX6urDWB3YfQzuuIqTWsbWhkfwt0K4FBZaV8o3wbb731Ftdff33c8fIS2jnNRVk7Ru/DEAnArLGVDCorpXtJCYPKSpk1tpJ7772XefPmxR1NJCX4DsPMpgNV7v5qlvHzgNPd/a7CJhPpPAP69mLppA++AGDhwoW4e0yJRA4UXMFw96o2V60FGtts0y/tx5eBT7cZH9cV2UQKqU+fPvTu3TvuGCIpoUxJjTGzSVnGFrp7rq8dKwHuzjYYvQ/jvEMJJxKHRx99lMWLF8cdQyQl9g4j6iiqcoyvb+f2m9sZvwW4pSPZROL0wAMPkEgkmDQp23MpkcKKvWCISGZPPfUUr7/+etwxRFJCmZISkTZ69epFz549444hkqIOQyRQ8+bNY+PGjcG8Bl9EHYZIoGbOnMljj2X8WDWRWKjDEAnU0qVLi+bjK+TIoA5DJFA9evSgR48ecccQSVGHIRKoqqoqNm3apDUMCYY6DJFAVVVV8bvf/S7uGCIpJS0tLXFniM1LL7105N55EZEszjjjjJJM1x/RBUNERPKnKSkREcmLCoaIiORFBUNERPKil9UGyMy+Avyru4/JMHYPcDawO7pqlLvvLGS+NnlyZb0auAZ4F7iz9TvZC83MegLzgI+SPG9j3b2hzTYLgY8A+4Amdz+/wBm7AfcDnwT2AhPSP9Y/lHMZZWkva1C/owBm9hngJ+4+tM31FwG3kjyvs9394RjifUCOrDcAE4DW391rvMDfsKWCEZjoj+1cIOM3DAJnAOe6+9bCpcosV1Yz6wdcD5wJHAM8b2ZL3X1vYVMC8A3gNXf/kZldSvLj7r/dZpuTgf/p7nG9CuTLwDHu/lkzGwz8DBgFwZ3LnFkjwfyOApjZTcDlwJ421/cApgGV0dhyM3uiva9M6ErZskbOAK5w95cKm+p9mpIKzwqSD3AHiJ7ZnQw8ZGbLzezKgiY7UNaswFnAcnffGz27fAP4RMGSfdDZwNPR5cXA8PRBMzsROB5YZGbPm9mFBc4HaRnd/UWSxaFVSOcScmQN8HcUkt/aOTrD9RXAG+6+3d3/BjwPnFPQZAfKlhWSBeP70e/o9wuYKUUdRkzM7CrghjZXj3f3R81saJabHQv8X2Aq0B2oNrP/dPc/d13SDmc9DkifhtgN9OmCeB+QJevmtCyZcnyI5LPke4APk3ymudLdt3Rl1jbanq/3zOwod383w1hBzmUOubLG8juai7s/ZmblGYZCO6+5sgIsAH4O7AIeN7MLCz01qYIRE3efBcw6yJslgHvcPQFgZn8kOY/cpX+MHcy6C0j/QurewI5OC5VFpqxm9tu0LJlyvAPMiB7wtpjZK4ABhSwYbc9XtyhPprGCnMsccmWN5Xe0g0I7r1mZWQkwvXUtyMyeBD4FFLRgaEqquJxC8tlv92j+9Wzg5ZgzZbMS+JyZHWNmfUi2/3+JKcty4ILo8vnAsjbjw4FfA5hZKXAaUOiPiU1ljNYFXksbC+lcQu6sxfQ7WgucbGYfNrMPkZyOeiHmTNkcB/zFzEqj4vF5oOBrGeowioCZTSI51/qEmT0CvEjy1Ty/cPf/ijfdB7XJei/JB+duwA/cvTmmWA8Ac83seeBvwJgo60+B37j7YjM718xeBPYDN8ewYPs4MMLMVgAlwPhAz2U+WUP/HR0DlLr7Q1HuP5A8r7PdfWO86T6oTdabgWqSr0x71t2fKnQefTSIiIjkRVNSIiKSFxUMERHJiwqGiIjkRQVDRETyooIhIiJ5UcEQEZG8qGCIiEhe/j/pXQtYKiekUQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_01.trace, varnames=['beta'], rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, beta3, beta2, beta1, alpha]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:32<00:00, 368.16draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_02:\n",
    "    alpha = pm.Normal('alpha', mu=0, sd=1)\n",
    "    beta1 = pm.Normal('beta1', mu=0, sd=1, shape=3)\n",
    "    beta2 = pm.Normal('beta2', mu=0, sd=1, shape=3)\n",
    "    beta3 = pm.Normal('beta3', mu=0, sd=1)\n",
    "\n",
    "    score_i = pm.Deterministic(\n",
    "        'score',\n",
    "        wines.score.mean()\n",
    "        + alpha\n",
    "        + beta1[0] * wines.wine_amer.values\n",
    "        + beta1[1] * wines.judge_amer.values\n",
    "        + beta1[2] * wines.flight_red.values\n",
    "        + beta2[0] * (wines.wine_amer.values * wines.judge_amer.values)\n",
    "        + beta2[1] * (wines.wine_amer.values * wines.flight_red.values)\n",
    "        + beta2[2] * (wines.judge_amer.values * wines.flight_red.values)\n",
    "        + beta3 * (wines.wine_amer.values * wines.judge_amer.values * wines.flight_red.values)\n",
    "    )\n",
    "\n",
    "    sigma = pm.HalfCauchy('sigma', beta=5)\n",
    "    score_obs = pm.Normal('score_obs', mu=score_i, sd=sigma,\n",
    "                          observed=wines.score.values)\n",
    "    \n",
    "    m_02.name = 'full interaction, v2'\n",
    "    m_02.trace = pm.sample(5000, tune=1000, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEcCAYAAAAP5CkrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8VPW57/FPpKkYY9Viqrs9YA4cfUi3bW098dhuK6EVvGwtltqWw0EBrdXa1ipeGq3V2mP3dheLyvGG3CvbYr1UpV4qYkzR2Bdsayt2x4dLGkBAiGxAwiQWJeePtaJpTEISZmatWev7fr14kZm1Zq3vrEzmmef3WzNT1NbWhoiISDbsF3UAERFJDhUVERHJGhUVERHJGhUVERHJGhUVERHJGhUVERHJmg9FHUCkv8xsf2AuMBR4C/iuu68ys68CNwPrw1WvB14CHgMOAC5091fM7ETgn9z937rZ/reBCcAeoBj4kbs/14+cFwFHAHcD17n7xWbWCAx399YO600Kr6vuxTbnAQvd/alulg8BPuPui/qady/7rQIucvdx2dyuJIc6FSlkFwDN7n4C8H3g9vD644Cr3L0q/FcLjCYoKhcD55tZEfAD4LauNmxm44BRwJfdvYqguNxrZof1N6y7v+HuF/f39n30JeCf8rQvkfeoU5FC9kngSQB3dzOrCK8/DvismV0KLAN+CDQTdCkHALuA8cBvOnYKnVwITHH33eH2/2pmx7r71rBLGBT++2fgKuCLwABgmrs/EHZBtwHbgHeAP5hZOUF3cUK4jxnhdZuBiR13bmbfDzO2hbeZ3lXIsHP4IfA3go5tIXATUA2UmFkd8FdgOlAEbAXOAz4L/Ft4u2eAKncfGW7zt8CPgWHAdwm6tDbgq532PRf4H+Exvc3d7+3mWEqKqFORQvYn4AwzKzKzE4BPmNkAYDFB53ISUApcRPDEeThBsbiH4Anyz2Y2w8yu6mLbHwcaOl7h7ls7XHzW3b8AnAD8d3c/ERgJ/MjMDgHuAv63u59M8KTelbvcfQTQSNB1AWBmnwS+CZxIUKzOMjPr4TgcCXwtzHKVu79LUFjuc/fHgJkEQ4NVwBMERRBgoLt/0d1vAAaa2ZFm9g/AYe7+MnA08M/hfftP4JQOGQ8iOL5jgVOBd3vIJymiTkUK2RygAlgKvAC85O7vmtkcd98OYGaPAl9z9z3ApeF1VxN0EdcSFJ+fmNnR7r6yw7bXAoOBHe1XmNkpwCvhRQ///xRwnJk9F14uBsqBwzts7wWCV/Qd/c3d/xD+XEcw1LY8vHwMQaFYEl4+FDiqwz47W+Hu7wDvmFlLF8srgDvDulQMrOp0HwBmA+cCbxPMUwFsAeabWTMwHHixfWV33xl2gvcAHwEWdJNNUkadihSySmBJ+Er6AaAhnCt5xcz+W7jOlwkm6QEws48B5u5LgRKCV9htwIGdtj0H+LGZfSi83dHALN5/Rb4n/P81oCbsAr4E/BpYA2zoMBxX2UX2D5vZseHPXwRe7bDMgb8AI8PtzuP9YtaVrj7Abw/v/307cG64rauA33a6DxAMm51B0MHdZ2YHAzcA44BvAS0Ew2cAhB3Nce7+VYIhwJ+3HytJNxUVKWSrgEvN7EXg/xLMgbQRPAk+bGa1BIVjZofbXAvcGP58J/A74B+AP3fcsLsvBP4APG9mvyd49T7B3bd0yrAIaDazpQTFq83ddxIMs/3SzJYQdB2dvQ18P9z2xwgKVvu+/0zQpTxvZv9B0KVs6P1hAWAFMCY84eA7YZbnCYbFPlCg3L05PAb1Yf63CDqsFwk6wRaCIcF2bwBHhHM2i4Gbw25JUq5In1IsIiLZok5FRESyRkVFRESyRkVFRESyRkVFRESyRkVFRESyJnHnlb/00ks6nU1EJMeOO+64oq6uT+IpxZHdofr6eioqKva+YkTing+UMRtmzpzJpk2buO6666KO0qO4H0eIf8YI83VZUCCBnYpI2l1wwQXU19dHHUNSSnMqIiKSNepURBKmqqqKTCbDsmXLoo4iKaSiIpIwkyZNYuPGjVHHkJRSURFJmEmTJmlORSKjORWRhNm9eze7d++OOoaklDoVkYQZNWqU5lQkMioqIgnzrW99iw0b+vr1KyLZoaIiElPrtmY4f/5yGpp2MbTsQGZPrGTIoJK93m7ChAmaU5HIqKiIZMHkucuo8aacbX/VlmZOmlrTx1s17HWNkVbG3MnH9y+USBciLSpm9oa7H9HNsirgIncfl99UEgejb6ll5ebmqGP0YO9P2IWgxpsor348wgT9O45HH17K05eNyHIWyQZ1KhJLcX7CyPXnLbUPe63a8n5R3a8IhpWVsnjK3o/L/fffz4YNG5gyZUrOMmZD3D9XS/onL0XFzI4BpgEDgMOA77h7XYflzwGvAcMJPqjsm+Gio8zsSeBjwCJ3/4mZjQCuJzgduhQY7+4r83E/RDrL9bBXuz1twRBY77qKUsCYvpd1NfQluZCvTuUfgcvdfYWZjQcmA3Wd1qlz94vM7GLgGuBhYCBwFkExWgf8JNzWBHffaGbXAF8Hfta+kcbGRlpaWnJ9f7rU2toa6wnSuOeDIOOIm55m7fa4v8+i8Ie/oh/6gr4exyMPKebuMYNzlOWD4v43E1W+njrMfBWVDcCPzawFOAh4q4t1ng3/rwPGhD+/6u5vA5jZOx22Nd3MmoFPAC903Eh5eXl2k/dB3Nv5uOeDIGNt9eioY/Qo18dx1LRa1jQ1s6etb8Ne7Qrls78K5fEY54xxzJevd9RPB65394nACrr+LP7jwv//CfhL+HNX340yE5js7pOAjd1sS6RgzZ5YybCyUgYUFTGsrJTZEyv7dPtLLrmECRMm5CidSM/y1aksAB4ws23A6wTzKp1NMrMpwC7gHOBTPWxrqZntAjYDH89BXpHIDBlU0qfOpLOxY8fGeshGki0vRcXdpxFM1He+vuPpxFe7+2sdLj8X/vu7dd093qe0iETszTffZNu2bVHHkJTSKcUiCXP22WcXxJyKJFMsioq7V0WdQSQpLr/8ctavXx91DEmpWBQVEcmeM888U3MqEhl9n4pIwrzxxhs0NeX+DZkiXVGnIpIw48aN05yKREZFRSRhqqurWbduXdQxJKVUVEQS5tRTT9WcikRGcyoiCbN+/Xo2bdoUdQxJKXUqIglzzjnnaE5FIqOiIpIw1157LWvXro06hqSUiopIwpx88smaU5HIaE5FJGEaGhr0jnqJjDoVkYQ577zzyGQyjB4d7++lkWRSURFJmBtuuEFzKhIZFRWRhBkxYoTmVCQymlMRSRh3569//WvUMSSl1KmIJMyFF15IJpPh9NNPjzqKpJCKikjC/Mu//AuNjY1Rx5CUUlERSZgvfOELHHrooVHHkJTSnIpIwrz66qusWrUq6hiSUupURArYuq0Zzp+/nIamXQwtO5DZEyv53ve+RyaT4Stf+UrU8SSFVFREsmzy3GXUeP6/eXHVlmZOmloDJ1wJQHn14znb10grY+7k43O2fSlcKioSS6NvqWXl5uaoY/SgIeoAkarxpiwVrXgcx6MPL+Xpy0ZEHSMRel1UzGwesNDdn+rFukOAz7j7or2sVwa8AHza3Vs7LXsOKAEmAcOA64B3gDnuPtPMfgucDBzS+bZS+OL8B15fX09FRUXUMQComlpD49bMe5fLB5Vw66hDaWhoYOzYsREm27s4HUfJnlx1Kl8ChgPdFhUzOwW4CTiih+2cC6wBHgMqgV3AC2b2mLufYWaN2Qoskk1RDYE1bs1w1sIMsD9Tlmn4S/Kvr0XlYjO7Mrzd+e6+2sy+D4wH2oCFwB1ANVBiZnXADuB6gjPNSoHx7r4S2EPQaby0l31WAKvdfRuAmT0PnAQ80NXKjY2NtLS09PFuZUdra2usPx4jrvkuenQ9a7fv7nBNPIZEelYIGXMnacNf7Y48pJi7xwx+73Jc/2baRZWvpw6zr0Wlzt1vMrPTgZ+b2bXAN4ETw+WLgd8RdCDD3f0xM7sYmODuG83sGuDrwM/cfTGAme1tnx8hKEztdgIHd7dyeXl5H+9S9sS9nY9rvtoOmeKasaM4ZRw1rZY1Tc3saYP9imBYWSmLp4yIVcbuKOO+i2O+vr5P5ffh/3WAAccARwJLwn+DgKM63WYDMD2ckxkJFPdxn28BB3W4fBCwvY/bEEmk2RMrGVZWyoCiIoaVlTJ7YiXLly9nxYoVUUeTlOprp3I8QUH5IvAq4MBfgNPcvc3MLgNeAap4v2DNBIa5+04zmw8U9XGf9cBRZvZRoJlg6OvmPm5DJJGGDCph8ZS/P6nh3K9dSSaT4Rvf+EZEqSTN+lpUTjCzZwnmT85z97VmtgR43sz2B5YRdCYrgB+Z2R+BBcBSM9sFbAY+3pcduvtuM5tCMKy2H8HZXxv6mFskNW6//XYaGuI1VyHp0eui4u6Turl+KjC109UvEwyPQTB539N2y3ux70X0cCaZiLzvmGOOYcCAAVHHkJSK+2d//dLMPtnVgvB9Kj2djiySSnV1dbz88stRx5CUiu076t29ai/Lz8hTFJGCcs0115DJZBg/fnzUUSSFYltURKR/ZsyYwZo1a6KOISmloiKSMGbGnj17oo4hKRX3ORUR6aPa2lqWL18edQxJKXUqIglz/fXXk8lkOPfcc6OOIimkoiKSMHPmzGH16tVRx5CUUlERSZihQ4fy9ttvRx1DUkpzKiIJ88wzz1BXVxd1DEkpFRWRhLnxxhuZMWNG1DEkpTT8JZIw9957L6tWrYo6hqSUiopIwgwePJjm5uaoY0hKafhLJGGeeuopli5dGnUMSSkVFZGEuemmm5g1a1bUMSSlNPwlkjALFy5k5cqVUceQlFJREUmYI444gm3btkUdQ1JKw18iCbNo0SJqamqijiEppaIikjC/+MUvmDdvXtQxJKU0/CWSMA8++KDmVCQy6lREEuawww7j0EMPjTqGpJQ6FZGEefjhh3n99depqKiIOsoHrNua4fz5y2lo2sUnPvIhFnzsSIYMKok6lmSRiopIwkyfPp1MJsMll1yy13Unz11GjTflIdUHrduxm5Om5ueEgpFWxtzJx+dlX2kX26JiZpOAnwK3uvs0M7sUOMLdq83sbOBG4BF3r44yp6TL6FtqWbm5GWiIOkr3TrgSgPLqxyMOEh813rQPxyN+v+ujDy/l6ctGRB2jS7EtKqH7gLvM7N+B44GHANz9QTMrBYZHGU7S5+nLRlBfXx/LoaWO4ppx1LRa1jQ1s6ft/euO+lgpsydWxnIYLK7HMc7iXlQABgLzgcWoiIjs1f3338+GDRv69GQY5TDYqi3NGgZLkKK2tra9rxWBcPhrePvw1t4ut2tsbGxraWnJb9hQa2srAwcOjGTfvRH3fBC/jBc9up6123dHHUOkR0ceUszdYwbnbX8VFRVF3S0rhE6lT8rLyyPbd9xb5bjng/hlrO0iS9wydpbJZHjttdf43Oc+F3WUbo2aVsvqLc20AfsVwbCyUhZPid8cQdx/13HMp/epiCRMSUkJBxxwQNQxejR7YiWDDy5mQFERw8qCORVJhsR1KiJpt2DBgj7PqeTbkEElzDhrcKwzSv8UTFFx93lRZxApBLNmzSKTyfDDH/4w6iiSQnEf/hpvZlM6Xxm+T0XvTxHpwuLFi/UlXRKZ2HYqYWcyr5tlDwIP5jOPSKEoLi6muLg46hiSUrEtKiLSP/PmzWPjxo2ar5BIxH34S0T6aN68eTzyyCNRx5CUUqcikjDPPfcc9fX1UceQlFKnIiIiWaOiIpIwM2fO5IEHHog6hqSUiopIwtx///08+eSTUceQlNKcikjCPPPMM5pTkcioUxERkaxRURFJmDvvvJNf/epXUceQlNLwl0jCLFq0iObm5qhjSEqpqIgkzJNPPqk5FYmMhr9ERCRrVFREEua2227j3nvvjTqGpJSGv0QSZsmSJZpTkcioqIgkzGOPPaY5FYmMhr9ERCRrVFREEubmm29mzpw5UceQlNLwl0jCvPjii+zcuTPqGJJSKioiCfPQQw9pTkUio6IiIn2ybmuG8+cvp6FpF0PLDmT2xEqGDCqJOpbEhIqKSIGbPHcZNd7UxZKGnO971ZZmTppasw9b6FvGkVbG3MnH78P+JNdiW1TMbBLwU2AWcBJB1iLg28CngBuBR9y9OqqMUphG31LLys37+j6O3D9hywfVeBPl1Y/nea/5+10ffXgpT182Im/7y4XYFpXQfcAw4HZ3f8TMTgH+1d3HmlkpMDzaeFKI9vWPtr6+noqKiiylyY1cZhw1rZbVTc20tQWXiwcUsWRKVZ+HwNJ+HJMq7kUF4HJgR/jzh4DWCLOIxFLXQ2D5eYW9+922fRgC63tGDYHFW+yLiru/CWBmBtwMnNXT+o2NjbS0tOQj2ge0trbG+qybuOeD7Ga86NH1rN2+Oyvb+iANf0Ul/0NgufldH3lIMXePGbxP24jqb7qn7i32RQXAzEYCdwLnuLv3tG55eXleMnUl7q1y3PNBdjPW5ui+xv04TpgwgR07drBo0aKc7WPUtFrWNDWzpw32K4JhZaUsntK3YcW4H0eIf8Y45ot9UQkLym3Aqe6+Nuo8InG3YMGCnL96nT2x8gOnFYtAARQV4Fbgw8D8YAQMd/cLo40kkm5DBpX0uTORdIh9UXH3z0SdQaSQXHfddbz55pvceeedUUeRFIr7B0qON7Mpna80s7MBvT9FpAvr169n06ZNUceQlIptp+Lu84B53Sx7EHgwn3lECsXcuXNjf5afJFfcOxURESkgKioiCXP11Vczbdq0qGNISqmoiCTM1q1b2bFjx95XFMmB2M6piEj/3HPPPZpTkcioUxERkaxRpyKSMFdccQVbt25l7ty5UUeRFFKnIpIwLS0tvP3221HHkJRSpyKSMHfccYfmVCQy6lRERCRr1KmIJMyll17Kf/3Xf/HLX/4y6iiSQupUREQka9SpiCTMrbfeqjkViYw6FRERyRp1KiIJ893vfpdt27Zx3333RR1FUkidikjCHHDAAey///5Rx5CUUqcikjA333yz5lQkMupUREQka9SpiCTMt7/9bbZv386vf/3rqKNICqmoiBSAdVsznD9/OQ1NuxhadiCzJ1YyZFBJl+sOGjSItra2PCcUCaioiPTB5LnLqPEmoCGyDKu2NHPS1Joe1jgRBkF59eM52f9IK2Pu5ONzsm0pfL0uKmY2D1jo7k/1Yt0hwGfcfVEP61wGjAsvPuHuN3Ra/hxQAkwChgHXAe8Ac9x9ppn9FjgZOMTdW3t7PyS7Rt9Sy8rNzVneanRP2LJ3Nd6UxYJVCL/r9zMefXgpT182IsIs8ZerTuVLwHCgy6JiZkOB/wP8L2AP8LyZ/cbdX+m06rnAGuAxoBLYBbxgZo+5+xlm1pib+NJb2f4Dq6+vp6KiIqvbzLYoMlZNraFxa+a9y+WDSnjuypFdrjt58mS2b9/Ob37zm3zF6xf9rpOpr0XlYjO7Mrzd+e6+2sy+D4wH2oCFwB1ANVBiZnXADuB6gjPNSsN1/wqc6u7vAphZMdBdt1EBrHb3beG6zwMnAQ/0MbvIPnl/6AuifoXduDXTfbdw+NlweO6Gv9ppGEy60teiUufuN5nZ6cDPzexa4JvAieHyxcDvgJuA4e7+mJldDExw941mdg3wdXf/GfCmmRUBU4GX3X1lN/v8CEFharcTOLi7gI2NjbS0tPTxbmVHa2trrN8fkKt8Fz26nrXbd2dxi4UwJCLZGQYrhN913zMeeUgxd48ZnIMsfy+q55yeure+FpXfh//XERSDY4AjgSXh9YcCR3W6zQZgupk1A58AXgAws4HAHIIicXEP+3wLOKjD5YOA7d2tXF5e3ou7kRtxb5Vzla82i9uM+zGEaDKOmlbLmqZm9rTBfkUwrKyUxVO6H3rUccyOuGeMY76+vvmxvdf9IvAq4MBfgJHuXgXMA14hmCdp3/ZMYLK7TwI2AkVhh/Io8Gd3v7B9GKwb9cBRZvZRM/swwdDXi33MLVLQZk+sZFhZKQOKihhWVsrsiZXdrjthwgSuuuqqPKYTeV9fO5UTzOxZgvmT89x9rZktIZho3x9YRtCZrAB+ZGZ/BBYAS81sF7AZ+DhwFjAC2N/MTgu3fbW7f6BYuPtuM5tCMKy2H8HZXxv6fE9FCtiQQSU9diYdmRlNTU17X1EkB3pdVMJOo6vrpxIMhXX0MmDhzwu72eTAPux7Ed2cSSYif+/HP/5xrOf2JNni/tlfvzSzT3a1IHyfyhF5ziMiIj2I7TvqwzmanpafkacoIgVl3LhxvPXWWzzxxBNRR5EUim1REZH+OfbYY9myZUvUMSSlVFREEqa6ulpzKhKZuM+piIhIAVGnIpIwX/va19i5cydPP/101FEkhVRURBLm85//PJs3b446hqSUiopIwlxxxRWaU5HIaE5FRESyRp2KSMJ85Stfobm5mWeffTbqKJJCKioiCfPlL39ZcyoSGRUVkYT5wQ9+oDkViYzmVEREJGvUqYgkzGmnnUZzczNLly6NOoqkkIqKSMKceeaZvPHGG1HHkJRSURFJmIsvvlhzKhIZzamIiEjWqFMRSZiTTz6ZXbt28eKLH/h2bpGcU1ERSZhvfvObbNq0KeoYklIqKiIJc8EFF2hORSKjORUREckadSoSe+u2Zjh//nIamnYxtOxAZk+sZMigkqhjxVZVVRWZTIZly5ZFHUVSSEVFejR57jJqvCnPe23odsmqLc2cNLUmj1m6NtJ2Mnfy8VHH6NKkSZPYuHFj1DEkpYra2tqiztAlM5sE/BS41d2nmdmlwBHuXm1mZwM3Ao+4e3Wnm+7THRp9Sy0rNzfvyyZEEuXow0t5+rIRWd9ufX09FRUVWd9uNsU9Y4T5irpbEPdO5T7gLjP7d+B44CEAd3/QzEqB4dne4b788egBuO+6yjhqWi2rm5ppf/1TPKCIJVOqIhsCi/tx3L17N/X19Xz605+OOoqkUNyLCsBAYD6wmBwUEdm7/A+BdT/8BbD73bZIh8AqP3EAD8S4qIwaNUpzKhKZuA9/DW8f3trb5XaNjY1tLS0t/d7vRY+uZ+323f2+vUgSHXlIMXePGZzVbba2tjJw4MCsbjPb4p4xqnwVFRUFO/zVZ+Xl5ft0+9p9eAUa92GRuOeD7jOOmlbLmqZm9rTBfkUwrKyUxVOyP87fG4V8HONEGfddHPPpfSpSEGZPrGRYWSkDiooYVlbK7ImVUUeKrUwmw7506yL7InGdiiTTkEElkXUmheb000/XnIpEpmCKirvPizqDSCH4zne+w4YNG6KOISkV9+Gv8WY2pfOV4ftUOr8/RUQIPlDytNNOizqGpFRsO5WwM5nXzbIHgQfzmUekUOzYsYOdO3dGHUNSKrZFRUT6Z8yYMZpTkcioqIgkzCWXXMLrr78edQxJKRUVkYQZO3asvk9FIhP3iXoR6aM333yTbdu2RR1DUkqdikjCnH322ZpTkcioqIgkzOWXX8769eujjiEppaIikjBnnnmm5lQkMppTEUmYN954g6amfH9bp0hAnYpIwowbN05zKhIZFRWRhKmurmbdunVRx5CUUlERSZhTTz1VcyoSGc2piCTM+vXr2bRpU9QxJKXUqYgkzDnnnKM5FYmMiopIwlx77bWsXbs26hiSUioqIglz8skna05FIqM5FZGEaWho0DvqJTLqVEQS5rzzziOTyTB69Oioo0gKqaiIJMwNN9ygORWJjIqKSMKMGDGiX3Mq67ZmOH/+chqadjG07EBmT6xkyKCSHCSUJFNREYmByXOXUePZ/ryuhn7fctWWZk6aWpPFLN3pPuNIK2Pu5OPzkEGySUVF8mr0LbWs3Ny8l7X6/2SYP4WQsbDVeBPl1Y9HHYNc/q6PPryUpy8bkbPtRyG2RcXMJgE/BWYBJxFkLQK+DXwKuBF4xN2ro8oofbe3P6D6+noqKirylKZ/4p6xrq6OxsZGxo8f36fbVU2toXFr5r3L5YNKeO7KkdmO9564H0cojIxxE9uiEroPGAbc7u6PmNkpwL+6+1gzKwWGRxtPJDuyP/x1MNe8sm+v8hu3ZvLQKWj4K2niXlQALgd2hD9/CGjtaeXGxkZaWlpyHqorra2tsX7TWZzzXfToetZu3x1eKoShpULIWNjSMPzVG0ceUszdYwZ3uSyqv+meurfYFxV3fxPAzAy4GTirp/XLy8vzkKprcW+V45yvNswV54zt4p6xqqqqX5/9NWpaLWuamtnTBvsVwbCyUhZPyd14f9yPI8Q/YxzzFcQ76s1sJPAIcI67e9R5ROJs6tSpXHHFFX2+3eyJlQwrK2VAURHDykqZPbEyB+kk6WLfqYQF5TbgVHfXO7pE9qKyspLS0tI+327IoJKcdiaSDoXQqdwKfBiYb2bPmdmMqAOJxNmf/vSn2M6dSfLFvlNx989EnUGkkFx66aVkMhnGjh0bdRRJobgXlfFmtsXdp3W80szOBqoJ5llEpINbb72VhgadnSbRiG1Rcfd5wLxulj0IPJjPPCKF4thjj2X//fePOoakVCHMqYhIHyxfvpwVK1ZEHUNSKradioj0z5VXXkkmk+Eb3/hG1FEkhVRURBLm9ttv15yKREZFRSRhjjnmGAYMGBB1DEkpzamIJExdXR0vv/xy1DEkpdSpiCTMNddcQyaT6fNH34tkg4qKSMLMmDGDNWvWRB1DUkpFRSRhzIw9e/ZEHUNSSnMqIglTW1vL8uXLo44hKaVORSRhrr/+ejKZDOeee27UUSSFVFREEmbOnDmsXr066hiSUioqIgkzdOhQ3n777ahjSEppTkUkYZ555hnq6uqijiEppaIikjA33ngjM2bou+wkGhr+EkmYe++9l1WrVkUdQ1JKRUUkYQYPHkxzc3PUMSSlNPwlkjBPPfUUS5cujTqGpJSKikjC3HTTTcyaNSvqGJJSGv4SSZiFCxeycuXKXq27bmuG8+cvp6FpF0PLDmT2xEqGDCrJcUJJMhUVkRybPHcZNd6U/x0/8XifVl+1pZmTptbkKEx3evdlYiOtjLmTj89xFsmGnBcVM5sHLHT3p3qx7hDgM+6+qId1vgtMAtqAm93911mKKlk0+pZaVm7u72RxIXxrYSFSGbclAAAHOUlEQVRkTI4ab6K8um9FMnvi/rvue76jDy/l6ctG5CBL/DqVLwHDgS6LipkdBnwH+CwwEPhPM3vA3dvyF1F6o78P2Pr6eioqKrKcJrvinrGqqopMJsOyZcv2uu6oabWsbmqmLfwLKh5QxJIpVXkZAov7cYT4Z4xjvnwVlYvN7Mpwf+e7+2oz+z4wnqDjWAjcAVQDJWZWB+wAric4maAUGO/uK83sWHd/x8zKgVYVFIm7vA9/nXAlQL9e2e9+ty3PQ2B7f5Wtoa/CUtTWltvn5HD46zV3v8nMTge+BVwL3AO0v5xdTNCBfB4Y7u7VZnYx8Ii7bzSza4Aid/9ZuM3vATcA0939ho77a2xsbGtpacnpfepOa2srAwcOjGTfvZGPfBc9up6123fndB8i0jdHHlLM3WMGZ217FRUVRd0ty1en8vvw/zpgKnAMcCSwJLz+UOCoTrfZAEw3s2bgE8AL7Qvc/XYzuwd40sxGuvt7L63Ky8tzcgd6I46taEf5yFe7j9uP+zGE+Gd8+OGHef3117nkkkt6tf6oabWsaWpmTxvsVwTDykpZPCU34+0dxf04QvwzxjFfvt6n0t67fhF4FXDgL8BId68C5gGvAHs6ZJoJTHb3ScBGoMgCD5tZEbAbeDu8jYiEpk+fzoIFC3q9/uyJlQwrK2VAURHDykqZPbEyh+kk6fLVqZxgZs8SzJ+c5+5rzWwJ8LyZ7Q8sI+hMVgA/MrM/AguApWa2C9gMfNzd3cz+DLwYbutJd6/N030QKQiPPvoo7t7r9YcMKslLZyLpkPOiEnYaXV0/lWAorKOXAQt/XtjN7W4gmE8RkS4cfPDBHHTQQVHHkJSK2ynFIrKP7r//fjZs2BC7sXZJB332l0jC3HXXXSxc2GWjL5Jz6lREEuaJJ57gtddeizqGpJQ6FZGEKSkp4YADDog6hqSUOhWRhFmwYIHmVCQy6lREEmbWrFk89NBDUceQlFKnIpIwixcvpr6+PuoYklLqVEQSpri4mOLi4qhjSEqpUxFJmHnz5rFx40bNqUgk1KmIJMy8efN45JFHoo4hKZXzj77Pt5deeilZd0hEJIaOO+64Lj/+PnFFRUREoqPhLxERyRoVFRERyRoVFRERyRqdUtxPZnYwwReJfQT4MDDF3V/stM4FwIXAO8CN7v7bvAcNcnwV+Lq7j+9i2W3AicDO8Kox7r4jn/nCHD1ljPQ4mtkBBL/rjxEcp4nu3tRpnUeBwwi+kbTF3U/LQ679gDuBzxB8C+q33H11h+WRP/56kTEWj78wy/8C/i38NtqO158JXEdwHOe4+8wI4rVn6S7jZcC3gPbH5YXel29qyyIVlf6bAixx91vNzIBfAZ9rX2hmRwCXAP8TGEjwLZeL3f3tfIYM/2hPAf7UzSrHAae4+5v5S/X3esoYk+P4HWCFu//EzMYB1wI/6LTOUcA/uns+z3w5Cxjo7p83sxOAXwBjIDbHrceMocgffwBmdhVwDrCr0/XFwC1AZbjsBTN7zN03xyVj6DjgXHd/Kb+pPkjDX/13CzAj/PlDQGun5ccDL7j72+Err9XAp/OYr10dwZPiB4SvIo8C7jGzF8zsvLwme1+3GYnHcTwReCr8+Ung5I4Lzexw4BBgkZk9b2Zn5DuXu/+BoIC0i8Nxgx4yxujxB7AGGNvF9RXAanff5u5/A54HTsprsvd1lxGConJ1+Pi7Oo+ZPkCdSi+Y2fnAZZ2unuzuy8NXhAuASzst/wjQsY3fCRwcQcb7zayqm5sdCPw/YBowAKgxs/9w91dilDEOx3Fzhwxd7f/DBK/AbwM+SvBqdpm7b8lVzlDnY/OumX3I3d/pYllOj1sPesqY18dfT9z9ITMr72JRXI5jTxkh+Pr1O4C3gN+Y2RlRDberqPSCu88GZne+3sw+RfDLvMLdazstfgvo+EXhBwHb851xLzLAbe6eATCzZwnGvnPyR93PjJEfRzN7uEOGrvb/BnB3+ES5xcxeBgzIdVHpfGz2CzN0tSynx60HPWXM6+Ovn+JyHLtlZkXAre1zUWb2OPBZIJKiouGvfjKzTwIPAOPd/ckuVlkGfNHMBoaT+hXAq/nM2AtHE7yqHhCOHZ8I/DHiTJ3F4Ti+AJwe/nwasLTT8pMJHguYWSlwDJCPjwl+L1c4X7Giw7I4HDfoOWMhPP7qgaPM7KNm9mGCoa8X93KbfPsI8KqZlYYF5ktAZHMr6lT6718JJkBvC+bp2eHuY8xsCsEY7GNmNp3gCWg/4Efu3nneJRKdMt4L/IHgrKVfuvtfok0XiNlxvAuYb2bPA38DxocZfw486O5PmtkpZvYHYA9wTZ4mnn8DjDKzOqAImByz49abjHF9/I0HSt39njDv7wiO4xx33xBtukCnjNcANQRn2C1x9yeiyqWPaRERkazR8JeIiGSNioqIiGSNioqIiGSNioqIiGSNioqIiGSNioqIiGSNioqIiGTN/wdZCFZHNdu8EAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_02.trace,\n",
    "              varnames=['alpha', 'beta1', 'beta2', 'beta3'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "INFO (theano.gof.compilelock): Refreshing lock /home/petrushev/.theano/compiledir_Linux-5.0-arch1-1-ARCH-x86_64-with-arch--3.7.3-64/lock_dir/lock\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, beta3, beta2, beta1, alpha]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:29<00:00, 400.44draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_03:\n",
    "    alpha = pm.Normal('alpha', mu=0, sd=1)\n",
    "    beta1 = pm.Normal('beta1', mu=0, sd=1, shape=2)\n",
    "    beta2 = pm.Normal('beta2', mu=0, sd=1, shape=3)\n",
    "    beta3 = pm.Normal('beta3', mu=0, sd=1)\n",
    "\n",
    "    score_i = pm.Deterministic(\n",
    "        'score',\n",
    "        wines.score.mean()\n",
    "        + alpha\n",
    "        + beta1[0] * wines.judge_amer.values\n",
    "        + beta1[1] * wines.flight_red.values\n",
    "        + beta2[0] * (wines.wine_amer.values * wines.judge_amer.values)\n",
    "        + beta2[1] * (wines.wine_amer.values * wines.flight_red.values)\n",
    "        + beta2[2] * (wines.judge_amer.values * wines.flight_red.values)\n",
    "        + beta3 * (wines.wine_amer.values * wines.judge_amer.values * wines.flight_red.values)\n",
    "    )\n",
    "\n",
    "    sigma = pm.HalfCauchy('sigma', beta=5)\n",
    "    score_obs = pm.Normal('score_obs', mu=score_i, sd=sigma,\n",
    "                          observed=wines.score.values)\n",
    "    \n",
    "    m_03.name = 'full interaction, v3'\n",
    "    m_03.trace = pm.sample(5000, tune=1000, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEcCAYAAAAP5CkrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8VfWZ7/FPoCjE4A1TmfaAKRx9TA+d2uHEY61yaQEto9JaZ+RwqID3Wx1E60FrdbzUMsUi0lZU5KZMD7ZqRSpaEDFFwReM1Sqd+HCJAeQaGW5hJzSanD/2ShvT3Nl7r73X+r5fL18k6/rdy5397N/v2Ze8+vp6REREUqFL2AFERCQ6VFRERCRlVFRERCRlVFRERCRlVFRERCRlVFRERCRlPhN2AJHOMrOjgblAP+AAcIO7bzCzbwMPAluDTe8G3gJeAHoA17j7u2Z2DvA1d/+3Fo5/NTAWqAO6AT9w99c6kfNaoDfwKHCXu19vZhXA6e5e02i78cGyye045jxgobu/3ML6vsCX3X1xR/O2cd4hwLXuPjqVx5Xo0EhFctlVQJW7nwV8D/h5sHwgcJu7Dwn+KwVGkCwq1wNXmFke8C/Aw80d2MxGA8OBb7j7EJLF5SkzO6mzYd19p7tf39n9O+jrwNcydC6Rv9BIRXLZF4GXANzdzaw4WD4Q+IqZTQTWAP8XqCI5SukBHALGAL9pPFJo4hpgkrvXBsf/wMzOcPc9wSihV/DfPwK3AecCXYFp7v7rYBT0MLAX+Bh408yKSI4uzgrO8ViwbBcwrvHJzex7Qcb6YJ8ZzYUMRg7/F/gzyRHbQmAKMBnIN7NVwAfADCAP2ANcDnwF+Ldgv1eAIe4+NDjmb4EfAv2BG0iO0uqBbzc591zgvwfX9GF3f6qFaykxopGK5LJ3gAvMLM/MzgI+b2ZdgWUkRy6DgALgWpIPnCeTLBaPk3yA/KOZPWZmtzVz7M8B5Y0XuPueRr++6u5nA2cBX3D3c4ChwA/M7HhgJvC/3X0YyQf15sx098FABclRFwBm9kXgUuAcksXqW2ZmrVyHU4DvBFluc/dPSBaWX7r7C8AsklODQ4AlJIsgQHd3P9fd7wG6m9kpZvZ3wEnu/jZwGvCPwW37T+C8Rhl7kry+FwPnA5+0kk9iRCMVyWVzgGJgJfAG8Ja7f2Jmc9x9H4CZLQK+4+51wMRg2e0kRxF3kiw+/2pmp7n7+kbH3gz0AfY3LDCz84B3g189+PdLwEAzey34vRtQBJzc6HhvkHxG39if3f3N4OdVJKfa1ga/DyBZKJYHv58AnNronE295+4fAx+bWXUz64uBR4K61A3Y0OQ2AMwGLgMOk+xTAewG5ptZFXA6sLphY3c/GIwEHweOBRa0kE1iRiMVyWUlwPLgmfSvgfKgV/Kumf23YJtvkGzSA2BmnwXM3VcC+SSfYdcDxzQ59hzgh2b2mWC/04An+Osz8rrg3/eBFcEo4OvAr4BNwLZG03ElzWQ/yszOCH4+F1jXaJ0DfwKGBsedx1+LWXOa+wC/Ov769+3AZcGxbgN+2+Q2QHLa7AKSI7hfmtlxwD3AaOBKoJrk9BkAwYhmoLt/m+QU4E8arpXEm4qK5LINwEQzWw3cR7IHUk/yQfA5MyslWThmNdrnTuD+4OdHgN8Bfwf8sfGB3X0h8Cbwupn9nuSz97HuvrtJhsVAlZmtJFm86t39IMlptifNbDnJUUdTh4HvBcf+LMmC1XDuP5IcpbxuZv9BcpSyrf2XBYD3gFHBCw6uC7K8TnJa7G8KlLtXBdegLMh/gOQIazXJkWA1ySnBBjuB3kHPZhnwYDBakpjL06cUi4hIqmikIiIiKaOiIiIiKaOiIiIiKaOiIiIiKaOiIiIiKRO515W/9dZbejmbiEiaDRw4MK+55VF8SXHablBZWRnFxcVtb5gFlDU9ciXrrFmz2LFjB3fddVfYUdolV65rg1zKm6aszRYUiOBIRUTgqquuoqysLOwYEkPqqYiISMpopCISQUOGDCGRSLBmzZqwo0jMqKiIRND48ePZvn172DEkhlRURCJo/Pjx6qlIKNRTEYmg2tpaamtrw44hMaSRikgEDR8+XD0VCYWKikgEXXnllWzb1tGvYBE5cioqIhE0duxY9VQkFCoqIllmy54EV8xfS3nlIfoVHsPscSX07ZXfoWMkEgmqq5v7unqR9FJRETlCE+auYYVXpuXYG3ZXMWjqis4f4FcvAjDUCpk74cwUpRJpmYqK5JwRD5WyfldViAnKQzx356zwSoomvxh2jL9x2skFLL15cNgxJIVCLSpmttPde7ewbghwrbuPzmwqyXZhPghl4oMEh08rZVNlFXX10CUP+hcWsGxSx29zLn3ooUSH3qcikkW27ElQ+0kddcFnbfc9MZ/Z40o6fJz9+/dz8ODBFKcTaVtGRipmNgCYBnQFTgKuc/dVjda/BrwPnE7yI5UvDVadamYvAZ8FFrv7v5rZYOBukgWxABjj7uszcTtEWpKuvkrFnsSR9VSee1H9FMmojHyfipldCvynu79nZmOAoe5+VcP0V1BU5rj7k2Z2Pcni8hzwC+AfSBajLe5+UrD+eXffbmZ3AHnu/qOGc1VUVNSn61UvNTU1dO/ePS3HTrUoZ7120VY279O7xaPilOO7Mf28wpy5v0K0/77ao7i4OPTvU9kG/NDMqoGewIFmtnk1+HcVMCr4eZ27HwYws48bHWuGmVUBnwfeaHyQoqKi1CZvJJfmqKOctTTE25Xu65qqfgpE+z4QtlzKm+msmeqpzADudvdxwHs0/61hA4N/vwb8Kfi5uWHULGCCu48HtrdwLJGcNHtcCf0LC+ial0f/woJO9VMAPvroI/bu3ZvidCJty9RIZQHwazPbC3xIsq/S1HgzmwQcAr4LfKmVY600s0PALuBzacgrEoq+vfI7PTJp7JJLLtFnf0koMlJU3H0ayUZ90+WNX058u7u/3+j314L/PrWtu09KT0qR6LjlllvYunVr2DEkhvTmR5EIuvDCC/XZXxKKrCgq7j4k7AwiUbJz504qKytzppks0ZEVRUVEUmv06NHqqUgoVFREImjy5Mls2bIl7BgSQyoqIhF0/vnnq6ciodBnf4lE0NatW9mxY0fYMSSGNFIRiaDvfve76qlIKFRURCLozjvvZPPmzWHHkBhSURGJoGHDhqmnIqFQT0UkgsrLy/WOegmFRioiEXT55ZeTSCQYMWJE2FEkZlRURCLonnvuUU9FQqGiIhJBgwcPVk9FQqGeikgEuTsffPBB2DEkhjRSEYmga665hkQiwciRI8OOIjGjoiISQQ888AAVFRVhx5AYUlERiaCzzz6bE044IewYEkPqqYhE0Lp169iwYUPYMSSGNFIRiaAbb7yRRCLBRRddFHYUiRkVFZEImjp1ql79JaFQURGJoJKSEgoKCsKOITGknopIBL3zzjt686OEQiMVkQiaOHEiiUSCiy++OOwoEjMqKiIRNH36dMrLy8OOITGkoiISQWeccQZHH3102DHabcfBWm6aVkp55SH6FR7D7HEl9O2VH3Ys6QQVFZEIWrt2LR988AHFxcWd2n/C3DWs8MoUp2qfDburGDR1RUbPOdQKmTvhzIyeM6pUVEQ6acRDpazfVRV2jFYcw21/eDHsEDlhhVdSNLmj1yqXpheTWU87uYClNw9O65naXVTMbB6w0N1fbse2fYEvu/viNrYrBN4A/t7da5qsew3IB8YD/YG7gI+BOe4+y8x+CwwDjm+6r0gmpPuP80isW7eO8vLynHnz46AfL+XDA7XU1UOXPOhfWMCySdl7fcvKyjo9Csy0TGdN10uKvw58rbUNzOw8YCnQu5XNLgM2AA8BI4DBwNVmdrK7XwDsTE1ckWgZMGAAp556atgx2mXLngQf19VTV5/8ve+J+cweVxJuKOm0jk5/XW9m3w/2u8LdN5rZ94AxQD2wEPgFMBnIN7NVwH7gbpIFrAAY4+7rgTqSI4232jhnMbDR3fcCmNnrwCDg1x3MLhIbq1atoqKiosPPUMPspTSo2JPIaE9F/ZTU6mhRWeXuU8xsJPATM7sTuBQ4J1i/DPgdMAU43d1fMLPrgbHuvt3M7gD+CfiRuy8DMLO2znksycLU4CBwXEsbV1RUUF1d3cGb1T41NTU584YyZU2PxlmvXbSVzftqQ07UmuO44131VNrSuX4K5EpPpc+xn+HxFB+ztScrHS0qvw/+XQVMBQYApwDLg+UnAE3H3NuAGWZWBXyeZA+lIw4APRv93hPY19LGRUVFHTx8+2keNT1yNWtpFmd2dzZt2pQTX9I1fFopG3dXUU9u9FMgd++zmdDRnkrDGPFcYB3gwJ+Aoe4+BJgHvEtyaqvh2LOACe4+HtgO5HXwnGXAqWZ2opkdRXLqa3UHjyESK2bGF77whbBjtMvscSX0Oa4bXfPy6F9YoH5KjuvoSOUsM3uVZP/kcnffbGbLgdfN7GhgDcmRyXvAD8zsD8ACYKWZHQJ2AZ/ryAndvdbMJpGcVutC8tVf2zqYWyRWSktL2bx5c048m+7bK5/HvtUnJ7JK29pdVIKRRnPLp5KcCmvsbaChWbKwjeMWtePci4FWX54sIn919913k0gkuOyyy8KOIjGT7W9+fNLMxrv7fzZdEbxPpbWXI4vE1pw5c9i4cWPYMSSGsraoBD2a1tZfkKEoIjmnX79+HD58OOwYEkP6PhWRCHrllVdYtWpV2DEkhlRURCLo/vvv57HHHgs7hsRQ1k5/iUjnPfXUU2zYsCHsGBJDKioiEdSnTx+qqrL5E5QlqjT9JRJBL7/8MitXrgw7hsSQiopIBE2ZMoUnnngi7BgSQ5r+EomghQsXsn79+rBjSAypqIhEUO/evdm7d2/YMSSGNP0lEkGLFy9mxYrMfs+7CKioiETST3/6U+bNmxd2DIkhTX+JRNAzzzyjnoqEQiMVkQg66aSTOOGEE8KOITGkoiISQc899xzLli0LO4bEkIqKSATNmDGDBQsWhB1DYkg9FZEIWrRoEe4edgyJIY1URCLouOOOo2fPnmHHkBjSSEUkgp5++mm2bdum732XjNNIRSSCZs6cycKFC8OOITGkkYpIBC1ZsoT3338/7BgSQxqpiERQfn4+PXr0CDuGxJBGKiIRtGDBAvVUJBQaqYhE0BNPPMGzzz4bdgyJIY1URCJo2bJllJWVZeRcW/YkuGL+WsorD9Gv8Bhmjyuhb6/8jJxbso+KikjETJi7hhVemfzll1szeu4Nu6sYNLWzH7lf3unzDrVC5k44s9P7S+qoqEhsjHiolPW7qlJwpM4/+El6rPBKiia/mOGzZv5+cNrJBSy9eXDGz9sRWVtUzGw8cC8w3d2nmdlEoLe7TzazS4D7gefdfXKYOSV3pOKPsaysLCea30OGDCGRSLBmzZq0n2v4tFI2VVZRVw9d8qB/YQHLJnXsWufKdW2Qa3kzKdsb9b8EZprZvwM3NCx092eAKaGlEslyr732GvPnz8/IuWaPK6F/YQFd8qBrlzw27a5i+LRStuxJZOT8kl2ydqTSSHdgPrAMOD3kLCJZ71M9lQxP0dR9Ug90treinkoUZH1Rcfe9wNJgOqxNFRUVVFdXpyVLTU1Nxl5Rc6SUtX2uXbSVzftqO7iXeirZJqo9lVOO78ajo/oc0THS8ffV2tRf1heVjioqKkrbsXNpHlVZ26e0g+fNles6bNgwDh06xOrVqzN2ziPpreTKdW2QS3kznTXbeyoi0gmvvPIKc+bMyeg5G3orXfPy6F9YwOxxJRk9v2SHyI1URCQcfXvld/hVXxI9OVNU3H1e2BlEcsUjjzzCzp07uffee8OOIjGT7dNfY8xsUtOFwftU9P4UkRYsXryYFSs6+852kc7L2pFKMDKZ18K6Z4BnMplHJJe89NJLOfPqP4mWbB+piIhIDlFREYmghx9+mKeeeirsGBJDWTv9JSKdt3z5cqqqUvHhmSIdo6IiEkEvvPCCeioSCk1/iYhIyqioiETQgw8+mPF31IuApr9EImn16tUcPHgw7BgSQyoqIhH07LPPqqciodD0l4iIpIyKikgETZkyhVmzZoUdQ2JI018iEfTOO+9w4MCBsGNIDKmoiETQwoUL1VORUGj6S0REUkZFRSSC7rvvPmbOnBl2DIkhTX+JRJC7s3///rBjSAypqIhE0IIFC9RTkVBo+ktERFJGRUUkgu666y5+9rOfhR1DYkjTXyIRtHXrVvbt2xd2DIkhFRWRCJo7d656KhIKTX+JiEjKqKiIRNDtt9/OtGnTwo4hMaTpL5EstmVPgivmr6W88hD9Co9h9rgS+vbKb3O/PXv26H0qEgoVFZEOuOuVHaydXx7KuTfsrmLQ1BXt2/jEUXAiFE1+MSXnHmqFzJ1wZkqOJdGmoiJHZMRDpazfVZWCI4XzQC3ts8IrU1agWpbZ+8BpJxew9ObBGT1nHLS7qJjZPGChu7/cjm37Al9298WtbHMzMDr4dYm739Nk/WtAPjAe6A/cBXwMzHH3WWb2W2AYcLy717T3dkhqpeKPsqysjOLi4hSkSb9MZx0+rZRNlVXU1UOXPOhfWMCySW1f81tvvZU9e/Ywd+7cDKQ8crl0H5DWpatR/3Xgay2tNLN+wP8BzgbOAkaY2d83s+llwAbgIWAEMBi42sxOdvcLgJ2pDi6STWaPK6F/YQFd8qBrlzw27a5i+LRStuxJtLpfdXU1hw8fzlBKkb/q6PTX9Wb2/WC/K9x9o5l9DxgD1AMLgV8Ak4F8M1sF7AfuJlnACoJtPwDOd/dPAMysG9DSaKMY2Ojue4NtXwcGAb/uYHaRIzJh7hpWeCVhTdXVfVIPtLO30nMk9ExdT6WBeivSlo4WlVXuPsXMRgI/MbM7gUuBc4L1y4DfAVOA0939BTO7Hhjr7tvN7A7gn9z9R8BHZpYHTAXedvf1LZzzWJKFqcFB4LiWAlZUVFBdXd3Bm9U+NTU1OfOGsrCyXrtoK5v31XZiT/VUckF6eyvZcR845fhuPDqqT6vbxP2xoLWpyo4Wld8H/64iWQwGAKcAy4PlJwCnNtlnGzDDzKqAzwNvAJhZd2AOySJxfSvnPAD0bPR7T6DFz58oKipqx83onFya9w0ra2knzqnr2raO9lYmTpzIf/3Xf/Hkk09mMGXn5dJ9AHIrb6azdrSn0jDuPRdYBzjwJ2Couw8B5gHvAnWNjj0LmODu44HtQF4wQlkE/NHdr2mYBmtBGXCqmZ1oZkeRnPpa3cHcIjmtobfSNS+P/oUFzB5XEnYkkWZ1dKRylpm9SrJ/crm7bzaz5cDrZnY0sIbkyOQ94Adm9gdgAbDSzA4Bu4DPAd8i2XQ/2sy+GRz7dnf/m2Lh7rVmNonktFoXkq/+2tbhWyqSw/r2ym/Xq74aTJ8+PWemZyRa2l1UgpFGc8unkpwKa+xtwIKfF7ZwyO4dOPdioMWXJ4uISHbI9s/+etLMvtjciuB9Kr0znEckJ9xwww3cd999YceQGMrad9QHPZrW1l+QoSgiOadHjx4kEq2/l0UkHbK2qIhI5z344IPqqUgosn36S0REcohGKiIRdPXVV7Nv3z5+9atfhR1FYkZFRSSCevXqRX19fdgxJIZUVEQi6Mc//rF6KhIK9VRERCRlNFIRiaAJEyawb98+fvOb34QdRWJGRUUkgvr06UOPHj3CjiExpKIiEkH33nuveioSCvVUREQkZTRSEYmgsWPHsn//fhYv1uewSmapqIhEkJlRWVkZdgyJIRUVkQj64Q9/qJ6KhEI9FRERSRmNVEQiaPTo0Rw4cIAlS5aEHUViRkVFJILOOOMMdu/eHXYMiSEVFZEImjx5snoqEgr1VEREJGU0UhGJoO985zscPHiQpUuXhh1FYkZFRSSCvvrVr7Jr166wY0gMqaiIRNCtt96qnoqEQj0VERFJGY1UJOts2ZPgivlrKa88RL/CY5g9roS+vfLDjpVTLrroIqqqqnj11VfDjiIxo6IizZowdw0rPJOfHVXe7NINu6sYNHVFBnN82lArZO6EM0M7f2d94xvfUE9FQpFXX18fdoZUS8sNGvFQKet3VaXj0CJZ77STC1h68+C0Hb+srIzi4uK0HT/VcilvmrLmtbQia0cqZjYeuBeY7u7TzGwi0NvdJ5vZJcD9wPPuPjkTeZbePFh3pDRpmnX4tFI2VVZRVw9d8qB/YQHLJqXvAa0jcvm6imRC1haVwC+BmWb278CZwLMA7v6MmRUAp4cZTlJvy54EtZ/UUReMN/uemM/scSXhhspB3/zmN6mqqmLlypVhR5GYyfaiAtAdmA8sQ0UkY7Klp1KxJxFKTyVXeykNLrzwQnbu3Bl2DImhrO2pBNNfpzdMb7X1e4OKior66urqTp/32kVb2byvttP7i+SCU47vxqOj+oQd4y9qamro3r172DHaLZfypiNrcXFx7vVUOquoqOiI9i9tZQ46l+aoczVrNvdTIHeva7bLpayQW3kznVVvfpSsMntcCf0LC+ial0f/wgL1Uzpp2LBhXH755WHHkBiK3EhFclvfXvlZNTLJVZdeeik7duwIO4bEUM4UFXefF3YGkVxx1VVX6bO/JBTZPv01xswmNV0YvE8lI+9PERGR9svakUowMpnXwrpngGcymUcklwwZMoREIsGaNWvCjiIxk7VFRUQ6b/z48Wzfvj3sGBJDKioiETR+/Hj1VCQU2d5TEZFOqK2tpbZWb+KVzNNIRSSChg8frp6KhEJFRSSCrrzySrZt2xZ2DIkhFRWRCBo7dqx6KhIK9VREIiiRSHAkH6wq0lkaqYhE0MiRI9VTkVCoqIhE0HXXXaeeioRCRUUkgi699FL1VCQU6qmIRND+/fs5ePBg2DEkhjRSEYmgUaNGqacioVBREYmgm266iQ8//DDsGBJDKioiEXTxxRerpyKhUE9FJII++ugj9u7dG3YMiSGNVEQi6JJLLlFPRUKhoiISQbfccgtbt24NO4bEkIqKSARdeOGF6qlIKNRTEYmgnTt3UllZGXYMiSGNVEQiaPTo0eqpSChUVEQiaPLkyWzZsiXsGBJDKioiEXT++eerpyKhUE9FJIK2bt3Kjh07Or3/lj0Jhk8rpf/tSxg+rZQtexIpTCdRppGKSJaYMHcNKzzFzfWlLx7xITbsrmLQ1BUpCNOW8ja3GGqFzJ1wZgaySGepqEjWGPFQKet3VYUdox3afvCT9FjhlRRNPvJCmRrZfT847eQClt48OOPnzdqiYmbjgXuBJ4BBJLPmAVcDXwLuB55398lhZZTUCuMPoKPKysooLi4OO0a7HEnW4dNK2VRZRV09dMmD/oUFLJuUvv8/uXRdIffyZlK291R+CfQHfu7uQ4AHgB+7+zPAlDCDiWSz8vLyI3pH/X2jBtC1Sx4AXbvkcd+oAamKJhGXtSOVRm4B9gc/fwaoCTGLSMqlpZfS4NUjnyqq/aSe0bPeTEGYtrQ+naR+Sm7I+qLi7h8BmJkBDwLfam37iooKqqur05KlpqYmZ16mmStZr120lc37asn2+elPy6Ws0ZFd/RTIlvvBKcd349FRfVpcn47Hgtam/rK+qACY2VDgEeC77u6tbVtUVJS2HLk0j5orWUuLi3MmK+TOdQX1VNIpl/JmOmu291QaCsrDwPnu/h9h5xHJBe7OBx980On9Z48roX9hAV3z8uhfWMDscSUpTCdRlgsjlenAUcD85AwY7u7XhBtJJLtdc801JBIJRo4c2an9+/bKT+vIRKIr64uKu3857AwiueaBBx6goqIi7BgSQ9k+/TXGzCY1XWhmlwB6f4pIC84++2y+8pWvhB1DYihrRyruPg+Y18K6Z4BnMplHJJesW7eO8vLynGkmS3RkbVERkc678cYbSSQSXHTRRWFHkZhRURGJoKlTpx7Rq79EOktFRSSCSkpKKCgoCDuGxFC2N+pFpBPeeeednPhEBYkejVREImjixIkkEgkuvvjisKNIzKioiETQ9OnTKS/Pjs+mknhRURGJoDPOOIOjjz467BgSQ+qpiETQ2rVree+998KOITGkkYpIBH3/+98nkUjwz//8z2FHkZhRURGJoJ///OfqqUgoVFREImjAgAF07do17BgSQ+qpiETQqlWrePvtt8OOITGkkYpIBN1xxx0kEgnGjBkTdhSJGRUVkQh67LHH2LRpU9gxJIZUVEQiyMyoq6sLO4bEkHoqIhFUWlrK2rVrw44hMaSRikgE3X333SQSCS677LKwo0jMqKiIRNCcOXPYuHFj2DEkhlRURCKoX79+HD58OOwYEkPqqYhE0CuvvMKqVavCjiExpKIiEkH3338/jz32WNgxJIY0/SUSQU899RQbNmwIO4bEkIqKSAT16dOHqqqqNrfbsifBFfPXUl55iH6FxzB7XAl9e+VnIKFElYqKSBpNmLuGFV4ZYoL2f1Lxht1VDJq6Io1Z2tJ61qFWyNwJZ2Yoi3SWioo0a8RDpazf1fYz3dTJpY9pz6Ws0bHCKyma/GLYMRrJpfvBp7OednIBS28enJYzpb2omNk8YKG7v9yObfsCX3b3xa1scwMwHqgHHnT3X6UoqjSSrjtcc8rKyiguLs7Y+Y5ErmTduXMn69evZ9CgQa1uN3xaKZsqq6irhy550L+wgGWTMvf/vkGuXNcGuZQ301mz7dVfXwe+1tJKMzsJuA44G/gG8FMzy8tQNpGc0bt3bwoLC9vcbva4EvoXFtA1L4++J+ZT+0kd/W9fwvBppWzZk8hAUomaTE1/XW9m3w/Od4W7bzSz7wFjSI44FgK/ACYD+Wa2CtgP3E2y8BUAY9x9vZmd4e4fm1kRUOPu9Rm6DSJHJJz+SvunaCoaFZFw+isdm05SjyU75dXXp/cxOZj+et/dp5jZSOBK4E7gcaBhnL2M5Ajkq8Dp7j7ZzK4Hnnf37WZ2B5Dn7j8KjnkjcA8ww93vaXy+ioqK+urq6rTclpqaGrp3756WY6daurJeu2grm/fVpvy4IpIZpxzfjUdH9TmiYxQXF7c4Q5Spkcrvg39XAVOBAcApwPJg+QnAqU322QbMMLMq4PPAGw0r3P3nZvY48JKZDXX3vzylKioqSssNAM2jApSm4Zi6rqn30UcfsX79es4+++zlb49VAAAFu0lEQVR27xNmfyVXrmuDXMob1Z5Kwxj1XGAd4MCfgKHuPgSYB7wL1DXKNAuY4O7jge1AniU9F/RRaoHDwT4i0shJJ53ECSec0KF9GvdX+hcWMHtcSZrSSZRlaqRylpm9SrJ/crm7bzaz5cDrZnY0sIbkyOQ94Adm9gdgAbDSzA4Bu4DPubub2R+B1cGxXnL30gzdBpGc8dxzz/Hhhx926Blq3175obzyS6Il7UUlGGk0t3wqyamwxt4GLPh5YQv73UOynyIiLZgxYwaJRIKbbrop7CgSM3rzo0gELVq0CHcPO4bEULa9T0VEUuC4446jZ8+eYceQGNJIRSSCnn76abZt25Yzr1CS6NBIRSSCZs6cycKFzbYlRdJKIxWRCFqyZAnvv/9+2DEkhjRSEYmg/Px8evToEXYMiSGNVEQiaMGCBeqpSCg0UhGJoCeeeIJnn3027BgSQxqpiETQsmXLKCsrCzuGxJBGKiIR1K1bN7p16xZ2DIkhjVREImjevHls375dPRXJOI1URCJo3rx5PP/882HHkBhK+5d0Zdpbb70VrRskIpKFBg4c2OwXdUWuqIiISHg0/SUiIimjoiIiIimjoiIiIimjlxS3wMyOI/mVxscCRwGT3H11k22uAq4BPgbud/ffZjzop/N8G/gndx/TzLqHgXOAg8GiUe6+P5P5mmojb1ZcWzPrQfJ+8FmS126cu1c22WYRcBJQC1S7+zcznLEL8AjwZeAwcKW7b2y0PiuuZZClrazZeD/9X8C/ufuQJssvBO4ieV3nuPusEOJ9SitZbwauBBruu9d4mr7FTUWlZZOA5e4+3cwM+H/APzSsNLPewE3A/wS6A6+b2TJ3PxxG2OCP8TzgnRY2GQic5+4fZS5Vy1rLm2XX9jrgPXf/VzMbDdwJ/EuTbU4F/oe7h/Wql28B3d39q2Z2FvBTYBRk3bVsNWsg2+6ntwHfBQ41Wd4NeAgoCda9YWYvuPuuzKf8S6ZmswYGApe5+1vpzqHpr5Y9BDwW/PwZoKbJ+jOBN9z9cPBMaiPw9xnM19Qqkg+AfyN4dngq8LiZvWFml2c0WfNazEt2XdtzgJeDn18ChjVeaWYnA8cDi83sdTO7IMP5oFFGd3+TZAFpkE3XElrJmqX3003Axc0sLwY2uvted/8z8DowKKPJ/lZLWSFZVG4P7qO3pzOERiqAmV0B3Nxk8QR3Xxs801sATGyy/lig8bD8IHBc+lImtZL1aTMb0sJuxwA/A6YBXYEVZvYf7v5u+pImdTJvNl3bXY2yNJfjKJLPth8GTiT5jHWNu+9OZ9Ymml6vT8zsM+7+cTPrMnItW9Fa1tDupy1x92fNrKiZVdl2XVvLCrAQ+AVwAPiNmV2QrmlQFRXA3WcDs5suN7Mvkfyfcau7lzZZfQBo/CXgPYF9aQsZaClrGxLAw+6eADCzV0nOaaf9j7WTebPm2prZc42yNJdjJ/Bo8KC428zeBgzIZFFper26BHmaW5eRa9mK1rKGdj/thGy7ri0yszxgekNvysxeBL4CpKWoaPqrBWb2ReDXwBh3f6mZTdYA55pZ96CpXwysy2TGDjiN5DPorsFc8DnAH0LO1JpsurZvACODn78JrGyyfhjJ+wlmVgAMADL98cB/yRj0Kd5rtC6briW0njWX7qdlwKlmdqKZHUVy6mt1G/uE5VhgnZkVBAXm60DaeisaqbTsxyQbmw8n+/Tsd/dRZjaJ5FzqC2Y2g+SDTBfgB+7etO8SqiZZnwLeJPkKpSfd/U/hpvtbWXptZwLzzex14M/AmCDrT4Bn3P0lMzvPzN4E6oA7Qmgy/wYYbmargDxgQpZey/Zkzer7qZmNAQrc/fEg9+9IXtc57r4t3HSf1iTrHcAKkq+4W+7uS9J1Xn1Mi4iIpIymv0REJGVUVEREJGVUVEREJGVUVEREJGVUVEREJGVUVEREJGVUVEREJGX+P2yV1VpNkEJ0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_03.trace,\n",
    "              varnames=['alpha', 'beta1', 'beta2', 'beta3'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "INFO (theano.gof.compilelock): Refreshing lock /home/petrushev/.theano/compiledir_Linux-5.0-arch1-1-ARCH-x86_64-with-arch--3.7.3-64/lock_dir/lock\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, beta2, beta1, alpha]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:24<00:00, 491.45draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_04:\n",
    "    alpha = pm.Normal('alpha', mu=0, sd=1)\n",
    "    beta1 = pm.Normal('beta1', mu=0, sd=1, shape=2)\n",
    "    beta2 = pm.Normal('beta2', mu=0, sd=1, shape=3)\n",
    "\n",
    "    score_i = pm.Deterministic(\n",
    "        'score',\n",
    "        wines.score.mean()\n",
    "        + alpha\n",
    "        + beta1[0] * wines.judge_amer.values\n",
    "        + beta1[1] * wines.flight_red.values\n",
    "        + beta2[0] * (wines.wine_amer.values * wines.judge_amer.values)\n",
    "        + beta2[1] * (wines.wine_amer.values * wines.flight_red.values)\n",
    "        + beta2[2] * (wines.judge_amer.values * wines.flight_red.values)\n",
    "    )\n",
    "\n",
    "    sigma = pm.HalfCauchy('sigma', beta=5)\n",
    "    score_obs = pm.Normal('score_obs', mu=score_i, sd=sigma,\n",
    "                          observed=wines.score.values)\n",
    "    \n",
    "    m_04.name = 'second level interaction'\n",
    "    m_04.trace = pm.sample(5000, tune=1000, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEcCAYAAAAP5CkrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XF8VPWZ7/FPQEQwqBRT3faKqVQf6aVbd7lxrVWBFqhlVVraLVwuCoiC0laRq71oXV1X17JFEamKikC0rMWKVqQiFRWzKPqCdbXFbnhAMBBBIGUBDQMWJfePmdSYJiETJvObc873/Xr5Iplz5pzvHCd58vyemaSorq4OERGRXOgQOoCIiMSHioqIiOSMioqIiOSMioqIiOSMioqIiOSMioqIiOTMEaEDiLSVmXUG5gGnAO8DP3D39Wb2HeAOoDqz683A68DTQBdggrv/3szOAb7m7v/azPHHA6OAg0An4Cfu/lIbcl4BnAjcD9zk7hPNrAo43d33N9hvTOa2Ka04ZjmwwN2XNrO9J/AVd1+cbd5DnLc/cIW7j8jlcSU+1KlIlF0O1Lr7WcCPgHsyt/cFfuzu/TP/VQCDSReVicA4MysCrgbuburAZjYCGAR8w937ky4uvzCz49sa1t23ufvEtt4/S18Hvpanc4n8mToVibIvAc8CuLubWe/M7X2BvzGzScAq4P8BtaS7lC7AXmAk8OuGnUIjE4DJ7n4gc/x3zOwMd9+Z6RJ6ZP77e+DHwLlAR2C6uz+e6YLuBnYBHwGvmVkp6e7irMw5Hsjcth0Y3fDkZvajTMa6zH1mNhUy0zn8P+BPpDu2BcBUYArQ1cxWAu8AM4EiYCdwKfA3wL9m7vc80N/dB2SO+RvgH4FewA9Id2l1wHcanXse8MXMNb3b3X/RzLWUBFGnIlH2JnCBmRWZ2VnA582sI7CMdOdyHlAMXEH6G+cJpIvFg6S/Qf7OzB4wsx83cezPARsb3uDuOxt8+qK7nw2cBXzB3c8BBgA/MbPjgFnA/3b3gaS/qTdllrv3A6pId10AmNmXgOHAOaSL1bfNzFq4DicD381k+bG7f0y6sDzq7k8Ds0kvDfYHlpAuggBHufu57n4LcJSZnWxmfwUc7+5vAKcBf595bP8FfLNBxm6kr+8w4Hzg4xbySYKoU5Eomwv0BlYArwCvu/vHZjbX3XcDmNki4LvufhCYlLntetJdxI2ki88/mdlp7r6uwbE3AScBe+pvMLNvAr/PfOqZf78M9DWzlzKfdwJKgRMaHO8V0j/RN/Qnd38t8/FK0kttqzOf9yFdKF7IfN4dOLXBORtb4+4fAR+Z2b4mtvcG7svUpU7A+kaPAWAOcAnwIek5FcAO4GEzqwVOB16t39ndP8h0gg8CxwDzm8kmCaNORaKsDHgh85P048DGzKzk92b2PzL7fIP0kB4AM/ssYO6+AuhK+ifsOuDoRseeC/yjmR2Rud9pwEN88hP5wcy/a4HlmS7g68CvgA3AlgbLcWVNZD/SzM7IfHwu8FaDbQ78ARiQOW45nxSzpjT1C/wO8snXtwOXZI71Y+A3jR4DpJfNLiDdwT1qZscCtwAjgMuAfaSXzwDIdDR93f07pJcAf1Z/rSTZVFQkytYDk8zsVeBW0jOQOtLfBJ80swrShWN2g/vcCNyW+fg+4LfAXwG/a3hgd18AvAa8bGb/Tvqn91HuvqNRhsVArZmtIF286tz9A9LLbI+Y2Quku47GPgR+lDn2Z0kXrPpz/450l/Kymf0H6S5lS+svCwBrgKGZFxxcmcnyMullsb8oUO5em7kGlZn875PusF4l3QnuI70kWG8bcGJmZrMMuCPTLUnCFem3FIuISK6oUxERkZxRURERkZxRURERkZxRURERkZxRURERkZyJ3evKX3/9db2cTUSknfXt27eoqdvj+JLirB5QZWUlvXv3PvSOBSSKmSGauaOYefbs2bz33nvcdNNNoaNkJYrXOoqZISe5mywoEMNORSTpLr/8ciorK0PHkITSTEVERHJGnYpIzPTv359UKsWqVatCR5EEUlERiZkxY8awdevW0DEkoVRURGJmzJgxmqlIMJqpiMTMgQMHOHDgQOgYklDqVERiZtCgQZqpSDAqKiIxc9lll7FlS7Z/fkUkN1RURGJm1KhRmqlIMJqpiMRMKpVi376m/lS9SPtTpyISM0OGDNFMRYJRURGJmSuvvFIzFQlGRUUksM07U4x7eDUba/ZySsnRzBldRs8eXdt8vOHDh2umIsGoqIi0wdh5q1juNTk/7vodtZw3bXmOjrYRgAFWwryxZ+bomCItU1GRgjX4rgrWba8NHYP6b85RtdxrKJ3yTOgYzTrthGKeu6Zf6BiSI0GLipltc/cTm9nWH7jC3UfkN5UUikL4RpOPv5cxaHoFG2pqOVgHHYqgV0kxyya3/bE/+eSTvPvuu1x11VU5TCnSOnpJsUhgtw7tQ8cO6b951LFDEbcO7XNYxxs2bBiDBg3KRTSRrOWlUzGzPsB0oCNwPHClu69ssP0lYC1wOum/KDY8s+lUM3sW+Cyw2N3/ycz6ATeTLojFwEh3X5ePxyFSr71mKgc+rmPE7NdydLSNmqdI3uXlzwmb2XDgv9x9jZmNBAa4++X1y1+ZojLX3R8xs4mki8uTwL3A35IuRpvd/fjM9qfcfauZ3QAUufu/1J+rqqqqLps3fu3fv5+jjjoqZ481H6KYGbLLfcWiajbt1i9FjLuTj+vE/UNPAqL5vI5iZjj83L179w7+54S3AP9oZvuAbsD7TezzYubflcDQzMdvufuHAGb2UYNjzTSzWuDzwCsND1JaWppVsCj+jekoZobsclcUyOOL4kxl8eLFVFdXM3HixBymbH9RfF5HMTO0b+58zVRmAje7+2hgDeklrsb6Zv79GvCHzMdNtVGzgbHuPgbY2syxRCJjzugyepUU07GoiF4lxcwZXXZYx7vwwgsZMGBAjtKJZCdfncp84HEz2wW8S3qu0tgYM5sM7AUuBr7cwrFWmNleYDvwuXbIK5I3PXt0PazOpLFt27ZRU1MTyZ+gJfryUlTcfTrpQX3j2xu+nPh6d1/b4POXMv99al93n9w+KUXiYcSIEfrdXxKM3vwoEjNTpkxh8+bNoWNIQhVEUXH3/qEziMTF+eefr9/9JcHozY8iMVNdXc17770XOoYkVEF0KiKSOxdffLFmKhKMiopIzNx4441s2rQpdAxJKBUVkZgZOHCgZioSjGYqIjGzceNGqqurQ8eQhFKnIhIzl156KalUisGDB4eOIgmkoiISM7fccotmKhKMiopIzPTr108zFQlGMxWRmHF33nnnndAxJKHUqYjEzIQJE0ilUgwZMiR0FEkgFRWRmLn99tupqqoKHUMSSkVFJGbOPvtsunfvHjqGJJRmKiIx89Zbb7F+/frQMSSh1KmIxMwPf/hDUqkUF110UegokkAqKiIxM23aNL36S4JRURGJmbKyMoqLi0PHkITSTEUkZt588029+VGCUaciEjOTJk0ilUoxbNiw0FEkgVRURGJmxowZbNy4MXQMSSgVFZGYOeOMM+jcuXPoGJJQmqmIxMzq1atZs2ZN6BiSUOpURGLmuuuuI5VK8f3vfz90FEkgFRWRmLnnnns0U5FgVFREYqZPnz507NgxdAxJKM1URGJm5cqVvPHGG6FjSEKpUxGJmRtuuIFUKsXIkSNDR5EEUlERiZkHHniADRs2hI4hCaWiIhIzZsbBgwdDx5CEUlERiZmKigo2bdpE7969Q0c5pM07U4x7eDUba/by+WOOYP5nT6Znj66hY8lhUFERiZmbb76ZVCrFJZdc0qr9x85bxXKvaedUh7Z5zwHOm7Y8yLkHWAnzxp4Z5Nxxo6IikoXBd1Wwbntt6BgtO+s6AEqnPBM4SHQs95rDuF7ReE/QaScU89w1/dr9PK0uKmZWDixw96Wt2Lcn8BV3X3yI/UqAV4C/dvf9jba9BHQFxgC9gJuAj4C57j7bzH4DDASOa3xfkfaSjy/KXKisrIzE8teg6RVsqKnlYB0UAV/8bDHLJkfjGkN0rnM+tdf7VL4OfK2lHczsm8BzwIkt7HYJsB64CxgM9APGm9kJ7n4BsC03cUXi4/nnn2flypWhY7TKnNFl9CoppkMRdOwAG3bUMmh6BZt3pkJHkzbKdvlropldl7nfOHd/28x+BIwE6oAFwL3AFKCrma0E9gA3ky5gxcBId18HHCTdabx+iHP2Bt52910AZvYycB7weJbZRRLhtttuI5VKMW7cuKzvG3K+crAu/e/6HbVBZiuaq+RGtkVlpbtPNbMhwM/M7EZgOHBOZvsy4LfAVOB0d3/azCYCo9x9q5ndAPwD8C/uvgzSL388hGNIF6Z6HwDHNrdzVVUV+/bta/UD2r9/f+T+Sl4UM0M0c9dnvmJRNZt2Hwgdp3U0U2mTts9VojFTaeikY47gwcO4f0tLftkWlX/P/LsSmAb0AU4GXsjc3h04tdF9tgAzzawW+DzpGUo23ge6Nfi8G7C7uZ1LS0uzOngU10SjmBmimbs+c0VEc0fFoOkVvL2jljqgQxH0KonGbCVq17lee+bOdqZS3xueC7wFOPAHYIC79wfKgd+TXtqqP/ZsYKy7jwG2kp7HZaMSONXMPmNmR5Je+no1y2OIJMbSpUtZsWJF6BhZmTO6jJOO7UTHoiJ6lRQzZ3RZ6EjSRtl2KmeZ2Yuk5yeXuvsmM3sBeNnMOgOrSHcma4CfmNl/AvOBFWa2F9gOfC6bE7r7ATObTHpZrQPpV39tyTK3SGJMnTqVVCrF+PHjQ0dptZ49uvLAt0+K5E/98mmtLiqZTqOp26eRXgpr6A2gfliy4BDHLW3FuRcDLb48WUTSFixYwLp160LHkIQq9F99/4iZfampDZn3qbT0cmSRRDrxxBMpKSkJHUMSqmDfUZ+Z0bS0/YI8RRGJlMWLF1NdXa2lJAmi0DsVEcnSnXfeSXl5eegYklAF26mISNssXLhQMxUJRp2KSMwcf/zxdO/ePXQMSSgVFZGYefLJJ1m2bFnoGJJQKioiMTNz5kzmz58fOoYklGYqIjGzaNEi3D10DEkodSoiMXPsscfSrVu3Q+8o0g7UqYjEzGOPPcaWLVv0PhUJQp2KSMzMmjWLBQta/O1IIu1GnYpIzCxZsoS1a9eGjiEJpU5FJGa6du1Kly5dQseQhFKnIhIz8+fP10xFglGnIhIzDz30EE888UToGJJQ6lREYmbZsmVUVlaGjiEJpU5FJGY6depEp06dQseQhFKnIhIz5eXlbN26VTMVCUKdikjMlJeX89RTT4WOIQmlTkUkZl566SXNVCQYdSoiIpIzKioiMTN79mwef/zx0DEkoVRURGLmscce49lnnw0dQxJKMxWRmHn++ec1U5Fg1KmIiEjOqKiIxMx9993HL3/5y9AxJKG0/CUSM4sXL6a2tjZ0DEkoFRWRmHn22Wc1U5FgtPwlIiI5o6IiEjN33303v/jFL/J6zs07UwyaXkGv65cwaHoFm3em8np+KRxa/hKJuLHzVrHcaxrc8kXgizw65ZkgedbvqOW8acvbeO+Nh3XuAVbCvLFnHtYx5PCoqEgsDL6rgnXb22s4fXjf6CR/lnsNpXkvpmGeH6edUMxz1/QLcu6WqKhILLTXF1dlZWUkf4V8vnMPml7BhppaDtZBhyLoVVLMssnZ/T+J4rWOYub2VrAzFTMbY2abzWxy5vNJZjY18/H3zGxt/eci8ok77riDuXPn5vWcc0aX0aukmI5FRfQqKWbO6LK8nl8KR6F3Ko8Cs8zs34AzgScA3H2hmRUDp4cMJ1KIXn31VT744IO8nW/zzhTjHl7Nxpq9nFJyNHNGl9GzR9e8nV8KS6EXFYCjgIeBZaiIiDTpU8P6Uy8FCDBbONwhPbRlPqHhfGEp+KLi7ruA58xsTGv2r6qqYt++fa0+/v79+yP3RrEoZoZwua9YVM2m3QcO4wga1BeyMMP5hsI+P04+rhP3Dz0pq/sc7tdiS3Okgi8q2SotLc1q/ygO2qKYGcLlrjiMc0bxWk+dOpUdO3Ywffr0vJwvF0N6iOa1jmJmaN/cBTuoF5G2efPNN1m7dm3ezqchvTQUu05FJOkWLFiQ12XGnj26tqkzkXiKTFFx9/LQGUREpGWFvvw1sv59Kg2Z2feAKQHyiBS8W2+9lVmzZoWOIQlVsJ1KpjMpb2bbQmBhPvOIRIW7s2fPntAxJKEKtqiISNvMnz8/ki85l3go9OUvERGJEBUVkZi56aab+PnPfx46hiSUlr9EYqa6uprdu3eHjiEJpaIiEjPz5s3TTEWC0fKXiIjkjIqKSMxcf/31efu9XyKNqaiIxMzOnTv1PhUJRjMVkZh58MEHNVORYNSpiIhIzqhTEYmZa6+9lp07dzJv3rzQUSSB1KmIxMy+ffv48MMPQ8eQhFKnIhIz9957r2YqEow6FRERyRl1KiIxM2nSJP77v/+bRx55JHQUSSB1KiIikjPqVERiZsaMGZqpSDDqVEREJGfUqYjEzA9+8AN27drFo48+GjqKJJA6FZGY6dKlC507dw4dQxJKnYpIzNxxxx2aqUgw6lRERCRn1KmIxMz48ePZvXs3v/rVr0JHkQRSURGJmR49elBXVxc6hiSUiopIzPz0pz/VTEWC0UxFRERyRp2KSMyMHTuW3bt38+tf/zp0FEkgFRWRiNm8M8W4h1ezsWYvp5QczZzRZfTs0fXP20866SS6dOkSMKEkmYqKSJbGzlvFcq8JHQOA9TtqOW/a8ka3/h0cA6VTnmmXcw6wEuaNPbNdji3Rp6IieXPFomo2PbwxdIw2iGLm9rPca9qtYBXStT7thGKeu6Zf6BiR0+qiYmblwAJ3X9qKfXsCX3H3xS3scw0wIvPpEne/pdH2l4CuwBigF3AT8BEw191nm9lvgIHAce6+v7WPQ8K5f+hJ9O7dO3SMrFRWVhZc5kHTK9hQU8vBOuhQBL1Kilk2+ZNvfqNGjWLPnj0sXtzsl19BKsRrLdlrr1d/fR34WnMbzewU4P8AZwNnAYPN7K+b2PUSYD1wFzAY6AeMN7MT3P0CYFuug4sUos07UwyaXkGv65dw4OOD9PxMVzoWFdGrpJg5o8s+ta+Z8YUvfCFQUkm6bJe/JprZdZn7jXP3t83sR8BIoA5YANwLTAG6mtlKYA9wM+kCVpzZ9x3gfHf/GMDMOgHNdRu9gbfdfVdm35eB84DHs8wuctgKYZ5StTP154+bnqn8LXSFp3O4RKU5irRWtkVlpbtPNbMhwM/M7EZgOHBOZvsy4LfAVOB0d3/azCYCo9x9q5ndAPyDu/8L8EczKwKmAW+4+7pmznkM6cJU7wPg2OYCVlVVsW/fvlY/oP3790fujWJRyHzFomo27T7QxJbCWTNvvShmzq32naM0VJjX+uTjOnH/0JP+4vYofC025XBzt7RMmW1R+ffMvytJF4M+wMnAC5nbuwOnNrrPFmCmmdUCnwdeATCzo4C5pIvExBbO+T7QrcHn3YDdze1cWlraiofxiSiu40Yhc0UT+aKQu7FCyXyoOUpDI0aM4P3332fJkiV5Tnl4CuVaZyOKmaF9c2c7U6nvf88F3gIc+AMwwN37A+XA74GDDY49Gxjr7mOArUBRpkNZBPzO3SfUL4M1oxI41cw+Y2ZHkl76ejXL3CKRNmd0Gb1KipudozR0xhlncPrpp+cxncgnsu1UzjKzF0nPTy51901m9gLwspl1BlaR7kzWAD8xs/8E5gMrzGwvsB34HPBt0kP3zmb2rcyxr3f3vygW7n7AzCaTXlbrQPrVX1uyfqQiEdazR9dmO5PGpkyZEsklGYmHVheVTKfR1O3TSC+FNfQGYJmPFzRzyKOyOPdiIFqvjxQRSaBC/4WSj5jZl5rakHmfyol5ziNS8L773e9y9dVXh44hCVWw76jPzGha2n5BnqKIRMpXv/pVtm/fHjqGJFTBFhURaZtrr71WMxUJptCXv0REJELUqYjEzEUXXURtbS0vvvhi6CiSQCoqIjHzjW98QzMVCUZFRSRmrr76as1UJBjNVEREJGfUqYjEzLe+9S1qa2tZsWJF6CiSQCoqIjFz4YUXsm2b/tSQhKGiIhIzEydO1ExFgtFMRUREckadikjMDBw4kL179/Lqq/oLEZJ/KioiMTN8+HDee++90DEkoVRURGLm8ssv10xFgtFMRUREckadikjM9O/fn1QqxapVq0JHkQRSURGJmTFjxrB169bQMSShVFREYmbMmDGaqUgwmqmIxMyBAwc4cOBA6BiSUOpURGJm0KBBmqlIMCoqIjFz2WWXsWXLltAxJKFUVERiZtSoUZqpSDCaqYjETCqVYt++faFjSEKpUxGJmSFDhmimIsGoqIjEzJVXXqmZigSjoiISM8OHD9dMRYLRTEUkZvbs2cMHH3wQOoYklDoVkZgZOnSoZioSjIqKSMxcddVVvPvuu6FjSEKpqEhB2LwzxbiHV7OxZi+nlBzNnNFl9OzRNXSsSBo2bJhmKhKMior82dh5q1juNe18lo2H3GP9jlrOm7a8nXN82gArYd7YM/N6zvbyxz/+kV27doWOIQmlotLI4LsqWLe9NnSMVjj0N2dpveVeQ+mUZ5rZGtFr/XRzjwdOO6GY567pl8cwkhQFW1TMbAzwz8AMd59uZpOAE919ipl9D7gNeMrdp+TyvFH4QqusrKR3796hY2StpdyDplewoaaWg3XQoQh6lRSzbHL4/xdRvNaLFy+murqaiRMnho4iCVToLyl+FJhlZv8G/KD+RndfCEwNlkpybs7oMnp+Jj1DOVgHBz4+yOadqcCpounCCy9kwIABoWNIQhVsp9LAUcDDwDLg9MBZYq/95yqtW0qq2pnK+1wF4jFb2bZtGzU1NZHrsCQeCr6ouPsu4LnMctghVVVVZfXL9Pbv39/sK2WuWFTNpt2F+seOIrrOX+Canq1E9FovaX6mAnDycZ24f+hJeQpzaC19LRaqKGaGw8/d0g8sBV9UslVaWprV/i2tmVcU6E96UVznh0PnLsS5ShSv9dKlS9m8eTPjx48PHSUrUbzWUcwM7Zu70GcqkiBzRpfRq6SYjkVF9CopZs7ostCRIun888/n3HPPDR1DEip2nYpEV88eXYN3JnFQXV3Ne++9F8mfoCX6IlNU3L08dAaRKLj44ov1u78kmEJf/hppZpMb35h5n0pO358iEhc33ngjEyZMCB1DEqpgO5VMZ1LezLaFwMJ85hGJioEDB0byFUkSD4XeqYhIljZu3Eh1dXXoGJJQBdupiEjbXHrppaRSKQYPHhw6iiSQiopIzNxyyy1s2rQpdAxJKBUVkZjp16+fZioSjGYqIjHj7rzzzjuhY0hCqVMRiZkJEyaQSqUYMmRI6CiSQCoqIjFz++23U1VVFTqGJJSKikjMnH322XTv3j10DEkozVREYuatt95i/fr1oWNIQqlTEYmZH/7wh6RSKS666KLQUSSBVFREYmbatGl69ZcEo6IiEjNlZWUUFxeHjiEJpZmKSMy8+eabevOjBKNORSRmJk2aRCqVYtiwYaGjSAKpqIjEzIwZM9i4cWPoGJJQKioiMXPGGWfQuXPn0DEkoTRTEYmZ1atXs2bNmtAxJKHUqYjEzHXXXUcqleL73/9+6CiSQCoqIjFzzz33aKYiwaioiMRMnz596NixY+gYklCaqYjEzMqVK3njjTdCx5CEUqciEjM33HADqVSKkSNHho4iCaSiIhIzDzzwABs2bAgdQxJKRUUkZsyMgwcPho4hCaWZikjMVFRUsHr16tAxJKHUqYjEzM0330wqleKSSy4JHUUSSEVFJGbmzp3L22+/HTqGJJSKikjMnHLKKXz44YehY0hCqaiIxMzzzz/Ppk2b6N27d86OuXlninEPr2ZjzV5OKTmaOaPL6Nmja86OL/GhoiJSwMbOW8Vyr2nDPU/k1inP5DwPwPodtZw3bXm7HBuy+/UyA6yEeWPPbKcs0hYqKlJwBt9VwbrttaFjNKDfo1WolnsNpe1UPFsv/PPjtBOKee6afqFjAAVcVMxsDPDPwEPAeaSzFgHjgS8DtwFPufuUUBmlfRTKFwdAZWVlTpeR8iXXuQdNr2BDTS0H66BDEfQqKWbZ5Nz+f4ritY5i5vZW6O9TeRToBdzj7v2B24GfuvtCYGrIYCKFaunSpaxYsSKnx7x1aB86digCoGOHIm4d2ienx5f4KNhOpYH/C+zJfHwEsD9gFpG8a9tc5fPc3k7LQgc+rmPE7Nfa5ditXUrSLKVwFXxRcfc/ApiZAXcA325p/6qqKvbt29fq4+/fv5/KysrDyphvUcwMbc99xaJqNu0+0A6JWiv8mrl8WmHMUuqFfX6cfFwn7h96Ulb3OdzvIS0t+RV8UQEwswHAfcDF7u4t7VtaWprVsaO4JhrFzND23BUBH2vSrnVzNFNpWhQzQ/vmLvSZSn1BuRs4393/I3QekUK3ePFili/P7Ut+54wuo1dJMR2LiuhVUsyc0WU5Pb7ERxQ6lRnAkcDD6RUw3N0nhI0kUrjuvPNOUqkUEydOzNkxe/bomvPOROKp4IuKu38ldAaRKFm4cCHr1q0LHUMSqtCXv0aa2eTGN5rZ9wC9P0WkCccffzzdu3cPHUMSqmA7FXcvB8qb2bYQWJjPPCJR8eSTT/Luu+9GcoAs0VfonYqIZGnmzJnMnz8/dAxJqILtVESkbRYtWsQhXnkv0m7UqYjEzLHHHku3bt1Cx5CEUqciEjOPPfYYW7Zs0UxFglCnIhIzs2bNYsGCBaFjSEKpUxGJmSVLlrB27drQMSSh1KmIxEzXrl3p0qVL6BiSUOpURGJm/vz5mqlIMOpURGLmoYce4oknnggdQxJKnYpIzCxbtiySf29H4kGdikjMdOrUiU6dOoWOIQmlTkUkZsrLy9m6datmKhKEOhWRmCkvL+epp54KHUMSqqiuri50hpx6/fXX4/WAREQKUN++fYuauj12RUVERMLR8peIiOSMioqIiOSMioqIiORM4l5SbGbHAvOBY4Ajgcnu/mqjfS4HJgAfAbe5+2/yHrQJZvYd4B/cfWQT2+4GzgE+yNw01N335DNfUw6RueCus5kKO21JAAAD0UlEQVR1If38+Czpazna3Wsa7bMIOB44AOxz92/lPWg6RwfgPuArwIfAZe7+doPthXh9D5W5IJ/H9czs74B/dff+jW6/ELiJ9LWe6+6zA8RrUguZrwEuA+qf3xM8B3/dLXFFBZgMvODuM8zMgF8Cf1u/0cxOBK4C/hdwFPCymS1z9w+DpP0k193AN4E3m9mlL/BNd/9j/lK1rKXMhXqdgSuBNe7+T2Y2ArgRuLrRPqcC/9PdQ7/K5dvAUe7+VTM7C7gTGAoFfX2bzZxRcM/jemb2Y+BiYG+j2zsBdwFlmW2vmNnT7r49/yk/rbnMGX2BS9z99VyeM4nLX3cBD2Q+PgLY32j7mcAr7v5h5iekt4G/zmO+5qwk/Q3vL2R++jsVeNDMXjGzS/OarHnNZqZwr/M5wNLMx88CAxtuNLMTgOOAxWb2spldkOd8Df05q7u/RrqA1Cv469s4cwE/j+ttAIY1cXtv4G133+XufwJeBs7La7LmNZcZ0kXl+szz+PpcnTDWnYqZjQOuaXTzWHdfnflJbj4wqdH2Y4CG7fYHwLHtl/LTWsj8mJn1b+ZuRwM/B6YDHYHlZvYf7v779kv6iTZmDnqdodnc2/kkV1OZjiT90/XdwGdI/1S6yt13tGfWZjS+hh+b2RHu/lET2/J+fZvRUuagz+NDcfcnzKy0iU2Feq1bygywALgXeB/4tZldkIsl0lgXFXefA8xpfLuZfZn0Bb3W3SsabX4faPgHvrsBu9stZCPNZT6EFHC3u6cAzOxF0mvWeflibGPmoNcZms5tZk82yNVUpm3A/ZlvgjvM7A3AgBBFpfE17JDJ1dS2vF/fZrSUOejz+DAU6rVulpkVATPq51Vm9gzwN8BhF5XELX+Z2ZeAx4GR7v5sE7usAs41s6MyQ/3ewFv5zNgGp5H+ibljZn33HOA/A2c6lEK9zq8AQzIffwtY0Wj7QNLPH8ysGOgDhPqVwH/OmplPrGmwreCvbxOZo/g8hvT//1PN7DNmdiTppa9XD3Gf0I4B3jKz4kyB+TqQk9lKrDuVZvyU9ODy7vScnj3uPtTMJpNeF33azGaS/mbSAfiJuzeeuxSERpl/AbxG+hVJj7j7H8Kma1oErvMs4GEzexn4EzASwMx+Bix092fN7Jtm9hpwELgh4FD518AgM1sJFAFjI3B9D5U5Es9jADMbCRS7+4OZx/Bb0td6rrtvCZuuaY0y3wAsJ/0qvBfcfUkuzqFf0yIiIjmTuOUvERFpPyoqIiKSMyoqIiKSMyoqIiKSMyoqIiKSMyoqIiKSMyoqIiKSM/8f9MavQQKzpDYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_04.trace,\n",
    "              varnames=['alpha', 'beta1', 'beta2'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "INFO (theano.gof.compilelock): Refreshing lock /home/petrushev/.theano/compiledir_Linux-5.0-arch1-1-ARCH-x86_64-with-arch--3.7.3-64/lock_dir/lock\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [sigma, beta2, beta1, alpha]\n",
      "Sampling 2 chains: 100%|██████████| 12000/12000 [00:21<00:00, 551.70draws/s]\n"
     ]
    }
   ],
   "source": [
    "with pm.Model() as m_05:\n",
    "    alpha = pm.Normal('alpha', mu=0, sd=1)\n",
    "    beta1 = pm.Normal('beta1', mu=0, sd=1, shape=2)\n",
    "    beta2 = pm.Normal('beta2', mu=0, sd=1, shape=2)\n",
    "\n",
    "    score_i = pm.Deterministic(\n",
    "        'score',\n",
    "        wines.score.mean()\n",
    "        + alpha\n",
    "        + beta1[0] * wines.judge_amer.values\n",
    "        + beta1[1] * wines.flight_red.values\n",
    "        + beta2[0] * (wines.wine_amer.values * wines.flight_red.values)\n",
    "        + beta2[1] * (wines.judge_amer.values * wines.flight_red.values)\n",
    "    )\n",
    "\n",
    "    sigma = pm.HalfCauchy('sigma', beta=5)\n",
    "    score_obs = pm.Normal('score_obs', mu=score_i, sd=sigma,\n",
    "                          observed=wines.score.values)\n",
    "    \n",
    "    m_05.name = 'second level interaction, v2'\n",
    "    m_05.trace = pm.sample(5000, tune=1000, chains=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEcCAYAAAD6GqKbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8VPW57/FPpKjEeMVUdntADhx9TI9t7WbHY60KtoiXrdJaKx4OSvB+q0W8bLRWa20tFeuFWtEiEJXdYr1LVWpUpGj0Bdtqxe74cIlACgqRcgsTLJacP9ZKTWMScpv8Ztb6vl8vXiSz1qz1/c1k8uT3e9YkBQ0NDYiIiHS3XUIHEBGRZFKBERGRrFCBERGRrFCBERGRrFCBERGRrFCBERGRrPhM6AAinWVmuwEzgUHAZuBSd19qZt8CbgNq4l1vBN4Angb6ABe6+9tmdhTwNXf/WSvHvwAYA+wAegPfd/eXO5HzIqAfcC9wg7tfYmYrgEPcfVuT/cri2ya245jlwGx3n9vK9gHAl919Tkfz7uS8w4CL3P3M7jyuJJNmMJLPzgfq3P0I4LvA3fHtQ4Br3H1Y/G8+MIKowFwCnGtmBcD3gLtaOrCZnQkcB3zD3YcRFZqHzGz/zoZ19w/c/ZLO3r+Dvg58rYfOJdIizWAkn30BeA7A3d3MSuLbhwBfMbPxwELgP4A6otlLH2ArMBp4oukMopkLgQnuvj0+/ntmdpi7r49nD33jf/8OXAMcDfQCbnf3R+LZ0V3ABuBj4HUzG0g06zgiPsd98W1rgbFNT25m340zNsT3mdJSyHhG8R/A34hmcrOBScBEoNDMKoH3gClAAbAeOAf4CvCz+H4vAMPc/dj4mL8DfgAMBi4lmr01AN9qdu6ZwP+KH9O73P2hVh5LSSnNYCSfvQWcbGYFZnYE8Hkz6wVUEM1ojgGKgIuIvokeQFQ4fkX0zfJPZnafmV3TwrE/B1Q3vcHd1zf59CV3PxI4Avif7n4UcCzwfTPbB5gK/F93H070Db4lU919KLCCaDYGgJl9ARgFHEVUuL5pZtbG43Ag8O04yzXu/neiIvNrd38amEa0fDgMeJaoIALs7u5Hu/tNwO5mdqCZ/Quwv7u/CRwM/Hs8tv8Gjm+ScU+ix/c04ATg723kk5TSDEby2QygBFgAvAq84e5/N7MZ7r4RwMyeAr7t7juA8fFt1xLNLq4nKkQ/NLOD3X1Jk2OvBPoDmxpvMLPjgbfjTz3+/4vAEDN7Of68NzAQOKDJ8V4l+km/qb+5++vxx5VEy3GL4s8PJSoaL8af7wsc1OSczS1294+Bj82svoXtJcA9cY3qDSxtNgaA6cDZwEdEfS2AdcADZlYHHAK81rizu2+JZ4i/AvYCZrWSTVJMMxjJZ6XAi/FP2I8A1XFv5W0z+x/xPt8gavADYGafBczdFwCFRD95NwB7NDv2DOAHZvaZ+H4HA/fzyU/qO+L/3wXmxbODrwO/BZYDq5ss2ZW2kH1XMzss/vho4J0m2xz4M3BsfNxyPilsLWnpFwru4JPXtwNnx8e6BvhdszFAtLR2MtHM7tdmtjdwE3AmcB5QT7TEBkA80xni7t8iWia8tfGxEmmkAiP5bCkw3sxeA24m6pk0EH1DfNzM5hMVkWlN7nM98OP443uA3wP/Avyp6YHdfTbwOvCKmf2B6Kf6Me6+rlmGOUCdmS0gKmQN7r6FaCnuQTN7kWg20txHwHfjY3+WqHg1nvtPRLOXV8zsv4hmL6vb/7AAsBgYGV+scHGc5RWipbNPFSt3r4sfg6o4/2aimddrRDPEeqJlw0YfAP3iHk8FcFs8ixL5hwL9NmUREckGzWBERCQrVGBERCQrVGBERCQrVGBERCQrVGBERCQrEnfd+htvvKHL4kREsmzIkCEFO9sniZcpN1RVVVFSUrLzPXNUvucHjSEXTJs2jffff58bbrghdJQuyffnARI5hp0WF0jgDEZEIueffz5VVVWhY0iKqQcjIiJZoRmMSEINGzaMTCbDwoULQ0eRlFKBEUmosrIy1qxZEzqGpJgKjEhClZWVqQcjQakHI5JQ27dvZ/v27aFjSIppBiOSUMcdd5x6MBKUCoxIQp133nmsXt3RPyMj0n1UYEQSasyYMerBSFDqwYgkVCaTob6+PnQMSTHNYEQS6qSTTlIPRoJSgRFJqIsvvlg9GAlKBUYkoUaNGqUejASlHoxIQm3atIktW7aEjiEpphmMSEKNHDlSPRgJSgVGJAesWp/h3AcWUV27lUHFezB9bCkD+hZ26ZiXX345f/nLX7opoUjHqcCIdNC4mQuZ57VZO/7SdXUcM3leNxxpN2Awt098hmOtmJnjDu+GY4q0nwqM5JwRd8xnydq60DFi1aEDdIt5XsvAic+EjtGigw8o4vkrhoaOIVmgAiM5J1e+2fTkn7k97vb5LK+tY0cD7FIAg4uLqJjQtcdBfw9GQgt6FZmZfdDGtmFmNrsn84iEMn1sKYOLi+hVUMDg4iKmjy3t8jGvvPJKysrKuh5OpJM0gxFJqFNOOUXvg5GgeqTAmNmhwO1AL2B/4GJ3r2yy/WXgXeAQoAAYFW86yMyeAz4LzHH3H5rZUOBGotlXETDa3Zf0xDhEst3gh+5s8jeqVpNfgihoaGjI+knMbBTw3+6+2MxGA8e6+/lm9oG794sLzAx3f9DMLiEqNI8DvwT+lagwrXL3/ePtT7r7GjO7Dihw9580nmvFihUNGzZsYPfdd8/6uLJl27ZteZ0f2jeGi56qYeVG/UEs+cSB+/Tm3pH9/+m2tLwecl3TMZSUlBS05z49tUS2GviBmdUDewKbW9jnpfj/SmBk/PE77v4RgJl93ORYU8ysDvg88GrTgwwcOJD6+voea85mQ082l7OlPWOYn+NjzPcm/9y5c1m1ahUXXHBBN6UMIy2vh1zXmTH0VJN/CnCju48FFhMtgzU3JP7/a8Cf449bml5NA8a5exmwppVjieSVbDT5TzjhBI4++uhuSCfSOT01g5kFPGJmG4C/EPVhmiszswnAVuAs4IttHGuBmW0F1gKfy0JekR41oG9hl2cszdXU1PD+++/n/U/Okr96pMC4++1ETf7mt/dr8um17v5uk89fjv/9077uPiE7KUWS5ayzztL7YCQoXaYsklDXX389K1euDB1DUiwnCoy7DwudQSRphg8frvfBSFD6ezAiCVVdXU1NTU3oGJJiOTGDEZHud84555DJZBgxYkToKJJSKjAiCXXTTTepByNBqcCIJNTQoUPVg5Gg1IMRSSh357333gsdQ1JMMxiRhLrwwgvJZDKcdNJJoaNISqnAiCTULbfcwooVK0LHkBRTgRFJqCOPPJJ99903dAxJMfVgRBLqnXfeYenSpaFjSIppBiOSUJdddhmZTIZTTz01dBRJKRUYkYSaPHmyriKToFRgRBKqtLSUoqKi0DEkxdSDEUmot956S2+0lKA0gxFJqPHjx5PJZDjttNNCR5GUUoERSag777yT6urq0DEkxVRgRBLqsMMOY7fddgsdQ1JMPRiRhFq0aBGLFy8OHUNSTDMYkYS6+uqryWQynHHGGaGjSEqpwIgk1N13360ejASlAiOSUIceeii9evUKHUNSTD0YkYSqrKzkzTffDB1DUkwzGJGEuu6668hkMowePTp0FEkpFRiRhLrvvvtYvnx56BiSYiowIgllZuzYsSN0DEkx9WBEEmr+/PksWrQodAxJMc1gRBLqxhtvJJPJcPbZZ4eOIimlAiOSUDNmzGDZsmWhY0iKqcCIJNSgQYP46KOPQseQFFMPRiShXnjhBSorK0PHkBRTgRFJqB//+Mfcd999oWNIimmJTCShHnroIZYuXRo6hqSYCoxIQvXv35+6urrQMSTFtEQmklBz585lwYIFoWNIiqnAiCTUpEmTuP/++0PHkBTTEplIQs2ePZslS5aEjiEppgIjklD9+vVjw4YNoWO026r1Gc59YBHVtVsZVLwH08eWMqBvYehY0gUqMCIJNWfOHGpqaigpKenU/cfNXMg8r+3mVO2zdF0dx0ye1+SW7P9lzmOtmJnjDs/6edJEBUakC0bcMZ8la3P1Sq1dgAO5deIzoYPkhXley8CsPlbZK5IHH1DE81cMzdrxO0sFRqQLcvFF3ejDDz9kyZIlHHnkkaGjtMtxt89neW0dOxpglwIYXFxExYShVFVVdXoWliuSMIbOaHeBMbNyYLa7z23HvgOAL7v7nJ3sVwy8CnzJ3bc12/YyUAiUAYOBG4CPgRnuPs3MfgcMB/Zpfl8Rgf3335/a2jBLXJ0xfWzpp3owkt+ydZny14GvtbWDmR0PPA/0a2O3s4GlwB3ACGAocIGZHeDuJwMfdE9ckeR5/PHHqaioCB2j3Qb0LWT62FIGFe9Bde1Wzn1gEavWZ0LHki7o6BLZJWZ2dXy/c919mZl9FxgNNACzgV8CE4FCM6sENgE3EhWzImC0uy8BdhDNQN7YyTlLgGXuvgHAzF4BjgEe6WB2kVSZMmUKmUyGyy+/vFP3D9nkh+aN/uz1L9Tcz56OFphKd59kZicBt5rZ9cAo4Kh4ewXwe2AScIi7P21mlwBj3H2NmV0HfAf4ibtXQPRnXXdiL6Ii1WgLsHdrO69YsYJt27ZRVVXVwaHljnzPD+kew0VP1bBy4/YsJOqgI64GyHLjOv9lv7nfKPtXwrXkwH16c+/I/l0+TtPXQ3v7SR0tMH+I/68EJgOHAgcCL8a37wsc1Ow+q4EpZlYHfJ6o59IRm4E9m3y+J7CxtZ0HDhxIfX19XjfUktAQTPMY5ufQuPPteWip0T/lxM/m1Rhakm/PQ0s6M4aO9mAa55FHA+8ADvwZONbdhwHlwNtEy1+Nx54GjHP3MmANUNDBc1YBB5nZfma2K9Hy2GsdPIZI6jz88MM899xzoWN0yPSxpQwuLqJXQQGDi4vU6M9zHZ3BHGFmLxH1W85x95Vm9iLwipntBiwkmrEsBr5vZn8EZgELzGwrsBb4XEdO6O7bzWwC0dLbLkRXka3uYG6R1Jk6dSqZTIYJEyaEjtJuA/oWUjHhny/9rloXKIx0WbsLTDwDaen2yUTLZU29CTQ2V2bv5LgD23HuOUCblzyLyD979tlneffdd0PHkBTL9d+m/KCZfaGlDfH7YNq6xFkk1QoLC+nTp0/oGJJiOftO/rin09b2k3soikhemjVrFqtXr8775rLkr1yfwYhIJ91///089thjoWNIiuXsDEZEuqaioiLv34sk+U0zGJGE6t27N7179w4dQ1JMMxiRhCovL2fNmjXqwUgwmsGIJFR5eTlPPvlk6BiSYprBiCTUyy+/rB6MBKUZjIiIZIUKjEhCTZs2jUce0V+1kHBUYEQSKh9/2aUki3owIgn1wgsvqAcjQWkGIyIiWaECI5JQ99xzD7/5zW9Cx5AU0xKZSELNmTOHurq60DEkxVRgRBLqueeeUw9GgtISmYiIZIUKjEhC3XXXXTz00EOhY0iKaYlMJKFefPFF9WAkKBUYkYR6+umn1YORoLREJiIiWaECI5JQt912GzNmzAgdQ1JMS2QiCfXaa6+xZcuW0DEkxVRgRBLqscceUw9GgtISmYiIZIUKjEhCTZo0iWnTpoWOISmmJTKRhHrrrbfYvHlz6BiSYiowIgk1e/Zs9WAkKC2RiYhIVqjAiCTUzTffzNSpU0PHkBTTEplIQrk7mzZtCh1DUkwFRiShZs2apR6MBKUlMhERyQoVGJGEuuGGG/jFL34ROoakmJbIRBKqpqaGjRs3ho4hKaYCI5JQM2fOVA9GgtISmYiIZIVmMCIJde211/Lhhx8G/31kq9ZnOPeBRVTXbmVQ8R5MH1vKgL6FQTNJz1CBEUmo9evXs2nTJsbNXMg8rw0dB4Cl6+o4ZvK8TtyzulvOf6wVM3Pc4d1yLNk5FRhJlRF3zGfJ2roO3KN7vrEFsd9I2A/IkeKSC+Z5LQMnPhPo7LnxtXTwAUU8f8XQHjmXCoykSkdeWFVVVZSUlGQxTfblwhiOu30+y2vr2NEAuxTA4OIiKiboeUiDnG3ym1mZma0yswnx5+PNbFL88elm9m7j5yLyaVdddRW33npr6BhMH1vK4OIiehUUMLi4iOljS0NHkh6S6zOYXwNTzew/gcOBxwDc/VEzKwIOCRlOJJfV19fz0UcfBTt/8+b+vKuGqbmfMrleYAB2Bx4AKlBBEWm3zL+NpTJoz+ETnW/uN+pa/0LN/TByvsC4+wbgeTMra8/+K1asYNu2bXn9BrN8zw/hx3DRUzWs3Li9G46UG41Z6Zqwzf1G4b6WDtynN/eO7N+lYzR9Tbe3n5TzBaajBg4cSH19fV431JLQEAw9hvndcO7QY+iq8ePH89e//pUHH3wwyPm72txvlO/PA6R3DDnb5BeR/KbmviRuBiMikTvvvDPoMuWAvoWdmrFIcuRNgXH38tAZRESk/XJ9iWx04/tgmjKz04GJAfKI5I1LL72Um2++OXQMSbGcncHEM5byVrY9Cjzak3lE8k2fPn3IZDKhY0iK5WyBEZGuue222/L+cnfJb7m+RCYiInlKMxiRhLrgggvYuHEjv/3tb0NHkZRSgRFJqL59+9LQ0BA6hqSYCoxIQv30pz9VD0aCUg9GRESyQjMYkYQaN24cGzdu5IknnggdRVJKBUYkofr370+fPn1Cx5AUU4ERSagf/ehH6sFIUOrBiIhIVmgGI5JQY8aMYdOmTcyZMyd0FEkpFRiRhDIzamtrQ8eQFFOBEUmoH/zgB+rBSFDqwYiISFZoBiOSUGeeeSabN2/m2WefDR1FUkoFRiShDjvsMNatWxc6hqSYCoxIQk2cOFE9GAlKPRgREckKzWBEEurb3/42W7Zs4fnnnw8dRVJKBUYkob761a+ydu3a0DEkxVRgRBLqqquuUg9GglIPRkREskIzGJGEOvXUU6mrq+Oll14KHUVSSgVGJKG+8Y1vqAcjQanAiCTU9773PfVgJCj1YEREJCs0gxFJqBNPPJG6ujoWLFgQOoqklAqMSEKdcsopfPDBB6FjSIqpwIgk1CWXXKIejASlHoyIiGSFZjAiCTV8+HC2bt3Ka6+9FjqKpJQKjEhCjRo1ivfffz90DEkxFRiRhDr//PPVg5Gg1IMREZGs0AxGJKGGDRtGJpNh4cKFoaNISqnASI9YtT7DuQ8sorp2K4OK92D62FIG9C0MHSvRysrKWLNmTegYkmIqMCkybuZC5nltD56xusVbl66r45jJ83owR+RYK2bmuMN7/LyhlJWVqQcjQanANDPijvksWVsXOgatfXOWzpvntQyc+EwH75WE5yEaw8EHFPH8FUMDZ5E0UYFpJhdegFVVVZSUlISO0SXNx3Dc7fNZXlvHjgbYpQAGFxdRMSH8Y92WfH8e1IOR0NpdYMysHJjt7nPbse8A4MvuPqeNfa4Azow/fdbdb2q2/WWgECgDBgM3AB8DM9x9mpn9DhgO7OPu29o7Dglj+tjST/VgJLvOO+88Vq9eHTqGpFi2LlP+OvC11jaa2SDg/wFHAkcAI8zsSy3sejawFLgDGAEMBS4wswPc/WRAv8kvTwzoW8j0saUMKt6D6tqtnPvAIlatz4SOlWhjxozh1FNPDR1DUqyjS2SXmNnV8f3OdfdlZvZdYDTQAMwGfglMBArNrBLYBNxIVMyK4n3fA05w978DmFlvoLVZSAmwzN03xPu+AhwDPNLB7BLruWZ/6/2LUI1+SE+zP5PJUF9fHzqGpFhHC0ylu08ys5OAW83semAUcFS8vQL4PTAJOMTdnzazS4Ax7r7GzK4DvuPuPwE+NLMCYDLwprsvaeWcexEVqUZbgL1bC7hixQq2bdvWpatnLnqqhpUbt3f6/t0jCc3l3NSxZn8CnofffnqsB+7Tm3tH9g8QpuO6+nrOBUkbQ3t7kx0tMH+I/68kKgyHAgcCL8a37wsc1Ow+q4EpZlYHfB54FcDMdgdmEBWMS9o452Zgzyaf7wlsbG3ngQMHUl9f36Xm7PzAjd18by5Dy2PIt0Z/vj8PDz/8MKtXr2bChAmho3RJvj8PkN4xdLQH07iucDTwDuDAn4Fj3X0YUA68DexocuxpwDh3LwPWAAXxzOUp4E/ufmHjUlkrqoCDzGw/M9uVaHlMvx42D00fW8rg4iJ6FRQwuLhIjf4sGzVqFCeeeGLoGJJiHZ3BHGFmLxH1W85x95Vm9iLwipntBiwkmrEsBr5vZn8EZgELzGwrsBb4HPBNoob9bmbW+Aq41t0/VTjcfbuZTSBaetuF6CoyXRqThwb0LczpGUvSbNq0iS1btoSOISnW7gITz0Baun0y0XJZU28CFn88u5VD7t6Bc88BWr3kWUQ+beTIkXofjASV679N+UEz+0JLG+L3wfTr4TwieePyyy9nzJgxoWNIiuXsO/njnk5b20/uoSgieem0007L+yuXJL/l+gxGRDrpww8/ZMOGDaFjSIrl7AxGRLrm9NNPVw9GglKBEUmoK6+8kpqamtAxJMVUYEQS6pRTTlEPRoJSD0YkoT744ANqa3vyD8yJ/DPNYEQS6swzz1QPRoJSgRFJqIkTJ7Jq1arQMSTFVGBEEuqEE05QD0aCUg9GJKFqamp4//33Q8eQFNMMRiShzjrrLPVgJCgVGJGEuv7661m5cmXoGJJiKjAiCTV8+HD1YCQo9WBEEqq6ulrv5JegNIMRSahzzjmHTCbDiBEjQkeRlFKBEUmom266ST0YCUoFRiShhg4dqh6MBKUejEhCuTvvvfde6BiSYprBiCTUhRdeSCaT4aSTTgodRVJKBUYkoW655RZWrFgROoakmAqMSEIdeeSR7LvvvqFjSIqpByOSUO+88w5Lly4NHUNSTDMYkYS67LLLyGQynHrqqaGjSEqpwIgk1OTJk3UVmQSlAiOSUKWlpRQVFYWOISmmHoxIQr311lt6o6UEpRmMSEKNHz+eTCbDaaedFjqKpJQKjEhC3XnnnVRXV4eOISmmAiOSUIcddhi77bZb6BiSYurBiCTUokWLWLx4cegYkmKawYgk1NVXX00mk+GMM84IHUVSSgVGJKHuvvtu9WAkKBUYkYQ69NBD6dWrV+gYkmLqwYgkVGVlJW+++WboGJJimsGIJNR1111HJpNh9OjRoaNISqnAiCTUfffdx/Lly0PHkBRTgRFJKDNjx44d3X7cVesznPvAIqprtzKoeA+mjy1lQN/Cbj+P5D8VGJEcN27mQuZ5bReOkL0ryZauq+OYyfOydvxPtG8Mx1oxM8cdnuUs0l4qMJJTRtwxnyVr60LHaEKX+eaTeV7LwInPhI7RivBfSwcfUMTzVwztsfOpwEhO6ckv/p2pqqqipKQkdIxOq66uZtmyZYwYMaJbj3vc7fNZXlvHjgbYpQAGFxdRMSF7z1u+Pw+QjDF0Rs5epmxmZWa2yswmxJ+PN7NJ8cenm9m7jZ+LyKcNGjSI/v37d/txp48tZXBxEb0KChhcXMT0saXdfg5JhpwtMLFfA1PN7D+BSxtvdPdHARUXkTa88MILVFZWdvtxB/QtZPrYUgYV70F17VbOfWARq9Znuv08kv/yYYlsd+ABoAI4JHAWkSA63+jvx81Z7kn0TKO/4/0LNfzDy/kC4+4bgOfNrKw9+69YsYJt27bl9V/yy/f80LkxXPRUDSs3bs9Sos4K35iVzsm9hn/ufC0duE9v7h3ZseXTpq/p9vaTcr7AdNTAgQOpr6/P64ZaEhqCnRnD/Bwbc1qfh/boyUa/nofc0Jkx5HoPRkQ6ae7cuSxYsCArx1ajX9ojcTMYEYlMmjSJTCbDBRdc0O3HHtC3MKuXJksy5E2Bcffy0BlE8sns2bNZsmRJ6BiSYrm+RDa68X0wTZnZ6cDEAHlE8ka/fv0oLi4OHUNSLGdnMPGMpbyVbY8Cj/ZkHpF8M2fOHGpqavK+uSz5K9dnMCLSST//+c8pLy8PHUNSLGdnMCLSNY8++qh6MBKUZjAiCbX//vuz7777ho4hKaYCI5JQjz/+OBUVFaFjSIqpwIgk1JQpU5g1a1boGJJi6sGIJNRTTz2Fu4eOISmmGYxIQu29997sueeeoWNIimkGI5JQDz/8MKtXr9b7YCQYzWBEEmrq1KnMnj07dAxJMc1gRBLq2Wef5d133w0dQ1JMMxiRhCosLKRPnz6hY0iKaQYjklCzZs1SD0aC0gxGJKHuv/9+HnvssdAxJMU0gxFJqIqKin/8DXWREDSDEUmo3r1707t379AxJMU0gxFJqPLyctasWaMejASjGYxIQpWXl/Pkk0+GjiEpVtDQ0BA6Q7d64403kjUgEZEcNGTIkIKd7ZO4AiMiIrlBS2QiIpIVKjAiIpIVKjAiIpIVeX+ZspntDcwC9gJ2BSa4+2vN9jkfuBD4GPixu/+ux4O2g5l9C/iOu49uYdtdwFHAlvimke6+qSfztcdOxpDTz4OZ9SH6Wvos0eM81t1rm+3zFLA/sB2od/cTezxoC8xsF+Ae4MvAR8B57r6syfacfuyhXWPIl9fA/wF+5u7Dmt1+CnAD0XMww92nBYjXLm2M4QrgPKDxdXGht/FX7fK+wAATgBfd/U4zM+A3wL82bjSzfsDlwL8BuwOvmFmFu38UJG0r4hfP8cBbrewyBDje3T/suVQd09YY8uR5uBhY7O4/NLMzgeuB7zXb5yDgf7t7rl0d801gd3f/qpkdAfwcGAl589hDG2OI5cNr4BrgLGBrs9t7A3cApfG2V83saXdf2/Mp29baGGJDgLPd/Y32HCsJS2R3APfFH38G2NZs++HAq+7+UfzTzjLgSz2Yr70qib7BfUr8k91BwK/M7FUzO6dHk7Vfq2MgP56Ho4C58cfPAcObbjSzA4B9gDlm9oqZndzD+dryj+zu/jpRMWmUD489tDGGPHoNLAdOa+H2EmCZu29w978BrwDH9Giy9mttDBAVmGvjr/9rd3agvJrBmNm5wBXNbh7n7ovin9JmAeObbd8LaDqN3gLsnb2UbWtjDA+b2bBW7rYH8AvgdqAXMM/M/svd385e0tZ1cgz58Dys5ZOMLeXblein6ruA/Yh+Cl3o7uuymbWdmj++fzezz7j7xy1sC/rYt6GtMeTUa6A17v6YmQ30IPVFAAAB2ElEQVRsYVO+PAdtjQFgNvBLYDPwhJmd3NZya14VGHefDkxvfruZfZFo4Fe5+/xmmzcDTf8w+Z7AxqyF3InWxrATGeAud88AmNlLROvUQV5cnRxDzj8PZvY4n2RsKd8HwL3xN7x1ZvYmYEAuFJjmj+8ucc6WtgV97NvQ1hhy6jXQCfnyHLTKzAqAOxv7Xmb2DPAVoNUCk/dLZGb2BeARYLS7P9fCLguBo81s9/iCgBLgnZ7M2A0OJvppuVe8lnsU8MfAmToqH56HV4GT4o9PBBY02z6c6GsNMysCDgVy5dcV/yN73L9Y3GRbPjz20PYY8v01UAUcZGb7mdmuRMtjr+3kPrlmL+AdMyuKi83XgTZ7MXk1g2nFT4kal3dFPX42uftIM5tAtOb5tJlNIfpmsQvwfXdv3qfJSc3G8BDwOtHVSw+6+5/DpmufPHsepgIPmNkrwN+A0QBmdivwqLs/Z2bHm9nrwA7guhxqOD8BHGdmlUABMC7PHnvY+Rjy7jVgZqOBInf/VTyW3xM9BzPcfXXYdO3TbAzXAfOIrvJ70d2fbeu++lUxIiKSFXm/RCYiIrlJBUZERLJCBUZERLJCBUZERLJCBUZERLJCBUZERLJCBUZERLLi/wPWVwXFvYJ/9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pm.forestplot(m_05.trace,\n",
    "              varnames=['alpha', 'beta1', 'beta2'],\n",
    "              rhat=False)\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/petrushev/etc/jup_env/lib/python3.7/site-packages/theano/tensor/subtensor.py:2197: 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",
      "  rval = inputs[0].__getitem__(inputs[1:])\n"
     ]
    },
    {
     "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>WAIC</th>\n",
       "      <th>pWAIC</th>\n",
       "      <th>dWAIC</th>\n",
       "      <th>weight</th>\n",
       "      <th>SE</th>\n",
       "      <th>dSE</th>\n",
       "      <th>var_warn</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>second level interaction, v2</th>\n",
       "      <td>864.67</td>\n",
       "      <td>4.65</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>17.19</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>second level interaction</th>\n",
       "      <td>866.17</td>\n",
       "      <td>5.31</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0</td>\n",
       "      <td>17.07</td>\n",
       "      <td>0.33</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>full interaction</th>\n",
       "      <td>866.45</td>\n",
       "      <td>6.59</td>\n",
       "      <td>1.78</td>\n",
       "      <td>0</td>\n",
       "      <td>17.46</td>\n",
       "      <td>2.29</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>full interaction, v3</th>\n",
       "      <td>866.64</td>\n",
       "      <td>5.63</td>\n",
       "      <td>1.97</td>\n",
       "      <td>0</td>\n",
       "      <td>17.16</td>\n",
       "      <td>0.37</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>full interaction, v2</th>\n",
       "      <td>867.83</td>\n",
       "      <td>6.15</td>\n",
       "      <td>3.16</td>\n",
       "      <td>0</td>\n",
       "      <td>17.16</td>\n",
       "      <td>0.42</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                WAIC pWAIC dWAIC weight     SE   dSE var_warn\n",
       "second level interaction, v2  864.67  4.65     0      1  17.19     0        0\n",
       "second level interaction      866.17  5.31   1.5      0  17.07  0.33        0\n",
       "full interaction              866.45  6.59  1.78      0  17.46  2.29        0\n",
       "full interaction, v3          866.64  5.63  1.97      0  17.16  0.37        0\n",
       "full interaction, v2          867.83  6.15  3.16      0  17.16  0.42        0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pm.compare({\n",
    "    m_01: m_01.trace,\n",
    "    m_02: m_02.trace,\n",
    "    m_03: m_03.trace,\n",
    "    m_04: m_04.trace,\n",
    "    m_05: m_05.trace,\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdgAAAE4CAYAAAD4neIbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X18VGed///XZJKQcBMg3BQKJITQXrYGqIW2tNxsaTdFv7TKuura2qpVqY9t3e96s98iq+t2bRetq2t11Z9bb7bW211XpQrWIoUKpC2lVDrEtlcbEgIUqNyE+wzJTOb3x3UCk2RCQjInZ2byfj4e88icmTnnfJJc1/nMdZ3rXCeUSCQQERGR9MoLOgAREZFcpAQrIiLiAyVYERERHyjBioiI+EAJVkRExAdKsCIiIj5QghUREfGBEuwgZ4x53Bjzfwdwf/9ojPl5H9f9hjHmkTSHJOIrY8z7jDHPDeD+yowxJ40xY/qw7hxjjCZHSJOQJpqQbGGM+QYw3Fr7waBjEclFxpg5wFZrbSjoWHJBftABZCpjzFTgFeBe71EM/AZYZq1tNcYMBf4VeA8wBNgIfNxau9sYcz3wn8AvgGVAHPihtfb/dbOvIuDfgMXAxcBh4AFr7Xe89xPedv4RuAj4EfBz4BvAJOBx4DZrbczb1gNJcf3Oi6vJi+t7wIvADd427wFWW2u/bIwpAb4K/BXQBqwB7rbWnjLGvAtYDkzDlZvHgQ9Za097rcpTwJuAa4AGb59Ppvhd7wPmWGtvTn7uvTcWOAhUWGt3GWNmAd8GZgLPAUeAE95n872/2fuBk8BD3vJ0b93Lga8Bs4FDwFestf+Z6u8vg0MvyvBh4ApgLvAqcBfwGeAvgUbgdmvtH71t3QLcD1QAFviktXaz994uXL17J7AJd9z4B2ttlff+h4BPAxOBl4G/s9ZuMcaUca7MjsMdf+6x1j59IccU79jV4G1jePtza+0h7/3VwPPW2vu848VDwHtxdevRTtu6BfiSF+uTQAiIJK2b8ljTm//HYKAu4vMbgivslwILgHcAf+29923gau9RBvwZWG2MKfDev9T7eTGuAH7cGDO3m/18ArgKV7FLgM8DXzfGDE/6zDuBWcAc4MPeZxYAVbhk+Q7vc1/CJbmrgem4A8n3k7YzDdiAS9SrO8XxTdwB401AJVAOrDTGTAF+CHzKWjvGi2MBcGvSuh/AfQEYAzzlbavPjDGFwK+B3wOjcAezJUkfWQFUA2/x4rkRCHvrDvfW+733e74b+Kwx5h3IoNTLMvwR4P8Bo4GjuC/N3wBKcV9KP+9taw7wU+BTuPL+JWCNMebipG1djqtLH+oUx43A13FJciTuy/JvvPL+HeB1XL0dDbwAPJi0+oUcU3rrPtyx4jLgSuD6pFinA/8DfBYYizte/FXSuj0dawY9tWB79qC19jTwknce5VLvm9t7gRusta8DGGM+jvsGfFXSuv9qrW0FNhpj6nEV5NkU+/g28F1cpZ4MNANFuIp90vvMf1hrTwAvG2P2A9+31h709v0SUG6MCeEOEjdZaw94730K2G+MGZe0v0ettWe89/F+DsElokXW2j97r92Oa7m/AVxurW0wxowGJuBahZOStvlba+0Wb72f4FrG/TEPd5D5vLU2Bqw3xvxv0vsfAP7JWrvb2+c/cC4BLwGarbVf8pZfNMZ8E9cieayfcUl26k0Z/p219gUAY8xmoMBau85bXgd83Pvch4GfJfXQ/K8x5iPA+3C9KAC/sNae8tZNjuN24MfW2j94730DeD5pu+2tv3Lc8SA5Puj9MaW3bgXutdbu9+L5Z9wXU3DHuI3W2l94y/9ljPmo97nzHmvaj02DnRJsz/6c9LwV1+ofDRQAu9rfsNZGvcRXBhwATrVXsOR1jTHvw3X1tHsbsBvX4rvO2+YO773kHoZDSc/juMrXrs377DhcQlzTaaBCFPdtGuC0tfZYit9zNK7Fvjvpd9oLZyvTnd5BJIr7Nj+0U3yd/079PYczEXjDS67tGoAp3vPJybGS9L/AHZzKjTHJf6M8b30ZnFrpuQz3po6BK1+LvC7ndvnAa0nL+7qJYyKu2xgAa20b8DSAMeZSXKuwHNd1fKJTfCmPKd3sp7cmAnuTlpPrSOc6BufqWU/HGiVYlGD76g3gDK4gtSehYlzXzYHzrWit/THw4+TXjDGP4wruRd753Utx5xaT9WY02iEvruustX/ytp2P676pA+afZzsHgRZcpWpvlV8JLPJ+xw8AVycl3ad7EU9P4kBh0nLyqMfXgYnGmEJrbYv3WvK3+T24LzM13vKUTutut9ae7U0wxoxH5X0wew89l+Hejvh8HfhG8vlPY0wF51qf59vWXlwdO8sY8yXcF+xf4c65/sh7/S7gzb2MqTtx7+f56lk55+pR5zq2sNP2puDODfd0rBF0DrZPvG+djwJfNMZc7A14+iruW2vNeVdObTSusMa9RNDetVnY/SrdxvUD4EvGmHHeeZ0HcOdEC3pYNw78BPi8MabUGDMKd/5nshdfHDhjjMk3xnwYd+7lguJL4VVgrjFmmvc3/DTnDkw1uL/nSmNMoTFmHu4g2e57wKeNMVOMMSOALya9twaYYoy5y4t3Mq7ba3k/45Xslc4y/AjwIWPMPGNMyCubL+K+jPbkR8BtxphrjTF5xpi7gTtwreWhuNNDGGNm4M7x9reOvQEcw3Vftw9amt3pd/lH4y7tKcWdk02OdYEx5h3GmLAx5r24XrZ+HWsGEyXYvvsksA3YCuzHffOr9s6PXKi/xw1UOuZtbzvum25VH7b1CVzi2oZrlc4F3mqtbe5lHLuAP+G6u+pxg4kewY3i3YlLektxA0b6El+yX+FGQ2/xYt4KnAbwuoaX4AZ1NeFGOv4qad2v4JJwrRfvK97rLdbao8BNwN/guq63AZtxo8FlcHqENJVha20NbpDSt3B19lHgH621vzrvim7dp4CP4b4gHsUlvrd5p20+CjxkjDmO+7L7XWB8p/ETFxprC/C3wEe87b4f+FnSR/4VN/p3O/ASrp60r9uIO2f877jxJe/E1dH2HqX+HGsGBV0HKwPKGPMvwBXW2n6N6DXGXAPUWWsPe8tVuFbEcFVwGcy87up6oLQ/l8x4lw2VWGtrk157HnjYWvtw/yPNfTonJQPGu3zG0HFAVF99APft/v24AVUrgD8oucpg5l3hUIU75ZRqMOOFmAg8boy5FtdSXeptu8v17ZKauohlIG3FdSN9Ow3b+ixuZOce3ECNIrzzTCKD2Cdwgyi/4J0n7TPvsrt/Bp4AjuPOz77bWruzv0EOFr3qIjbGfJJzF0xvxZ0ruBR3jmAU7jzYBzoNIRcRERm0emzBGmOuBu7EjbibgetWvgc3wmy5tfZNuJPjn/UxThERkazSmy7iJuBj1tpT1toEbiDJTGBk+2wkuJbs3/gUo4iISNbpcZCTtfY1vBlKvGs0P4Y7h3ZJ0sf20eniaYBt27aFcJMvnEhHsCKD1Ahg3+zZs9M25F91UyQtzls3ez2K2LtDwxrchNQbgZs7fSTVCfWL6TgNl4j0zdkZttJEdVMkPbqtm71KsMaYK3DJ9YvW2v/wro+amPSRzvNZtjsBMGPGDMLh8AVF3F+1tbVUVfV3HgSRc4IoU/F4nB07dkD6W5qqm5IzMrVu9phgvVlEfoe7L+gvAbx7np42xvyFdx72Q7gEnFI4HB7wShwKhQZ8n5LbcrFMqW5KLsjUMtWbFuzHcfco/Zwx5nPea2twtzn6jjFmJO4ODLf5E6KIiEj26c0gp88An+nm7WvTG46IiEhu0ExOIiIiPlCCFRER8YESrIiIiA9yL8HG47ByJRXLlsHKlW5ZRERkgOXe7eoefBDuv58R0ShEIhAKwYoVQUclIiKDTO61YNevh2jUPY9G4UndulBERAZe7iXYG26AoiL3vKgIbrwx2HhERGRQyr0u4uXLIRTixKpVjFi6FO69N+iIRERkEMq9BBsOw4oVNCxZwsyZM4OORkREBqnc6yIWERHJAEqwIiIiPlCCFRER8YESrIiIiA+UYEVERHygBCvSg5YzzXzvXZUc+asr+d67Kmk50xx0SNmppQWqq7l8/nyornbLIjlMCVakBz98XxW3PVbP9fVxbnusnh/ePiPokLLTkiWwbh35J07AunVw881BRyTiKyVYkR5M376b4ph7XhyD6X9sDDagbLVtW8fl558PJg6RAaIEK9KDuivKaPamZGnOh7q3lAcbULaaPbvj8pw5wcQhMkCUYEV6cMePa/nJ0kqemhbmJ0srueNHO4IOKTutWQPV1cRGjHDnYFevDjoiEV/l3lSJImlWOKSYD/+8jkgkwvWafrPvCgth7VpeikQ0jakMCmrBioiI+EAJVkRExAdKsCIiIj5QghUREfGBEqyIiIgPlGBFRER8oAQrIiLiAyVYERERHyjBioiI+EAJVkRExAdKsCIiIj5QghUREfGBEqyIiIgPen03HWNMCVAD3GKt3WWMuQn4Mi5J/xH4sLW2xZ8wRUREskuvWrDGmGtxydUkvfw94L3W2iqgGHh/+sMTERHJTr3tIr4LuBvYl/RaGBhhjAkDQ4DmNMfWN/E4rFxJxbJlsHKlWxbph5ZYC9WPVjN/9XyqH62mJaaOmr5oOXGMveOLueyKWewdX0zLiWNBhyTZrrkZKiupuvJKqKx0yxmkVwnWWnuntXZTp5fvBp7CJd1xwP+mN7Q+evBBuP9+Rjz3HNx/P3zpS0FHJFluyU+WsK5hHSdaT7CuYR03//TmoEPKSn+unMCkg1EKEjDpYJQ3pk8MOiTJdlVVUF9PXjwO9fUwY0bQEXXQ63OwyYwxE4AvAlVAA/Dv3uOeVJ+vra0lFAr1NcYLUvHYY4yIRt1CNMqJVatoWLJkQPYtuWnL3i0dlp/d8yyRSGRA9p1IJHzd/kDWzcsORWnfUwiYcLB5wP6OkpuqGhs7tBLbdu2iNoPqZp8SLLAAqLXW7gQwxnwH+J/uPlxVVUU4HO7jri7QO94BkQhEo1BUxIilS5k5c+bA7Fty0jXbr2Fdw7qzy3OnzB2wMhWPx9m+fbtv2x/Iurl3bBGTDrokmwAOjCtW3ZT+KS93LVdP3tSpGVU3+3qZTi0w1xjT3sfzdmBbH7eVXsuXw+c+x4mrr4bPfQ7uvTfoiCTLrbltDdXTqhlRMILqadWsvnV10CFlpfE7D7B3fDGtIdg7vpiL6vYHHZJku9paqKykLRx252B37Ag6og761IK11r5sjPkM8KQxJgbUA8vSGllfhcOwYgUNS5bo27GkRWF+IWvvWEskElGZ6ofCESOZ8sZp/R0lfYqLoa6O2gwtUxeUYK21U5Oe/wD4QboDEhERyQWayUlERMQHSrAiIiI+UIIVERHxQc4l2HhbnJWbVrJs8zJWblpJvE0zOYlkgni0mfo5lZRfeyX1cyqJRzNr1h2RdOvrdbAZ68GaB7l/4/1EY1EiGyOECLFiwYqgwxIZ9BrnV1GxrZ4QULKtnvoFM6jcWhd0WCK+ybkW7PqG9URjbianaCzKkw1PBhyRiACMfXl3h5mcxr3UGGQ4Ir7LuQR7Q8UNFOUXAVCUX8SNFTcGHJGIABy6rIz2yeUSwMHLy4MMR8R3OddFvHzeckKEWLVjFUtnLOXeeZrJSSQTlG+upX7BDMa+tItDl09l6qbMmnVHJN1yLsGG88KsWLCCJSM1k5NIJgkXFVO5tU4zOcmgkXNdxCIiIplACVZERMQHSrAiIiI+UIIVERHxgRKsiIiID3IvwcbjsHIlFcuWwcqVblmkHzT9Zno0nzzGJlPM5GtmsckU03zyWNAhSZY7dvgAdaUhLrtiFnWlIY4dPhB0SB3k3GU6PPgg3H8/I6JRiEQgFIIVmipR+k7Tb6bH87MnMP/VKCFg/qtRNs2eyEJ7OuiwJIsdvGQilU1uZrDKJqi7dCIjDyd6XG+g5F4Ldv16iLqpEolG4UlNlSj9o+k30+PNu6Mdpkqs2q3J/qV/yo/SoUxNbQoymq5yL8HecAMUuakSKSqCGzVVovSPpt9Mjz+VFXWYKrG2rDjIcCQHNI6iQ5naNTrIaLrKvS7i5cshFOLEqlWMWLoU7tVUidI/mn4zPeZsO8Cm2ROp2t1MbVkxV23bH3RIkuXGvbafuksnMrXJJdfxr2ZWmcq9BBsOw4oVNCzRVImSHpp+Mz2Kh49koT1NJBJhof6OkgYjx0xg5OFExk6/mXtdxCIiIhlACVZERMQHSrAiIiI+UIIVERHxgRKsSA80k5OI9EXujSIWSTPN5CQifaEWrEgPNJOTiPSFEqxIDzSTk4j0hbqIRXqgmZxEpC+UYEV6oJmcRKQv1EUsIiLiAyVYERERHyjBioiI+KDX52CNMSVADXCLtXaXMeZa4KvACCACfMBa2+JPmCIiItmlVy1YL5nWAMZbLgF+CdxlrX2z97FlvkQoIiKShXrbRXwXcDewz1uuBp6x1ka85b/DJdzAxaPN1M+ppPzaK6mfU0k82hx0SJLl4sePcaqkmKpZszhVUkz8+LGgQ8pKLbEWqh+tZv7q+VQ/Wk1LTB1e0j8H99RxKj/EjFmzOJUf4uCeuqBD6qBXCdZae6e1dlPSS9OBE8aYXxhjIsDngSY/ArxQjfOrqNhWz8jTcSq21bNrwYygQ5IsF508gaEnouQBQ09EiU6eGHRIWWnJT5awrmEdJ1pPsK5hHTf/9OagQ5IsN7TiEobGIQQMjUNxxSVBh9RBX6+DzQf+DzAX2AV8D/g0cF+qD9fW1hIKhfq4qwtT/nIj7XsKAWNf2kUkEjnfKiLnVXUi2qFMFZ9oHrAylUgkfN3+QNbNLXu3dFh+ds+zqpvSLzO85Aru57A4GVU3+5pgDwDPWWt3Ahhj/gf4WHcfrqqqIhwO93FXF6b+snJKttUTAhLAocunanIA6ZdTI4oY6iXZBNA8onjAylQ8Hmf79u2+bX8g6+Y1269hXcO6s8tzp8xV3ZR+ORXmbAs24S1nUt3s62U6a4G3GGPKveX/A7zQx22lVfnmWurnVHJsaJj6OZVM3bQj6JAkyxXtPcDpEcW0AadHFFO0d3/QIWWlNbetoXpaNSMKRlA9rZrVt64OOiTJcqcbXuNU+FxybW54LeiQOuhTgrXW7sGNGv61MeYVYDzwhXQG1lfhomIqt9bR+MwLVG6tI1xUHHRIkuXCJSMZdvw0tS++yLDjpwmXjAw6pKxUmF/I2jvWsvnmzay9Yy2F+YVBhyRZbtyU6QyPJdjx4osMjyUYN2V60CF1cEFdxNbaqUnP1wBr0h2QiIhILtBMTiIiIj5QghUREfGBEqyIiIgPlGBFRCQrxVtb2LCsmuh75rNhWTXx1syaHUw3XBcRkay08e4lzH1kHcUxaN65jo2hm1n08NqgwzpLLVgREclKI5/eRnHMPS+Owcia54MNqBMlWBERyUrHrptNs9cP25wPx+bNCTagTtRFLCIiWWnht9awMXQzwzY+y6mFc1n4zcyaHUwJVkREslK4oJBFD68lEolwdQbOa60uYhERER8owYqIiPhAXcQiMjASCaipYfSGDXD8OMybBwN0L1qRICjBioj/Ghth8WJoaODicBgeeAAqKuCJJ6C8vOf1RbKQuohFxF+JhEuuO3dCSwvh5mZoaXHLb32re18kB+Vcgm05eoRjw/KpmjWLY8PyaTl6JOiQJMsdOXKAtdNCTL5mFmunhThy5EDQIWWXmhrYtQtisY6vx2Lw2mtw773w+99DXZ1LvCK9dGz/bk7lh5gxaxan8kMc27876JA6yLku4uZJ4yk5HScElJyOc2zyRRSebA06LMliz8+ZSHUDhIDqBvj9VRO5aadaXb1WVwcFBXDmTNf34nH48pfdAyAvD6ZMgWnTXBfytGkdn48bp/O2clb+lHKGxl3dHBqH02Xl0Jo5dTPnEuwIL7mC+6OXnIqd7+MiPZqznw5las6+IKPJQtOnd98yLSiAlSth2DCor3ePhgbYtg02bOj6+WHDuk++U6fC0KG+/iqSWdqTK3hJNsMO9zmXYE8MDZ9twSaA48PyGRV0UJLVnp/I2RZsAnj+Yrgp6KCyybx5UFFB26uWvKTGRVsI8ior4VOf6toqTSSgqelcwm1Pvu3Lv/1t1y5ngIkTOybf5AR88cWuhSw543T4XJJNAKfzYVjQQSXJuQRb/PqfOTb5IkpOxTg+LJ+he98IOiTJcnOe38/vr5rInH0uuc7Zuj/okLJLKARPPEHoppuI171KIuFeyrvEuFHEqbp8QyEoLXWPOSnml43FYO/ejsm3/flrr8HTT3ddp7DQtXJTJd+KChg5Mu2/uvgrtqeR02XlDI255Brb3Rh0SB3kXIItHFVK4clWIpEIMzNw6izJPqWlE7hpZ4JIJMJNKlN9U15O6JVXCNfUsGfDBqYsWtS/62Dz812ynDoVFi3q+v7Jk6lbvvX18NRT8LvfdV2ntDR18p02zZ0XLijoW6zim5ETy6A1kbHH+5xLsCKSoUIhmD+fppISpvh9MBw+HGbMcI/O2trgjTdSJ9/6eng+xS3PwuFzg69SnQMeO1aDr6QLJVgRGVzy8ty52okTXSu6s2jUTYyRKvlu3Qrr13ddZ/jw7pPv1KlQXOz7ryWZRwlWRCRZUREY4x6dJRJw5Ejq5FtfD7/5jbv0qLOLL+6afNuXJ07U4KscpQQrIgMikUhQs6eGDY0bOD7qOPOmzCOUbd2qoRCMGeMeV13V9f1YDPbsSZ18X3kFNm/uus6QIS7Rpkq+FRVQUuL/75WlMr1M5VyCjZ86ydFLy7jsQBOHJ4xm1Ku7CQ8bHnRYIoNa49FGFv9oMQ1HGwgT5oEXH6BiVAVP3P4E5aNyaC7i/PxziTGV48dd4k01AGv9enj88a7rjB3b/aVHU6a4fQ5CjUcbuelHN/Ha4ddIkCD0QohLxlzC2tvXZkyZyrn/zNFLyyjd10QIKN3XxBFTzpi9h4MOS2TQSiQSLP7RYnYe2Uksce7a1Z1HdvLWH7+Vl+5+KaNaHb4qKYFZs9yjs7Y22L+/+9HPW7d2XScchrKy7kc/l5bm5OCr9jLVnlwBEiR47fBrGVWmci7Bjtzf1GFmj5H7NBexSJBq9tSw6+iuDskVIJaIYQ9Zpn19GkMLNANTF+O9x1yAUQxpKeHiI61MPtTClEOtTD7cwqTDrUw5tI/JmxsZ9uSTXTZxoiiP18cUsGdsIXvHFvD6mEL2jC1g79hC9pUW0FKQned+T7eepvFo49nk2i5Bgvqmemr21DC/bH5A0Z2TWwm2sRHCeSRibWdn9iCc517XLbFEAlF3pI6CcAFn4inmIgZa4i1KsL1wpjCPhglDaJgwpOubiQSjT8aZfNgl4MmHWs49P9zKDZEThFNM0XtgVD57x7jke/bnWJeED5XkZ2zrtyXe/U0hCvMKqTtSpwSbVt4tscLkEaINcC3YMHnullgvvZSxhUUkl00vnd7tAbEgXMB/v+u/M+JgmNNaW88Nvkrqdp7gPebsTNHTV1TUcaarzt3Qw4Mb27J592ZufPTGlOWqpa2F6aXTA4iqq9xJsN4tsUKd5icNxWJgLcyf76ZcSy4kU6e6ycNFxDfzpsyjYlRFl3Ow+aF8po2exrwpKa5FlfQqKDh33Evl2LHU007W17tbCaa6WcO4cd0n38mTfR18lS1lKncS7PluiQWwZUvq+Ukvuqj7QjJpkhtEICJ9FgqFeOL2J86OIs4nnxgxpo2exhO3P5ERg1EGvZEj4Yor3KOztjbYty918m1ocMfWzvLz3Wm57u58NHp0v3oUz5apH97ExEg9lx7O49UxbRyYVZlRZSp3EmxPt8Rau9b9c7u7QPyZZ1Kv1z45eKpCMkr36RHpjfJR5bx8z8vumsUXN7Bo1qKMu2ZRupGX51qkkyfDwoVd3z99GnbtSp18n37atYA7Gzmy++Nqebm7NrgH5cfg5W+GSNSHaA1DQTxEaBqEbklAhty3IXcSrHdLLHbu7Hgbq/x8909buNB9Yyorg+uv77r+qVPnCkfnQrJxo7vrR2ejR3dfSMrK3N07RARwrY75ZfMpOVrCzLLMm5hd+mjoULj8cvfoLJGAgwdTH1fr6+HFF10LOVko5JJ5quPqtGmu1xFg8WJCO3cSisUY0uqtu3NnRo256XWCNcaUADXALdbaXUmv3wO821p7fdqjuxDeLbESixeTqK+nNRyiIJ4gVFlJqLtbYiUbNgyqqtyjs0Ti3OTgqQrJCy+4zyTLy3MXgXdXSMaNy4gCID3L9NliRDJWKATjx7vH3Lld329pgd27U1/7u2OHa9x0Vlzskuzu3V2Tcyzm1q+pceNuAtarBGuMuRZ4GDCdXr8cWAHUpT+0C9c4Ehbfk2BiJMElh+G1MQn2z2zjiZIE/bpIJxSCCRPc47rrur5/5oy7FKi75PvUU13XGTas+9lZpk513wolcINmBiKRIBQWutN707sZ9Xv0aOrku3171+SavM26uuxJsMBdwN3AD9tfMMYMAf4T+CfgA+kP7cKcnS2mqR47OcZTk93r+U31/s/sMWQIXHqpe3QNDJqaup+d5be/7dil3W7ChO5nZ7n4Yk0OPgDOzUBURyxxbgL3nUfqMmq2mGxx5MBu2qaWU3UGDg2BvF2NlE4oCzosyWSjRsFb3uIeyTZvhhtvTD3upqWl+4Q9wHqVYK21dwKYjneX+ALwfaAh/WFduKyZLWYEMMt7MIRw/BLGH21lytkLxN0MLZMPHWNyZAtjUox8bskP8XppwdmLwveOaf/pLhQ/WayRz+nQ3WwxsUQ8o2aLyRZtU8sZc8Zdnz7mDByuKIfmFLMfiPTEG3PT9qolL6kItYUgb9q01LchDECfBjkZY6qBMmvtJ40x1/f0+draWt+/6W9o3EAe3bfqTkVPkRfPzFZfwzBoGBaGsmKg430jh55pY/KRGGWHXRKeciTGlMMxphxp5apXT7PgpVNdttc0NI89YwrYU5rvHmMK2D3G/dw/Kp9YWK2u3jjV2vVv2y6ffDa8uIGSo/7e6STR+dx+mg1E3WxX5SVXcD9LoxCJRAZk35J7Ch56iLK3v41hredea86H3Q89ROuOHb7vvzd1s6+jiG/8eeWLAAAVyElEQVQF3myM2Q4MByYYY/7XWvuuVB+uqqoi7PP1pMdHHeeBFx9I+V5BuIBf3vrL3GtttLWdG3yV1O08ur6e0Q0NzIy8nnrwVfvk4KkGYI0dq8FXnvPNFhMjxqJZi3wfDRuPx9m+fbtv2x+Iutnu0BDOtmATwJEimDlTo4mlj2bO5Hu3VHD7Yw0MicOZMPzs7dP48OLFA7L73tTNPiVYa+2H2p97Ldj7ukuuAyVbZvZIq7w8d7PmiRNTd4lEo27wVarrfrdudbfH6mz48O4n3pg61Y3gGyTOlamO52DzQ+HcLVM+ytvVyOGKckqjLrnmNTQGHZJkuTt+8id+ePsMKl/Yxc4rp3LHj/xvuV6InLkOVrPFpFBUBMa4R2eJBBw5kjr5NjTA6tUQj3dd7+KLux/9PHFiTg2+UplKr9IJZdCcIBKJqOUqaVE4pJgP/7yOSCTC9RlYpi4owVprp6Z47Sng+vSE0z+aLeYChEIwZox7XHVV1/djMTc5eKrka627zqyzIUPOzXzVOflWVLh7YWYZlSkR6aucacG202wxaZKf75JiRQXccEPX948fd8k21bW/69fD4493XWfs2O4n3pgyxdfJwftDZUpE+iIzj2iS+UpKYNYs9+isrQ327+/+7hxbt3ZdJxw+N/gq1Tng0lINvhKRrKIEK+mXl+fuRDRpUurZVJqbzw2+6px8t2yBJ5/sus6IEd1PvFFe7s43i4hkECVYGXjFxfCmN7lHZ4kEHD7c/V2Pfv3r1IOvJk1KnXwrKtysWP0ZfJVIQE0NozdscF3j8+apNS0iPVKClcwSCrlztWPHwtVXd32/tdUNvkqVfF96CTZt6rpOUVHHwVbJybeiwrWOu9PYCIsXQ0MDF4fD8MADbp0nnnAtZxGRbijBSnYpKDiXIFM5dqz72w7+/vep5y4dNy71ed+pU+Ftbzt7C8Sz0zFk2C2xRCQzKcFKbhk5Eq64wj06a2uDffu6v+3gli2920eG3RJLRDKTEqwMHnl57kbOkyfDwoVd3z99GnbtOpdwf/Mbd8lRqnO+GXRLLBHJTEqwIu2GDoXLL3cPcLfI+sMfUifYDLollohkptyZ104k3bxbYnWZACM/352jzZBbYolIZlKCFelOKORGC1dWQmEh8aFDXdfw9OnudQ1wEpHzUBexyPmUl8PLL0NNDfs2bGDKokW6DlZEekUJVqQnoRDMn09TSQlTMvCOHSKSmdRFLCIi4oPcS7DxOKxcScWyZbByZeoRoCIXoPnkMTaZYiZfM4tNppjmk8eCDikrxaPN1M+ppPzaK6mfU0k82hx0SJLtTp6E0lJmzJrlbghy8mTQEXWQe13EDz4I99/PiGgUIhHXvbdiRdBRSRZ7fvYE5r8aJQTMfzXKptkTWWhPBx1W1mmcX0XFtnpCQMm2euoXzKBya13QYUk2KyuDpiZCAE1NbszE4cNBR3VW7rVg16+HaNQ9j0ZT35lF5AK8ebdLrgAhoGq3Wl59Mfbl3R3+juNeagwyHMkFTU0dl48cCSaObuRegr3hhnO3LisqghtvDDYeyXp/Kisi4T1PALVlxUGGk7UOXVbW4e948HLdLEH6afTojsulpcHE0Y3c6yJevhxCIU6sWsWIpUvh3nuDjkiy3JxtB9g0eyJVu5upLSvmqm37gw4pK5VvrqV+wQzGvrSLQ5dPZeqmHUGHJNlu924oLydx5Aih0lJ396sMknsJNhyGFStoWLKEmbqkQtKgePhIFtrTRCIRFqpM9Vm4qJjKrXVEIhHVTUmP4cPh8GF2ZGiZyr0uYhERkQygBCsiIuIDJVgREREfKMGKiIj4QAlWRETEB0qwIiIiPlCCFRER8YESrIiIiA+UYEVERHygBCsiIuIDJVgREREfKMGKiIj4QAlWRETEB72+m44xpgSoAW6x1u4yxtwGfNp7eyfwIWttU7cbEBERGUR61YI1xlyLS67GW54EfAm40Vo7E3gJuM+nGEVERLJOb7uI7wLuBvZ5y23AR621B73lPwJlaY6tTw4eP0jev+Qx61ezyPuXPA4eP9jzSiLnET91ksOTSrnsLbM4PKmU+KmTQYeUleJtcVZuWsmyzctYuWkl8bZ40CFJlos3HaFlSD4zZs2iZUg+8aYjQYfUQa+6iK21dwIYY9qX9wNrvNeGAiuAr/sT4oW56KsXkSABQIIEF331Itr+uS3gqCSbHb20jNJ9TYSA0n1NHDHljNl7OOiwss6DNQ9y/8b7icaiRDZGCBFixYIVQYclWSw+YTwFLXFCQEFLnNYJFxE+0xp0WGf1+hxsKsaYMcAq4AVr7Q+6+1xtbS2hUKg/u+q19uSavByJRAZk35KbLjvgkitACCjZf2TAylQikej5Q/0wkHXzsR2PEY1FAYjGoqzasYolI5cMyL4lN83wkivgJdlYRtXNPidYY0w58ATwGOcGO6VUVVVFOBzu664uSOhXoQ5JNkSImTNnDsi+JTcdnjD6bAs2ARyfWDpgZSoej7N9+3bftj+QdfMdx95BZGOEaCxKUX4RS2csVd2UfmkpDJ9twSaA1sL8jKqbfUqwxpgiXHL9trX2ob5swy9vfOKNs93EIUK88Yk3gg5JstyoV3dzxJRTsv8IxyeWMso2Bh1SVlo+bzkhQqzasYqlM5Zy77x7gw5Jslz4wJ9pnXARBS0xWgvzCR/IrON9X1uwdwDTgQ8aYz7ovfbH9nO1QRpXMo62f24jEono27GkRXjYcMbsPawy1U/hvDArFqxgycgl+jtKWoRHlxI+05qxdfOCEqy1dqr39DveQ0RERFLQTE4iIiI+UIIVERHxgRKsiIiID5RgRUREfKAEKyIi4gMlWBERER8owYqIiPhACVZERMQHSrAiIiI+UIIVERHxgRKsiIiID5RgRUREfKAEKyIi4gMlWBERER/kXIKNt7awYVk10ffMZ8OyauKtLUGHJNnu5EkoLWXGrFlQWuqW5cI1N0NlJVVXXgmVlW5ZpB9ORk9S+sVSZv1qFqVfLOVkNLPqZs4l2I13L2HuI+u42p5g7iPr2HjPzUGHJNmurAyamggBNDVBeXnQEWWnqiqorycvHof6epgxI+iIJMuVPVRG05kmAJrONFH+tcyqmzmXYEc+vY3imHteHIORNc8HG5Bkv6amjstHjgQTR7bbvbvjcmNjMHFIzmhPru2ORDOrbuZcgj123Wya893z5nw4Nm9OsAFJ9hs9uuNyaWkwcWS7srKOy+oJkH4aPaRj3Swtyqy6mXMJduG31vDsndU8Z0bw7J3VLPzm6qBDkmy3ezeUlpIAl1zV8uqb2lqorKQtHHbnYHfsCDoiyXK7P777bFItLSql8e8zq27mBx1AuoULCln08FoikQhXz5wZdDiSC4YPh8OH2RGJMFNlqu+Ki6Gujlr9HSVNhhcN5/Dyw0QytEzlXAtWREQkEyjBioiI+EAJVkRExAdKsCIiIj5QghUREfGBEqyIiIgPlGBFRER8oAQrIiLiAyVYERERHyjBioiI+EAJVkRExAdKsCIiIj5QghUREfFBr++mY4wpAWqAW6y1u4wxfwF8DRgKPAn8nbU25k+YIiIi2aVXLVhjzLW45Gq85SHAo8B7vNeGAh/2KUYREZGs09su4ruAu4F93vLVQL219lVrbQL4HvA3PsR3weKtLWxYVk30PfPZsKyaeGtL0CFJlms5cYy944u57IpZ7B1fTMuJY0GHlJVUNyXdWmItVD9azfzV86l+tJqWWGaVqV4lWGvtndbaTUkvTQJeT1reB0xOZ2B9tfHuJcx9ZB1X2xPMfWQdG++5OeiQJMv9uXICkw5GKUjApINR3pg+MeiQspLqpqTbkp8sYV3DOk60nmBdwzpu/mlmlalen4PtJFVibuvuw7W1tYRCoT7u6sIM27SFYu9McHEMhm18lkgkMiD7ltx02aEo7aU3BEw42DxgZSqRSPi6fdVNyWZb9m7psPzsnoErU72pm31NsHuB5K/xE73XUqqqqiIcDvdxVxdmw4JraN65juIYNOfDqYVzuXrmzAHZt+SmvWOLmHTQJdkEcGBcMTMHqEzF43G2b9/u2/ZVNyWbXbP9GtY1rDu7PHfK3Iyqm329TGcLcKkx5k3e8p3Amj5uK60WfmsNz95ZzXNmBM/eWc3Cb64OOiTJcuN3HmDv+GJaQ7B3fDEX1e0POqSspLop6bbmtjVUT6tmRMEIqqdVs/rWzCpTfWrBWmvPGGNuB35qjCkGngO+mdbI+ihcUMiih9cSiUT07VjSonDESKa8cZpIJDJg345zkeqmpFthfiFr71ibsXXzghKstXZq0vM/AG9Jd0AiIiK5QDM5iYiI+EAJVkRExAdKsCIiIj5QghUREfGBEqyIiIgPlGBFRER8oAQrIiLiAyVYERERHyjBioiI+EAJVkRExAdKsCIiIj5QghUREfGBEqyIiIgPlGBFRER8oAQr0oN4awsbllUTfc98NiyrJt7aEnRIWSneFmflppUs27yMlZtWEm+LBx2SZLlML1N9uuG6yGCy8e4lzH1kHcUxaN65jo2hm1n08Nqgw8o6D9Y8yP0b7ycaixLZGCFEiBULVgQdlmSxTC9TasGK9GDk09sojrnnxTEYWfN8sAFlqfUN64nGogBEY1GebHgy4Igk22V6mVKCFenBsetm0+z19TTnw7F5c4INKEvdUHEDRflFABTlF3FjxY0BRyTZLtPLlLqIRXqw8Ftr2Bi6mWEbn+XUwrks/ObqoEPKSsvnLSdEiFU7VrF0xlLunXdv0CFJlsv0MqUEK9KDcEEhix5eSyQS4eqZM4MOJ2uF88KsWLCCJSOXMFN/R0mDTC9T6iIWERHxgRKsiIiID5RgRUREfKAEKyIi4gMlWBERER8owYqIiPhACVZERMQHSrAiIiI+UIIVERHxgRKsiIiID5RgRUREfKAEKyIi4gMlWBERER/0+246xpjbgfZbyD9urf2H/m5TREQk2/WrBWuMGQp8HVgEzAIWGGP+Mh2BiYiIZLP+dhGHvW0U41rD+UBzf4MSySTxtjgrN61k2eZlrNy0knhbPOiQRASgpQWqq7l8/nyornbLGaRfXcTW2hPGmH8CXgFOA38Anu78uXh84A9IiUQikP1K7vny01/mKzVf4Uz8DHU1deSTz6eu+9SA7NvvMqy6KVntne+EZ54hBMSfeQb++q9h1aoB2XVvynAokUj0eQfGmJnAD4DFwDHgR8Bz1tp/A9i2bdskYG+fdyAi7SbPnj379XRtTHVTJG26rZv9HeS0GHjSWvtnAGPMI8DdwL957+8DJgMn+rkfkcFsBK4upZPqpkj/nbdu9rcFexPwFeBa4BTw/wF/ttZ+rs8bFRERyQH9GuRkrV2L6xbeBkSAIcAX0xCXiIhIVutXCzYXGWM+Dzxvrf110LFI5jLGbLfWXhF0HIOJ6qb0RibVTSVYERERH/R7JqdUjDHXA5/FjSy+DDda8b3AdcADuK7peuCj1to3jDG7gEeBm4AxwMestU+k2O67gU/hrrsdCtxlrd1gjHkK+CNwIzAMWA58CKgCvmat/YoxZhjwH8AVuOt3/8Na+11jzAeBDwKlwFpgLPCUtfYRY8z/xQ3aagPWAX8PvBn4hrefi4CvWmu/aoy5D5gEVAJTgd9aaz/Wrz+k9IkxJh83HmAGMB54DfgE8DPA4srAK8Bq4AO4MvdOa+2fjDGzgYdw/9+jwN3W2le8MtaEK88fBJ6x1oaMMaOA7+LKWgvwaWvtb40xHwPejyunAO+11tb2tqz7RXVTdTNIg61u+jkX8VzgE9bay3GTT/wd8DDujzUTqMFVhnbHrLVzgXuBlZ03ZozJA/4WeLu1dpb3mU8mfSTsbfdR3HngdwELcAcTgM8AtdbaK3EHk7u9y4wAyoDZydM8ev/Mv8cN4KoCKoC/AD4CfMFae5W3/PmkGK4A3gbMBP7KGDOjl38rSa/rgLhXniqBQtz/ZQaubBhcZZxqrb0O+AnwEWNMAfB94A6vnPwTboxBuz9Za99krX026bXPA3ustW8C3gl8wRgz0nt+vbW2CvgF8NGkdc5b1geA6qbqZlAGVd30pQXrqbXW7vaeb/d+PmetrfeeP8y5OYwB1iR9trTzxqy1bcaYpcAtxhgDXA8kX+m72vvZCDxrrT0NNHrfYsB9KxlmjHm/t1yCq2zgzuu0dtrl9cCvrbVN3vItAMaYTcBbjTErvPWHJ63zpLX2DHDGGFOX6vcQ/1lrNxpjDhtj7gHeBFyK+z/tt9ZuAzDG7AWe9FZpxB2kDTAdWOWKGADjjDGF3vMuk6jgpgm93dtvHW7KUIwx7wX+xhhzKfBWztUB6KGsDwDVTdXNQAy2uulngo0mPU8A7wB2J70WAgpSfD7hvYcxJvkXnw88D/wQ2IgbtZzczZM8R1YsRTxh4H3W2he8bY/HdZPdipuFqrOYFwve5yd4+/gOcBxYhevWeG+K36HD7yEDyxjzdlx351eB/8J1LTbSsYxA13ISBna2D5AwxoSASdbaFq9S96acGO9zG3Hdno8DB4C3JK3TpawPMNVN1c1ADLa6OZC3q9sCXGOMmeYt34WbWrFb1tor2h+4bzptwL8C63HdCuEL2P964G+NMSFjzDjgBVxXRHc2Am8zxozwusC+j/tGVA18xlr7GK4bCmPMhcQh/vtL4GfW2v/CVaCF9K6svAKUGmMWesu3Ab/pYZ2n8A7kxpipwAZcN9hr1tp/x5X7Cy2rA011UwbKoKqbA5lg38BV3F8aY/4E3EDHvu+evIgbLPEKrgKeBKZ632R641+AImAHsAl4wFq7vbsPW2v/CHwN1/WwA3jJWvsL4D5gszHmBVz3QgMwrbvtSCC+A9xqjPkj8EvgGVx5Oy+vC/HdwL8ZY9pbYe/rYbX7gGnGmBdxLacP4LqZQsaYl4BngV1kdhlR3ZSBMqjqpi7TERER8cFAtmBFREQGDSVYERERHyjBioiI+EAJVkRExAdKsCIiIj5Qgh1kjDFzjDHfDToOEelIdTP36DIdERERH/g5VaIExLsQ+3ZrbcQY8z1gqLX2VmPMdOC/gRPW2uu9u1A8i5t4/WLgPmvtD7q7u0kgv4xIDlHdHFzURZyb1uCmJAOYA1zjPX8b0Plm1cOstfNwd5j4d++1893dRET6TnVzEFGCzU1rgL/05pbdDew3xlyCq8QnUnwWOt494iZgmTehew0winN3NxGRvlPdHESUYHPTM7j7ZC7GTXD9FK4CT8fNFZssCmCtTT4Z3353k/bJ3OcCP/c5ZpHBQHVzEFGCzUHW2jhu0vRP4irxemA55+6x2JMLvbuJiPSC6ubgogSbu9YAY3B3Oqnxnvd0e6d2F3R3ExG5IKqbg4Qu0xEREfGBWrAiIiI+UIIVERHxgRKsiIiID5RgRUREfKAEKyIi4gMlWBERER8owYqIiPhACVZERMQHSrAiIiI++P8BexqzxaSw9LYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 480x320 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmp = agg_prediction(m_05.trace['score'], wines,\n",
    "                     ['wine_amer', 'judge_amer', 'flight_red'], summarize())\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(6, 4), dpi=80, sharey=True)\n",
    "\n",
    "for j_amer, title in zip((0, 1), ('non-american', 'american')):\n",
    "    ax = axes[j_amer]\n",
    "\n",
    "    for flight, c in zip((0, 1), ('g', 'r')):\n",
    "        sub = tmp[(tmp.judge_amer == j_amer) & (tmp.flight_red == flight)]\n",
    "        data_sub = wines[(wines.judge_amer == j_amer) & (wines.flight_red == flight)]\n",
    "\n",
    "        ax.scatter(data_sub.wine_amer, data_sub.score,\n",
    "                   c=c, marker='.')\n",
    "\n",
    "        ax.plot(sub.wine_amer, sub.mean_, c=c, marker='o')\n",
    "        \n",
    "        ax.set(xticks=[0, 1], xticklabels=['non-american', 'american'],\n",
    "               xlabel='wine', title=f'{title} judge',\n",
    "               xlim=(-.2, 1.2), ylim=(8, 20))\n",
    "        \n",
    "fig.tight_layout()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* American judges give better scores\n",
    "* Non-american red wines score better than non-american white wines\n",
    "* American white wines score better than American red wines\n",
    "* White whines score about the same regardless of origin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOgAAAE4CAYAAABR4L/jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuclGX9//HXMMuqgIBYhqKAin5AlpMbp0B+nkhN8VB+VfKQmpliZmqZlkeyqO8vNTMtD+VZs9IipcJQPLAKJQrLqnxKBARBMxREwJZd5vvHdaPDNsAwyzDX7ryfjwcPZu657pnP7t7vua+57+u+JpXJZBCROLUpdQEisnEKqEjEFFCRiCmgIhFTQEUipoCKREwBFYlYRb4NzawjUAOMcfcFZjYcuAHYEagFvuTu9U3W6QbcB3QFlgInufu/tlbxIq1dXnvQJIw1gCX3OwKPAGe7e9+k2VdyrHozcJe79yEE9cZmVyxSRvLt4p4NjAOWJPdHA8+7e21y/3xCYD9iZm2Bg4AHkkX3AEcmy0UkD3l1cd39DAAzW7+oF7DSzB4G9gGmARc1WW1nYKW7r02eo8HM3gc+ycdBB2DmzJkpYDdgZWE/hkiLtSOwpLq6OueY27w/g+ZY73PAMGAB8EvgUuDqrDYb2zuvy7FsN2BxgbWItHS7A2/meqDQgL4F/M3d5wGY2W+ArzVp8w7Q0cwqkr1nBeHdYlmO51sJ0K9fP9LpdIElbd5z1bsw/LX/kAIywHO9tmfEzLeL9nrSsqzr1IkUfLR9ZIA2K1YU7fUaGxuZM2cObKLnWGhAHwe+Z2Y93H0hYW/6YnYDd19rZk8DJwN3J/8/tb7Lm0s6nS5qQD897XWeq96VqjfWUNd9BwZPm1fU15OW5Z0ZT7NL1ac/Cui/6l6ga4m3j4LOg7r7IsJR2z+a2VxgF2ACgJndYWZHJ03HAaea2cvAWYSDSSWzQ4dOjPLVLJ4xm1G+mh06dCplORKZrn2raZPJMGf2bNpkMnTtW13qkrZsD+ruPbNuTwIm5WhzVtbtRcChzahPpKxpJJFIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxBRQkYgpoCIRq8i3oZl1BGqAMe6+wMwmAGOB5UmT29395ibrGPB34PVk0dvufljzyxYpD3kF1MyGA7cBlrV4CHCcu7+0iVWHAHe6+wWFlyhSvvLt4p4NjAOWAJhZChgEXGlmtWZ2o5ltl2O9wcD+ZjbTzJ4ws75bpWqRImhcW8/Ur4zmwxNGMvUro2lcW1/qkvILqLuf4e7PZi3aGXgOuBDYH/gE8N0cq64Gfunu1cB1wEQza9u8kkWK45lxRzLsrikM8ZUMu2sKz5x3VKlLyv8zaDZ3/zfwUfVm9mPgl8CVTdpdmnX7T8nn1j5Aba7nraurI5VKFVLSFlm7di21tTlLkDLW/tkZ7NAQbu/QAO2fmV7U7SSTyWy2TUEBNbN9gKHufl+yqA3QkKPdt4BfuPvKTbVbr6qqinQ6XUhJW6S2tpb+/fsX/XWkZZl6wFDWzJvCDg2wpgJWjRrGkCJuJ42NjcyaNWuTbQoKKFAPXG9mTwFvAucDv8/R7lDgQ+AmMzsYSANzC3xNkaIadcsknkkdRftnprNq1DBG3fxYqUsq7Dyouy8kfP6cDDjQSPiMiZmNN7NzkqbnAJ83szrgh8CJ7r6u2VWLFEG6bSUH3fY42/9mGgfd9jjptpWlLmnL9qDu3jPr9v3A/TnaXJl1ez5wUDPqEylrGkkkEjEFVCRiCqhIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxBRQkYgpoCIRU0BFIqaAikRMARWJmAIqEjEFVCRiCqhIxCrybWhmHYEaYIy7LzCzCcBYYHnS5HZ3vznHOvcC+wCrgC+6+z+3SuUiZSCvgJrZcOA2wLIWDwGOc/eXNrHqeOBFdz/GzA4G7gY+U2ixIuUm3y7u2cA4YAmAmaWAQcCVZlZrZjea2XY51hsD3AXg7k8CXc2se7OrLtAHb7/JqrYp+g0YwKq2KT54+81SlSIRmv/y32hIhe2jIZVi/st/K3VJ+QXU3c9w92ezFu0MPAdcCOwPfAL4bo5VuwHZKVgC7F5Yqc2X2n132jVACsL/u5esFInQHlVDSRO2j3Ryv9Ty/gyazd3/DRy1/r6Z/Rj4JXBlk6a53gDWbex56+rqSKVShZSUl35JOOHjkNbW1hbt9aRl6ceG20ea4m4fmUxms20KCqiZ7QMMdff7kkVtgIYcTRcDuwKLkvu7JstyqqqqIp1OF1JSXlZV8NEeNAOsroD+/fsX7fWkZWmAj/agGaCR4m4fjY2NzJo1a5NtCj3NUg9cb2a7J59Hzwd+n6PdJOBMADM7EPjA3Tca0GLLLF7M6oqPw5lZXLJSJEKL6mbQyMfhXFQ3o8QVFRhQd19I+Pw5GXDCz3MdgJmNN7NzkqZXAP3NrC55/LRmV9wMHT7VjfZrM8yZPZv2azN0+FS3UpYjkdmz7xAqMmH7qMhk2LPvkFKXtGVdXHfvmXX7fuD+HG2uzLq9HPhCM+oTKWsaSSQSMQVUJGIKqEjEFFCRiCmgIhFTQEUipoCKREwBFYmYAioSMQVUJGIKqEjEFFCRiCmgIhFTQEUipoCKRKygKU9arEwGamrYaepUeP99GDECijgHkkhzlU9AFy6Eww6D+fPZLZ2Ga6+FPfeEyZOhR49SVyeSU3l0cTOZEM5586C+nvSaNVBfH+4ffnh4XCRC5RHQmhpYsAAamkw82NAAr78eHheJUHkE9LXXoG3b3I+1bRseF4lQeQS0V6/Qpc1l1Sp4551tW49InsojoCNGhANCFU2OibVpE47iXnIJHHdcOJAkEpHyCGgqFY7W7r03VFbS2K4dVFbCvvvC88/DCSfAH/4AffrAhAkb39uKbGPlEVAIp1JefRWeeIIll14KTzwBr7wCQ4fCQw/B44/DHnvAd74D/fvDlCmlrlikjAIKYU86ciTvHXMMjBy54SCF0aOhtha+/314441w/8QT4U19RaGUTnkFdHO22y7sQV95BY49Fn7zG+jdG667DtauLXV1UoYU0Fx69oTf/x4eewx22QW++U0YNAiefrrUlUmZUUA35cgjoa4OrroqnCs98EA49VR4661SVyZlQgHdnB12gKuvDkE94gi47z4wg5tu+u+RSSJbmQKar169YNKk0PXt3Bm+/nUYPDicphEpEgV0S6RS4eDRK6/AZZfByy/DZz4DX/6yRiNJUSighWjfHn7wg3Ba5pBD4Fe/Ct3eW2+FxsZSVyetiALaHL17w1//Cr/+dfises45MHw4vPBCqSuTVkIBba5UKgxomDsXLr4YXnwRhgyBc8+Fd98tdXXSwimgW8uOO8KPfwyzZoVRSr/4Rej23nknrFtX6uqkhco7oGbW0czmmFnPJsvPM7OnNrKOmdn7ZjYr+Te5eeW2AFVVYUDDPfeEq2XOPBMOOABmzy51ZdIC5RVQMxsO1ADWZPl+wGWbWHUIcKe7D0z+HVZwpS1JKhUGNLjD+efD9Omw//7wjW/AihWlrk5akHz3oGcD44Al6xeY2XbArcAVm1hvMLC/mc00syfMrG/BlW4Fjas+YFm3LvQZNIBl3brQuOqD4r5g587w05+Gg0ZDhsCNN4YDS/ffr3mQIlTfUM/oe0Yz8rGRjL5nNPUNpb/sMK+AuvsZ7v5sk8UTgF8B8zex6mrgl+5eDVwHTDSzjcw9UnzL9+1OlyXv0XYddFnyHsttG83mN2hQmPfojjvCoPtTToGDDw7nUyUaRz5wJFPmT2Hl2pVMmT+Fox48qtQlFTbtppmNBrq7+0VmduDG2rn7pVm3/2RmE4A+QG2u9nV1daSKOE9tn7feY/2zp4COS9+ltjZnKcUxeDDpRx6h60030eXhh2HAAN455RT+9dWvsq5du21Xh+Q0Y/GMDe5PXzS9qNtHJo9eVKHz4o4F+prZLKAD0NXMfufux2c3MrNvAb9w95XJojbARgewVlVVkU6nCyxp85Z13YkuS0JIM8D7u3ahf//+RXu9jRo1Cv72Nxg3jl3uuotd/vpXuOEGOP54TaRdQkNnDWXK/I8v1B+2x7Cibh+NjY3MmjVrk20KOs3i7me6ex93HwicBbzQNJyJQ4HTAczsYCANzC3kNbeGzv94g3e7dWFtG3i3Wxc6ewnnIBoyBGbMgFtuCROXnXBCmLv3H/8oXU1lbtIXJzF6r9Hs2HZHRu81msfGPlbqkrb+eVAzG29m5yR3zwE+b2Z1wA+BE929ZCcF0+07sPPiZbz60mx2XryMdPsOpSolKSgdBjS4w+mnh1FJ/frB5ZfD6tWlra0MVVZU8vipjzPtqGk8furjVFZUlrokUvn0g4tt5syZHYEVAwcOLGoXd73a2trSdG03p6YGxo0LY3x79AhHfY8+Wt3ebWxbbR9ZXdxO1dXV7+dqo5FEMRkxAmbOhJ/8JAwTPPZYGDMmzH4vZUkBjU1FBVxwQej2fvGL4RrUvn1h/Hj48MNSVyfbmAIaq113DQMannwyTLp91VVhGOFf/lLqymQbUkBjd9BBYQD+//5vmAvpiCPgC18IU4NKq6eAtgSVlfCtb4WJt48/Hh55JMyC/8Mfahb8Vk4BbUn22AN++9vQze3WLUy7MmBA6AZLq6SAtkSHHQZz5oRvCV+wIEy7MnasZsFvhRTQlmq77eC73w3d3qOPDtOu9O4N11+vWfBbEQW0pevZEyZOhEcfhU9+Mky7sv/+8GzTi4+kJVJAW4ujjgrTgF55ZRjPO2oUfOlL8Pbbpa5MmkEBbU122AGuuSYE9fDDw7QrZvCzn2k60BZKAW2NevWCP/0JHn4YOnYM064MHhymXpEWRQFtrVIp+Pznw0Gkb387HPUdPhy+8hX4979LXZ3kqawCmslkmPbGNCYunMi0N6bldUV7i9e+fRjQUFsbplm5447Q7b3tNk0H2kSM20fZXG62cPlCDrvvMOYvn0+aNI00smfnPZl8ymR6dN5GcxOVWiYDDz0EF10ES5eGi8ZvuQWqq0tdWcmVYvvQ5WaJTCbDYfcdxrx351HfWM+axjXUN9Yz7915HH7/4VG8U24TqRScdFKYBf/CC8OlbYMHw3nnwXvvlbq6kol5+yh0TqIWpWZRDQuWL6Ahs+F0SA2ZBvzfzl4/3Yt2bcts0q5esM+3e3DFQ0upvuUWlt1zK9cd9ykmDu1UdheIr167moXLF5JhwyA2ZBp4/b3XqVlUw8juI0tSW1nsQV979zXapjc+22d9Y3kOOP9nt+057cKeXHbabmRS8IN7l3DPDQvY983yuu50U3//yjaVvPbua9uwmg2VxR60V5deG/0jtE235aHjHyrZO2QUzgNuXA6XX071z3/O73+0MJyaueaacJqmlZv2xjQOueeQnNtI/bp6enXpVYKqgrLYg47YYwR7dt6TitSG70cVqQr22mkvRuwxokSVRaRz5zCg4e9/h09/Oky70rs3PPhgq58FP+btoywCmkqlmHzKZPbusjeV6UrapdtRma6k1869mHzK5KJOlt3i7L8/PPcc3H47/Oc/YdqVQw5p1bPgx7x9lM1pFghH62oW1TB19lQOGnAQI/YYoXBuyrJl4ZrT228PcyVddBFccQV0KPF0pUWyrbcPnWZpIpVKMbL7SI7pcQwju49UODdn553DgIbp06F//zDtSp8+YQhhBG/sW1uM20dZBVQKNHRo+KqKm2+GlSvDtCtHHAH//GepK2v1FFDJTzodJtX+xz/CZWyTJ4dZBq+4QrPgF5ECKltml13grrvgmWfCmN5rrw3z9j76aKkra5UUUCnMAQfAiy+GKVaWLQvTrhx9NMzf1NfFypZSQKVwFRVhTO/cuWHSskcfhf32g+99T7PgbyUKqDTfbrvBAw/AE0+EOZKuvDJ8S9vkyaWurMVTQGXrOfhgmD07XH+6ZEmYduX442HRolJX1mIpoLJ1VVaGGRxefTV8RcXDD4chgz/6kWbBL4ACKsXRvTv87nfw5z+HLvCll8LAgTB1aqkra1EUUCmuww8P8yGNHx+O8B58cBjfu3RpqStrERRQKb7ttw8DGl55Jczf++CD4RzqT34CDQ2bX7+MKaCy7ey5ZzgV88c/hnG+F14Yrp6ZNq3UlUUr74CaWUczm2NmPZssP8/MntrEOhPN7BUz+7uZ7dO8cqVVGDMmTK59+eXhm8QPOABOPx3+9a9SVxadvAJqZsOBGsCaLN8PuGwTq44HXnT3/YBvA3cXWKe0Nu3ahQENc+bAZz8Ld98dur233KJZ8LPkuwc9GxgHLFm/wMy2A24FrtjEemOAuwDc/Umgq5l1L6jSreCthc6qihT9BgxgVUWKtxZ6qUqR9fbdN3zf6e9+F64zPe+8MB3ojBnbvJQVS9/YYPtYsbT032KeV0Dd/Qx3b/p1WROAXwGbGnzZDcj+0solwO5bVOFWtOPevWnXCCmgXSN02Lt3qUqRbKlUOGf66qtwySVhku3hw+Hss8M4322kYo8eG2wfFd1LP19yQZOGmdlooLu7X2RmB26iaa43gI1OZ15XV1fUi2T7Jb98CP+3b4Ta2tqivZ4U4OST2W7YMLr94Ad0uP12Gn77W9664ALePfZYaFPcY5pNt492DcXdPvKZzaTQWf3GAn3NbBbQgdB1/Z27H9+k3WJgV2D9WK9dk2U5VVVVFXXKk1VpPnqHzCT3+/fvX7TXkwL17w/HHgsPPkjFxRez+zXXsPvkyeHz6aBBRXvZptvH6oribh9ZU55sVEFvSe5+prv3cfeBwFnACznCCTAJOBMg2dN+4O4bDWixrZw3l1Xpj8P5wby5pSpFNieVCgMa3OEb3/h4tsHzz4fly4vykg2LFrK64uNwNryxsCivsyW2ep/BzMab2TnJ3SuA/mZWB1wHnLa1X29LdO1hdGjIMGf2bDo0ZOjawza/kpRWx45www3hayqGDw9Tg5qF7z7dyvMiddq1O+3Xhu2j/doMnXYt2fHMj5TVrH7r1dbWqmvbEq1bF4J5ySXwzjvh/OnNN4dL27aibbV9aFY/aV3atAkDGtzD/EjTpoXPpBdfDO/n3L5bPAVUWp6ddgp7zr//PXx14vXXh0vafv3rVjcdqAIqLVd1NTz/PNx6a5gFf+xYOPTQMAVLK6GASsvWpk0Y0OAOX/4yPPlkOE1z2WWwalWpq2s2BVRah098Au64I+xRq6rCtCt9+sAjj7Tobq8CKq3LsGHhs+lNN4UDR1/4Anzuc/Ba6b7jszkUUGl90mn42tdCt/fUU8Ng/KoquOoqWLOm1NVtEQVUWq9PfSqcN336adhnnzDtSt++8Nhjpa4sbwqotH6jRoVZ8K+7LgxwGDMGjjkGFiwodWWbpYBKeWjbNny/6dy5cOKJYdqV/faD738/nKKBcDBp2jR2mjgxDIKI4OCSAirlpVu3MKBhypQwNejll4ehgvfeG476HnIIu02YEL5VvE8fWFjaAfMKqJSnQw4JF4ZPmACLF8Npp4WvVqyvJ71mTZhke968MG1oCfekCqiUr8rKMKH2vfeGAQ9Ng9jQAK+/DjU1pakPBVQkfGt4u3a5H6usLOk5VAVUpFevjX9vTH19eLxEFFCRESPCpNoVTWYAqqiAvfYKj5eIAiqSSoXvMt17b6ispLFdu9C17dUrLC/iRHabU+ikYSKtS48eYdrPmhqWTJ3KHgcdFPacJQwnKKAiH0ulYORI3uvYkT0imRJHXVyRiCmgIhFTQEUipoCKREwBFYmYAioSMQVUJGIKqEjEFFCRiCmgIhFTQEUipoCKREwBFYmYAioSMQVUJGIKqEjE8r5g28w6AjXAGHdfYGYXAF9NHp4EXOLumSbrHAbcDyxOFr3k7mc0v2yR8pBXQM1sOHAbYMn9vsA4YBDwH+AZYDTweJNVhwDfd/cbtlbBIuUk3y7u2YRALgFw95eBKndfDXQGOgHLc6w3GDjSzF4ys4lmtvtWqFmkbOQVUHc/w92fbbJsrZmdC8wHlgKzcqz6HvBjdx8E/AV4oJn1ipSVVGYLvnfCzBYAB7r7gqxlFcDdwEJ3/85m1l8O9HD3FdnLZ86c2RFYUVFRQWobzKK2du1a2rZtW/TXkZZpW20fmUyGhoYGgE7V1dXv52pT0Kx+ZtYT6Oru0929wcweAM5t0mZ74EJ3n5C1OAWs3djzVlVVkU6nCylpi9TW1tI/klnbJD7bavtobGxk1qxcHc+PFXqa5RPAfWa2o5m1AU4gHCj6iLt/CJxuZmMAzOwMYHryuVVE8lBQQN39BeBGYDrhs+f7wA0AZnaHmR2dNB0LXG5mLwOnAmc1u2KRMrJFXVx375l1+ybgphxtzsq6/SIwtBn1iZQ1jSQSiZgCKhIxBVQkYgqoSMQUUJGIKaAiEVNARSKmgIpETAEViZgCKhIxBVQkYgqoSMQUUJGIKaAiEVNARSKmgIpETAEViZgCKhIxBVQkYgqoSMQUUJGIKaAiEVNARSKmgIpETAEViZgCKhIxBVQkYgqoSMQUUJGIKaAiEVNARSKmgIpETAEViZgCKhIxBVQkYhX5NjSzjkANMMbdF5jZBcBXk4cnAZe4e6bJOt2A+4CuwFLgJHf/11apXKQM5LUHNbPhhHBacr8vMA74NNAP+AwwOseqNwN3uXsfQlBv3Ao1i5SNfLu4ZxMCuQTA3V8Gqtx9NdAZ6AQsz17BzNoCBwEPJIvuAY5MlpdE46oPWNatC30GDWBZty40rvqgVKVIhNa89Sb1bVL0GzCA+jYp1rz1ZqlLyq+L6+5nAJhZ9rK1ZnYu8CNgBjCryWo7AyvdfW3SvsHM3gc+SRL0phobG7e0/i3y7oDe7LSintQO7em0op53B/ahy9wFRX1NaUH2Mdq0a886wp6rcR+jcfmKor1cPtt73p9Bc3H3n5vZ7cDdwNXAd7Ie3tjeeV2OZTsCzJkzpznlbN5DE3mjyaI3ZjV9X5Gy9dTT/71s22wfOwLv53qgoICaWU+gq7tPT/aMDwDnNmn2DtDRzCqSNhVJIctyPOUSYHdgZSH1iLRgO7KRHiUUvgf9BHCfmQ0CVgEnAM9kN0i6wE8DJxP2sCcDT63v8marrq7OAKXv8Itsezn3nOsVdB7U3V8gHJGdTvjs+T5wA4CZ3WFmRydNxwGnmtnLwFnA+YW8nki5SmUymc23EgDMbDzwgrv/sdS1SPGZ2Sx3H1jKGhRQkYg16yhuIczsQOByYAXQB1gMnEQY7HAtodv9OvBVd3/bzBYQzqF+lnDq5mvuPjnH8/4PcDGwA9AOONvdp5rZU8BLwCFAe+DbwJlAFXCju19nZu2Bm4CBQBq4yd3vMLPTgdOBLsDjhM/eT7n7XWb2dUIXfh0wBbgA6Av8LHmdTwE3uPsNZnY10A3YG+gJ/Mndv9asX2SZSw46/pwwUGYX4J/AhcCvASf8LecCjwFfImw7n3f3l82sGvgJ4e+0HBjn7nOTbeU9wnZ5OvC8u6fMrDNwB2GbqQcudfc/mdnXgNMI2xuEkXJ1+W6z+SjVWNxhwIXuvh+whvDZ9DbCL7A/YdTSz7Lar3D3YcAlwA+aPpmZtSEcRT7a3QckbS7KapJOnvce4IfA8cABhDcKgO8Cde6+P+GNYpyZ9U8e6w5Uu/s3s16vmhDI4YQ/2p7A/yN8zp7g7oOT++OzahgIHAH0B44zs355/q4kt88Ajcl2sTdQSfj99iP8jY0QtJ7u/hnCgJmzkoEyvwJOTf7eVxBGua33srv3dvfpWcvGA4vcvTfweWCCmXVKbh/o7lXAw3w89BU2s83ma5vvQRN17r7+lOT6E01/c/fXk9u3AZdltZ+U1bZL0ydz93VmdiwwxsJoigOB7LPAjyX/LwSmJyOgFibvjBDe6dqb2WnJ/Y6EIEH4zNn0yPOBwB/d/b3k/hgAM3sWONzMLkvW75C1zhPu/h/gP2b2Wq6fQ/Ln7s+Y2TIzOw/oDexL+H0vdfeZAGa2GHgiWWUh4Y3UgF7AH7IG3nzSzCqT28/leLmDgFOS130NGJA8/0nAiWa2L3A4Gw7W2eQ2m69SBfTDrNsZ4BjYYAxBCmibo30meQwzy/5ljAReAO4lnO6pBbK7kPVZtxty1JMGTnb3F5Pn3oXQBR8LrM7RviGphaR91+Q1bicc0f4Doat1Uo6fYYOfQwqTnCm4lnD24E7Cx4+FbPi3hv/+e6eBeesP/phZCujm7vVJYPP5e1vS7hnCR6M/A28Bg7LW+a9tthCxXG42AxhqZnsl988Gcgzr+Ji7D1z/j/DuuQ74PvAkoauT3oLXfxI418xSZvZJ4EVC92hjngGOMLMdk+71rwjvsqOB77r7REIXFzPbkjokf4cCv3b3OwnhGEV+f/O5QBczG5Xc/yLw6GbWeYrkzTYZpDOV0MX+p7tfT9h+t3Sby0ssAX2bEMpHknOmB7Nhf35zZhMOBM0lhOsDoGfy7piPa4DtgTnAs8C17r7RMV7u/hLhPPBzyTqvuPvDhOGO08zsRUKXZz6w18aeR5rldmCsmb0EPAI8T9huNin5mPE/wP83s/U9rZM3s9rVwF5mNpvQO/oSoQubMrNXCOMBFlCEv7VOs4hELJY9qIjkoICKREwBFYmYAioSMQVUJGIKaJkxs0+b2R2lrkPyo9MsIhEr1VA/KaLkBPwp7l5rZr8E2rn7WDPrBTxEmMztwOTqjemECwd2A65297s3dnVPSX6YMqcubus0iTAUDsLcxUOT20cATS82b+/uIwhXZlyfLNvU1T2yDSmgrdMk4NBkbPMbwFIz24cQ0KYTs+W66uKzwFeSCxJqCHMfK6AloIC2Ts8TrlM9jDCw+ylCOHsRxipn+xCgydd2rL+6Z/3FCMOA3xaBqZ5ZAAAAZ0lEQVS5ZslBAW2F3L2RMOj/IkJAnyTMJPHEptbLsqVX90iRKKCt1yTCdBuzCd3Undn8ZVXrbdHVPVI8Os0iEjHtQUUipoCKREwBFYmYAioSMQVUJGIKqEjEFFCRiCmgIhFTQEUi9n8M4kee1gO/hAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 240x320 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tmp = agg_prediction(m_05.trace['score'], wines,\n",
    "                     ['wine_amer', 'flight_red'], summarize())\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(3, 4), dpi=80, sharey=True)\n",
    "\n",
    "\n",
    "for flight, c in zip((0, 1), ('g', 'r')):\n",
    "    sub = tmp[(tmp.flight_red == flight)]\n",
    "    data_sub = wines[(wines.flight_red == flight)]\n",
    "\n",
    "    ax.scatter(data_sub.wine_amer, data_sub.score,\n",
    "               c=c, marker='.')\n",
    "\n",
    "    ax.plot(sub.wine_amer, sub.mean_, c=c, marker='o')\n",
    "\n",
    "    ax.set(xticks=[0, 1], xticklabels=['non-american', 'american'],\n",
    "           xlabel='wine',\n",
    "           xlim=(-.2, 1.2), ylim=(13, 16))\n",
    "        \n",
    "fig.tight_layout()\n",
    "None"
   ]
  }
 ],
 "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"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}