diff --git a/exercises/Chapter5_part1.ipynb b/exercises/Chapter5_part1.ipynb new file mode 100644 index 0000000..df949e0 --- /dev/null +++ b/exercises/Chapter5_part1.ipynb @@ -0,0 +1,692 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Chapter 2 Exercises" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING (theano.configdefaults): install mkl with `conda install mkl-service`: No module named 'mkl'\n" + ] + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n", + "from scipy import stats\n", + "import arviz as az\n", + "import pymc3 as pm\n", + "np.random.seed(123)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 1\n", + "***" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'y')" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFz9JREFUeJzt3X+s3fV93/HnG2PIzRrlktgN+IKL6RjLMqtzekeTWppYEkqKKuy4qUb+WKBKZGUrqvaLySxTGkXb7BRpUrpkTd0kKkgdkDHiuIPMS+JGWVeRcYkhDjAXBxVxr1lwICaL4qU2vPfH+V44vpz743vvOefz/Z7zfEhH95zv+fict7+2vq/7/X5+fCMzkSRppc4rXYAkqV0MDklSLQaHJKkWg0OSVIvBIUmqxeCQJNVicEiSajE4JEm1GBySpFrOL13AIGzYsCEvv/zy0mVIUms8/PDDP8jMjStpO5LBcfnllzMzM1O6DElqjYh4eqVtvVQlSarF4JAk1WJwSJJqMTgkSbUYHJKkWgwOSVItxYbjRsRlwJ3AxcDLwP7M/NSCNgF8Crge+Alwc2Z+e1A1HTgyx+2HjnHi1Gk2TU5w63VXsXPb1KC+TpJaqeQ8jrPAP8/Mb0fEG4CHI+Krmfl4V5tfBa6sHr8E/H71s+8OHJnjtvuOcvrMSwDMnTrNbfcdBTA8JKlLsUtVmfns/NlDZv5f4Alg4RF6B3BndjwITEbEJYOo5/ZDx14JjXmnz7zE7YeODeLrJKm1GtHHERGXA9uAby14awp4puv1LK8Nl/nP2B0RMxExc/Lkydo1nDh1utZ2SRpXxYMjIn4G+C/AP8nMHy18u8cfyV6fk5n7M3M6M6c3blzRcivn2DQ5UWu7JI2rosEREevphMYfZ+Z9PZrMApd1vb4UODGIWm697iom1q87Z9vE+nXcet1Vg/g6SWqtYsFRjZj6PPBEZv77RZodBD4YHe8AXszMZwdRz85tU+zdtZWpyQkCmJqcYO+urXaMS9ICJUdVbQf+IXA0Ih6ptv0rYDNAZn4WeIDOUNzjdIbj/uYgC9q5bWpoQeHQX0ltVSw4MvPP6N2H0d0mgd8aTkXD49BfSW1WvHN8HDn0V1KbGRwFOPRXUpsZHAU49FdSmxkcq3TgyBzb9x1my5772b7vMAeOzK34zzr0V1KbjeQ9xwdtrZ3b820cVSWpjQyOVViqc3ulB/9hDv2VpH7yUtUq2LktaZwZHKtg57akcWZwrIKd25LGmX0cq2DntqRxZnCskp3bksaVl6okSbUYHJKkWgwOSVIt9nE0mPfskNREBkdDec8OSU3lpaqG8p4dkprK4GgolzWR1FQGRwMdODLHedH7rrouayKptKLBERFfiIjnIuK7i7x/TUS8GBGPVI+PDbvGYZvv23gp8zXvuayJpCYo3Tn+R8CngTuXaPM/MvPXhlNOeb36NgDWRbB311Y7xiUVV/SMIzO/CbxQsoamWawP4+VMQ0NSI7Shj+OdEfFoRHwlIt5WuphBc8l2SU3X9OD4NvBzmfkLwH8ADizWMCJ2R8RMRMycPHlyaAX2m0u2S2q6RgdHZv4oM39cPX8AWB8RGxZpuz8zpzNzeuPGjUOts592bpti766tTE1OEMDU5IR9G5IapXTn+JIi4mLg+5mZEXE1naB7vnBZA+eS7ZKarGhwRMRdwDXAhoiYBX4HWA+QmZ8F3g/8o4g4C5wGbszsMU5VkjQ0RYMjMz+wzPufpjNcV5LUEI3u45AkNY/BIUmqxeCQJNVicEiSajE4JEm1GBySpFoMDklSLQaHJKkWg0OSVIvBIUmqxeCQJNVicEiSajE4JEm1GBySpFoMDklSLQaHJKkWg0OSVIvBIUmqpeitY9V+B47McfuhY5w4dZpNkxPcet1V7Nw2VbosSQNkcGjVDhyZ47b7jnL6zEsAzJ06zW33HQUwPKQRVvRSVUR8ISKei4jvLvJ+RMTvRcTxiPhORLx92DVqcbcfOvZKaMw7feYlbj90rFBFkoahdB/HHwHvXeL9XwWurB67gd8fQk1aoROnTtfaLmk0FA2OzPwm8MISTXYAd2bHg8BkRFwynOq0nE2TE7W2SxoNpc84ljMFPNP1erbapga49bqrmFi/7pxtE+vXcet1VxWqSNIwNL1zPHpsy54NI3bTuZzF5s2bB1lTY5Qe0TT/XY6qksZL04NjFris6/WlwIleDTNzP7AfYHp6ume4jJKmjGjauW3KoJDGTNMvVR0EPliNrnoH8GJmPlu6qCZwRJOkUoqecUTEXcA1wIaImAV+B1gPkJmfBR4ArgeOAz8BfrNMpc3jiCZJpRQNjsz8wDLvJ/BbQyqnVTZNTjDXIyQc0SRp0Jp+qUqLcESTpFKa3jmuRTiiSVIpBkeLOaJJUgkGhxbV73kipeedSOoPg0M99XueSFPmnUhaOzvH1VO/54k470QaHQaHeur3PBHnnUijw+BQT/1e+daVdKXRYXCop37PE3HeiTQ67BxXT/2eJ+K8E2l0RGdVj9EyPT2dMzMzpcuQpNaIiIczc3olbb1UJUmqxeCQJNVicEiSajE4JEm1GBySpFoMDklSLc7jUCu4sq7UHAaHGs+VdaVm8VKVGs+VdaVmKRocEfHeiDgWEccjYk+P92+OiJMR8Uj1+HCJOlWWK+tKzVLsUlVErAM+A1wLzAIPRcTBzHx8QdN7MvOWoReoxtg0OcFcj5BwZV2pjJJnHFcDxzPzqcz8K+BuYEfBetRQrqwrNUvJ4JgCnul6PVttW+jXI+I7EXFvRFw2nNLUJDu3TbF311amJicIYGpygr27ttoxLhVSclRV9Ni2cKnePwHuysyfRsRHgDuAd/X8sIjdwG6AzZs397POsdWkIbA7t00ZFFJDlDzjmAW6zyAuBU50N8jM5zPzp9XLPwR+cbEPy8z9mTmdmdMbN27se7HjZn4I7Nyp0ySvDoE9cGSudGmSCisZHA8BV0bEloi4ALgRONjdICIu6Xp5A/DEEOsba20bAnvgyBzb9x1my5772b7vsAEnDVCxS1WZeTYibgEOAeuAL2TmYxHxCWAmMw8Cvx0RNwBngReAm0vVO27aNATWCYLScBWdOZ6ZDwAPLNj2sa7ntwG3DbsutWsI7FJnRwaH1H/OHFdPbRoC26azI2kUGBzqqU1DYBc7C2ri2ZE0ClzkUItqyxDYW6+76pw+Dmju2ZE0CgwOtd58uDVlzok06pYNjmrk0x9n5g+HUI+0Km05O5JGwUr6OC6mswDhF6vVbHvN+JYkjYllgyMz/zVwJfB5OvMonoyIfxcRPz/g2iRJDbSiUVWZmcD/qR5ngYuAeyPidwdYmySpgVbSx/HbwE3AD4DPAbdm5pmIOA94EviXgy1RktQkKxlVtQHYlZlPd2/MzJcj4tcGU5YkqamWDY7uJUB6vOeigxp7TVp+XhoG53FIa+ACixpHLjkirUHblp+X+sHgkNbABRY1jgwOaQ1cYFHjyD4OqabuzvA3Tqxn/brgzEv5yvsusKhRZ3BINSzsDD91+gzrzwsuev16Tv3kjKOqNBYMDqmGXp3hZ15OXn/B+Rz52K8UqkoaLvs4pBrsDJcMDqkWO8OlwsFRLdN+LCKOR8SeHu9fGBH3VO9/KyIuH36V0qvadC92aVCK9XFExDrgM8C1wCyde34czMzHu5p9CPhhZv71iLgR+CTwD4ZfrdpkkEuAeLdBqWzn+NXA8cx8CiAi7gZ2AN3BsQP4ePX8XuDTERHVMu/SawxjCRDvNqhxV/JS1RTwTNfr2WpbzzaZeRZ4EXjzUKpTK7kEiDR4JYOj1y1oF55JrKRNp2HE7oiYiYiZkydPrrk4tZOjnqTBKxkcs8BlXa8vBU4s1iYizgfeCLzQ68Myc39mTmfm9MaNGwdQrtrAUU/S4JUMjoeAKyNiS0RcANwIHFzQ5iCduw8CvB84bP+GluKoJ2nwigVH1WdxC3AIeAL4YmY+FhGfiIgbqmafB94cEceBfwa8Zsiu1G3ntin27trK1OQEAVz0+vVceP55/NN7HmH7vsMcODJXukSp9WIUf4Gfnp7OmZmZ0mWosIUjrKBz9rF311ZHRUkLRMTDmTm9krbOHNfIcoSVNBgGh0aWI6ykwTA4NLIcYSUNhsGhkdWWEVYHjsyxfd9htuy53w58tYL349DIasO6UsNYIkXqN4NDI63p60ot1YHf5Lo13rxUJRVkB77ayOCQCrIDX21kcEgFtaUDX+pmH4dUUBs68KWFDA6psKZ34EsLealKklSLwSFJqsXgkCTVYnBIkmoxOCRJtTiqStJQHDgy57DjEWFwSCOuCQdsF3McLV6qkkbY/AF77tRpklcP2MNeut27MY4Wg0MaYU05YLuY42gpEhwR8aaI+GpEPFn9vGiRdi9FxCPV4+Cw65TarikHbBdzHC2lzjj2AF/PzCuBr1evezmdmX+netwwvPKk0dCUA7aLOY6WUsGxA7ijen4HsLNQHdJIa8oBe+e2Kfbu2srU5AQBTE1OsHfXVjvGWyoyc/hfGnEqMye7Xv8wM19zuSoizgKPAGeBfZl5YInP3A3sBti8efMvPv300/0vXGqhJoyqUvNFxMOZOb2itoMKjoj4GnBxj7c+CtyxwuDYlJknIuIK4DDw7sz83nLfPT09nTMzM2uoXpLGS53gGNg8jsx8z2LvRcT3I+KSzHw2Ii4BnlvkM05UP5+KiG8A24Blg0NS83jmMzpK9XEcBG6qnt8EfHlhg4i4KCIurJ5vALYDjw+tQkl905T5JOqPUsGxD7g2Ip4Erq1eExHTEfG5qs1bgZmIeBT4Uzp9HAaH1EJNmU+i/iiy5EhmPg+8u8f2GeDD1fM/B7YOuTRJA7CW+SRe4moeZ45LGrjVzifxElczGRySBm6180m8xNVMro4raeDmLy3VveTUlCVTdC6DQ9JQ7Nw2VbtvYtPkBHM9QqKfS6bYh1Kfl6qkMXPgyBzb9x1my5772b7vcKP7Cwa9ZIp9KKtjcEhjpG0HykGvcWUfyup4qUoaI0sdKJt6eWY1l7hWyj6U1fGMQxojHijP1ZRl59vG4JAaahB9ER4oz9WUZefbxuCQGmhQfREeKM/lfUJWxz4OqYEG1Rex2vkUo2yQfSijyuCQGmiQfREeKLVWBofUQMOY+Kb+G5fJhAaH1EC3XncVt9139JzLVePcF9E0vQICOOffbL5fChi58DA4pAbqV19Ev38DHpffqJcyP3BhYUC8bv15rZsjs1oGh9RQa+2LWOwAN//ZpT+vrRYbuLBw27xRnCPjcFxpRPV7OQ2X5+ioGwSj2C9lcEgjqt8js5x13rFYEExOrB+bOTIGhzSi+j1L3FnnHYtNovz4DW8bm8mE9nFII6rfI7Mc6dWx3MCFUQyKhYoER0T8BvBx4K3A1Zk5s0i79wKfAtYBn8vMfUMrUmq5fs8Sd9b5q8Z9EmVk5vC/NOKtwMvAHwD/oldwRMQ64C+Aa4FZ4CHgA5n5+HKfPz09nTMzPbNIktRDRDycmdMraVvkjCMznwCIiKWaXQ0cz8ynqrZ3AzuAZYNDkjQ4Te7jmAKe6Xo9C/zSYo0jYjewG2Dz5s2DrUySGqLEpMyBBUdEfA24uMdbH83ML6/kI3psW/S6WmbuB/ZD51LVioqUpBYrNSlzYMGRme9Z40fMApd1vb4UOLHGz5SkkVHqVsBNnsfxEHBlRGyJiAuAG4GDhWuSpMYoNSmzSHBExPsiYhZ4J3B/RByqtm+KiAcAMvMscAtwCHgC+GJmPlaiXklqolKTMosER2Z+KTMvzcwLM/MtmXldtf1EZl7f1e6BzPwbmfnzmflvS9QqSU1V6lbATR5VJUlaQqlJmQaHJLVIr+G3/3PPu4Zag8EhSQMwiPkVTbknSpNHVUlSK80f4OdOnSZ59QB/4Mjcmj63KfdEMTgkqc8GdYBvyj1RDA5J6rNBHeCbck8Ug0OS+mxQB/hSw28XMjgktd6BI3Ns33eYLXvuZ/u+w2vuS1irQR3gd26basRdBh1VJanVmjLSqNsg51c04SZSBoekViu10N9ymnCAHxQvVUlqtaaMNBonBoekVmvKSKNxYnBIarWmjDQaJ/ZxSGq1Ugv9jTODQ1Lr1e2ILnGf7lFicEgaK00cvts29nFIGitNWSiwzQwOSWPF4btrZ3BIGisO3127IsEREb8REY9FxMsRMb1Eu7+MiKMR8UhEzAyzRkmjyeG7a1eqc/y7wC7gD1bQ9u9n5g8GXI+kMeHw3bUrEhyZ+QRARJT4ekljbpTXkRqGpvdxJPDfI+LhiNhduhhJ0gDPOCLia8DFPd76aGZ+eYUfsz0zT0TEzwJfjYj/nZnfXOT7dgO7ATZv3ryqmiVJyxtYcGTme/rwGSeqn89FxJeAq4GewZGZ+4H9ANPT07nW75Yk9dbYS1UR8dci4g3zz4FfodOpLkkqqNRw3PdFxCzwTuD+iDhUbd8UEQ9Uzd4C/FlEPAr8L+D+zPxvJeqVJL2q1KiqLwFf6rH9BHB99fwp4BeGXJokaRmROXrdARFxEni6dB1dNgBtmovStnrBmoelbTW3rV4oV/PPZebGlTQcyeBomoiYycxFZ8g3TdvqBWselrbV3LZ6oR01N7ZzXJLUTAaHJKkWg2M49pcuoKa21QvWPCxtq7lt9UILaraPQ5JUi2cckqRaDI4BaNv9RmrU+96IOBYRxyNizzBr7FHLmyLiqxHxZPXzokXavVTt30ci4uCw66xqWHK/RcSFEXFP9f63IuLy4Vd5Tj3L1XtzRJzs2q8fLlFnVz1fiIjnIqLnyhLR8XvV3+c7EfH2YdfYo6blar4mIl7s2scfG3aNS8pMH31+AG8FrgK+AUwv0e4vgQ1tqBdYB3wPuAK4AHgU+FsFa/5dYE/1fA/wyUXa/bjwvl12vwH/GPhs9fxG4J6G13sz8OmS+3VBPX8PeDvw3UXevx74ChDAO4BvtaDma4D/WrrOxR6ecQxAZj6RmcdK17FSK6z3auB4Zj6VmX8F3A3sGHx1i9oB3FE9vwPYWbCWpaxkv3X/Xe4F3h3lblbTtH/nZWVnxewXlmiyA7gzOx4EJiPikuFU19sKam40g6OsNt1vZAp4puv1bLWtlLdk5rMA1c+fXaTd6yJiJiIejIgS4bKS/fZKm8w8C7wIvHko1b3WSv+df7267HNvRFw2nNJWrWn/d1fqnRHxaER8JSLeVrqYbqVuHdt6w77fyFr1od5evwEPdEjeUjXX+JjN1T6+AjgcEUcz83v9qXBFVrLfhr5vl7CSWv4EuCszfxoRH6FztvSugVe2ek3avyv1bTpLgPw4Iq4HDgBXFq7pFQbHKuWQ7zfSh+9aa72zQPdvlpcCJ9b4mUtaquaI+H5EXJKZz1aXHZ5b5DPm9/FTEfENYBuda/jDspL9Nt9mNiLOB95IucsYy9abmc93vfxD4JNDqGsthv5/d60y80ddzx+IiP8YERsysxHrbnmpqpAW3m/kIeDKiNgSERfQ6cQtMkqpchC4qXp+E/Cas6aIuCgiLqyebwC2A48PrcKOley37r/L+4HDWfWQFrBsvQv6B24AnhhifatxEPhgNbrqHcCL85c5myoiLp7v54qIq+kcq59f+k8NUene+VF8AO+j81vOT4HvA4eq7ZuAB6rnV9AZsfIo8BidS0aNrbd6fT3wF3R+Yy9Wb1XLm4GvA09WP99UbZ8GPlc9/2XgaLWPjwIfKlTra/Yb8Anghur564D/DBync++ZKwrv2+Xq3Vv9n30U+FPgbxau9y7gWeBM9f/4Q8BHgI9U7wfwmervc5QlRjo2qOZbuvbxg8Avl665++HMcUlSLV6qkiTVYnBIkmoxOCRJtRgckqRaDA5JUi0GhySpFoNDklSLwSENWET83WpBwNdVKwY8FhF/u3Rd0mo5AVAagoj4N3RmiE8As5m5t3BJ0qoZHNIQVOs+PQT8PzrLR7xUuCRp1bxUJQ3Hm4CfAd5A58xDai3POKQhqO53fjewBbgkM28pXJK0at6PQxqwiPggcDYz/1NErAP+PCLelZmHS9cmrYZnHJKkWuzjkCTVYnBIkmoxOCRJtRgckqRaDA5JUi0GhySpFoNDklSLwSFJquX/A4zyttikCzDzAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "dummy_data = np.loadtxt('../code/data/dummy.csv')\n", + "x_1 = dummy_data[:, 0]\n", + "y_1 = dummy_data[:, 1]\n", + "\n", + "# Set order of polynomial fit for Question 1\n", + "order_1 = 5\n", + "x_1p = np.vstack([x_1**i for i in range(1, order_1+1)])\n", + "x_1s = (x_1p - x_1p.mean(axis=1, keepdims=True)) / \\\n", + " x_1p.std(axis=1, keepdims=True)\n", + "y_1s = (y_1 - y_1.mean()) / y_1.std()\n", + "plt.scatter(x_1s[0], y_1s)\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "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: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [00:32<00:00, 154.94draws/s]\n", + "There was 1 divergence after tuning. Increase `target_accept` or reparameterize.\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (2 chains in 2 jobs)\n", + "NUTS: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [02:05<00:00, 39.74draws/s]\n", + "There were 5 divergences after tuning. Increase `target_accept` or reparameterize.\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (2 chains in 2 jobs)\n", + "NUTS: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [00:07<00:00, 641.02draws/s]\n", + "The acceptance probability does not match the target. It is 0.8850791171163658, but should be close to 0.8. Try to increase the number of tuning steps.\n" + ] + } + ], + "source": [ + "trace_1 = []\n", + "for sd in [1, 100, np.array([10, 0.1, 0.1, 0.1, 0.1])]:\n", + " with pm.Model() as model_p:\n", + " α = pm.Normal('α', mu=0, sd=1)\n", + " β = pm.Normal('β', mu=0, sd=sd, shape=order_1)\n", + " ϵ = pm.HalfNormal('ϵ', 5)\n", + "\n", + " μ = α + pm.math.dot(β, x_1s)\n", + "\n", + " y_pred = pm.Normal('y_pred', mu=μ, sd=ϵ, observed=y_1s)\n", + "\n", + " trace_p = pm.sample(2000)\n", + " trace_1.append([sd, trace_p])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4m9WZ+P3v0S7Zsi1b3ncnIfu+QcmwhVD2rZRAaYG2DJ2WaYeZa3jpvL1aWqZD6W86v/f9tWXolYEWOk2hb9jCQENZAgQo2ZyE7Hvs2I7jfbdlWdJ5/3hs4yS2Y1uSZSv357p02ZYfneeW4tw6Os8591Faa4QQQsQPU6wDEEIIEVmS2IUQIs5IYhdCiDgjiV0IIeKMJHYhhIgzktiFECLOSGIXQog4I4ldCCHijCR2IYSIM5ZYnNTr9eqioqJYnFoIISat0tLSeq11+vmOi0liLyoqYvv27bE4tRBCTFpKqfKRHCdDMUIIEWcksQshRJyRxC6EEHEmJmPsYnz09PRQWVmJz+eLdShCiFFwOBzk5eVhtVrH9HhJ7HGssrISt9tNUVERSqlYhyOEGAGtNQ0NDVRWVlJcXDymNmQoJo75fD7S0tIkqQsxiSilSEtLC+uTtiT2OCdJXYjJJ9z/t5Mvsfs74cRH0FYT60iEEGJCmnyJXSko/wQajsY6EjHOioqKqK+vD/uYkfrggw+48cYbw2qjqKiIuXPnsmDBApYsWTLmdq644opBF/U1NDRw5ZVXkpiYyN///d+HE6qII5Pv4qnVCQleaD4JXBrraEQcCQaDmM3mMT8+EAhgsZz7X+r999/H6/WGE9qQHA4H//qv/8revXvZu3dvVM4hJp8R99iVUr9VStUqpfYOuO/flVIHlVK7lVKvKqVSohPmWVIKobUSQsFxOZ0Ym7KyMmbMmMEDDzzAnDlzuOeee3j33Xe59NJLmTZtGlu3bgWgsbGRW2+9lXnz5nHxxReze/duwOiNXnPNNSxcuJBvfetbaK372/7DH/7AsmXLWLBgAd/61rcIBof/W3jhhReYO3cuc+bM4dFHH+2/PzExkR/96EcsX76cTz/9lLfeeosZM2awYsUKXnnllf7jOjo6+MY3vsHSpUtZuHAh69evB+C5557jy1/+MjfddBPXXHPNmF6ndevWMWfOHObPn89ll10GQFdXF3fddRfz5s1j9erVdHV1DfrYhIQEVqxYgcPhGNO5RXwaTY/9OeDXwO8H3PcO8C9a64BS6ufAvwCPDvLYyErOh8rt0HYaknOjfrp48MGhWurauiPaZrrbzhXTM4Y95ujRo6xbt441a9awdOlS/vjHP/Lxxx/z+uuv88QTT/Daa6/x2GOPsXDhQl577TU2btzIvffey65du/jJT37CihUr+NGPfsSbb77JmjVrADhw4AB/+tOf+OSTT7BarXznO99h7dq13HvvvYPGcOrUKR599FFKS0vxeDxcc801vPbaa9x66610dHQwZ84cHn/8cXw+H9OmTWPjxo1MnTqV1atX97fxb//2b1x11VX89re/pbm5mWXLlnH11VcD8Omnn7J7925SU1PPObdSimuuuQalFN/61rd48MEHzznm8ccf5y9/+Qu5ubk0NzcD8PTTT+Nyudi9eze7d+9m0aJFI/tHEYJR9Ni11puAxrPue1trHej9cTOQF8HYhpaSb3xtqRiX04mxKy4uZu7cuZhMJmbPns3KlStRSjF37lzKysoA+Pjjj/na174GwFVXXUVDQwMtLS1s2rSJr371qwDccMMNeDweAN577z1KS0tZunQpCxYs4L333uP48eNDxrBt2zauuOIK0tPTsVgs3HPPPWzatAkAs9nMl770JQAOHjxIcXEx06ZNQynVf26At99+myeffJIFCxZwxRVX4PP5OHnyJACrVq0aNKkDfPLJJ+zYsYMNGzbw1FNP9Z93oEsvvZT777+f//qv/+r/5DHwuc+bN4958+aN7AUXgsiOsX8D+NNQv1RKPQg8CFBQUBDemWwJ4EozxtkLLg6vrQvE+XrW0WK32/u/N5lM/T+bTCYCAaNPMHCIpU/fdK/Bpn1prbnvvvv42c9+NqIYBmu/j8PhOGNcfahpZlprXn75ZaZPn37G/Vu2bCEhIWHI9nNycgDIyMjgtttuY+vWrf3DLX1+85vfsGXLFt58800WLFjArl27hozl1Vdf5Sc/+QkAzzzzTFgXZEX8isisGKXUD4AAsHaoY7TWa7TWS7TWS9LTz1tO+PxSCoweeygUflsipi677DLWrjX+dD744AO8Xi9JSUln3L9hwwaampoAWLlyJS+99BK1tbWAMUZfXj50NdPly5fz4YcfUl9fTzAY5IUXXuDyyy8/57gZM2Zw4sQJjh07Bhjj8n2++MUv8qtf/ar/TWLnzp3nfV4dHR20tbX1f//2228zZ86cc447duwYy5cv5/HHH8fr9VJRUXHGc9+7d2//dYfbbruNXbt2sWvXLknqYkhh99iVUvcBNwIr9XBdo0hLyYdTO6GjFtxZ43ZaEXk//vGP+frXv868efNwuVw8//zzADz22GPcfffdLFq0iMsvv7z/k96sWbP46U9/yjXXXEMoFMJqtfLUU09RWFg4aPvZ2dn87Gc/48orr0RrzfXXX88tt9xyznEOh4M1a9Zwww034PV6WbFiRf9Mkx/+8Ic8/PDDzJs3D601RUVFvPHGG8M+r5qaGm677TbAmDHzla98hWuvvfac4x555BGOHDmC1pqVK1cyf/58pk+f3v+aLFiwgGXLlg15nqKiIlpbW/H7/bz22mu8/fbbzJo1a9jYRHxTo8nFSqki4A2t9Zzen68F/jdwuda6bqTtLFmyRIe90YavFT59CqZeDflLw2srTh04cICZM2fGOgwhxBgM9v9XKVWqtT7vR7XRTHd8AfgUmK6UqlRKfRNjlowbeEcptUsp9ZvRhR4GRxI4U6B5RBuKCCHEBWPEQzFa67sHufvZCMYyeikFUH8EtDZWpAohhJiEJQUGSs6Hni7oiMwSciGEiAeTO7H3zWdvPhnbOIQQYgKZ3IndkQKOZGgui3UkQggxYUzuxK4UeAqNHrvMZxdCCGCyJ3YATxH0+KBd6rPHOynbO7qyvaWlpcydO5epU6fyve99r39xVWNjI6tWrWLatGmsWrWqf+GXiB+TP7Gn9C5KaSqLaRhi8jtfhcjz6SuRcLb333+fXbt2DZqYw9VXtvcXv/jFOb/79re/zZo1azhy5AhHjhzhrbfeAuDJJ59k5cqVHDlyhJUrV/Lkk09GPC4RW5M/sdsTe+uzy3z2iUbK9o5MNMr2VldX09rayiWXXIJSinvvvZfXXnsNgPXr13PfffcBcN999/XfL+LH5NtoYzCeIqjeBcEAmOPjKUXckXcjP1yVmAnTrh72ECnbG5uyvVVVVeTlfV5sNS8vj6qqKsAodZCdnQ0Y5Rb6au6I+DH5e+xgJPZgAFqrYh2JOIuU7Y1N2d7hKmaK+Bcf3dvkfGOGTFOZMUtGnOs8PetokbK9sSnbm5eXR2VlZf/PlZWV/bFkZmZSXV1NdnY21dXVZGTEpqSziJ746LFbHeDOlnH2SUrK9ka+bG92djZut5vNmzejteb3v/99f0XLm2++ub+C5vPPPz9opUsxucVHjx2MnvrJLRDoBov9/MeLCUPK9kanbO/TTz/N/fffT1dXF9dddx3XXXcdAN///ve58847efbZZykoKGDdunXDPg8x+YyqbG+kRKRs79maymDXCzD3y+CdGtm2Jykp2yvE5DUuZXsnvKQ8Y0ZM04lYRyKEEDEVP4ndbDEWKzVKYhdCXNjiJ7EDpJZAZwN0yRJpIcSFK/4SO0Dj0HOahRAi3o1ma7zfKqVqlVJ7B9yXqpR6Ryl1pPerJzphjpDTY2yXJ8MxQogL2Gh67M8BZ8/V+j7wntZ6GvBe78+xo5TRa28qg1B4BZ2EEGKyGnFi11pvAhrPuvsW4Pne758Hbo1QXGOXWgLBHmipiHUkQggRE+GOsWdqrasBer/Gfm1ySgGYzDLOHoekHvu5az+ee+450tPTeeCBB4Cx1Wcfi7feeovp06czderUIcv+btq0iUWLFmGxWHjppZdG1O5IYjx48CCXXHIJdrt90HLFgzlx4gTLly9n2rRprF69Gr/ff84xw712Qxlpbftrr72WlJSUc/6e7rnnHlJTU0f8+ozUuF08VUo9qJTarpTaXldXF70TWeyQnCeJXYzaZKzHDrB69WqeeeYZYGz12UcrGAzy0EMPsWHDBvbv388LL7zA/v37zzmuoKCA5557jq985SsjbnskMaampvLLX/6Sf/7nfx5xu48++ij/+I//yJEjR/B4PDz77LPnHDPcazeUkda2f+SRR/jv//7vc+5fu3YtN99884jPN1LhlhSoUUpla62rlVLZwJD1P7XWa4A1YKw8DfO8w0stgWPvg68VHElRPdVk8XHVx9R3RaYn28fr9LIid8WQvy8rK+Paa69lxYoVbN68mfnz5/P1r3+dxx57jNraWtauXcuyZctobGzkG9/4BsePH8flcrFmzRrmzZtHQ0MDd999N3V1dSxbtuyceuy//OUv8fv9LF++nP/8z/88o5jX2V544QWeeOIJtNbccMMN/PznPweMeuz/9E//xF/+8hf+4z/+g/b2dh5++GG8Xu8ZpXI7Ojr47ne/y549ewgEAvz4xz/mlltu4bnnnuPNN9/E5/PR0dHBxo0bR/06rlu3jp/85CeYzWaSk5PZtGkTXV1dfP3rX2f//v3MnDlzyHrsZ+urz3706NEz7h9Ynx3or8/eV2ZgNLZu3crUqVMpKTFmod11112sX7+eWbNmnXFcUVERYBR8G4mRxpiRkUFGRgZvvvnmiNrVWrNx40b++Mc/AkYN+h//+Md8+9vfPuO4oV674axfv54PPvigv90rrrii/29roJUrV/YfNx7C7bG/DtzX+/19wPow24sMmfY4YRw9epR/+Id/YPfu3Rw8eLC/HvsvfvELnnjiCYD+euy7d+/miSee6K+r3lePfefOndx88839ZXIH1mPftWsXZrO5v2DWYPrqsW/cuJFdu3axbdu2/s0l+uqxb9myhSVLlvC3f/u3/M///A8fffQRp0+f7m+jrx77tm3beP/993nkkUfo6OgAjHrszz///KBJva8e++LFi/vryZ+trx77Z599xuuvvw6cWY/9Bz/4AaWlpaN96c8wXH32sbSVn58fkbaiFeNADQ0NpKSkYLFYItouTNza9iPusSulXgCuALxKqUrgMeBJ4P9TSn0TOAl8ORpBjlpCOtjd0HgMchbEOpoJYbiedTT11WMHhq3H/vLLLwPn1mPv28VoqHrsYOw2NFzp2YH12IH+euy33nrrkPXYAb761a/2J+O3336b119/vf9j+mjqsefk5FBbW8uqVauYMWPGOWV7++qx33nnndx+++2AMT79ve99DxhbPfazRbI+e7RqvU+2dieyESd2rfXdQ/xqZYRiiRylIG0q1OyVXZViTOqxR7ce+0gNV599tPLy8qio+HzWWThtRSvGgbxeL83NzQQCASwWS8TahYlb2z6+Vp4O5J1mTHuUGu0TntRjj3w99rMNV599tJYuXcqRI0c4ceIEfr+fF198MSIXACMZ40BKKa688sr+mSeRrEE/YWvba63H/bZ48WIddYEerT/8d60Pboj+uSao/fv3x/T8J06c0LNnz+7/+b777tPr1q0753cNDQ365ptv1nPnztXLly/Xn332mdZa6/r6er1q1Sq9cOFC/fDDD+uCggJdV1entdb6xRdf1PPnz9dz587VixYt0p9++qnWWuvCwsL+YwZau3atnjNnjp49e7Z+5JFH+u9PSEg447gNGzbo6dOn60svvVQ/+uij+oYbbtBaa93Z2akffPDB/jb67v/d736nH3rooUGf/7Fjx/S8efP0vHnz9KxZs/RPf/rTQY+77bbb+tv93ve+p0OhkO7s7NSrV6/Wc+fO1V/72tf0JZdcordt23bOYwc7f2FhofZ4PDohIUHn5ubqffv2aa213rZtm549e7YuKSnRDz30kA6FQlprrZ9++mn99NNPDxrbUN588009bdo0XVJScsbz+uEPf6jXr1+vtdZ669atOjc3V7tcLp2amqpnzZp13nZHEmN1dbXOzc3VbrdbJycn69zcXN3S0jJsu8eOHdNLly7VU6ZM0XfccYf2+Xxaa63Xr1+vf/jDH/YfN9RrN5T6+np91VVX6alTp+qrrrpKNzQ09D+Pb37zm/3HrVixQnu9Xu1wOHRubq5+6623+n838P/FQIP9/wW26xHk2Pipxz6YvS9DazVc8pAxPHOBkXrs8e+5555j+/bt/PrXv451KGKM7r//fm688UbuuOOOM+6XeuxDSZsG3W3QXhPrSISICqfTyYYNG/oXKInJ5Z577uHDDz/E4XBEtN34vqqYNsXoqdcfAXdWrKOJCa113M8AuJCtXr2a1atXxzoMMUZDTdMNdyQlvnvstgRIyoWGI7GOJCYcDgcNDQ1h/5EIIcaP1pqGhoawevHx3WMHY9rj8Q8uyFWofdPHolrCQQgRcQ6H44zFWqMV/4ndO81I7A1HIHdxrKMZV1arleLi4liHIYQYZ/E9FAPgSgNXqjHOLoQQF4D4T+xKgfciaCoHf2esoxFCiKiL/8QOkD4DdOiCvYgqhLiwXBiJ3Z0FjmSoOxTrSIQQIuoujMSuFKRPN/ZC7fHFOhohhIiqCyOxgzEcEwpCw8iL6AshxGR04ST2pByjRnvdwVhHIoQQUXXhJPa+4ZjGExDojnU0QggRNRdOYgcjsYcC0HAs1pEIIUTUXFiJPSnPqB9TdyDWkQghRNREJLErpf5RKbVPKbVXKfWCUiqyNSgjxWSCjJnQcFxmxwgh4lbYiV0plQt8D1iitZ4DmIG7wm03ajJmGcMx9TKnXQgRnyI1FGMBnEopC+ACTkWo3chLygGnB2r2xzoSIYSIirATu9a6CvgFcBKoBlq01m+H227UKAWZs4xNrrvbYh2NEEJEXCSGYjzALUAxkAMkKKW+OshxDyqltiultse8PnjGbNAaamVOuxAi/kRiKOZq4ITWuk5r3QO8Anzh7IO01mu01ku01kvS09MjcNowJKSBOxNq9sY2DiGEiIJIJPaTwMVKKZcyNtdcCUz8+YQZs6HtNHQ2xjoSIYSIqEiMsW8BXgJ2AHt621wTbrtRlzHTGG+v2RfrSIQQIqIiMitGa/2Y1nqG1nqO1vprWuuJv2bfkQQpBcZwjGz2LISIIxfWytOzZc2DrmZoPhnrSIQQImIu7MSePh0sNji9O9aRCCFExFzYid1sNVai1h2Uio9CiLhxYSd2MIZjggGonfgTeYQQYiQksSflQIJXhmOEEHFDErtSRq+9pQo6GmIdjRBChE0SO0DmbFAmOP1ZrCMRQoiwSWIHsCdC2hQ4vccYbxdCiElMEnufnIXg75Q67UKISU8Se5/UEqNOe9WOWEcihBBhkcTeRymj195SCe21sY5GCCHGTBL7QFlzwWSBUztjHYkQQoyZJPaBbC6j6uPpPbISVQgxaUliP1vuIgj2yCYcQohJSxL72dzZ4M4yLqJKOV8hxCQkif1sSkHeEuioh8bjsY5GCCFGTRL7YDJmGYuWKrfFOhIhhBg1SeyDMZkhdzE0npCpj0KISSciiV0plaKUekkpdVApdUApdUkk2o2pnIVgtkivXQgx6USqx/5/gLe01jOA+cDkL25udRpVH2v2QXd7rKMRQogRCzuxK6WSgMuAZwG01n6tdXO47U4IeUtBh6CqNNaRCCHEiEWix14C1AG/U0rtVEo9o5RKOPsgpdSDSqntSqntdXV1ETjtOHClQtpUOLVDFiwJISaNSCR2C7AIeFprvRDoAL5/9kFa6zVa6yVa6yXp6ekROO04KbgEenxSZkAIMWlEIrFXApVa6y29P7+EkeijptXXgx6vxUPJuZBaDBVbjBWpQggxwYWd2LXWp4EKpdT03rtWAvvDbXcoDe3d/P6vZeypaonWKc5V+AWjVnu17LAkhJj4IjUr5rvAWqXUbmAB8ESE2j1HaoKN7GQnHx2pp9U3Tj3olAJIyYeTmyEUHJ9zCiHEGEUksWutd/WOn8/TWt+qtW6KRLuDUUpx9cxMtNa8f7B2/IZkCi6B7jaj8qMQQkxgk3LlabLLyhemejle18GhmrbxOWlqiVEc7OSn0msXQkxokzKxAyzISyE72cEHh+ro9I/DBtRKQfFl0NUM1buifz4hhBijSZvYTSbF1bMy8QdCvHdgnIZkUksgOQ/K/yozZIQQE9akTewA3kQ7X5iSxtHadg5Uj8OQjFJQcrlRYkBWowohJqhJndgBFhV4yPU4ef9QLS1d49CLTimAtCnGWHuPL/rnE0KIUZr0id1kUnxxVhYAb+87PT5DMsWXGUm9Ysv5jxVCiHE26RM7GLNkLr8oncqmLkrLozbT8nPuLMiYAZVbwdca/fMJIcQoxEViB5idk8S0zEQ+OdpAdUtX9E9YcoWxJ+qJTdE/lxBCjELcJPa+hUuJDgtv7q7G1xPlueZOj1HW9/QeaK2O7rmEEGIU4iaxAzisZm6Ym01Hd5C399dEf7y98Atgc8HRd43euxBCTABxldgBspIdrJjm5VhtOzsrorzfh8UOxZdDSyXUHYzuuYQQYoTiLrEDLCpIYUpGIh8drqeisTO6J8uaB4kZcGwjBPzRPZcQQoxAXCZ2pRRfnJ1JisvKn/dUR7cKpMkE01YZs2PKP4neeYQQYoTiMrED2C1mbpqfQyCkeeOzagLBUPROllIA2fOgYiu0T5Jt/4QQcStuEzsYtdu/ODuLmlYf70a7nkzJlcaY++ENciFVCBFTcZ3YAaZmJHLJlDQOVLey9URj9E5kc8GUq6ClSnZaEkLEVNwndoDlxanMzE7ir8caOHg6iitFs+YawzLH3zc25RBCiBi4IBK7sXgpg1yPk3f21VDVHKWVqUrB9OsgGIBDb8mQjBAiJiKW2JVSZqXUTqXUG5FqM5IsZhM3zcvB7bDw+q5T1Ld3R+dErlSj3EDDUdlGTwgRE5Hssf8DcCCC7UWc02bmtoV5mE3w2s6q6JX5zVtibH599F0pEiaEGHcRSexKqTzgBuCZSLQXTckuK7ctzMMfDPHqjsrobKunFEy/HnQQDsksGSHE+IpUj/3/Bf4vYMjJ4kqpB5VS25VS2+vqYjvXO91t55YFubT5Ary6syo6BcNcqcYsmcbjULkt8u0LIcQQwk7sSqkbgVqt9bB7xWmt12itl2itl6Snp4d72rDlpji5aX4ODe1+XtkRpeSeswi80+D4B1IBUggxbiLRY78UuFkpVQa8CFyllPpDBNqNuiJvAjfOy6a+vZvXdlbRHYhwclcKZtwAtgTYvx4CUbpgK4QQA4Sd2LXW/6K1ztNaFwF3ARu11l8NO7JxUpKeyPVzs6lp7ebVaPTcrU6YeTP4muGwTIEUQkTfBTGP/XymZiRyw7xsatu6eak0ChdUU/Kh6G+gZj9U7Yhs20IIcZaIJnat9Qda6xsj2eZ4mZqRyM3zc2ju9LNueyVtka4IWfgFY7z96LvQfDKybQshxADSYx+gyJvArQtzae8O8KdtFTREchGTUjDjRmNLvX2vgq8lcm0LIcQAktjPkudx8eXFeYS05k/bK6hsiuBGHVYHzPkShAKw9xUIRrFOvBDigiWJfRAZSQ5WLy0gwWbhlR1VHK6JYEGvhDTjYmp7DRx4HUJRrBMvhLggSWIfQrLTyuql+WQlOXhzdzWfHmuIXD137zSYejXUHTa21BNCiAiSxD4Mh9XM7YtymZWTxObjDby5pxp/IEI97LwlkLfUWJVaIStThRCRY4l1ABOdxWzimlmZeBPtfHSkjqbOCm6cm40nwRZ+41Ougu4WOPaesVFH5uzw2xRCXPCkxz4CSikWF3q4dUEu7b4Af9x6kqO1ERh3N5mM8faUAjjwhjE0I4QQYZLEPgpF3gS+sryA1AQb//NZNR8ergt/k2yz1Zgp486C/a9Bw7HIBCuEuGBJYh+lZKeVLy/OY0F+CjvKm/jT9goaO/zhNWqxw7w7wZVmTINsPB6ZYIUQFyRJ7GNgMZu4ckYGN83Poc0X4I9bytld2RzerBmrE+bfBS4P7HkZ6o9GLmAhxAVFEnsYpmYk8tWLC8lOdvLegVpe3VlFazilCGwJsOAeSPDC3peh9mDkghVCXDAksYcp0W7h9kW5XDUjg+oWH//9aTl7KlvG3nu3OmH+3ZCUY4y5n9oZ2YCFEHFPEnsEKKWYn5/CV5cXkuG28+6BGtaVVo691ozVAfNWg6cYDr0Fxz+Ucr9CiBGTxB5ByS4rdyzOY9WsTBra/azdcpJPjtaPbVGTxQZz74Ds+VD+Vzj4BoSisMuTECLuyAKlCFNKMSc3mZL0BDYdrmfriUYOVLeyYpqX6ZlulFIjb8xkhunXgSMZTmyCrmaYfRvYE6P3BIQQk5702KPEZbNw7Zws7lyaj9NmZsOe06zbXkl1S9foGlIKii6F2bdC+2kofQ5aT0UlZiFEfJDEHmW5KU7uXlrA1TMzaer08+LWCt7YfYqm0c59z5gJC+8FZYKda42dmGTcXQgxiLATu1IqXyn1vlLqgFJqn1LqHyIRWDwxmRRz85K5/9IiLi5Jo7yhk99/Ws47+2to6RrF9Eh3Jiy+3yhBcPgvxqyZHl/U4hZCTE4q3FK0SqlsIFtrvUMp5QZKgVu11vuHesySJUv09u3bwzrvZNbRHWBrWSN7Ko1dlGbnJLGkKJVkp3VkDWgNJzcb4+52N8y6GZLzohixEGIiUEqVaq2XnPe4iNUY//zE64Ffa63fGeqYcBN7aXkTm483cHFJGosLPWNuJ9bafD1sK2tkb1UrWsOMbDdLi1JJHWnlyJZK2P86dLdC/jIougzMcj1ciHgVk8SulCoCNgFztNatQx0XTmIvLW/inmc24w+EsFlMrH3g4qgm9/F4E2nz9VBa3sTeqhYCIU2xN4FFBR7yPM7zz6IJdBubdZzaZaxYnX49JOdGJU4hRGyNNLFHrHunlEoEXgYeHiypK6UeBB4EKCgoGPPZfK4sAAAcHElEQVR5Nh9vwB8IEdLQEwix+XhD1BLueL2JuB1WrpiewbLiVHZVNLO7soWXSivJTHKwID+FizITsZiHuBxisRtTIr0XwaENsPO/IXsBlFxurGIVQlxwIjIrRillxUjqa7XWrwx2jNZ6jdZ6idZ6SXp6+pjPdXFJGjaLCbMCq8XExSVpwx5fWt7EU+8fpbS8adTnGuxNJJpcNgtfmOLlmyuKWTkzg55giL/sO82zH5/gk6P1w19oTZsCy/7W2JmpehdsXQPVu2XmjBAXoLB77MoYK3gWOKC1/t/hhzS8xYUe1j5w8YiGR8Ltcfe9ifQEQiN6E4kUq9nEvLwU5uYmU9HYxa7KZraVNbKtrJHCNBdzc5Mp9iZiNp01TGOxG3upZs4xZs0cfBOqSo37UvLHJXYhROxFYlbMCuAjYA/Qt3b+/9Za/3mox4R78TSkQ/gCPnxBH76Aj+5gN/6gH3/QT0AHCISM2/pdVby84yQhrTEpuH1RHjfNz8WECZPJhFmZMSszVpMVi8mC1WTFZrZhM9uwm+3YzXb2VXWy9UTzsG8i4zEO3+rrYW9VC/tPtdLmC+CymZmRncSs7CTS3fZzH6A11OyD4x9AdxukX2RcXE0c+6clIURsxWxWzEiEk9jb/G38Yf8f0Awft0mZqG7pZt32SoJBMJsUq5fmk+txEtIhgjo44gqMNrMNp8WJ0+LEZXHhsrpwWpycrA+y6VAbGz5rJBCwYzPbWPvAJVG9mBsKacoaOthf3crxug6CIU26287MbDcXZbpxO86aMhnwQ+VWqNgCwR5jX9XCS8GVGrUYhRDREbeJPRAKsKNmBw6LA6fFicPi6O9dW01WrCYrZpMZkzIuHwzXm9ZaE9ABeoI9/V/9QT/+kJ/uYDfdge7+TwVdgS58QR+dPZ10Bjo5Xt/EKzsrCQY/f4tR2sJV0wu5YU4xSbYk3DY3SbYk42ZPwm4epGcdhi5/kIOnWzl0uo3qFh9KGStdp2e5mZqRiMs2YKTN3wkVm6GyFHQQMmZB4ReMmTRCiEkhbhP7RPGrjYf4f97bg1bdKFM3yuzDZvPz3avz8CaFaPO30R08s2yv3WwnxZ5Csj2ZZHsyKfaU/pvVPMLFSUNo6vBz8HQbh2vaaOzwY1KKXI+TaRmJTM1IJMHem+S726Biq1HnPRSAtKmQv9xY4DSaAmVCiHEniT3K+i7M9gRCmM0m7licx5cW5Z3xqaA72E2bv43W7lZa/C20dH9+6+jpOGM4KcGagMfhIcWeQqojtf+r0zKCuewDaK2pb/dzuKaNo7XtNHb4UQqykx1MSU9kSnoingSb0YOv3AandhhlCZKyIW8peKfLIichJihJ7OMgnIumPaEeI+F3t9DU3USzr9n42t2MP/h5gTC72U6qI7X/5nF4SHWk4rK6znsOrTUNHX6O1rZzrK6d2lbjE0Rqgo1ibwLF3gRy3RZMtXuNJN/ZCDaXUQM+ewE4U0b3ggghokoS+ySltaajp4MmXxON3Y00+5pp9DXS6Gs8Y2jHaXGS5kwj1ZFKmiONNGcaHocHq2noIZ2Wrh6O17Vzor6DyqYugiGN3WqiMDWBwlQnxZY6Euo+g4bejbRTCo0k771IevFCTACS2OOM1prOQCeNXY00+Br6k32jr5FAKACAQpFkT8Lr9PYn+zRnGm7ruRt8dAeCVDR2cqK+k7L6Dtq7jTa8bjtT3EGmBI/jbT+EubvVmB+fPsOYUZNSIGPxQsSIJPYLhNaaVn8rDV0N1HfV0+hrpL6rnlb/51UdbGYb7R0uqhstLC8o5G9KSkh1pGI2mfvbqGvvprzBSPLVLT6CIY3VBFNtDUzTZWR2l5NgCaHsSZAxA9JnGhtuS5IXYtxIYp9AYlGN0h/09yf5LeXl/GLjNkKqDbM5xO0L88jzJBjDOM40vE5v/81utuMPhKhs6uRko3FraPdjCvWQHahgKhVkB0+R4jDjSExBpc8wFj8l5YFJ9m0RIprGvQiYGNx4V6PsYzPbyErIIishiw/3OOhuNFbsmi1d2Lu9LMhwUddZR0VrBYcaD/U/LsmW1J/ki7K8LC1JJxiwUtnko6IpjR2NF9HZ2UFqRxk5zSfJq9xEsv0j3EnJOLOmobwXgafY2IxbCBETktijbDyqUZ7vE8HnNW/AqtxcP2Mhi7M/P66zp5P6rnrqu+qp66pjZ1UFB2u3k+dxkpNirLj1Or14PV4uz/FiV1k0t+VT1dzFpoZmbC1leE6Xk1mxBY/9U9wuOwkZxSRmT0elTQFXmgzZCDGOJLFHWbQLiY3kE8H5Cqe5rC4KrAUUJBVQWt7Ec2+F8AfTsdk7ePz2QrxJ3dR31bOrbld/GQab2UaaI43pU73YVT5+3wzq2s0cO30SW8txPIdOknRgD26HFVeSB3f2NJJypmPyFBpTKoUQUSOJPcpGU41yLEb6iWBxoWdE5/68PQs9vmTq6jO5c+5UwCjn0ORroq6rrr+Hf6DxQP+sHLMyk1qUisOcQ1vPdJraQqi6WuzNFSTXbMW26xMSHVYcnhwSs6bgyZ2G1VMgwzZCRJgk9nEw0qQ6FpH+RDBcexaThXRXOumuzytEhnSI5u5m6jrraOhqYFtlGXuqS8lKsZCT4kTlKZwFHhqCBag2P5bmFpLr60k59T7mHe/hstuwpeaRmFlCau5UnN5CCLO8ghAXOpkVEwciPetmrO19PiwUxGbz8++ri8nw+I0efmc9nYFOAALBEDpgw9YZxNbaQWJrCxn+dhyYsdus2D05uDKK8eRMJTmzCGV1hP2chIgHMivmAhGNqZRj/YTx+TCOosdv52SNm5tmTe3//cfHqvjw2DGKM4OkpnZT11VHm0fTFXJywtcNnV04OvwktJSTVrOf5D023GYbtuRsXOlFJGUVk5Y7BaszKSLPU4h4JYl9EovVVMqhDDeMU1rexAPP7T4j1mtnefAFfP3j9X2zcpp8TZzo9tPR2gSdHbg6j5N8eB8phyx4lJ3khCyc3kLcmUWkZReTmJots26EGEAS+yQ2nht7j8RwF4qHitVhcZDnziPPndd/bE+whwZfA3Wddf3Jvra9jiOtjfjam9CdJ0isPEBKuYlkbcdjTcKbOhV3RhGerBJSswsxWyNb+16IyUQS+yQWqz1ZhzPUMM5oYrWarf2Lq/oEQ0GaupuMRN9ZR21nHSfryzjdUou/oxndtAl37UaS99pJUQ5SE/PI8E7HmzWFtJwSEpJSKT3ZPO4rgIUAo7PS0dNBimN8KqZG5OKpUupa4P8AZuAZrfWTwx0vF08jJxblCsYq0rH21cmp66yj3ldPRWMFFbVHaG45TWtLI/7OFtLMVtItTqw6kR2VLur9mTSqbH7x9WtZOiUzAs9KiKFprTnecpyPqz7GbrazevrqUe2vcLZxqxWjlDIDh4FVQCWwDbhba71/qMdIYhfR9MmxKr659j0Cupk0ex1XXxSkpqmS1pZGHPixaBO5djcl6flkeKaQkzGdooK5pKSko6TejYiQZl8zH1V9REVbBWmONC7Pv/yMT6FjMZ6zYpYBR7XWx3tP/CJwCzBkYhcCovdpY9fJLvxdaYR0GvW+KeR4pnPTQjf3//49zLoOr72G/IIAxwN1HD71HlS9g3mHIsmSiCchl0xPMXlZM5lSMIfkRNn0W4xOV6CL0ppS9tbvxWKycGnupcz1zu3fh3k8RCKx5wIVA36uBJZHoF0RxwbO6LEMsbXgWA02nr+40MPa+687440kGArS2NXIsar9lFcfpKbpBFXtVRxvPQjlG2ALJFhTSE3MITO1mPysGUzNnYXHlRrWx2kRn3qCPext2EtpTSk9wR5mps1kWdayEe12FmmRGIr5MvBFrfUDvT9/DVimtf7uWcc9CDwIUFBQsLi8vDys84rJ7an3j/Ifbx8i1PvnpwC7NXJTNsf6aUBrTUNbHccrdlNx+jC1zSdo6qii0deOLxjEbrWQ4PKQlphDhqeIgsyLKM6aTnqCN+wNycXk1BPqYV/9PnbW7qQr0EVBUgGXZF9CmjPykxnGcyimEsgf8HMecOrsg7TWa4A1YIyxR+C8YhLr61V394TQgCayUzbHushKKYU3KQPv7KtZNvtqALafaODvnn2bFFVFqq2GSwp7CDSVc6D+ILuP/hmUGWV3k5yQQUZyAXnpUynMnEJmYgZJtiTp3ceprkAXe+v3srd+L12BLvLd+SzNWhr2OHokRCKxbwOmKaWKgSrgLuArEWhXxLG+Oe+v7Khk3fYKY8emCTJl82xbyppoDCZRr5MwB2Zyef50vnN5MS31p6g8dZiq2oPUN5fT1HiaU/UnOHLsXUImM8qWiN2VQpo7l1xvCQUZJWQlZox4M3IxMTX6Gtlbv5eDjQcJhAIUJhWyKGMR2YnZsQ6tX9iJXWsdUEr9PfAXjOmOv9Va7ws7MhH3+nrVty/Km9BTNgcbs1cmMykZ+aRk5DOHlQDoYA+t9dXU1ZzgdE3vME7LaZobd/JZ+Ra2Ko22JmB2JOJO8JKRkkeet4Q8T07/xuQ2c2wrXU6m6bPjKRAKUNZSxr6GfVS1V2FSJqZ5prEgfUFUhlzCJUXAhBiBMY/ZBwN0NJ2mqaaC+voT1DeW0dheSVNPO61000oPPRYbypZAp7bTEXQzJbOYZVNnkuNOJ9WZOm4Jf6KVqIg1rTU1nTUcbjrM4abD+IN+3DY3s9NmMyN1RmwuikoRMCEiZ8xj9mYLid48Er155HOJcafW+NubaKmvoq2ukrrGExw7fYI9lWWYzKc4Wb2XxrK3MNldmGwuzI4EkhK8pCflkp1aQF5KFqlODx67B6fFGbEx/M3HG/qvefh7Yl+iIhb6kvnx5uMcbT5Ke087FpOF4uRiZnhmkOvOHddpi2MliV2I8aYUNncq6e5U0ovnUgJsf/8orxzaR4puxatauSPDzkxvB41t1TQ319HU2EgdezlBD91mK8rqxGRz4XQmk5aYiTcpl2xPPtlJXjzOFJLtyVhNo5ul43HZ6Pv8Hur9eaQm8xBOV6CLqrYqytvKOdl6kq5AFyZlIt+dz/Ls5RQlF2E3T67aQ5LYhZgALi5J41cWG/WBVFosXuwzZlPa6efi2Wlcn5uAr7WWtsYa2htP09hyisbWKlq6amlqPUVrTRmH6GGXChA02VEWByabg0R7CskJ6UbiT84mMyWXLLeXZEfSoImqqdOPSUFIg0kZP4/EeAzhRPKNo7Onk9Mdp6nuqKaqvYr6rnoA7GY7BUkFFCYVUphUOOmS+UCS2IUIUySSzsDKmB6Xjcff2HdWoszDkZpHOlDc9yCt0d2tdLQ00N5UR2vLaRpaTtHUcZqWzgZaO6pprS/nED3sVkEAQiYbmO3YrS7cdg9JDg/JrlQ8iemk2hJItDXQ7ndjs9hHPEMp2lVGx/rGobWmo6eDRl8jdV11xq2zjjZ/G2Bs5ZiZkMmyrGXkufPIcGVMimGWkZDELkQYItlb7RvHf+r9oyNLlEqhHMkkOpJJzCzhnNnTAT/+jmbaWxtoaamjofUUTW01NHfW09bdTHvrKRqbjlGhAwSVMQizKi+IL6BxWh1s+HgdH25OxGVz47IlkmA3zuV2ppDk9OBOTCXJlcLiHDs2s6YnqKIyZXW4N45gKEhHoIN2fztt/jZa/a20drfS4m9hT3UVJxpayPM4yUlxkmRLIt2Vji1QRFmNlSunTmFZcfp5zj45SWIXIgyR6q0O7PWHW475zE8QGaQmZ5CaP/Pznv4AuseHr6OV+ubTNLXV0thRR0tHI22+Zjr8LXR0t9HYUUlVm49g0M9gc+hs2sSXCxQ9ARNJdheb//oGpVscWK1ObGYHNqsTq9WJxerEZnVisbmw2hxYbS4sNgdWmxOT1YbVbMJk0kAIjSagA/QEe0hJbcLuPkpQ+zFbggRcLbxwcBtdgS66A93os6JKsCbQ1GblxU966PF7seDmd19bySUlWZSWN/HwOuON+HebmuJ25o8kdiHCEIma+IP1+ofasGQsbQ33eGV14ExxkJ+Sccby8bNprfH1dNLS3kBTWwOtnY20dbXQ7muj09dKV08Hvp5Ouns6aQv46Am24u+sIxjsBh2EUGDQN4UBkRBSZkLKQkhZ0WYLmKzGxuZmG39TYqGl24I3KZU2XzeBgIsEWzpuqwu3LZFkh5tkm5tUVwoJNhu//7QMX8tRY1Wzgh3l7VxSMvE2p4kWSexChGG4XaNGarBk89CVUyPWViQSl1IKpy0BZ2oCWakFI36c1kbP29/Tja+7k+7uDrq72unxd9DT4yPQ3Umw20co4EP3+NB+H6ZgN8rvw9zThaWnGxVQoCFk1gTb6gm01tNlSqDLlIDP7MZnSeKkNZkuSzI+Sxshk4XTLT5MJgiFjNibO/39K5zNJgVB42u+x0VNq49EuwWXzRy18g/jPWtIErsQYRrrHPc+kdwJa6LtqqWUwqqsWO1WEuyJo28gFAR/R++tHbpbobsNfC1oXwvBzmZCvjqCIU0wpAkENX5TEr70NFZd6mRHs4O0rAIyvSl09QRJdlpZvTSf43Ud5HmcfHy0jhe3nSTP4yTP48LtsJDksJLstJLi6rvZ8LhsxhvCCA1M5MC4L/ySxC5EjEWi1x+NtiYEkxkcScbtLIreBBbwg68ZOhugox466qCjngJdxSqXBvZAZzIk50F2HiwoImhP4a/HG3jg+W34AxqrWfH962aS4XbQ5uvhWF07nf7g52EoRbLTgtdtJz3RTrrbTkaSg0T7uSn07OGwLy3KG/fhH0nsQkwA4fb6o9VWJEVtOMJig8QM4zZQoBvaa6GtGloqoakMaowyVmZHEq2VDvJDQcrIJBiy0NUT5IZ5nxfy8vUEaenqoanTT2O7n/oOP3Vt3Rypae8/xu2wkJXsIDvZSZ7HSXqi/ZzhMA3j/ilKErsQIurGMi007DcCix1S8o1b/jLQGrqaoOkENJUx33GQWy3H6A6ZqTRlcYXHC8FC44It4LCacVjNZCY5zmi2OxCkvt1PTauP0y3GrS/Z9yVwi0kRCBoVS7+0yNhERsbYhRBxZbQXdaOymlUpcKUat9zF5M0KcnnxXsoP7uROZy35TRvhrx9D+kzIng9JOcZjzmK3mMlNcZKb4uy/r83XQ1VzF1VNXZQ3dHLrwlwqm7qYmZ1Epz/A1IxEvnPFlHGrzS+JXQgRdaO9qDsu0xJNZubOmc/cOfON6TPN5cZQTe1+qP4M3FmQuwgyZoN5+FTpdliZkWVlRpZxLaC5009ZQyfH69rZUd7M9rImMpMcfGX5yGcUhUMSuxAi6kZ7UXfcZ/eYTJBabNymrYKavVC1Aw7+GY5/CPnLIWehMZ4/AikuGwtcNhbkp+DrCXKivoNAcPxKpEs9diHEhBTzipFaGxdcT242vlodkH8x5C3pH4cfb1KPXQgxqcV8do9Sn/fiW6qg/BM4/gFUlULxZZA5x+jpT0ATMyohhJhIknNh3p2w4CtgS4CDb8KO56HtdKwjG1RYiV0p9e9KqYNKqd1KqVeVUimRCkwIISYcTyEsvh9m3mSsgi19Do6+ayySmkDC7bG/A8zRWs8DDgP/En5IQggxgSkFWXNg2YPGtMiKbbD9WWMR1AQRVmLXWr+ttQ70/rgZyAs/JCGEmASsTph+HSy8x7jQunMtnPjImDoZY5EcY/8GsGGoXyqlHlRKbVdKba+rq4vgaYUQIoZSCmDJNyBjJpR9DJ/9Ebrbz/+4KDrvdEel1Ltw7uYswA+01ut7j/kBsAS4XY9g/qRMdxRCxKXTe+HwBrA4Yc7txurVCIrYdEet9dXnOdF9wI3AypEkdSGEiFtZcyAhHfa+DDv/ABddC9nzxj2McGfFXAs8Ctyste6MTEhCCDGJuTONmTPJ+ca0yLJPjDH4cRTuGPuvATfwjlJql1LqNxGISQghJjeby5j3njkbTmyCI2+P60XVsFaeaq2nRioQIYSIKyazMd/dnggntxg7QM26bVxWq0pJASGEiBalYMpVYHNDT8e4lSCQxC6EENGWv3RcTye1YoQQIs5IYhdCiDgjiV0IIeKMJHYhhIgzktiFECLOSGIXQog4I4ldCCHijCR2IYSIM+ct2xuVkypVB5SP+4mH5wXqYx3EKEy2eGHyxTzZ4oXJF/NkixdiG3Oh1jr9fAfFJLFPREqp7SOpczxRTLZ4YfLFPNnihckX82SLFyZHzDIUI4QQcUYSuxBCxBlJ7J9bE+sARmmyxQuTL+bJFi9MvpgnW7wwCWKWMXYhhIgz0mMXQog4c8EmdqXUl5VS+5RSIaXUkFe4lVJlSqk9vVv/bR/PGM+KY6TxXquUOqSUOqqU+v54xjhILKlKqXeUUkd6v3qGOC7Y+/ruUkq9HoM4h33NlFJ2pdSfen+/RSlVNN4xnhXP+eK9XylVN+A1fSAWcQ6I57dKqVql1N4hfq+UUr/sfT67lVKLxjvGQWI6X8xXKKVaBrzGPxrvGIeltb4gb8BMYDrwAbBkmOPKAO9kiBcwA8eAEsAGfAbMimHM/wv4fu/33wd+PsRx7TGM8byvGfAd4De9398F/GmCx3s/8OtYxThIzJcBi4C9Q/z+emADoICLgS2TIOYrgDdiHedQtwu2x661PqC1PhTrOEZqhPEuA45qrY9rrf3Ai8At0Y9uSLcAz/d+/zxwawxjGcpIXrOBz+MlYKVSSo1jjANNtH/j89JabwIahznkFuD32rAZSFFKZY9PdIMbQcwT2gWb2EdBA28rpUqVUg/GOpjzyAUqBvxc2XtfrGRqrasBer9mDHGcQym1XSm1WSk13sl/JK9Z/zFa6wDQAqSNS3TnGum/8Zd6hzVeUkrlj09oYzbR/m5H6hKl1GdKqQ1KqdmxDmaguN7zVCn1LpA1yK9+oLVeP8JmLtVan1JKZQDvKKUO9r6bR1wE4h2sFxnVaU/DxTyKZgp6X+MSYKNSao/W+lhkIjyvkbxm4/66DmMksfwP8ILWulsp9XcYnzauinpkYzeRXt+R2oGxvL9dKXU98BowLcYx9YvrxK61vjoCbZzq/VqrlHoV46NwVBJ7BOKtBAb2zvKAU2G2OazhYlZK1SilsrXW1b0frWuHaKPvNT6ulPoAWIgxjjweRvKa9R1TqZSyAMnE7mP6eePVWjcM+PG/gJ+PQ1zhGPe/23BprVsHfP9npdR/KqW8WusJUfdGhmKGoZRKUEq5+74HrgEGvUo+QWwDpimlipVSNowLfeM+y2SA14H7er+/DzjnU4dSyqOUsvd+7wUuBfaPW4Qje80GPo87gI269wpaDJw33rPGp28GDoxjfGPxOnBv7+yYi4GWviG8iUopldV3nUUptQwjlzYM/6hxFOurt7G6Abdh9BS6gRrgL7335wB/7v2+BGPWwWfAPowhkQkbb+/P1wOHMXq8MYu3N5Y04D3gSO/X1N77lwDP9H7/BWBP72u8B/hmDOI85zUDHgdu7v3eAawDjgJbgZIYv67ni/dnvX+vnwHvAzNiHO8LQDXQ0/s3/E3g74C/6/29Ap7qfT57GGaW2gSK+e8HvMabgS/EOuaBN1l5KoQQcUaGYoQQIs5IYhdCiDgjiV0IIeKMJHYhhIgzktiFECLOSGIXQog4I4ldCCHijCR2IYSIM/8/TVy53vspFqoAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "x_new = np.linspace(x_1s[0].min(), x_1s[0].max(), 100)\n", + "\n", + "for i, (sd, trace) in enumerate(trace_1[:3]):\n", + "\n", + " α_p_post = trace['α'].mean()\n", + " β_p_post = trace['β'].mean(axis=0)\n", + " x_new_p = np.vstack([x_new**i for i in range(1, order_1+1)])\n", + " \n", + " y_p_post = α_p_post + np.dot(β_p_post, x_new_p) \n", + " plt.plot(x_new_p[0], y_p_post, f'C{i}', label=f'model order {order_1} sd-{sd}', alpha=.5)\n", + "\n", + "plt.scatter(x_1s[0], y_1s, c='C0', marker='.')\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Questions for Osvaldo:\n", + "Are these plots right? At a high level seems to make sense that the very restrice SD \"Green line\" only fits intercept and nothing else, and that the widest SD has the largest slope because it has the largest coefficients. But could use your double check" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "How do we increase dummy data to 500 points as that data is provided in a csv and not generated?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercise 3\n" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "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: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [00:07<00:00, 689.37draws/s]\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (2 chains in 2 jobs)\n", + "NUTS: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [00:07<00:00, 643.79draws/s]\n", + "The acceptance probability does not match the target. It is 0.9126823488818671, but should be close to 0.8. Try to increase the number of tuning steps.\n", + "Auto-assigning NUTS sampler...\n", + "Initializing NUTS using jitter+adapt_diag...\n", + "Multiprocess sampling (2 chains in 2 jobs)\n", + "NUTS: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [00:15<00:00, 325.47draws/s]\n", + "The acceptance probability does not match the target. It is 0.893638190087912, but should be close to 0.8. Try to increase the number of tuning steps.\n" + ] + } + ], + "source": [ + "dummy_data = np.loadtxt('../code/data/dummy.csv')\n", + "x_1 = dummy_data[:, 0]\n", + "y_1 = dummy_data[:, 1]\n", + "y_1s = (y_1 - y_1.mean()) / y_1.std()\n", + "\n", + "traces_3 = []\n", + "for order_3 in [1, 2, 3]:\n", + " x_1p = np.vstack([x_1**i for i in range(1, order_3+1)])\n", + " x_1s = (x_1p - x_1p.mean(axis=1, keepdims=True)) / \\\n", + " x_1p.std(axis=1, keepdims=True)\n", + " \n", + " with pm.Model() as model_p:\n", + " α = pm.Normal('α', mu=0, sd=1)\n", + " β = pm.Normal('β', mu=0, sd=1, shape=order_3)\n", + " ϵ = pm.HalfNormal('ϵ', 5)\n", + "\n", + " μ = α + pm.math.dot(β, x_1s)\n", + " y_pred = pm.Normal('y_pred', mu=μ, sd=ϵ, observed=y_1s)\n", + "\n", + " trace_3 = pm.sample(2000)\n", + " traces_3.append([order_3, x_1s, trace_3])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Question for Osvaldo:\n", + "I don't get why order 2 model has lowest WAIC. It looks like the worst by far in posterior predictive check of mean values.\n", + "\n", + "See plots below" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 33, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XdwXfd16Pvv7zTgoPdCooMVBEmwgEWkKHZRomS12LItl1xbI+faeZO8vDvPztxxSSaJkzvjN/H185339BInjiMrji2bkij23hvYALABBAGi9w4c4Byc3/tjAyQlgkQ7DQfrM3MG7WDvhSNq4Ye1114/pbVGCCFE8DD5OwAhhBCeJYldCCGCjCR2IYQIMpLYhRAiyEhiF0KIICOJXQghgowkdiGECDKS2IUQIshIYhdCiCBj8cdJExISdFZWlj9OLYQQ01ZRUVGL1jpxrOf5JbFnZWVx6dIlf5xaCCGmLaVU1XieJ6UYIYQIMpLYhRAiyEhiF0KIIOOXGrsQwj+cTic1NTU4HA5/hyKeIjQ0lLS0NKxW66S+XxK7EDNITU0NkZGRZGVloZTydzhiFFprWltbqampITs7e1LHkFKMEDOIw+EgPj5eknoAU0oRHx8/pb+qJLELMcNIUg98U/1v5LHErpQyK6WuKKV2e+qYo3G5XZypPUP3YLc3TyOEENOWJ1fsfwbc9ODxRtXr7KW0tZQDlQcYcg95+3RCiACWlZVFS0vLlJ8zXseOHeOll16a0jG+8Y1vkJSURH5+vkdiGo1HErtSKg3YCfyTJ473NNEh0WzO2ExjXyNn6s54+3RCiBlsaGhqi0eXy/XY5/74j/+Yffv2Tem4Y/HUiv0fgf8TcHvoeE+VG5PL0sSlFLcUc6f9ji9OKYTwgMrKShYsWMDbb79Nfn4+b731FocOHWLdunXMnTuXCxcuANDW1sarr77KkiVLWLNmDdevXwegtbWV7du3s2zZMr71rW+htX5w7H//939n1apVFBQU8K1vfWvMpPz++++zePFi8vPz+e53v/vg8xEREfzgBz9g9erVnD17ln379rFgwQLWr1/P73//+wfP6+3t5Rvf+AaFhYUsW7aMDz/8EIB//dd/5fOf/zwvv/wy27dvf+y8GzZsIC4ubvIv4jhMud1RKfUS0KS1LlJKbXzK894B3gHIyMiY6mlZk7qGpr4mjlUfIz40nnh7/JSPKcSMUnYIeho9e8yIZJi79alPKS8v57e//S3vvvsuhYWF/PrXv+bUqVN89NFH/N3f/R27du3ihz/8IcuWLWPXrl0cOXKEr33ta1y9epW/+qu/Yv369fzgBz/gk08+4d133wXg5s2b/OY3v+H06dNYrVa+/e1v89577/G1r31t1Bjq6ur47ne/S1FREbGxsWzfvp1du3bx6quv0tvbS35+Pn/913+Nw+Fg7ty5HDlyhDlz5vDmm28+OMbf/u3fsnnzZn7xi1/Q0dHBqlWr2LrV+NnPnj3L9evXvZ7An8QTK/Z1wOeUUpXAfwCblVL//tknaa3f1Vqv1FqvTEwcczjZmMwmM9uztmMz2dhXuY/BocEpH1MI4X3Z2dksXrwYk8nEokWL2LJlC0opFi9eTGVlJQCnTp3iq1/9KgCbN2+mtbWVzs5OTpw4wVe+8hUAdu7cSWxsLACHDx+mqKiIwsJCCgoKOHz4MBUVFU+M4eLFi2zcuJHExEQsFgtvvfUWJ06cAMBsNvPGG28AcOvWLbKzs5k7dy5KqQfnBjhw4AB///d/T0FBARs3bsThcHD//n0Atm3b5rekDh5YsWut/xL4S4DhFft/01p/5anf5CHh1nC2ZW3jo7sfceT+EZ7Pel5auYQYrzFW1t4SEhLy4H2TyfTgY5PJ9KAm/WiJZcTI/9uj/T+utebrX/86P/7xj8cVw2jHHxEaGorZbH7svKMd44MPPmD+/Pmf+vz58+cJDw8fVxzeMu372GdHzGZt6loqOiu42nzV3+EIITxgw4YNvPfee4DRiZKQkEBUVNSnPr93717a29sB2LJlC7/73e9oamoCjBp9VdWTJ9yuXr2a48eP09LSwtDQEO+//z7PPffcY89bsGAB9+7d4+7du4BRlx/x/PPP87Of/ezBL4krV6544Cf3DI8mdq31Ma311HqBJmFp4lJyY3I5V3eO6u5qX59eCOFhP/rRj7h06RJLlizhe9/7Hr/85S8B+OEPf8iJEydYvnw5Bw4ceHC9Li8vj7/5m79h+/btLFmyhG3btlFfX//E46empvLjH/+YTZs2sXTpUpYvX84rr7zy2PNCQ0N599132blzJ+vXryczM/PB177//e/jdDpZsmQJ+fn5fP/73x/Xz/alL32JtWvXcvv2bdLS0vjnf/7nibw046Ke9ieJt6xcuVJ7eqMN55CT35X9jn5XP38074+IskV59PhCBIObN2+ycOFCf4chxmG0/1ZKqSKt9cqxvnfal2JGWM1WXsh6Abd2s79yP063098hCSGEXwRNYgeICY1hS8YWmvuaOVFz4qkXSIQQIlgFVWIHyI7OZmXySm633aa0tdTf4QghhM8FXWIHKEwpJCMqg5O1J6nvefIFFCGECEZBmdiVUmzL3EaULYr9lfvpGezxd0hCCOEzQZnYAULMIbyQ/QJOt5N9lftwuR8fxiOEEMEoaBM7QFxoHFsyttDU1yQXU4UIQtNtbG91dTWbNm1i4cKFLFq0iJ/+9KceieuzgjqxA+TE5LAieQW32m5R0lLi73CEENOIp8f2WiwWfvKTn3Dz5k3OnTvHz3/+c27cuDGlc4wm6BM7wKqUVWRFZXGq7hS1PbX+DkeIGWumj+1NTU1l+fLlAERGRrJw4UJqaz2fk6Y8BGw6UEqxNXMrH5R9wP7K/XJnqhDAqdpTtPR7pkQxIsGewPrZ65/6HBnba6isrOTKlSusXr16oi/zmGbEih3AZrY9uDN1b8VenENyZ6oQ/iBje6Gnp4c33niDf/zHfyQqyvOLzBmxYh8RExrD9sztfFLxCYerD/N8poz5FTPXWCtrb5npY3udTidvvPEGb731Fq+//vq44p2oGbNiH5ERlcHaWWup6KjgUqNnB5EJITwjWMf2aq355je/ycKFC/mLv/iL8bwUkzLjEjsYY37nx83nYsNF7nbc9Xc4QojPCNaxvadPn+ZXv/oVR44coaCggIKCAvbs2TPRl2dMQTO2d6Jcbhcfln9Iq6OV1+a8RmLY1LfrEyLQydje6UPG9k6CxWThhewXCDGHsPfeXvqcff4OSQghPGLGJnaAMGsYL2a/iGPIwd57e2XsgBAiKMzoxA6QGJbIlowtNPY1crz6uIwdEEJMezM+sQPkxuRSmFLI7fbbXG667O9whBBiSmZUH/vTrExeScdAB+frzxMbEktOTI6/QxJCiEmRFfswpRSb0jeRHJbMofuHaO5r9ndIQggxKZLYHzHSKRNqDmXPvT2yQYcQAW66je11OBysWrWKpUuXsmjRIn74wx96JK7PksT+GWHWMF7MeRGn28mee3tkpowQM5inx/aGhIRw5MgRrl27xtWrV9m3bx/nzp2b0jlGI4l9FAn2BLZlbqO1v5WDVQdxa7e/QxIiKMz0sb1KKSIiIgBjZozT6fTKvCq5ePoEmVGZrJ+9npO1Jzlbd5Z1s9f5OyQhPKrn5ElczZ4d22tJTCDi2Wef+pyZPrZ3aGiIFStWUF5ezne+853AHNurlApVSl1QSl1TSpUqpf7KE4EFgsWJi1mcsJhrzddk9yUhPGSmj+01m81cvXqVmpoaLly4QEmJ53OLJ1bsA8BmrXWPUsoKnFJK7dVae75w5AfrZq+ja7CLkzUnibRFkhmVOfY3CTENjLWy9paZPrZ3RExMDBs3bmTfvn3k5+ePK+7xmvKKXRtG2kesw4+guX3TpExsz9xOvD2eA5UHPL7jjBDiccE6tre5uZmOjg4A+vv7OXToEAsWLBjz+ybKIxdPlVJmpdRVoAk4qLU+74njBgqr2crOnJ3YzDY+qfhE2iCF8LJgHdtbX1/Ppk2bWLJkCYWFhWzbtm1K7ZNP4tGxvUqpGOAPwP+mtS75zNfeAd4ByMjIWPG036aBqqW/hV3lu4i0RvLa3NewmW3+DkmICZGxvdNHwIzt1Vp3AMeAHaN87V2t9Uqt9crExOk5+zzBnsDzmc/TNtDG/sr9DLmn1uMqhBDe4ImumMThlTpKKTuwFbg11eMGqvSodDambaS6u5rjNTINUggReDzRFZMK/FIpZcb4RfGfWuvdHjhuwFoYv5DuwW4uNV4i0hZJYUqhv0MSYty01rKJe4Cb6oJxyolda30dWDbV40w3hSmFdA92c7HhIhHWCBbGS91SBL7Q0FBaW1uJj4+X5B6gtNa0trYSGho66WPInaeTpJRiY/pG+lx9HKs5Rrg1nIyoDH+HJcRTpaWlUVNTQ3OzTC8NZKGhoaSlpU36+yWxT4HZZOb5rOfZVb6L/ZX7eWXOKySFJfk7LCGeyGq1kp2d7e8whJfJELApsplt7MzZSagllE8qPqFzoNPfIQkhZrhpmdjdvb3+DuFTwq3hvJTzEm7tZnfFbvqcff4OSQgxg027xD7U00Pbr35F1969uNra/B3OA7GhsezM2UnPYI/McRdC+NW0S+zKZsO+bDmDVfdp//X7dB08yFBnYJQ/UsJT2J61nea+ZvZV7pMbmIQQfjHtErvJZiN89Srivv417AUFDJaX0/bee/QcPx4QJZrs6GyeS3+O6u5qjlYflRuYhBA+N227Ykx2OxHr12EvWErfxYv0l5TguHkLe0EB9uXLMNn8N8clLz6Pflc/5+vPY7fYeWbWM9IzLITwmWmb2EeYIyKI3LQJe0EBfefP03fxIo7SEsIKCwldtAj1yFxlX1qetJw+Zx/Xmq9ht9hZnrzcL3EIIWaeaZ/YR1hiY4nasQPnskZ6T5+h5/gJ+q9eI/yZtdhyc32+YlZKsX72ehxDDs7VnyPUEkpefJ5PYxBCzExBk9hHWJOTiX7tVZxVVfScOUPX3n1YU1MIX7cOa2qqT2NRSrE5fTMOl4Pj1ccJNYeSE5Pj0xiEEDPPtLt4Oh5KKWxZWcR+8YtEbN7EUFc3Hb/7gK59+33eQWM2mdmRtYOksCQOVB2gurvap+cXQsw8QZnYRyiTCfuiRcR95S3CCgsZrLxH23vv0XvmDO7BQZ/FMbIDU0xIDPvu7aOxt9Fn5xZCzDxBndhHKJuN8DWrif3KVwiZO5e+osu0/+pX9JeWot1un8QQagnl5dyXsVvs7K7YTZsjcG6uEkIElxmR2EeYIyKI2raNmC98AXNMDD1HjtLxn/+Js7bWJ+cPt4bzcu7LmJWZj+9+LHNlhBBeMaMS+whrchLRr79O5PPbcTscdPz+D3Tt28dQd7fXzx0dEs3LuS/jcrv4+O7H9Dr9f1OVECK4zMjEDsYF1tB584h76y3CVq1isLKS9vfeo/fCBbTTu3Ne4u3xvJTzEv2ufj66+xH9rn6vnk8IMbPM2MQ+QlmthK9eRexbb2HLyqLv/AXaf/1rBu7e9eo4gOTwZF7MeZGugS52V+xmYGjAa+cSQswsMz6xjzBHRhK1YwfRr70KVitde/bS9fHHuNrbvXbO2RGz2ZG9g9b+VvZUyERIIYRnSGL/DFtaGrFvvknEs+tx1jfQ/v779J45g/ZSe2RmVCbbMrfR0NvA3sq9uNwur5xHCDFzSGIfhTKbsRcUEPeVtx60R7Z5sTyTG5PL5ozN1HbXsr9yv4z7FUJMiST2pzCFhxvtka+/hikk5EF5Zqijw+Pnmh83nw1pG6jqquJg1UHc2jf99UKI4COJfRyss2cT89nyzPkLaJdnyyaLEhaxbvY6KjorOFR1SJK7EGJSgm4ImLcokwl7QQG2OXPoPXWavgsXGLh9m4iNz2HLyPDYeZYmLsWt3ZytO4vZZGZz+maZ5S6EmBBZsU+QOSKCqB3PE/3K50ApOj/8iK79Bzy6e9OypGUUphRyu+02x6qPyS5MQogJkRX7JNkyMoj90hfpK7pMX9ElBquqCF+7ltD8RR5ZYa9MXolbuylqLMKkTGxI2yArdyHEuEhinwJlsRC+ehUh8+bSc+w4PceO4bh1k8hNm7AkJEzt2EqxKmUVbu3mStMVlFI8O/tZSe5CiDFNuRSjlEpXSh1VSt1USpUqpf7ME4FNJ5bYWKJffYXIbVtxd3bS/pvfGL3vUxxNoJRiTeoaChILKGkp4WTtSSnLCCHG5IkVuwv4P7TWl5VSkUCRUuqg1vqGB449qqKqds5VtLImJ54VmbHeOs2EKKUIXbAAW2YmvWfO0Fd0mYHyciI2bpzSxVWlFGtnrUWjudZ8DYWx5Z6s3IUQTzLlxK61rgfqh9/vVkrdBGYDXknsRVXtvPVP5xh0ubFZTLz39pqASe4AJrudyC1bCJm/gJ6jR+n88CNC5s8jYv16TGFhkzqmUopnZj0DwLXmawCS3IUQT+TRrhilVBawDDg/ytfeUUpdUkpdam5unvQ5zlW0Muhy49bgdLk5V9E66WN5ky1tNrFf+iJhhYUMlJfT/utf47h1a9KllJHkvjRxKcUtxZyqPSVlGSHEqDyW2JVSEcAHwJ9rrbs++3Wt9bta65Va65WJiYmTPs+anHhsFhNmBVaLiTU58VOI2ruUxWLs3PTmm5hjYug+eIiujz6a9L6rI8m9ILGA4pZiqbkLIUalPJEYlFJWYDewX2v9f431/JUrV+pLly5N+ny+rLF76lxaaxzFxfSeOQtowlatxl6wFGWa+O9WrTVn689ytekqefF5PJf2nJRlhJgBlFJFWuuVYz1vyjV2ZWSUfwZujiepe8KKzFif1NU9Wc9XSmFfsgRbdjY9x0/Qe/o0A2VlRG7ZPOHWSKUUa1PXYlImLjdexq3dbEzfiEnJ/WZCCM+UYtYBXwU2K6WuDj9e9MBx/c4b9XxzZCRRO180tuXr6TZaI8+dm/DcGaUUq1NWszJ5JbfabnH4/mGZLSOEADzTFXMKCNg6wFRKKSP1fKfL7dF6/si2fLb0dHpOnaLv4iUGyu8SuWUz1tTUCR1nVeoqzCYz5+vPM6SH2JaxDbPJ7JE4hRDTk0dq7BM11Rr7eHmilOKLev5gVRXdR4/i7unFvmQx4WvWoGy2CR3jatNVztSdISsqi+1Z27GY5KZiIYLNeGvsQV2U9UQpZUVmLN/ZNMerNX1bZiaxX/4y9sX59F8vpv399xmsrp7QMQqSCtiQtoHKrkr23tuL0y3b7AkxUwV1Yp9OrZEmm42I554j5vXXwGSmc9eHdB85gntg/Jtc5yfkszljMzXdNey+u5vBIe9s5yeECGxBXYqBwBw/MBbtdNJ38SJ9l69gCgsjYtNGQrKzx/39Ze1lHLp/iAR7Ai/lvITdYvditEIIXxlvKSboE/t05mxsoufIYVwtrYTMm0fEhmcx2ceXpCs7K9lfuZ/okGhezn2ZcGu4l6MVQnib1NiDgDU5iZgvfIGwVasYKC+j/de/ZqC8fFzfmxWdxc6cnXQPdrOrfBddg4/dDCyECFLTM7EPem63Ik8oqmrn50fLKapq9/ixldlM+OpVxL75JqbwCLr27qNr717cfX1jfm9aZBov576Mw+XgD2V/oM3R5vH4hBCBZ/qVYgZ74ezPITIVUhZD4gKwhno2wAnw5bRJPTRE/5Ur9F64gMlmI3zDBkLmzh1znEBrfyu7K3bjcrvYmbOTlPAUr8QnhPCu4C3FKBNkPQvOfri9F878DEp3QetdcPv2zsuiqnb+8dAdBpy+mTapzGbCVq40Vu9RUXTvP0DXnj1j7rcab4/n1TmvEmIO4eO7H1PdPbFWSiHE9DL9VuwjtIbuemgogaZScDrAFg7JeZCyBCKSPBPsEzy6Undr4zekzeq7+fDa7ab/6lX6zp8Hi4WIDc8RMu/pq/deZy+77+6mfaCdrRlbmRM7x+txCiE8x2dDwPxGKYiaZTzmbIHWcmgohpoiqL5oJPaUxZCUByERHj/9ozc/mRSsm5PAn2+d57OWSmUyEbZ8uTFU7PBhug8cYKCsjIiNGzFHjN4BE24N59W5r7KnYg8Hqw7iGHKQn5Dvk3iFEL4zfVfsTzLYC023oOE6dDcYpZu4HEjJh/i5YPbM77KRFfvIHBl/7uSk3W76r12j79y54dX7BkLmzXvi6t3pdnKw8iCVXZWsTF5JYUqhjP0VYhqQPnaA3hZjFd9YCgPdYAmBpIXGSj5qtrHqn4JAu/nJ1d5Oz+HDOOsbsGVnP3X17tZujlUf41bbLfLi89iQtkHG/goR4CSxP8rtho4qI8m33IYhF9hjjVV88iLj/SDxqdW72ULEhmcJmT9/1BW51poLDRcoaiwiKyqLbVnbsJqsfohaCDEektifxDUAzbehsQTaq4zPxaQ/bJ20hPgnLg/71Oo9K4uITZueuHovaSnhZM1JksOTeSH7BRlBIESAksQ+Hv0dRpmmsQT62oz6e8I8SM6H2GyYxLZ1gWQiq/eKjgoOVh0kwhbBSzkvER0S7YeIhRBPI4l9IrSGrjojwTfdMFonQyKMMk3yYoiY/ObbgWC8tff6nnr23NuDSZnYmbOTpDDvtowKISZGEvtkDbmM1snGEuOmJ+2GyGQjwSfnGb3y09B4O2faHe3srthNv6uf7ZnbyYrO8k/AQojHSGL3hMFeaLwBjcXQ3Wi0TsbnGqWa+Dkea530pfGs3vucfXxy7xNa+lp4Nu1Z6XUXIkBIYve0nmYjwTeWwkCPMZ8mKc9I8lGzptw66UvGXavX6Dv/5NW7c8jJwSqj131p4lLWzlor7ZBC+Jkkdm9xu6Gjcrh18o5RugmLMxJ88iKwx/g7wnFztbfTfegQroZGbDnZRG7ciCn84erdrd2crj1NcUsx2dHZbM3YitUs7ZBC+Iskdl9wDUDzLSPJdwwP1orNNJJ84vxp0Tr5+Or98Zkz15uvc7r2NAlhCbyY/aJs2iGEn0hi97X+dqNM01BivG+2QMJ84yaomKyAb50ca/Ve2VnJwaqD2Mw2Xsx+kcSw6d0pJMR0JIndX7SGrlpjFd9001jVh0QaZZqUxRCe4O8In2isiZEt/S3sqdiDY8jB1syt5ETn+DliIWYWSeyBYMgFrWXGKr6tYrh1MsUYK5y0EGxh/o5wVK62NroPHx519d7n7GPPvT009zWzOnU1y5KWyQAxIXxEEnugGegxbn5qKIaepoetkymLjdZJk9nfEX7K01bvTreTo/ePUt5RzrzYeWxM34jFNP1aP4WYbnya2JVSvwBeApq01mM2Pc/IxP6oniZjrHDjDaNX3hoKSYuMenxkqldaJyc7ifJJtXetNUWNRVxouEBSWBIvZL8gF1WF8DJfJ/YNQA/wb5LYJ8DthvZ7w62TZeB2GTX45HzjLtdQz8xrmeq+rI+v3h/2vVd0VHD4/mFsZhs7snaQHJ7skZiFEI/z6Q5KWusTSqksTxxrRjENl2Pic435NCOtkxXH4N5xiMk0VvEJ88Fim/RpHt3taWRf1okk9sd3azrIQFk5ERs3khOTQ1RIFHvv7WVX+S42pm9kftz8SccqhJg6j9XYhxP77iet2JVS7wDvAGRkZKyoqqryyHmngwmXQfraHk6d7O8As9UYKZySbyT7CZZqPLnb06f63h+ZGOkYcnCg8gC1PbVyp6oQXuLzi6djJfZHzaRSzJTKIFpDZ42xim++Ca5BCI16OHUyPH5CcXhytydXezs9R47grKvHlpVJxKZNEGbnTN0ZiluKmR0xm22Z2wizBmbnjxDTkST2APHzo+X85MBt3BrMCv5i+3y+s2nOxA805DTq8I0jrZPamFGTkm/MrLH6fnMM7XbjuH6d3nPnwGQmYv06QhYu5Hb7bY5XH8dusbMje4eM/xXCQ8ab2OVvZS9bkxOPzWLCrMBqMbEmZ/yr7E8xW40Lqku+AGu/A7mbwe2EOwfgzM+g5IPhC7BDnv0BnkKZTNgLCoj94hexxMfTffgIXR9/zFzrbF6b+xoAfyj7Azdbb/osJiGE57pi3gc2AglAI/BDrfU/P+n5M2nFDl7c9Fpr6Gk0VvGNpTDYZ6zckxcZnTWRKT6bOqm1xlFSQu/pMwCEr1uHe342h6sPU9NdQ158Hutnr5d+dyGmQG5QmkGKqto5f7eJDQk95JvvP1y5hycYN0Al5Rm1+Qkcb7K/iIa6uug+cgRndQ3W2bMJ3/Qcl/pvc6XpCklhSWzP2k6UbfyxCCEeksQ+Q4x6cXZWqDGnprEEOmuNVXts1sOpk08ZvTvVnncYXr3fuEHvqdOg3YStXk19ViRHa46hUGzN3EpmVOYUf3IhZh6psc8Qo/WoY7XD7OWw/Guw+luQ+YzRQnnzYzj9U7j1CbRXGaWc8RxvgpRS2BctIvatL2OdnUbvqdPEHbrC6/FbiLRFsqdiD+frz+PWbk+8BEKIz5CC5zQ3cnF2pEf9sYuzYXGQvQGynoWO+8YqvvkW1F837mwdmToZFje+402AOSKCqJd2MnCnjN6TJ3Dv2suOFcu4lBRHUWMRDb0NbM3cKqMIhPAwKcUEgQnXxIecxu5PDcXQXmms3KNnG6WapIUU1Tk8frHX3ddHz4mTDJSVYUmIp2F5JicHSrCZbGzJ3EJ6ZLpHziNEMJMauxgfR9fDqZO9LWCyQMIc4waouGyPT50cqLhHz7FjuPv6cOblcDyplQ5XNyuSV7AyZaXcrSrEU0hiFxMz0jrZUAJNw62TtrCHUycjkj3WOukeGKD39BkcpaUQFUHpAjs3QlpJDU9lW+Y2ImwRHjmPEMFGEruYPPeQcXdrQzG0lhsfRyQaq/jkPGNHKA8YrKml58gRhjo7aU2P4nRaD4SEsCljk+zOJMQoJLELz3D2D5dqSqCrbrh1Mtu44Jow96mtk+OhnU76Ll6k78oVBixQlKOpTVDkJy5m7ay1WE1TO74QwUQSu/C83lZoLDbucnV0GaOEExcapZrodFBq0jc3OZua6Dl6jMHGBu7Huricq4iOS2Vr5lYS7IG7T6wQviSJXXiP1tBRZazim28ZXTb2GG6Tydc+6aXZFTaFDT2u0XfhPB0DnVxJH6ItO441aWtZkrBE9lYVM55PN9r6r8h8AAAZWklEQVQQM8zInayxWTB3O7TchoYSOq8c48u6mRpTAneGMrhYnjGJDT2WEZKbg/XYMdZVVnC3oY4LeQe5n36fzRmbpeddiHGQxC6mxmIz6u0pi7GGr+N8xQfM0/fYZrnCyz3tUFpmfD0229gxahzM0dFEfe5zhNwpw3byBHEXayitOcdvuurZkLWZObGTGHssxAwipZhpzmuTI6cSz90Wnk1xscRSbbROOh1gCzc6alKWQMT457O7HQ56T5+m/fplyl31VOUnkJ63ivWz1xNqCfXiTyJE4JEa+wzgiYFdXuceMlomG4qprSiltq2H5FkZZOatNqZOhoyvZ32wppauo0eorblBWewA/SsX8Ny8HWREZXj5BxAicEiNfQaY6ibVPmEyQ+J8ivqS+OYlC1lDVSy+V8V/7W9mVsxR4+7WlMUQPxfMT/7naEubTfyXvkTYlStEnz3O3X3XOF5ZQ0bhJtamrcNmnvxm30IEG0ns05gnB3Z527mKVrpcVq7qORQ755BlT+Kb6T1G62TpLrCEQNJCY15NdNqod7kqi4WwwkLS5swh6vgxYm+c437l7/j9mpusW/qyzJsRYpgk9mlsRWYs7729JqBq7E/y2V9CBQvmQGYsZD9ntE42lhiPuqtgjzV645MXGe9/hiU2lthXXiUsbxFxh/dy99A1zty5R8pzz7M2a4Os3sWMJzV24TNjXuh1DUDzbSPBd9w3+uVj0h9MncQS8ti3uAcG6Dp7moqz+6lztdK/fD6F6z9PVnSW938gIXxMLp6K6c3RadwA1VhibBJiskDiPCPJj9I66Wxqom7/R1SUX6Q91krMxs2szX8Bu8Xupx9ACM+TxC6Cg9bGjJrGEmNmjdNhdNIkjbROJj7yVE1vyXXKD/6euvYqHPPSydv6Bean5MtdqyIoSGIXwWfIBW13h6dO3gXthsjkh1MnbcZdqe7+fuqPH6D8/AG6TAPY1hSy+tk3iQmN8fMPIMTUSGIXwW2w19iwu6EYuhtAmSA+1yjVxM8Bs4XBxgbu7H6f2nvXccRHMLBsK43uBTyTmxTQF5qFeBJJ7GLm6Gl+OHVyoAesoUapJjkfHZlKe/EVTv/hV1wqu0t5TCLFcc/xL29/TpK7mHbkBiUxc0QkQsRmyN4IHZXGKr7hOtReRoXFEZecT+OqL3Gj9hPmtxWT2fUBhz5pYME3vk14qGc2DREikMgGkyJ4mEwQlwN5r3A57av8pnsxld0K7p3gxcFPyE0d5M6sNXTYZpNdUczRn32P4qsHcWu3vyMXwqNkxS6CTlFVO2/9yxUGXQqbJZ3/+Op2CrJr+XbIeRoa6pkVl0qILY8b125S/dtfUXf+GIt2foW0tIX+Dl0Ij/DIil0ptUMpdVspVa6U+p4njinEZH12hs7pWhdkrSNjx//Oqtf/jLQFq0iMcPLsilnMzY7BVHmH6//rbzj5u5/S2dns7/CFmLIpr9iVUmbg58A2oAa4qJT6SGt9Y6rHFmIynjhDRyljDk10GszZiqm1jLnJJWSk3aSk5BZNJ/dy+vJpkp97kcWb3sRmk5ubxPQ05a4YpdRa4Eda6+eHP/5LAK31j5/0PdIVI8C7s+QndOyBHmi6QfvN45ScO0N3YwemmFiytnyOeRvewmSRDbVFYPBZu6NS6o+AHVrrt4c//iqwWmv9p5953jvAOwAZGRkrqqqqpnReMb0F7Cz5niaqz/6eO8f2MNDRhS0+igUbd5K28mWITB116qQQvjLexO6JGvto/9If+22htX5Xa71Sa70yMTFxlG8RM8los+QDQkQS6dv+hE0/+oCcV7/OkNPC9Q/+gxM//Q6t+/8Bqs4Yc2yECGCe6IqpAR4dhJ0G1HnguCKIPVoHN5tN1Hb0U1TVHhirdsBktbJg238hZ8OblB75LQ0n93B+33Fir10mP38+kbPzjNHCCfONfV+FCCCeKMVYgDvAFqAWuAh8WWtd+qTvkRq7AKMc8/vLNfz2UjUutw6sksxn9HS1UnLoP2krOoNy9pKQHsuiOWmEhUdC4gIjycdkSqlGeJXP7jzVWruUUn8K7AfMwC+eltSFGLEiM5ZzFa243Dqwt/cDIqLiWfP6f6X92ZcpPfSfNJde43hTBclz01jkLCakoRhCo4zNQZIXQ3jg7mYlgp9HblDSWu8B9njiWGJmmS7b+z3aZbP+S39BU20ZNw/+lobbt2mqbCF18SIWZkRju38Oqs5C1CxjFZ+UB1ZpmxS+JUPAhN95s+3RE57WwVNXWcLtQx/Qf68Cc1gYqYXPsCA3E1vbHWM4mck8PHVysfHWZPbzTyOmMxkCJqaNFZmxAZnQR4zWwTMS76ysfGa9nU9N+VXuHP49NccPUXcxjNRnNrNwyVZs7eXG1MnmO8bKPXmRMVo4MkXq8cJrJLELMYbxlIvS5hSQNqeA6rLLlB/5kNpDu6k7c5iUNRtZuPa/ENrXZIwWrrsKNZcgPMFI8MmLjNq8EB4kpRghxmGi5aLa8quUHdlFX9U9lD2EhML15K17mfAQOzTfNPZz7awxVu2xWUaST5gnrZPiqWSjDSECQH1FCeXHPqTn7h20zUrsskIWPPsKsbEpxibdjSVGknd0gtkKSQuNJB+TIaUa8RhJ7EIEkJb7dyg7/hEdt0vQJkV4/hLmPvsSs1LnGht2d9wf3rD7Jgw5ITTaKNOkLIawOH+HLwKEJHYhAlBXUy23j39Ia0kR7iEX1twcste/QM6clZiUyUjqLXeMVXz7PSPpR882VvFJC6V1coaTxC5EABvobOfO6U9oLDrNoKMPUhJJXb2RBcu2YB9J3gPdRkdNQzH0tgy3Ts6BlCUQl+2V1slAbz2d6SSxCzENuBz9VF46Qu25ozTU19FsshC/chWbn3+V5Jg040laQ0+jkeAbS8HZD7YwSFpk3AQVkeyRenzATtwUD0gfuxDTgCXUzpz1O+mcvYZ/+p/vM7f1KkkfHoAr54lbmkfG6i3MyVqGNTLF6H3P3QxtFUaSr7sMNReNzbyTF0NyHoRMfnPup/Xri+lFErsQAeBcZTt3Q7MpS8kmaaCFWbH1JN+p4l7Jz7mTGk/C8jXMW7qRxPAkSJhrPJz90HTDqMffPQIVRyE2e3jq5Dyjy2YCpst4BzE2KcUIEQBGyiAjSfW9t9dQEG+l9vIp6i6doKOtDqfdimneHDILNzInbSmhltCHB+hrGy7VlICjy+iHT1xoJPno9HGXaqTGHtikxi7ENPOkpKrdbrrLb1N18Qit5SX0OftwpMQSvXgZOYvXkxmbbXTUwHDrZJWxim++ZXTZ2GMejjIYo3VSEntgk8QuRBBydXbSeOUs9VfP0NZay6AVBrNnkbRkFbm5K0kJT0GNrM5dg8Otk8VGstfa2Mg7Jd9YzVtDP3VsuXga+OTiqRBByBIdzeyNO5i1YTuOykpqrpyk5dY1Om9/wPmoPejcDFKXrGZO2hLiQ+NRKflGInd0PWydvL0Pyg5BwnDrZGw2mExTungqK/3AIoldiGlImUzYc3KYm5NDrsNB961S6q6cpq20jLYr/8aJhEhUbgap+avITc4jwZ6AylwLGWugu8FI8E03oOkW2MIhOY/1Ken8bBIXT2Wl/3Su9nYcJaWg3URs2OCTc0piF2KaM4WGEl2wguiCFQx1dNB5o5j64vO0X6mi7UIptYlRqOx0kheuIDtlASmRKZiiUmHOFmi9a0ydrL3MUvdF9q+P4GL/bHLzV7FsnMlZ2iQfp10uBisr6S8pwVldAyZFyLx5aK0flsq8SBK7EAFosqUNc0wMcc88S+za9Qy1tNB1u5SGkot0XKmm83wp5xIicaYnkTB/KRmzF5Iem0Fo4jwY7IOmm2Q2FpPZdRtqyqAve3jq5Nyntk76ok1yOpR6tNa4mpsZuHULx+3baMcApsgIwtesJjQvD1N4uM9ikYunQgQYT5c2tNYMtbbSe+cWDTeK6GyoomOgk75IK47UWCKy55KSnU96TCaJYYmYHp06OdANlpCHUyej00ZtnfRm4g30Uo+rvZ2BO2UMlJUx1N4OZhMhOTmELlyINT0dZTJ57Fxy8VSIacrTpQ2lFJaEBKIT1hP9zHqGOjtxVFTQcvsabZV36Lp7kQZOUZkUxVBqAnG5ecxKncfs1C8T3deBaio1En3dVbDHGhdjkxcZ7w/z5i5YgVbqGVmZD1bcY6DiLkOtbaAU1lmzsBcUEDInF1No6NgH8iJJ7EJMkadXq94ubZijowlftozwZctIHxzEWV1Nd0UZLeUldJXW0XVlD3dD9lKaGIVKTSIuawEp87Yxa6CXmLZKVOUpuHcSYtIfTp20hHg0xkcFwh2xbocDZ20tg1VVDFZW4e7tfZjMNzyLLTcXc0SEz+N6EinFCDEF3ioT+KOmrLXG3dnJwP37dFaW0VZ5m+7OJroHe3DYYCA+Ap0UR1RyKgl2E8n9bSQ5B7CZQyBxnpHkh1snPc3Xr4d7YABXfT3OujoGa2pwNTWD1iibDVtGOrasLGyZmZjCwrwey6OkFCOED3irTOCPDb6VUphjYgiLiSFsyRJStGaovR1nXR2d98tpq7pDz61GeopraRxyUB1pxxlhIiRMEVVbRlz4EZJiUoiftQJr6jJjOJmHePL1+OwvCe1y4Wprw9XYhKu5CVdjI67WNuOGLpPCmpJCWGEhtvQ0LMnJKLPnxyV7miR2IaYgEMoE3qKUwhIXhyUuDnt+PimAu68PZ2MjvfXVtNdU0F1/n97qFvoG+6kZrKOKCly241gjQgiLTyIqfRExOYXEps4hNiLx4egDH/hsAncPDnKl9D7f+5eT2Pu7ueTq489XJDBLDYDbqFyo0BCsycmE5eZinTULa3IyyjqxYWqBQEoxQkyRJ8oEniw1+Lxs0deHs6WFnsZa2uvv0d1YRX9tOY6Oegac/WgAqx13ZAzWuGTs0fHYo+IIi4onPCqOyMgEwiNiMYfaUTYbympFWSwoiwVMpse6SrTW4Hajh4bA6UQ7nbgHB9GDTrSjH3e/g1uVTfxk12WsAw4i9CBvL01glt3Exco2ztxtwY2i3xrKusJ5vLBhEZb4eCxJSZiiojzSZ+7WbnqdvXQNdtE50EnnQCc2s40VySumdFwpxQjhI1MtE3iyTu+P1kBTWBghGRmEZGQQz9oHn3f39zNQU0Z7+Wm6qkvo7myjv68GR3cjXUNmTK6HCVspsJls2MzGw2qyYTNbsZpsWM02rBYbFmXGoswwjrVofWUbqR1d9FlC6LOEcDM0gbnrFpCyRHNw913aTaGYbFa+8uYawifw+mitGXQP0u/sp9/VT5+rj15n74NHj7OHnsEeup3dPLpoVkqRFpE25cQ+XlNK7EqpzwM/AhYCq7TWsgwXYoI8WacPpNZAk92Ofe4S7HOXMMvtho5Koze+5TZul5MeWzid4bPotsXQ43TS19tBd18nDkcP/Y4ecLlQegDc/SitAQUmhc1sw2INxWoJwWK1YQm1Y7HZsYSGYraHYbaH0dXpYvcnd3ANacxmMzs353J/VhjhKL6XkM31mg7yZ0cSGl7H9eZq3NqNy+3CpV243C6cQ04G3YM43U4GhwZxuBzG2yEHbu1+/GdVJsKt4URYI0gOT2aObQ6RtkiirFFEhUQRaYv0aRlqqiv2EuB14P/1QCxCzEierNMHbM3fZIK4HOPhGsDUfIuohhKiOsqNr8dkQE4+JC4ASwha6wcr4n7n8FtXP44hBw6Xg4GhAQaHBukdcuAccuJ0D+J09+J0NaC7NCh4Y30/Ne39pMXaaXAV0XD/YTiJydDogsaaT4eplMKiLMYvD5PxNsQcQpw9jlBzKCHmEOwW+4NHmDWMMEsYdovdJ6MCxssjNXal1DHgv413xS41diE+bTrX2Kekv8OYOtlYYmwWYrZAwnzjJqiYrEm1Tg65hxjSxkNrbbxFP1YaUShMyoRJmTArs/HwwgbhnuTTeezjSexKqXeAdwAyMjJWVFVVTfm8QoggoTV01RqlmqYb4Bow9m9NXgQpiyE8wd8RTomnftl6LLErpQ4BKaN86b9rrT8cfs4xZMUuhPCEIRe0lhlJvq0CtNvYyDtlMSTlgc23NwVNlScvaHusK0ZrvXVSEQghxGSYLcaYgqSFMNgLjTeM0cJlB6H8MMTnGkk+Ltd4boDzxwXtwH9VhBAzly0c0guNR0/T8IbdpdBSZmztl5RnjDKImjXuDbt9zR8XtKdUY1dKvQb8DEgEOoCrWuvnx/o+KcUIISbN7Yb2e8YF1+Y74HZBWPzDqZOh0f6O8LGaesDV2L1BErsQwiOcDmi+ZST5jmpj1R6TYZRqEuaDxfbUb/dGB5E3bxKTO0+FEMHPGgqzCoxHf7txwbWxBG7uBvN+SJw/PHUy67FSjbcScCDcJCaJXQgRHOyxkP0sZK2HzmqjFt90w0j2oVFGmSZ5MYQbNW5vJeBAuElMErsQwid8duPUSDkmJgPmbDUutDaWwP1zUHUWolIheTFrM1K9koBXZMby3ttr/HqTmNTYhRBeFxD7lg70DN/lWgw9zWAyU+FO4WxfKgvylrEiO/BvgpIauxAiYARC3ZmQCMhYbTy6G6GxmJzGUnKohdpScOYZF10jUwK2dXK8JLELIbwuEOrOnxKZbDxyNht3tzYWQ/01qC0yxhckj7RORvk3zkmSUowQwicmWmP3+TAzpwOabxo3QXXWGqv22CwjySfMG7N10hekj10IMW35vSbf12ZccG0oAUcnmK3GSOGUfIjJ9FupRmrsQohpy+81+bA4yN4AWc8arZMNJQ9X86FRxio+ZbHxvAAkiV0IEXACpib/aOvk3G3QcsdI8vfPQtUZY0ZNymJjYJnV7p8YRyGlGCFEQAroDUMGuo3WyYZi6G0Bkxni5wxPncwxPvYCqbELIYS3aQ09jcMbhJTCYJ8xLz5pkVGPj0j2aD1eauxCCOFtShl975EpkLvJaJ1sKIa6y1Bz0WidTFlstE6GRPosLEnsQgjhCSYzJMw1Hs7+h3Nq7h6FimNGy+Si13zSUSOJXQghPM1qh9krjEdfm7GK126ftUlKYhdCCG8Ki4Oc53x6SpNPzyaEEMLrJLELIUSQkcQuhBBBRhK7EEIEGUnsQggRZCSxCyFEkJHELoQQQUYSuxBCBBm/DAFTSjUDVT4/8dMlAC3+DmICplu8MP1inm7xwvSLebrFC/6NOVNrnTjWk/yS2AORUurSeKamBYrpFi9Mv5inW7ww/WKebvHC9IhZSjFCCBFkJLELIUSQkcT+0Lv+DmCCplu8MP1inm7xwvSLebrFC9MgZqmxCyFEkJEVuxBCBJkZm9iVUp9XSpUqpdxKqSde4VZKVSqlipVSV5VSftuodQLx7lBK3VZKlSulvufLGEeJJU4pdVApVTb8dtQdiZVSQ8Ov71Wl1Ed+iPOpr5lSKkQp9Zvhr59XSmX5OsbPxDNWvH+slGp+5DV92x9xPhLPL5RSTUqpkid8XSml/ufwz3NdKbXc1zGOEtNYMW9USnU+8hr/wNcxPpXWekY+gIXAfOAYsPIpz6sEEqZDvIAZuAvkADbgGpDnx5j/B/C94fe/B/zDE57X48cYx3zNgG8D/8/w+18EfhPg8f4x8H/7K8ZRYt4ALAdKnvD1F4G9gALWAOenQcwbgd3+jvNJjxm7Ytda39Ra3/Z3HOM1znhXAeVa6wqt9SDwH8Ar3o/uiV4Bfjn8/i+BV/0Yy5OM5zV79Of4HbBFKR/tcfa4QPtvPCat9Qmg7SlPeQX4N204B8QopVJ9E93oxhFzQJuxiX0CNHBAKVWklHrH38GMYTZQ/cjHNcOf85dkrXU9wPDbpCc8L1QpdUkpdU4p5evkP57X7MFztNYuoBOI90l0jxvvf+M3hssav1NKpfsmtEkLtH+347VWKXVNKbVXKbXI38E8Kqj3PFVKHQJSRvnSf9dafzjOw6zTWtcppZKAg0qpW8O/zT3OA/GOtor0atvT02KewGEyhl/jHOCIUqpYa33XMxGOaTyvmc9f16cYTywfA+9rrQeUUn+C8dfGZq9HNnmB9PqO12WM2/t7lFIvAruAuX6O6YGgTuxa660eOEbd8NsmpdQfMP4U9kpi90C8NcCjq7M0oG6Kx3yqp8WslGpUSqVqreuH/7RuesIxRl7jCqXUMWAZRh3ZF8bzmo08p0YpZQGi8d+f6WPGq7VufeTD/w/4Bx/ENRU+/3c7VVrrrkfe36OU+l9KqQStdUDMvZFSzFMopcKVUpEj7wPbgVGvkgeIi8BcpVS2UsqGcaHP510mj/gI+Prw+18HHvurQykVq5QKGX4/AVgH3PBZhON7zR79Of4IOKKHr6D5wZjxfqY+/Tngpg/jm4yPgK8Nd8esATpHSniBSimVMnKdRSm1CiOXtj79u3zI31dv/fUAXsNYKQwAjcD+4c/PAvYMv5+D0XVwDSjFKIkEbLzDH78I3MFY8fot3uFY4oHDQNnw27jhz68E/mn4/WeA4uHXuBj4ph/ifOw1A/4a+Nzw+6HAb4Fy4AKQ4+fXdax4fzz87/UacBRY4Od43wfqAefwv+FvAn8C/Mnw1xXw8+Gfp5indKkFUMx/+shrfA54xt8xP/qQO0+FECLISClGCCGCjCR2IYQIMpLYhRAiyEhiF0KIICOJXQghgowkdiGECDKS2IUQIshIYhdCiCDz/wNz3A6NPPSu8gAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# When I use np.linspace things go wrong?\n", + "x_new = np.linspace(x_1s[0].min(), x_1s[0].max(), 100)\n", + "\n", + "for order_3, x_1s, trace in traces_3:\n", + " α_p_post = trace['α'].mean()\n", + " β_p_post = trace['β'].mean(axis=0)\n", + "\n", + " x_new_order = np.vstack([x_new**i for i in range(1, order_3+1)])\n", + "\n", + " y_post = α_p_post + np.dot(β_p_post, x_new_order) \n", + " plt.plot(x_new, y_post, f'C{order_3}', label=f'model order {order_3}', alpha=.5)\n", + "\n", + "plt.scatter(x_1s[0], y_1s, c='C0', marker='.')\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd43NWZ8P3vmaLem9Usq7jITZZkueFugwFTbHpYSmgBNuyzT97svm/ItVdCNm92SZ4ned5kd9lNHGAhCQGHTgwY496bbLk39W7JapZkS5py3j9+spFtyWqj0Uhzf65Ll2Y0vzm/M4OZe37nPuc+SmuNEEII72Ma7g4IIYQYHhIAhBDCS0kAEEIILyUBQAghvJQEACGE8FISAIQQwktJABBCCC8lAUAIIbyUBAAhhPBSluHuwM1ERUXp5OTk4e6GEEKMGLm5uRe01tF9OdajA0BycjIHDx4c7m4IIcSIoZQq6euxMgQkhBBeSgKAEEJ4KQkAQgjhpTw6ByCEGB42m43y8nLa2tqGuyuiB35+fiQmJmK1WgfchgQAIcQNysvLCQ4OJjk5GaXUcHdHXEdrTV1dHeXl5aSkpAy4HRkCEkLcoK2tjcjISPnw91BKKSIjIwd9hSYBQAjRLfnw92yu+O8zagNAY1sje6v2IlteCiFE90ZtAChsKuTQ+UPsqdwjQUAIL5ecnMyFCxcGfUxfbd26lbvvvntQbTzzzDPExMQwbdo0l/SpO6M2AGTFZDEtcip5tXnkns8d7u4IIUYxh8MxqOfb7fYb/vbUU0+xfv36QbXbm0EHAKXUWKXUFqXUKaXUCaXU/+zmGKWU+jelVL5S6qhSKnuw5+1NR1Ex0w/UMTEwhf3V+zlWe2yoTymEcJHi4mLS09N57rnnmDZtGo899hgbN25k/vz5TJgwgf379wNQX1/P6tWrycjIYO7cuRw9ehSAuro6VqxYQVZWFi+88MI1owB/+tOfmD17NpmZmbzwwgu9fni/++67TJ8+nWnTpvGDH/zg6t+DgoL48Y9/zJw5c9izZw/r168nPT2dBQsW8NFHH109rrW1lWeeeYZZs2aRlZXFp59+CsBbb73FQw89xD333MOKFStuOO+iRYuIiIgY+JvYB66YBmoH/kFrfUgpFQzkKqW+1lqf7HLMncCEzp85wH91/h4yJyoaqdpzkvjqeFKWxLOjYgc+Zh8mRUwaytMKMfqc2wgt513bZtAYmHDrTQ/Jz8/n/fffZ82aNcyaNYs///nP7Ny5k88++4x//dd/5ZNPPuGVV14hKyuLTz75hM2bN/Pkk0+Sl5fHP//zP7NgwQJ+/OMf8/nnn7NmzRoATp06xdq1a9m1axdWq5Xvfve7vPPOOzz55JPd9qGyspIf/OAH5ObmEh4ezooVK/jkk09YvXo1ra2tTJs2jZ/+9Ke0tbUxYcIENm/ezPjx43nkkUeutvEv//IvLFu2jDfffJPGxkZmz57Nrbcar33Pnj0cPXp0yD/oezLoAKC1rgKqOm83K6VOAQlA1wCwCviDNsLwXqVUmFIqrvO5Lpdb0sBjG6qJupzIkk1HecAxGdvSaDaXbcbH7ENK6MDnzQoh3CMlJYXp06cDMHXqVJYvX45SiunTp1NcXAzAzp07+fDDDwFYtmwZdXV1NDU1sX379qvfwu+66y7Cw8MB2LRpE7m5ucyaNQuAy5cvExMT02MfDhw4wJIlS4iONoprPvbYY2zfvp3Vq1djNpt54IEHADh9+jQpKSlMmDABgMcff/xq0NmwYQOfffYZv/zlLwFjim1paSkAt91227B9+IOLF4IppZKBLGDfdQ8lAGVd7pd3/m1IAsDewjo67E7Kg6LZnJDJ3LIy5h3yYUdGCBuKN3B32t0kBCUMxamFGH16+aY+VHx9fa/eNplMV++bTKarY+bdTfC4Mj2yu2mSWmu+/e1v8+qrr/apDzebQOLn54fZbL7hvN218eGHHzJp0rWjD/v27SMwMLBP/RgqLksCK6WCgA+B72mtL17/cDdP6fadVUo9r5Q6qJQ6WFtbO6C+zE2NxMdiwqygPjSa+AfvQ11qY35eO2F2H74o/ILzrS6+pBVCuN2iRYt45513AGPmTVRUFCEhIdf8/csvv6ShoQGA5cuX88EHH1BTUwMYOYSSkp6rJ8+ZM4dt27Zx4cIFHA4H7777LosXL77huPT0dIqKiigoKACMvMEVt99+O//+7/9+NZgcPnzYBa/cNVwSAJRSVowP/3e01h91c0g5MLbL/USgsru2tNZrtNY5WuucK5dd/TVzXDjvPDeX76+YxDvPzSV79hRCV6/C3OFg0RE7QZdhXeE66tvqB9S+EMIz/OQnP+HgwYNkZGTw8ssv8/bbbwPwyiuvsH37drKzs9mwYQNJSUkATJkyhZ/97GesWLGCjIwMbrvtNqqqeh6IiIuL49VXX2Xp0qXMmDGD7OxsVq1adcNxfn5+rFmzhrvuuosFCxYwbty4q4/96Ec/wmazkZGRwbRp0/jRj37Up9f26KOPMm/ePM6cOUNiYiJvvPFGf96aPlGDnSOvjOuet4F6rfX3ejjmLuDvgJUYyd9/01rP7q3tnJwc7coNYWw1NTR9+int2saWaeAIDWT1+NWE+oa67BxCjAanTp1i8uTJw90N0Yvu/jsppXK11jl9eb4rrgDmA08Ay5RSeZ0/K5VSLyqlXuw85gugEMgHfg981wXn7TdrTAxh992Hn9mXxcccqIYm1hWuo9XWOhzdEUKIYeWKWUA76X6Mv+sxGnhpsOdyBUtUFKH33w8ff8Kio41snVrDOrWOVeNX4WfxG3C7uSUN7C2sY25qJDPHhbuwx0IIMTRG7Urgm7GEhxN2/30EB4Sz8Kid1spSPi/8HJvDNqD2cksaeOz1vfxqwxkee30vuSUNLu6xEEK4nlcGAABzWBhh999PeGgsC47ZaSrN54uiL7A5+x8Erkw7dWqw2Z3sLawbgh4LIYRreW0AADCHhBB6/31ERSZxy1Eb9QWnWF+0HrvzxrocN9N12qnVYmJuauQQ9VgIIVzH63cEMwcFEXbfasr++z2CvjrJqYb9KKW4I/kOLKa+vT1Xpp1KDkAIMZJ49RXAFYcvdPBEZRSbqgJp+aiC07t383XJ1zicfa/wN3NcOC8tHS8f/kJ4oJFWDrqsrIylS5cyefJkpk6dym9+8xuX9Ot6EgAwxvBbtIVNiTOpsiSRdLCVC4f28nXp1zi1c7i7J4TwcK4uB22xWPjVr37FqVOn2Lt3L6+99honT57s4dkDJwGAb8bwnRYLu1NymDjjFjLPOajbu5NNpZskCAjhZt5eDjouLo7sbKNqfnBwMJMnT6aiomIQ72j3vD4HADeO4c9IDKH56zDI285x22a2LFQsHbcMk5J4KbzPzoqdXLjsmqGRK6L8o1iQsOCmx0g5aENxcTGHDx9mzhzXV9CXANBp5rjwa8bvg1fcRrKPD859X3Nqywa2LTOxJGmpbJQthJtIOWhoaWnhgQce4Ne//jUhISH9fAd7JwGgB8pkImjpEtJ8fNA71nF2w+dsX6FYlLREgoDwKr19Ux8q3l4O2maz8cADD/DYY49x//3396m//SVjGjehlCJw/i1MWHYfE2ot1H7+KTtLt8km80J4iNFaDlprzbPPPsvkyZP5/ve/35e3YkAkAPRCKUXgnNlMuvNRUht9qPn0I/aU7JAgIIQHGK3loHft2sUf//hHNm/eTGZmJpmZmXzxxRf9fXt6Nehy0EPJ1eWgB+vyiROc+PQtSv0vE3//I8xJWiDDQWJUknLQI4MnlIP2Gv5TpzLtwe8wti2A6vff5WDxruHukhBCDNjoDQBtTVBzClx8heM3cSIZ3/pb4m2BVK79I7mFO/r0vNySBl7bki+VQoUQHmP0BoCKQ3DiEzj8J7jo2r3nfVNSmPHE/2SMCqXivT9yuODmQUDKRQshPNHoDQApi2HSnXC5HnLfglN/hbbr96ofON/ERLKf+j7RljDK332bI2e393islIsWQnii0RsATCaIz4Q5L0LSXGM4aP/voHgnDHDjl+v5jIkl++l/JMIvkrL33uboya3dHhce4INJKUxIuWghhOdwSQBQSr2plKpRSh3v4fElSqmmLnsG/9gV5+0Tiy+kLYXZ34GINCjaAft+B+dPuCQ/4Bsdw6xnfkBYcAxla98mL++rax7PLWngp+tO4NQak0nx47unSsVQIYRHcNUVwFvAHb0cs0Nrndn581MXnbfv/MNh2v2Q9Rj4BMDJz+DQH6Bp8AWWfCIimPPcDwmNSqDiw3fJ27/u6mNdh3+01jRc6hj0+YQQ/TPSykG3tbUxe/ZsZsyYwdSpU3nllVdc0q/ruSQAaK23A/WuaGvIhSXBzKch/S5ov2gEgZOfGrOGBsEaHMKcZ39IWEIKlZ+9z6EdHwCyW5gQ3sDV5aB9fX3ZvHkzR44cIS8vj/Xr17N3795BnaM77swBzFNKHVFKfamUmurG895IKYjLgNkvwLhboPYs7F8DRdvBPvBv6NaAQOY+/QNC09KpXv9Xcjf8ieykMN55bi7fXzGJd56bK8M/QvSBt5eDVkoRFBQEGDWBbDbbkCw6dVcxuEPAOK11i1JqJfAJMKG7A5VSzwPPA1eXbw8Ziw+kLjaSxYVboXgXVB2B1CUwZpoRKPrJ7OvH3Cf+kb1rf8P5bV9zqK2drLuflg9+MWK17NiBvda15aAt0VEELVx402O8vRy0w+Fg5syZ5Ofn89JLLw1JOWi3XAForS9qrVs6b38BWJVSUT0cu0ZrnaO1zrlSgnXI+YXClFWQ/QT4hsCpdcbU0cayATVntliZ9+j3Cc3K4fy+7Rx65zc47f3baF4Ib3elHLTJZLppOegnnngCuLEc9OOPPw70XA46MzOTTZs2UVhY2GMfupaDtlgsV8tBAz2Wg1ZKXT03GOWgf/7zn5OZmcmSJUv6XA7abDaTl5dHeXk5+/fv5/jxbufYDIpbrgCUUrHAea21VkrNxgg8njcZPjQRsp80ZggVbjUWkcWkQ+pS8A/rV1Mmk4l597/EHusaavbvYf0rz7Lih69h6bysE2Kk6O2b+lDx9nLQV4SFhbFkyRLWr1/PtGnT+tTvvnLVNNB3gT3AJKVUuVLqWaXUi0qpFzsPeRA4rpQ6Avwb8C3tqVXolILYaTDnBUhZCHX5sP/3RkCwt/erKZPJxC33vnD1/oZXX8Jx+bKLOyyE9xqt5aBra2tpbGwEjE1rNm7cSHp6eq/P6y+XXAForR/t5fH/AP7DFedyG7MVkhdAbAYUbYOSPVB1FFIWGX8z9S12KqW482dv8dkvXsTa3Mb+N35Ozrf/AWuw63f3EcLb/OQnP+Hpp58mIyODgICAa8pBP/roo2RnZ7N48eJuy0E7nU6sViuvvfbaNeWbu+paDlprzcqVK3stBx0VFcWCBQuuDtn86Ec/4nvf+x4ZGRlorUlOTmbdunU3tNFVVVUV3/72t3E4HDidTh5++OFBTSvtiZSD7quLlZC/0Vg3EBQD42+F8O7/0XRHa03u4S84/+kHhIfGkvPtf8Ansts0iBDDTspBjwxSDtpdQuIh6wkjWWxvh7w/w/EP4VLflj8opcjJvouER56krrWWA6//nMuV5UPcaSGE6JkEgP5QCsZMMcpKpC6G+iI48DoUbAZbW5+ayJyylJTHnqXOcZGD//1LWgsLhrjTQgjRPQkAA2G2GgvI5rwAY6ZC2X6j0FzFIXA6e336tNR5THjyb6n36SD3T/8fF08ec0OnhegfTx4eFq757yMBYDB8g42SEjOfgoBIOPsV5L5pXBn0YnJiFlOe/B9cCDNx+C//ScMB1y/zFmKg/Pz8qKurkyDgobTW1NXV4efnN6h2JAnsKlpD7Rko3AKXGyFqAqQtg4DuF3lcUVRXwKH3/5OIqlYylj1ExKJlss+wGHY2m43y8nLa2vo2tCncz8/Pj8TERKxW6zV/708SWAKAqznsUHEQSnYZtxNmQvJ8sPr3+JTSphL2ffxbIovqmDLrTsbccQ/K4q4qHUKI0aQ/AUA+ZVzNbDE2oBkzDYp3QPkB4wdg4feN/QmukxQ6DssD/4Md63/PsX2fY2u+SMLqhzH59xw0hBBisCQHMFR8g4wtKXOe+eZvO/4P1HU/6yc+OJ5b7/47Gm5J58SZHRS/8waOzpWAQggxFCQADLXgMbDk5W9yAUf/Yvy03lhdMdI/kjuXv0DL8hxOVeaR/8c12Kqr3dxhIYS3kByAOznsUJHbmR+wQXyWUW7CJ+Caw9rsbXx15EPUV9tJ9Y1jwuon8E1LG6ZOCyFGElkJ7KnMFkiaY6wfiM+EykOw77dQdgCc32xK4WfxY2XWw/jcdyfnVC0n/vJ7Lh3OG8aOCyFGIwkAw8EnECbeDjnPGiUm8jcaK4ov5F/dqN5qsrIi/V4i7nuAwtA2jn/+By5u24buw0IzIYToCwkAwykoGjIegekPGfePvQ9H10JLLQAmZWJR8lLGrX6UwkQrx7Z9QMMXn6NttmHstBBitJBpoMNNKYgaDxEpUHnYmDp68A2Iy4SUhSifQGbFzyHgniDyNq/Fnvsl6S0Xib73PkwBAb23L4QQPZArAE9hMkNijrFRfcJMY2/ifb+F0n3gsDM1ciq3rHiaklljOX52J1Xv/Ql75yYXQggxEBIAPI1PAEy4DWY9C6FjjUqjB16H2rOkhCSzYvHTVC6axPGqw1S8+za2ysrh7rEQYoSSAOCpAqMg42Hjx2Q29h7I+zOxTsVdc5/kwvJMjrecpfQvf6Dt7FmXnTa3pIHXtuSTWyJXF0KMdq7aE/hNpVSNUqrbbeuV4d+UUvlKqaNKqWxXnNcrRKYZs4UmroDWWsj9byKK97Jq2mpabp/HMSoo/vRdLh06NOjKjbklDTz2+l5+teEMj72+V4KAEKOcq64A3gLuuMnjdwITOn+eB/7LRef1DiaTkReY84KRJ6g+RtDhP7EqOAl9+0KOBTdR+PVHtGzdOqhponsL6+iwO3FqsNmd7C2sc+GLEEJ4GpcEAK31duBmeyOuAv6gDXuBMKVUnCvO7VWs/sZexLO/A2Hj8CvZzT0NNQTPncKxWBvn9qyn6fPP0R0dA2p+bmokPhYTZgVWi4m5qZEufgFCCE/irmmgCUBZl/vlnX+rctP5R5eACJj+INQXYSnYxIr6GnYmB3DCtwXb0a2Mb2kh7J57MQcF9qvZmePCeee5uewtrGNuaiQzx4UP0QsQQngCdwWA7nY46XbAWin1PMYwEUlJSUPZJ4+SW9LQ/w/eiBQIewZT9REWFm4jIKyRY8kt2Ap2MekvlwlftQpLZP++xc8cFy4f/EJ4CXcFgHJgbJf7iUC38xe11muANWAUgxv6rg2/K8nXDrsTH4uJd56b2/cPYZMJ4rNQ0ZPJKd1NgOkr9jiL6Di3mclrW4m89z58EhOH9gUIIUYkd00D/Qx4snM20FygSWstwz+dXJJ8tfpB2jKmLPgByyYvp2iyjaPlG6j503/Sdvq06zsthBjxXHIFoJR6F1gCRCmlyoFXACuA1vq3wBfASiAfuAQ87YrzjhZXkq82u3PwydeACFJyXuCu+Jl8EbgG28E8pr1VQfTtjxCw5E7Zb1gIcZXsB+AhBpQD6EX9pQusy1tD0NZ9ZNSaGJMxm6CHX0AF3nyjeiHEyCWbwourWjpa+PzcJzh3biXzdB1xMVGErH4E04RFYPEZ7u4JIVxMNoQRVwX5BLE6/SH8l9/N/oVTKb7YRsM7b2Pf+GuoPnZ1/wEhhPeRctBewNfsy92pd7PNJ5iDYblc3ldO6q5zhDU24DthkrG4LFRmCgnhbSQAeIkjZc2cKkgiKsrKscWHuHSwismFirCOs/g3VaLGTIHUJeAf1qf2hiJnIYRwLwkAXuD6dQb/+9HFFC7Yy6XcUjJqQ4iwBBGkzqAunIOxsyFp3k3zA4NatyCE8BiSA/AC168zKD0fzOr0B6mfN4l98Zepqmyk6XwCzuAUKNltbERTdaTH/IAUjRNidJAA4AW6K/IWHRDNg5MeQs2czt5JirKyszQcbcU+7h7wC4XTX0Duf0NjaZ/aE0KMPDIN1Ev0NGZvc9rYVLqJ8sKjTM9rIsU/kdAVK/ANaofCLdB2EaInQdpS8A/vtT0hxPCSdQCiX7TW7K/eT17RHiYeqmGSI5rQefMJyM5ElR+A0t3GcFBiDoybDxbf4e6yEKIH/QkAkgQWKKWYEzeHMN8wtvluovVIJdN3bSX8wgWCbr0VU1wGFG4zNqivPgYpiyB2hlGITggxYkkAEFdNiphEmG8Y632/ZPfpErJP5hLd2EjoypWYJ99t7EqWvxHOrIeKXGP9QHjycHdbCDFA8hVOXGNM4BgenPQQvjOms2u6heLqU9T/ZS0dpaUQEgdZj8PU+8DeAXnvwrEP4NLNNoMTQngqyQGIbtmddnZW7ORsySHSD9cxkRhC5i/CPyvTqCjqsEP5ASjZBdoJCdkwboFRlloIMWwkByAGzWKysDhxMVH+Uez030rz4Sqmb9tIWG0twcuWoqxWGDcPYqdD0XYoPwjVxyFlIcRlSX5AiBFAAoDokVKKaVHTiPCL4Cuf9ew9VUrW8f04GhoIuWsl5uBg8A2C9JXf5AfOboCKQzB+OUSkDvdLEELchHxNE72KD4o3Fo1lTmXnZEVJxQka1q7FVlHxzUHBYyDzb2DaA+B0wJG1cPR9aJVVwkJ4KgkAok+CfYJZPX41Y6fO4cCsUM60FlH30Ydczsvjah5JKYieCLOeg7Rl0FQKB16Hc1+D7fLVtnJLGnhtSz65JQ3D9GqEECBDQKIfrCYry5OWc8Q/in1+O2g+cp6pWzYSWl1N0LJlmHw6C8iZLZA0B2KnQdEOY8ro+eOQvJBcWzKPvXlACskJ4QFccgWglLpDKXVGKZWvlHq5m8efUkrVKqXyOn+ec8V5hfsppciMyWTl5NVUzk1hV0ILVccP0PiX97HXXzcd1CcQJt0BOc9AUCyc+5rWXb8jwVGBU2spJCfEMBt0AFBKmYHXgDuBKcCjSqkp3Ry6Vmud2fnz+mDPK4bX2OCxPDTpYXyyMtiVYaXw/Cnq166l7ezZGw8OioEZ34LpD5EcGcB9lt08YN7JGEuzFJITYhi54gpgNpCvtS7UWncA7wGrXNCu8HDBPsGsSlvFxKkLyL0limPOMi588VdaduxAOxzXHqwURI0n6fa/59Z7HuPBiWY+mHmKme37oaN1eF6AEF7OFTmABKCsy/1yYE43xz2glFoEnAX+L611WTfHiBHGbDKzIGEBcYFxbPHfROPRUjL2bSP6/HmC77gDc1DQtU8wmUmftRxmzIPinVB5GM6fgOQFxlRSk3l4XogQXsgVVwCqm79dv7z4r0Cy1joD2Ai83WNjSj2vlDqolDpYW1vrgu4JGPqZN2lhaTyY/jD22dPZNUlTVJRHw3vv0VFe3v0TfAJg4gqY9ayxH3H+Jtj/e6g9KxvVC+Emgy4FoZSaB/xEa3175/0fAmitX+3heDNQr7UO7a1tKQXhGu7cwtHmtLGzfCf5xYeYcKiWiaZYwuYvxD872ygh0ZO6AijYDK0XjJxB0lwYM3VI+ijEaNafUhCuuAI4AExQSqUopXyAbwGfXdehuC537wVOueC8oo/cuYWj1WRladJSFk67i3OLxnEw4DzVWzdw8fMvcLa39/zEyDRjttCEFdBSAyc/gy2vQn3RkPVVCG836ByA1tqulPo74CvADLyptT6hlPopcFBr/Rnw90qpewE7UA88Ndjzir67soWjze502xaO6RHpRPlH8ZXfenafPEfGyb0k1NcReuedWKKju3+SyQyJMznSEkz7nt+TGO5P/JH3jMfm/70xrVQI4TJSDdRLDNcWjh2ODraUbaE8P49JeXWk+Y8lfNmt+E2e3GM/rwxXzbbk839m1hEf6m88aDLDov/bmFEkhOiWu4eAxAgwc1w4Ly0d7/ZVtz5mH1aMW8GsGXdycuFYDlNK9Zef0bx5C9pmu+H4rsNVB+zj+Tj0SQiONR50OmDrz6G52q2vQYjRSkpBiCGnlCIjOoOYgBg2BHzF3iNnmH5gC3FVVYTceQeWiIirx3Y7XDXuaXDY4NDbRpI49y2jDHXKIvANHr4XJsQIJ0NAwq0u2y+zqXQTtWePMvFYA6kBSYQtWYbflClXZwnddLjK1mZsUl9+EJQJkubB2Nlgtg7DqxHC8/RnCEgCgHA7rTVHao9woHA7cbmlpLeFEzM1h6ClSzD5+vatkUv1ULjFWDfgFwKpSyFmsuQHhNeTACBGhAuXL/B18QaceSeYXOxgbHw6YXfciXXMmL430lBibETTUgOhCcZG9SHxQ9dpITycBAAxYticNvZU7uHc6T0kH6piou9YIhcu+2bv4b5wOqH6qLE1ZUerUYY6ZbFxZSCEl5EAIEacoqYituVvIOTAWdKbAokKiCb84UewjonpeyP2dijdA2UHjAIlY+caK4olPyC8iAQAMSK12lrZXLIJ/Zd1WFrbmRo1lbAJUwm9+67+NXS50cgP1Jw2ZgmlLjHKSkh+QHgBWQcgRqRAayB3p91D7T1zaEmO4sSFE5w49BW1//4f3a4Z6JF/GEy9D7IeN1YPn/orHPoDNPVQmE4ILyVXAMIjfVX8FaWlx0n46igAWWOyCL/9TvwmTepfQ1pD9TEo2gbtLTBminFF4NdrLUIhRiQZAhKjQkNbA++e+jPJH+wHIDE4gZScZQQtX973BPEV9g4o2wul+4z7Y2cbawgsPtccNlwlM4Rwlf4EAFkJLDxWuF843816iTVmK375lZBbRNOuj5nW2krI8lsxB/WjOJzFx1g5HDcDCrdCyW5j5lDKYmNVsVJuLZsthCeQACA83vMZz1OaXMq6lL8SXHAe+9ENpFYUMebWlfhNnNi/xvxCYcoqY/ex/E1w+nOoyIXxy9lb2HFD2WwJAGI0kyEgMaIUNRWx+9RXBO4+xrj2YJIzFxG6ZCkmf//+N6Y11JyEgi3Q3kyhSuRbm/yoswdglSsAMUJJDkCMam32NnaV76B6z1biztaROmYK8Xfci09y8sAadNigbB+U7qGi4RIHHeMZm7mc7NRYl/ZbCHeQACBOEZ7OAAAdf0lEQVS8QsnFEnYdXUfg7uOMtYeQOudWQhYuxuTj0/uTO12T9B1jNmYLVR839ixOWQyxGWCS2dJi5JAAILxGu6OdPWU7qdr+NTGFDaQlTCdh5X1YExJ6fW6PSd+LlUZ+oKnc2J94/HIITx76FyOEC8hCMOE1fM2+LElezoL7XuLC0ukcrT3G0T/8hqbt29B2+02f2+NeySHxxiKyqauN8hJ578KxD4wKpEKMIi4JAEqpO5RSZ5RS+Uqpl7t53Fcptbbz8X1KqWRXnFeIK8YGj+X+BS8Q9ugj5MdoDm58h+I/rsFWU9Pjc65sPmNW3LhXslJGeenZ34HUxdBQDAdeN64MbG1D/4KEcINBDwEppczAWeA2oBw4ADyqtT7Z5ZjvAhla6xeVUt8C7tNaP9Jb2zIEJAaioqWCPfs+wnfPEeLMEUxYdh8hs+aguhnL7/PCr/Zmo9po9TGw+HWuKciU/IDwOG7NASil5gE/0Vrf3nn/hwBa61e7HPNV5zF7lFIWoBqI1r2cXAKAGOjKXJvTxv7inVR+vY7wioukjs8h6Z6HsYQPclpnc7Wx/0BjGQRGGfmBiNTBtSmEC7l7JXACUNblfjkwp6djtNZ2pVQTEAlccMH5xSg1mJW5VpOV+alLqX58Mrt3vMfJffu58LsCJt3+MCGZM/tfSuKK4FjIfAwunIWCzXBkLUSOh7RlEBjZ+/OF8CCuuH7t7v+k67/Z9+UY40ClnldKHVRKHaytrR1058TI1WOSth9iA2O5Z8VLRDz2N5QFtXPg499S+O4bOJqbB94xpSB6Esz6DqQthaZSIz9wbiPYLg+8XSHczBUBoBwY2+V+IlDZ0zGdQ0ChQLdTKrTWa7TWOVrrnOjoaBd0T4xUN03S9oPVZOWW8cuZ/9QPaZ09hdNndnHwP/9fmo4fYVBDoGaLseHM7BcgLgMqDsK+3xob1jsdA29XCDdxRQ7AgpEEXg5UYCSB/0ZrfaLLMS8B07skge/XWj/cW9uSAxCurs5pd9rJzd9O5RcfE1R3meSM+aStfBhTYD8Ky/WkpcbIDzSUQECkkR+ITBt8u0L0g9sXgimlVgK/BszAm1rrf1FK/RQ4qLX+TCnlB/wRyML45v8trXVhb+1KABBDpablPPs3/hGVe5zw4GjS73yU8GlZA88NXKE11OUb00UvNxgJ4vHLjYSxEG4gK4HFqOTqqwGH08Hhs9uo+PJj/OpaiJuUzeR7Hsca5oICcE6HUWW0eKdRayg+C5IXGCUmhBhCEgDEqDOUtfob2xo5sPVdOvbsJ8ASwIRb7ydhztJu1w30W8clKN4BlXnG5vTJC4xS1Cbz4NsWohtSCkKMOq6YEdSTML8wbr39RdK/831aowM5uu4P7P2vn9JUVjD4xn0CYOLtkPOMUWIif5MxY+jCOWO4SIhhJAFAjAiumhHUE6UUaYnTWfHcz4i6ZxV1jZXs/d3PyPtoDR2tg5gyekVQNGQ8AhkPA8qoLXTkPSNxLMQwkSEgMWK4c7/expYL5H31DpfyDmP1DSD51lWkzb4NkyuGhZwOqDxsDA3Z241tKlMWgY8LZiIJryc5ACFcpLTkGKfW/QlHZTX+8YlMuecJxiSl93h8v4KU7TIU7zKSxWYLjJsPCTnGbSEGSAKAEC7kcDo4te9Lyjevw3m5jfDMHDLueIzAoGs/4AecqG6tM8pK1OWDf5hRViJqorHiWIh+kiSwEC5kNpmZNu9uFn3vF0Rmz6Eh7yDbfvMyR/d8ht3xzZ4DA05UB0ZCxkMw4xEwWeD4R5D3Z6PwnBBDSAKAEH0UEBjKnPv/llnP/5DAsGjK133Ixt/+E8UlRwAXJKojUiHnWZi4AlprIfctOP05tLe4/sUIgQwBCTEgTqeT4v2bKNz0CR1trfjOyCDr9scpqre6JlFta4OSzvyAMsG4WyBxtuQHRK8kByCEm9hamjn11Vqq8nbT4W8hasltZM2+Bz+Ln2tOcKneyA9cOGcEgvBkYyqp5AdEDyQACOFmF0sLOPnXP9FQWYgtIYqkW1cxbfwtWEwu+sbeUGzsTXzFzKcgJM41bYtRRQKAEMNAO51U7dtGwda/0nypHvvEcUxc/gCT4qYPvsgcGDWFtv/ym/ux0yBlMfiFDL5tMWpIABBiGDkvXaJk6xeUHNhEC+3orKlkLLyfpLBk1wQCezuU7IbyA8ZQUNI8GDvHqDUkvJ4EACE8gO3CBQo3fET56YO0BCgsc3PInLmShOAE15zgcgMUbIHaM+AbbOxOFjNF8gNeTgKAEB5Ca01bUREFGz6kquIMzRF++M7JIStjBfFB8a45SWOpsRFN83mj4Nz4WyHURUFGjDgSAITwMNrhoOXEMYq3rqO6tpiLYwLxnZ1D1pRlJAS54MNaa6g+BoVboaMVxkyB1CXgFzr4trvhzrpMon8kAAjhobTNRvORw5TsXE91XSmNCcH4zs5h5sQlrgkE9g4o3QNl+0Fh5AbGzgWLz+Db7jSUezOIwetPAJBVJUK4kbJaCcmZzdTpM0g6dJDS3V9T/elWdo49jP+sHGamLbxpIOj1m7fFB1IXGxVGC7caxeaqjhizhWKnuyQ/0F3JCwkAI9OgAoBSKgJYCyQDxcDDWuuGbo5zAMc675Zqre8dzHmFGOlMvr6EzpvP1BlZjD2wj/L9W6n6eDM7kg/jNzObnOT5JAQlXDNrqF/fvP3DYOpqSMwx8gOnPzdWFY+/FcLGDqrvV0pe2OzOIdmbQbjPoIaAlFL/C6jXWv9cKfUyEK61/kE3x7VorYP6274MAQlv4WhupnnfHsoP7aCyrYYLaZEEZGaSk3TL1UDw2pZ8frXhDE4NZgXfXzGJl5aO771xreH8CeOKoL0ZYtIhdakRJAZIcgCey205AKXUGWCJ1rpKKRUHbNVaT+rmOAkAQvSBvaGB5j27qTy2lwrbBWonRBKQkUlO4hzO1wXy+Bv7rn7z7vfYu70DyvZB2V7QwNhZxhoCi++QvR7hfu4MAI1a67Au9xu01jf8i1RK2YE8wA78XGv9SV/alwAgvJWtpobmPbupOp1LubOOmonRWNIn0tERTXNDCvPSogb+zbvtonE1cP6EsQtZyiKIzQBX7HYmhp1LA4BSaiMQ281D/wS83ccAEK+1rlRKpQKbgeVa62533FZKPQ88D5CUlDSzpKSkL69DiFHJVlHBxd27qS48RnFrGe0RgVQvSueOtJWkhKQMbmVxUwUUbDJ+B8UY+YHwca7rvBgWHjcEdN1z3gLWaa0/6K19uQIQwlhMZisp4cKnH5NXkwdA7Zw0rGlpZMZmMTF84sCLzmkNNaegcItxZRA90cgPBES48BUId3JnAPjfQF2XJHCE1vr/ue6YcOCS1rpdKRUF7AFWaa1P9ta+BAAhvqHtdi787ndop5P6tnoqTE1UpIbiSE1kasx0pkZNJdA6wI3lHTZj7UDpHtBOY/ZQ0i1g7b6stSSBPZc7A0Ak8BcgCSgFHtJa1yulcoAXtdbPKaVuAX4HODF2IPu11vqNvrQvAUCIG2mt6SgspPXAAeorCqjmIiWJvrSkRpMak870qOmMCRgzsOGh9mYo2m6sKrb6Q/JCiMu8Jj8gC8E8m6wEFsILaK3pKC7m8uE8mksLqOmooyjWRH1aFBHRY5keNZ20sLSBDQ9drDLyA41lEBQNacshIgVg4NNRhVvISmAhvIBSCt+UFHxTUgisqSE0L4+EM6ep21VHSUQdO1IL2RMTw5TIKf0fHgqJg8zHjEqjhVvgyHsQNQHSlg1qIZgMHXkWuQIQYhRxNDdz+ehR2k6coPFiLZVBHRSN9eFyfARpEeOZHjWd2MDuJvXdrFG7sfdA6W5wOiAhm0NqKntKW/v1QS5DR+4hVwBCeClzcDBB8+cTMGsWgSdPEnHkCCnnajlfWEN+wgU+TjpNdEhc/4aHzBYYN8+oJVS8A8oPkq0PkB0cC2O/3ee+SQ0hzyMBQIhRyOTjg39mJn4ZGQQWFhJ0+DCJ5ZVcKG2kMO4iW5PK2RMcxuTIyaRHpBPq24ey0b5BMOlOiM+Gg29CczUcfAPSlkFkWq9Pd0cNIRli6h8ZAhJiBOvPB56tqorLeXm0FxRw0dZMabSiMMFCe3gAcYFxTI6cTFpoGta+bC2pNVw4BwWbjZ3JItOMRHHgzT/Uh/IDWoaYDDIEJIQX6O8HnjUuDmtcHI6mJvyPHCH05CkmHmnlQmALhWPy2Rpbxg4/f1LDUpkcMZm4wLiep5IqZSwai0g1qoyW7IQDr0NCNoybDz4B3T5t5rjwIftQliGm/pMAIMQINdAPPHNoKEGLFhEwezbtZ88ScPIkcSW1tBZepirGTkFMHmciTxPiG8qkiEmkR6QT7BPcQ2MWSJoDsdOgeKcRDM4fN9YPxGeByeziV90zKVPdfzIEJMQIdeUKYMDVQbuw1dTQdvIk7WfPYW+7RKPVRmmcmaIYhTPAl4TgBCaFTyI1LBWr6SZDRC21xvqB+iIIiPwmP+CmjeolByALwYTwOEP1weTqdrXNRnthIW0nT2ErL6fd2UFthJn8MU5qo6z4WP0YHzae9Ij0nlcbaw11BUYguFRvLCBLW24sKBNDTgKAEB5kpCYnHY2NtJ0+TdvJUzhaW2gx26iM9eFcjJ3LQcYew5kxmYwPG0+0f/SNwcDpgIpDxtRRR4cxJJS8wChBLYaMJIGF8CAjNTlpDgsjcO5cAmbPpqOkBL9TpwguKmJ8uaKg8RR1Ub4cn2kjryaPYJ9g0kLTSA1L/ebKwGQ2Np0ZMxVKdhnB4PwJI0mcmOPW/IDongQAIYbYSE9OKpPpaskJZ2srbWfOMnGXGRxg3+egMUxTEdXMifAD5NXmEWgNJC0sjbTQNGIDY1E+ATDhNuMKoGCz8VN52MgPRE1wW35A3EiGgIRwg9GWnNRa46ivpz2/gPaCfBx19di1naZQCxWRiqIIOx0BVgKtgaSEppAWlkZcYBwmZerMD2yG1gvGBjRpyyF4zHC/pBuM1P9mkgMQQriVvaGBjoIC2vMLsNfW4tAOGoNNVEWZKQzvoD3IB3+LP6mhqaSEphAXMAZr9XEjP2BvM7akTFlkrDZ2s+4+6Edq3gYkByDEqOTqb6SubM8SHo4lJ4eAnBwcTU20FxTgW1BAZPF5JhcpmgLtVEVfojDiMCeCTmAxWUgISmDc+IUkNVYSUn0Mak8Zm9AkzjLWF7hBTx/0IzVv018SAIQYAVz9jXQov+GaQ0MJyM4mIDsbR3MzHQUF+BUUElFaRXqxotXPRm2Ek7KQc+wIL0RbzYRHRDGu+QJJZ9cRV5GLefxyiE4f8vxATx/0Iz1v01cSAIQYAVz9jdRd33DNwcH4Z2bin5mJo6WVjqIi/EtKCC0vJ7VK0+7soCHERHVEB6dDzOQF+OPTdIy4/cdICEshbsJdRMdMNXIHQ6CnD/qZ48J557m5IzIH0B8SAIQYAVz9jXQ4vuGagwLxnz4N/+nT0A4H9upqOkpLCSotI660hgytuKgc1IdPpNpUzcGGPDpqjuATmkBs4lwSwtOID4wnOiDaZQHhZh/0Q1m3yFMMdk/gh4CfAJOB2VrrbjO2Sqk7gN8AZuB1rfXP+9K+JIGF+IYn5wAGy3npEh1l5djKSukoLcPZ2kqH7RItFwtppJaaEKiLi6V9XBqmoGBiA2NJCEogPiieaP9ozLKm4Cp3bgo/GWOz998B/9hdAFBKmYGzwG1AOXAAeFRrfbK39iUACOF9tNY4W1qwVVVhP38eW2kB9rMHsTWfp1nBxeg4aqNDqQ9RtEUG0REWcHXT+nEh41ietBw/i98wv4rh47ZZQFrrU50nvNlhs4F8rXVh57HvAauAXgOAEML7KKUwBwdjDg6GiROBhWj749jz8wg79CX2ynJSKxtor4uhubiNgsYTXZ69jz8vPoZfTCxRYfHEBMQQExBDtH903/Y58DLuyAEkAGVd7pcDc9xwXiHEKKEsFqzpOVgnZsP5Y1C4DUdTA5HmeCJOBoLFl4sdF6lsqSD2YAPtjmqazQcoDDBxJsQPW0gAAVGxhMUmER0xlujAGKL8o/q2JeYo1uurV0ptBLrbRfqftNaf9uEc3V0e9DjupJR6HngeICkpqQ/NCyG8hskEcTMgOh1z6R7MZQeIzrQDHURmPkN8STWm4GAcTRdxNNRzqbaaizUVXDrfQOvJU7TaDlJi0eQH+2ML8ScgKpbQMWMJj0okLCKeiJAxBFi738zGHdydl+k1AGitbx3kOcqBsV3uJwKVNznfGmANGDmAQZ5bCDEaWXwhdQnEZcLe/wLAlPff+ANk/OBqTiAYiNEaZ+slHI0N2OvqaKmtpOl8Gc01FVyqLKHVdpJmp4NSQFtMqKAg/ELCCAiNIig0iqCQKIJDown0D8VkNoPZjDKZjHOYzCiz6Zu/dXnsyn1Mpt6GyQHjw/+J3++mw+7EYrW4ZfWxO65/DgATlFIpQAXwLeBv3HBeIcRo5x8GS38Ixz4w9igG2PYLyH4CQhOBzpxCUCDmoEB8EhMJYAYxnU93trVhr6+npaGGpoYqLjacp7WxlstN9dQXnqLuchvKqa+242Oy4mP2wWrywcfc+WPywWq+8ncrqrtBD7MJZTJ3/u4SOJQJbbejbTaqz1Zzf34t+8dMpiAiyS2rjwcVAJRS9wH/DkQDnyul8rTWtyul4jGme67UWtuVUn8HfIUxDfRNrfWJmzQrhBD9M/1BY/+BvHfgYiUc+iPETDauEvzDenyayc8Pn/h4IuLjiSDzmse01rR0tNDYXMvFphqaLzdyqb2FpvYWLne0cKm9FafDjnJqlLaBsx2TEwJMvvib/PA3+RJg8sNP+RBgNn77KR98lBWlNdrhAKdGWS0oq5X4mGROtpyhwS/UbWszpBicEGJ0sXdA2V4o3WfcHzsbkuaBxcelp9Fa0+Zoo6WjhRZbC622VuN3h/H7yt/sTvs1z1MoAqwBBFmDCPQJxM/sh4/ZBz+zH4U17RSf92HphNQBf/uXYnBCCO9l8TEqi8bNgMKtULIbqo5A6mKj6qiL6gsppfC3+ONv8Sea7re71FrT7mg3AkLHN0HiSnBoaGug3d5Om6MNp3YCcEfmIqZFuWdhngQAIcTo5BcKU1ZBQg7kb4TTX0BFLoy/FcLcM8NQKYWfxQ8/ix9R/lE9Hqe1xq7tdDg63Do1dWgqLAkhhKcITYDsJ2HKvWC7DIffgeMfweWGYe1WbkkDr23JJ7ekAaUUVpOxgY6v2ddtfZArACHE6KeUsTdx1EQo2wele6Au39ibeNx8Y1rpTQxFHSZP2HBGAoAQwnuYrZC8wMgFFG03EsXVx4ycQeyMq+sHuhqKD2tP2XBGhoCEEN7HLwQm3w0znwL/CDizHnLfhIbiGw7t7sN6sK6U4zYrhnXDGbkCEEJ4r5A4yHocak9DwRbIexeiJkDaMgiIAIZm7wRP2XBG1gEIITzKsO1T4LBD+QEo2QXaCQkzjfyA1c+j9k7ojawDEEKMSMOaHDVbYNw8iJ1u5AfKD3TmBxYyc2yWx3/wD4TkAIQQHmMoxtv7zTcI0lfCzKchMBrOboCDb0B9ofv7MsQkAAghPIanJEcBCB4DmX8D0x4w6gwdWQtH34fWYQhKQ0RyAEIIj+KR4+0Ou7GKuGSncTsh25hOavUf7p7dQHIAQogRa+a4cM/54L/CbIGkORA7DYp2GMHg/HFIXgjxWTBCN6WXISAhhOgrn0CYdAfkPANBsXDuazjwBlzIBw8eTemJBAAhxIjWtaaO2wTFwIxvwfSHAA3H3oeja6Gl1n19cAEZAhJCjFjDOm1UKYgaDxEpUHEIincYs4Xis4yhIZ/h21u4r+QKQAgxYnnEtFGTGcbOgjkvQnw2VObBvt9C2X5j9pAHkwAghBixPGraqE8ATFwBs56FkATI3wT7fw+1Zz02PzCoaaBKqYeAnwCTgdla627nbCqlioFmwAHY+zpFSaaBCiF645HTRgHqCqBgM7RegPBxxkY0QTG9P2+Q3DkN9DhwP/C7Phy7VGt9YZDnE0KIa3jktFGAyDQITzaGhIq3w8E3jW0qkxcaq409wKACgNb6FBjbngkhhLiOyQyJM2HMFCjeZawfqDlpFJlLyDHWFwxn99x0Hg1sUErlKqWed9M5hRDCM1j9YcKtMOs5CBtnlJ4+8HuoOT2s+YFew49SaiMQ281D/6S1/rSP55mvta5USsUAXyulTmutt/dwvueB5wGSktyzcbMQQrhFYCRMf9AoLJe/CU58DGFjjfxAcHcfs0PLJbWAlFJbgX/sKQl83bE/AVq01r/s7VhJAgshRi2nE6ryjNLT9jajDHXKIvANHlSz/UkCD/kQkFIqUCkVfOU2sAIjeSyEEN7LZDKKys15ERJnwfkTsO93UO2+j8dBBQCl1H1KqXJgHvC5Uuqrzr/HK6W+6DxsDLBTKXUE2A98rrVeP5jzCiHEqGH1g/HLjfxARMrVrSjdQcpBCyHEKOJRQ0BCCCE8kwQAIYTwUhIAhBDCS0kAEEIILyUBQAghvJQEACGE8FISAIQQwktJABBCCC/l0QvBlFK1QMlw96OLKGCk7Wkw0vo80voLI6/PI62/MPL6PJz9Hae1ju7LgR4dADyNUupgX1fYeYqR1ueR1l8YeX0eaf2FkdfnkdJfGQISQggvJQFACCG8lASA/lkz3B0YgJHW55HWXxh5fR5p/YWR1+cR0V/JAQghhJeSKwAhhPBSEgBuQin1kFLqhFLKqZTqMaOvlCpWSh1TSuUppYZ1A4N+9PkOpdQZpVS+Uupld/bxun5EKKW+Vkqd6/wd3sNxjs73N08p9Zm7+9nZh5u+Z0opX6XU2s7H9ymlkt3fy2v601t/n1JK1XZ5X58bjn526c+bSqkapVS3W2Ipw791vp6jSqlsd/fxuv701t8lSqmmLu/vj93dx15preWnhx9gMjAJ2Ark3OS4YiBquPvb1z4DZqAASAV8gCPAlGHq7/8CXu68/TLwix6Oaxnm97XX9wz4LvDbztvfAtZ6eH+fAv5jON/X6/qzCMgGjvfw+ErgS0ABc4F9Ht7fJcC64X5fb/YjVwA3obU+pbU+M9z96I8+9nk2kK+1LtRadwDvAauGvnfdWgW83Xn7bWD1MPWjN315z7q+lg+A5Uop5cY+duVJ/437RGu9Hai/ySGrgD9ow14gTCkV557e3agP/fV4EgBcQwMblFK5Sqnnh7szfZAAlHW5X975t+EwRmtdBdD5O6aH4/yUUgeVUnuVUsMRJPrynl09RmttB5qASLf07kZ9/W/8QOdwygdKqbHu6dqAedK/276ap5Q6opT6Uik1dbg7cz3LcHdguCmlNgKx3Tz0T1rrT/vYzHytdaVSKgb4Wil1uvPbwZBwQZ+7+1Y6ZNPBbtbffjST1PkepwKblVLHtNYFrulhn/TlPXPr+9qLvvTlr8C7Wut2pdSLGFcvy4a8ZwPnSe9vXxzCKMvQopRaCXwCTBjmPl3D6wOA1vpWF7RR2fm7Rin1Mcbl95AFABf0uRzo+m0vEagcZJs9ull/lVLnlVJxWuuqzsv5mh7auPIeFyqltgJZGGPc7tKX9+zKMeVKKQsQyvANEfTaX611XZe7vwd+4YZ+DYZb/90Oltb6YpfbXyil/lMpFaW19piaRjIENEhKqUClVPCV28AKoNtZAR7kADBBKZWilPLBSFgOy8yazvN+u/P2t4EbrmCUUuFKKd/O21HAfOCk23po6Mt71vW1PAhs1p3ZwGHQa3+vGz+/Fzjlxv4NxGfAk52zgeYCTVeGDz2RUir2Sg5IKTUb4/O27ubPcrPhzkJ78g9wH8a3jnbgPPBV59/jgS86b6dizLA4ApzAGIbx6D533l8JnMX4Fj1sfcYYI98EnOv8HdH59xzg9c7btwDHOt/jY8Czw9TXG94z4KfAvZ23/YD3gXxgP5A6zP8Weuvvq53/Zo8AW4D0Ye7vu0AVYOv8N/ws8CLwYufjCnit8/Uc4yYz8zykv3/X5f3dC9wynP3t7kdWAgshhJeSISAhhPBSEgCEEMJLSQAQQggvJQFACCG8lAQAIYTwUhIAhBDCS0kAEEIILyUBQAghvNT/D9ENARfp4G/VAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# When I keep the original x_1S things look right?\n", + "for order_3, x_1s, trace in traces_3:\n", + " α_p_post = trace['α'].mean()\n", + " β_p_post = trace['β'].mean(axis=0)\n", + "\n", + " y_post = α_p_post + np.dot(β_p_post, x_1s) \n", + " plt.plot(x_1s[0], y_post, f'C{order_3}', label=f'model order {order_3}', alpha=.5)\n", + "\n", + "plt.scatter(x_1s[0], y_1s, c='C0', marker='.')\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
waicpwaicdwaicweightsedsewarning
Model order 29.046742.59952014.8008300
Model order 310.50583.108771.459081.83187e-154.727320.9593150
Model order 128.49562.3561719.44898.88178e-165.325585.118090
\n", + "
" + ], + "text/plain": [ + " waic pwaic dwaic weight se dse \\\n", + "Model order 2 9.04674 2.59952 0 1 4.80083 0 \n", + "Model order 3 10.5058 3.10877 1.45908 1.83187e-15 4.72732 0.959315 \n", + "Model order 1 28.4956 2.35617 19.4489 8.88178e-16 5.32558 5.11809 \n", + "\n", + " warning \n", + "Model order 2 0 \n", + "Model order 3 0 \n", + "Model order 1 0 " + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "traces_d = {f\"Model order {order}\":trace for order, trace in traces_3}\n", + "az.compare(traces_d)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercise 4" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First fit the linear and quadratic model again" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'y')" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFz9JREFUeJzt3X+s3fV93/HnG2PIzRrlktgN+IKL6RjLMqtzekeTWppYEkqKKuy4qUb+WKBKZGUrqvaLySxTGkXb7BRpUrpkTd0kKkgdkDHiuIPMS+JGWVeRcYkhDjAXBxVxr1lwICaL4qU2vPfH+V44vpz743vvOefz/Z7zfEhH95zv+fict7+2vq/7/X5+fCMzkSRppc4rXYAkqV0MDklSLQaHJKkWg0OSVIvBIUmqxeCQJNVicEiSajE4JEm1GBySpFrOL13AIGzYsCEvv/zy0mVIUms8/PDDP8jMjStpO5LBcfnllzMzM1O6DElqjYh4eqVtvVQlSarF4JAk1WJwSJJqMTgkSbUYHJKkWgwOSVItxYbjRsRlwJ3AxcDLwP7M/NSCNgF8Crge+Alwc2Z+e1A1HTgyx+2HjnHi1Gk2TU5w63VXsXPb1KC+TpJaqeQ8jrPAP8/Mb0fEG4CHI+Krmfl4V5tfBa6sHr8E/H71s+8OHJnjtvuOcvrMSwDMnTrNbfcdBTA8JKlLsUtVmfns/NlDZv5f4Alg4RF6B3BndjwITEbEJYOo5/ZDx14JjXmnz7zE7YeODeLrJKm1GtHHERGXA9uAby14awp4puv1LK8Nl/nP2B0RMxExc/Lkydo1nDh1utZ2SRpXxYMjIn4G+C/AP8nMHy18u8cfyV6fk5n7M3M6M6c3blzRcivn2DQ5UWu7JI2rosEREevphMYfZ+Z9PZrMApd1vb4UODGIWm697iom1q87Z9vE+nXcet1Vg/g6SWqtYsFRjZj6PPBEZv77RZodBD4YHe8AXszMZwdRz85tU+zdtZWpyQkCmJqcYO+urXaMS9ICJUdVbQf+IXA0Ih6ptv0rYDNAZn4WeIDOUNzjdIbj/uYgC9q5bWpoQeHQX0ltVSw4MvPP6N2H0d0mgd8aTkXD49BfSW1WvHN8HDn0V1KbGRwFOPRXUpsZHAU49FdSmxkcq3TgyBzb9x1my5772b7vMAeOzK34zzr0V1KbjeQ9xwdtrZ3b820cVSWpjQyOVViqc3ulB/9hDv2VpH7yUtUq2LktaZwZHKtg57akcWZwrIKd25LGmX0cq2DntqRxZnCskp3bksaVl6okSbUYHJKkWgwOSVIt9nE0mPfskNREBkdDec8OSU3lpaqG8p4dkprK4GgolzWR1FQGRwMdODLHedH7rrouayKptKLBERFfiIjnIuK7i7x/TUS8GBGPVI+PDbvGYZvv23gp8zXvuayJpCYo3Tn+R8CngTuXaPM/MvPXhlNOeb36NgDWRbB311Y7xiUVV/SMIzO/CbxQsoamWawP4+VMQ0NSI7Shj+OdEfFoRHwlIt5WuphBc8l2SU3X9OD4NvBzmfkLwH8ADizWMCJ2R8RMRMycPHlyaAX2m0u2S2q6RgdHZv4oM39cPX8AWB8RGxZpuz8zpzNzeuPGjUOts592bpti766tTE1OEMDU5IR9G5IapXTn+JIi4mLg+5mZEXE1naB7vnBZA+eS7ZKarGhwRMRdwDXAhoiYBX4HWA+QmZ8F3g/8o4g4C5wGbszsMU5VkjQ0RYMjMz+wzPufpjNcV5LUEI3u45AkNY/BIUmqxeCQJNVicEiSajE4JEm1GBySpFoMDklSLQaHJKkWg0OSVIvBIUmqxeCQJNVicEiSajE4JEm1GBySpFoMDklSLQaHJKkWg0OSVIvBIUmqpeitY9V+B47McfuhY5w4dZpNkxPcet1V7Nw2VbosSQNkcGjVDhyZ47b7jnL6zEsAzJ06zW33HQUwPKQRVvRSVUR8ISKei4jvLvJ+RMTvRcTxiPhORLx92DVqcbcfOvZKaMw7feYlbj90rFBFkoahdB/HHwHvXeL9XwWurB67gd8fQk1aoROnTtfaLmk0FA2OzPwm8MISTXYAd2bHg8BkRFwynOq0nE2TE7W2SxoNpc84ljMFPNP1erbapga49bqrmFi/7pxtE+vXcet1VxWqSNIwNL1zPHpsy54NI3bTuZzF5s2bB1lTY5Qe0TT/XY6qksZL04NjFris6/WlwIleDTNzP7AfYHp6ume4jJKmjGjauW3KoJDGTNMvVR0EPliNrnoH8GJmPlu6qCZwRJOkUoqecUTEXcA1wIaImAV+B1gPkJmfBR4ArgeOAz8BfrNMpc3jiCZJpRQNjsz8wDLvJ/BbQyqnVTZNTjDXIyQc0SRp0Jp+qUqLcESTpFKa3jmuRTiiSVIpBkeLOaJJUgkGhxbV73kipeedSOoPg0M99XueSFPmnUhaOzvH1VO/54k470QaHQaHeur3PBHnnUijw+BQT/1e+daVdKXRYXCop37PE3HeiTQ67BxXT/2eJ+K8E2l0RGdVj9EyPT2dMzMzpcuQpNaIiIczc3olbb1UJUmqxeCQJNVicEiSajE4JEm1GBySpFoMDklSLc7jUCu4sq7UHAaHGs+VdaVm8VKVGs+VdaVmKRocEfHeiDgWEccjYk+P92+OiJMR8Uj1+HCJOlWWK+tKzVLsUlVErAM+A1wLzAIPRcTBzHx8QdN7MvOWoReoxtg0OcFcj5BwZV2pjJJnHFcDxzPzqcz8K+BuYEfBetRQrqwrNUvJ4JgCnul6PVttW+jXI+I7EXFvRFw2nNLUJDu3TbF311amJicIYGpygr27ttoxLhVSclRV9Ni2cKnePwHuysyfRsRHgDuAd/X8sIjdwG6AzZs397POsdWkIbA7t00ZFFJDlDzjmAW6zyAuBU50N8jM5zPzp9XLPwR+cbEPy8z9mTmdmdMbN27se7HjZn4I7Nyp0ySvDoE9cGSudGmSCisZHA8BV0bEloi4ALgRONjdICIu6Xp5A/DEEOsba20bAnvgyBzb9x1my5772b7vsAEnDVCxS1WZeTYibgEOAeuAL2TmYxHxCWAmMw8Cvx0RNwBngReAm0vVO27aNATWCYLScBWdOZ6ZDwAPLNj2sa7ntwG3DbsutWsI7FJnRwaH1H/OHFdPbRoC26azI2kUGBzqqU1DYBc7C2ri2ZE0ClzkUItqyxDYW6+76pw+Dmju2ZE0CgwOtd58uDVlzok06pYNjmrk0x9n5g+HUI+0Km05O5JGwUr6OC6mswDhF6vVbHvN+JYkjYllgyMz/zVwJfB5OvMonoyIfxcRPz/g2iRJDbSiUVWZmcD/qR5ngYuAeyPidwdYmySpgVbSx/HbwE3AD4DPAbdm5pmIOA94EviXgy1RktQkKxlVtQHYlZlPd2/MzJcj4tcGU5YkqamWDY7uJUB6vOeigxp7TVp+XhoG53FIa+ACixpHLjkirUHblp+X+sHgkNbABRY1jgwOaQ1cYFHjyD4OqabuzvA3Tqxn/brgzEv5yvsusKhRZ3BINSzsDD91+gzrzwsuev16Tv3kjKOqNBYMDqmGXp3hZ15OXn/B+Rz52K8UqkoaLvs4pBrsDJcMDqkWO8OlwsFRLdN+LCKOR8SeHu9fGBH3VO9/KyIuH36V0qvadC92aVCK9XFExDrgM8C1wCyde34czMzHu5p9CPhhZv71iLgR+CTwD4ZfrdpkkEuAeLdBqWzn+NXA8cx8CiAi7gZ2AN3BsQP4ePX8XuDTERHVMu/SawxjCRDvNqhxV/JS1RTwTNfr2WpbzzaZeRZ4EXjzUKpTK7kEiDR4JYOj1y1oF55JrKRNp2HE7oiYiYiZkydPrrk4tZOjnqTBKxkcs8BlXa8vBU4s1iYizgfeCLzQ68Myc39mTmfm9MaNGwdQrtrAUU/S4JUMjoeAKyNiS0RcANwIHFzQ5iCduw8CvB84bP+GluKoJ2nwigVH1WdxC3AIeAL4YmY+FhGfiIgbqmafB94cEceBfwa8Zsiu1G3ntin27trK1OQEAVz0+vVceP55/NN7HmH7vsMcODJXukSp9WIUf4Gfnp7OmZmZ0mWosIUjrKBz9rF311ZHRUkLRMTDmTm9krbOHNfIcoSVNBgGh0aWI6ykwTA4NLIcYSUNhsGhkdWWEVYHjsyxfd9htuy53w58tYL349DIasO6UsNYIkXqN4NDI63p60ot1YHf5Lo13rxUJRVkB77ayOCQCrIDX21kcEgFtaUDX+pmH4dUUBs68KWFDA6psKZ34EsLealKklSLwSFJqsXgkCTVYnBIkmoxOCRJtTiqStJQHDgy57DjEWFwSCOuCQdsF3McLV6qkkbY/AF77tRpklcP2MNeut27MY4Wg0MaYU05YLuY42gpEhwR8aaI+GpEPFn9vGiRdi9FxCPV4+Cw65TarikHbBdzHC2lzjj2AF/PzCuBr1evezmdmX+netwwvPKk0dCUA7aLOY6WUsGxA7ijen4HsLNQHdJIa8oBe+e2Kfbu2srU5AQBTE1OsHfXVjvGWyoyc/hfGnEqMye7Xv8wM19zuSoizgKPAGeBfZl5YInP3A3sBti8efMvPv300/0vXGqhJoyqUvNFxMOZOb2itoMKjoj4GnBxj7c+CtyxwuDYlJknIuIK4DDw7sz83nLfPT09nTMzM2uoXpLGS53gGNg8jsx8z2LvRcT3I+KSzHw2Ii4BnlvkM05UP5+KiG8A24Blg0NS83jmMzpK9XEcBG6qnt8EfHlhg4i4KCIurJ5vALYDjw+tQkl905T5JOqPUsGxD7g2Ip4Erq1eExHTEfG5qs1bgZmIeBT4Uzp9HAaH1EJNmU+i/iiy5EhmPg+8u8f2GeDD1fM/B7YOuTRJA7CW+SRe4moeZ45LGrjVzifxElczGRySBm6180m8xNVMro4raeDmLy3VveTUlCVTdC6DQ9JQ7Nw2VbtvYtPkBHM9QqKfS6bYh1Kfl6qkMXPgyBzb9x1my5772b7vcKP7Cwa9ZIp9KKtjcEhjpG0HykGvcWUfyup4qUoaI0sdKJt6eWY1l7hWyj6U1fGMQxojHijP1ZRl59vG4JAaahB9ER4oz9WUZefbxuCQGmhQfREeKM/lfUJWxz4OqYEG1Rex2vkUo2yQfSijyuCQGmiQfREeKLVWBofUQMOY+Kb+G5fJhAaH1EC3XncVt9139JzLVePcF9E0vQICOOffbL5fChi58DA4pAbqV19Ev38DHpffqJcyP3BhYUC8bv15rZsjs1oGh9RQa+2LWOwAN//ZpT+vrRYbuLBw27xRnCPjcFxpRPV7OQ2X5+ioGwSj2C9lcEgjqt8js5x13rFYEExOrB+bOTIGhzSi+j1L3FnnHYtNovz4DW8bm8mE9nFII6rfI7Mc6dWx3MCFUQyKhYoER0T8BvBx4K3A1Zk5s0i79wKfAtYBn8vMfUMrUmq5fs8Sd9b5q8Z9EmVk5vC/NOKtwMvAHwD/oldwRMQ64C+Aa4FZ4CHgA5n5+HKfPz09nTMzPbNIktRDRDycmdMraVvkjCMznwCIiKWaXQ0cz8ynqrZ3AzuAZYNDkjQ4Te7jmAKe6Xo9C/zSYo0jYjewG2Dz5s2DrUySGqLEpMyBBUdEfA24uMdbH83ML6/kI3psW/S6WmbuB/ZD51LVioqUpBYrNSlzYMGRme9Z40fMApd1vb4UOLHGz5SkkVHqVsBNnsfxEHBlRGyJiAuAG4GDhWuSpMYoNSmzSHBExPsiYhZ4J3B/RByqtm+KiAcAMvMscAtwCHgC+GJmPlaiXklqolKTMosER2Z+KTMvzcwLM/MtmXldtf1EZl7f1e6BzPwbmfnzmflvS9QqSU1V6lbATR5VJUlaQqlJmQaHJLVIr+G3/3PPu4Zag8EhSQMwiPkVTbknSpNHVUlSK80f4OdOnSZ59QB/4Mjcmj63KfdEMTgkqc8GdYBvyj1RDA5J6rNBHeCbck8Ug0OS+mxQB/hSw28XMjgktd6BI3Ns33eYLXvuZ/u+w2vuS1irQR3gd26basRdBh1VJanVmjLSqNsg51c04SZSBoekViu10N9ymnCAHxQvVUlqtaaMNBonBoekVmvKSKNxYnBIarWmjDQaJ/ZxSGq1Ugv9jTODQ1Lr1e2ILnGf7lFicEgaK00cvts29nFIGitNWSiwzQwOSWPF4btrZ3BIGisO3127IsEREb8REY9FxMsRMb1Eu7+MiKMR8UhEzAyzRkmjyeG7a1eqc/y7wC7gD1bQ9u9n5g8GXI+kMeHw3bUrEhyZ+QRARJT4ekljbpTXkRqGpvdxJPDfI+LhiNhduhhJ0gDPOCLia8DFPd76aGZ+eYUfsz0zT0TEzwJfjYj/nZnfXOT7dgO7ATZv3ryqmiVJyxtYcGTme/rwGSeqn89FxJeAq4GewZGZ+4H9ANPT07nW75Yk9dbYS1UR8dci4g3zz4FfodOpLkkqqNRw3PdFxCzwTuD+iDhUbd8UEQ9Uzd4C/FlEPAr8L+D+zPxvJeqVJL2q1KiqLwFf6rH9BHB99fwp4BeGXJokaRmROXrdARFxEni6dB1dNgBtmovStnrBmoelbTW3rV4oV/PPZebGlTQcyeBomoiYycxFZ8g3TdvqBWselrbV3LZ6oR01N7ZzXJLUTAaHJKkWg2M49pcuoKa21QvWPCxtq7lt9UILaraPQ5JUi2cckqRaDI4BaNv9RmrU+96IOBYRxyNizzBr7FHLmyLiqxHxZPXzokXavVTt30ci4uCw66xqWHK/RcSFEXFP9f63IuLy4Vd5Tj3L1XtzRJzs2q8fLlFnVz1fiIjnIqLnyhLR8XvV3+c7EfH2YdfYo6blar4mIl7s2scfG3aNS8pMH31+AG8FrgK+AUwv0e4vgQ1tqBdYB3wPuAK4AHgU+FsFa/5dYE/1fA/wyUXa/bjwvl12vwH/GPhs9fxG4J6G13sz8OmS+3VBPX8PeDvw3UXevx74ChDAO4BvtaDma4D/WrrOxR6ecQxAZj6RmcdK17FSK6z3auB4Zj6VmX8F3A3sGHx1i9oB3FE9vwPYWbCWpaxkv3X/Xe4F3h3lblbTtH/nZWVnxewXlmiyA7gzOx4EJiPikuFU19sKam40g6OsNt1vZAp4puv1bLWtlLdk5rMA1c+fXaTd6yJiJiIejIgS4bKS/fZKm8w8C7wIvHko1b3WSv+df7267HNvRFw2nNJWrWn/d1fqnRHxaER8JSLeVrqYbqVuHdt6w77fyFr1od5evwEPdEjeUjXX+JjN1T6+AjgcEUcz83v9qXBFVrLfhr5vl7CSWv4EuCszfxoRH6FztvSugVe2ek3avyv1bTpLgPw4Iq4HDgBXFq7pFQbHKuWQ7zfSh+9aa72zQPdvlpcCJ9b4mUtaquaI+H5EXJKZz1aXHZ5b5DPm9/FTEfENYBuda/jDspL9Nt9mNiLOB95IucsYy9abmc93vfxD4JNDqGsthv5/d60y80ddzx+IiP8YERsysxHrbnmpqpAW3m/kIeDKiNgSERfQ6cQtMkqpchC4qXp+E/Cas6aIuCgiLqyebwC2A48PrcKOley37r/L+4HDWfWQFrBsvQv6B24AnhhifatxEPhgNbrqHcCL85c5myoiLp7v54qIq+kcq59f+k8NUene+VF8AO+j81vOT4HvA4eq7ZuAB6rnV9AZsfIo8BidS0aNrbd6fT3wF3R+Yy9Wb1XLm4GvA09WP99UbZ8GPlc9/2XgaLWPjwIfKlTra/Yb8Anghur564D/DBync++ZKwrv2+Xq3Vv9n30U+FPgbxau9y7gWeBM9f/4Q8BHgI9U7wfwmervc5QlRjo2qOZbuvbxg8Avl665++HMcUlSLV6qkiTVYnBIkmoxOCRJtRgckqRaDA5JUi0GhySpFoNDklSLwSENWET83WpBwNdVKwY8FhF/u3Rd0mo5AVAagoj4N3RmiE8As5m5t3BJ0qoZHNIQVOs+PQT8PzrLR7xUuCRp1bxUJQ3Hm4CfAd5A58xDai3POKQhqO53fjewBbgkM28pXJK0at6PQxqwiPggcDYz/1NErAP+PCLelZmHS9cmrYZnHJKkWuzjkCTVYnBIkmoxOCRJtRgckqRaDA5JUi0GhySpFoNDklSLwSFJquX/A4zyttikCzDzAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "dummy_data = np.loadtxt('../code/data/dummy.csv')\n", + "x_1 = dummy_data[:, 0]\n", + "y_1 = dummy_data[:, 1]\n", + "\n", + "# Set order of polynomial fit for Question 1\n", + "order_1 = 5\n", + "x_1p = np.vstack([x_1**i for i in range(1, order_1+1)])\n", + "x_1s = (x_1p - x_1p.mean(axis=1, keepdims=True)) / \\\n", + " x_1p.std(axis=1, keepdims=True)\n", + "y_1s = (y_1 - y_1.mean()) / y_1.std()\n", + "plt.scatter(x_1s[0], y_1s)\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "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: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [00:04<00:00, 1204.72draws/s]\n" + ] + } + ], + "source": [ + "# Fit Linear Model\n", + "with pm.Model() as model_l:\n", + " α = pm.Normal('α', mu=0, sd=1)\n", + " β = pm.Normal('β', mu=0, sd=10)\n", + " ϵ = pm.HalfNormal('ϵ', 5)\n", + "\n", + " μ = α + β * x_1s[0]\n", + "\n", + " y_pred = pm.Normal('y_pred', mu=μ, sd=ϵ, observed=y_1s)\n", + "\n", + " trace_l = pm.sample(2000)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "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: [ϵ, β, α]\n", + "Sampling 2 chains: 100%|██████████| 5000/5000 [00:08<00:00, 599.00draws/s]\n" + ] + } + ], + "source": [ + "# Fit quadratic model\n", + "order_4 = 2\n", + "x_1p = np.vstack([x_1**i for i in range(1, order_4+1)])\n", + "x_1s = (x_1p - x_1p.mean(axis=1, keepdims=True)) / \\\n", + " x_1p.std(axis=1, keepdims=True)\n", + "y_1s = (y_1 - y_1.mean()) / y_1.std()\n", + "\n", + "with pm.Model() as model_p:\n", + " α = pm.Normal('α', mu=0, sd=1)\n", + " β = pm.Normal('β', mu=0, sd=10, shape=order_4)\n", + " ϵ = pm.HalfNormal('ϵ', 5)\n", + " μ = α + pm.math.dot(β, x_1s)\n", + "\n", + " y_pred = pm.Normal('y_pred', mu=μ, sd=ϵ, observed=y_1s)\n", + "\n", + " trace_p = pm.sample(2000)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 2/2 [00:00<00:00, 27.55it/s]\n", + "100%|██████████| 2/2 [00:00<00:00, 13.82it/s]\n" + ] + } + ], + "source": [ + "# Generate posterior predictive draws\n", + "posterior_predictive_draws = 2\n", + "y_l = pm.sample_posterior_predictive(trace_l, posterior_predictive_draws,\n", + " model=model_l)['y_pred']\n", + "\n", + "y_p = pm.sample_posterior_predictive(trace_p, posterior_predictive_draws,\n", + " model=model_p)['y_pred']" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'Linear Posterior Predictive Fit')" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X+cXHV97/HXOyGBRS0BEgj5gUHMTYtABRf8Ab0XBQxyFSJiBL1VrJTSFih6G0mqReTWSyStXim2itGqVcGAGFPFGxGKCopmA0L4YSogXnYXIYCBIgv59bl/nLMwO5nZnV/nx+y+n4/HPHbmzJkznxnC+cz5/vh8FRGYmZk1alLRAZiZWXdx4jAzs6Y4cZiZWVOcOMzMrClOHGZm1hQnDjMza4oThzVE0h9J2lh0HGUm6TOS/rboOBoh6QxJN1c8flrSy1o4zrskfa+z0Y35nndLOibP97SRnDhsBEkPSjquentE/CgiFhQRUzVJF0namp7sNkv6saTXtnnMmySd2c4xIuLsiPhf7RyjIp55kiL9jE+n/12WduLYtUTEiyPigQZj2qXidV+NiDd2Oh5JX5S0peLzPy3pHel7viIibkr3u0jSVzr9/jY6Jw4rtcqTVJWvR8SLgRnAzcC1kpRfZCNJmtzGa+t9RoBp6ec8HbhQ0glNvr6bXZomtOHb14sOyBJOHNYQScdI6q94/KCkv5Z0p6QnJX1d0m4Vz79Z0s8rrggOrXhuqaT7Jf2npHskvbXiuTMk3SLpk5KeAC4aLa6I2Ap8CZgJ7C1pkqQPS/q1pEclfVnSHumxd5P0FUmPp3Gtk7SvpI8BfwRcnv6yvTzd//clXS/pCUkbJS2uiPOLkv5Z0nWSfge8Pt32dxX7/Kmk+9LXr5E0q+K5kPSXkn4J/HKs7z8ifgLcDRxc7/VjxLt3GsNTkn4GHFh5/PR4L0/v90j6h/Q7fFLSzZJ6gB+mu29Ov6fXVjZ5pU11f1913G9J+kB6f5akb0jaJOlXks4b63PXMnxVnCbRvwHekcZzRyvHsxZEhG++PX8DHgSOq7H9GKC/ar+fAbOAvYB7gbPT5w4HHgVeDUwG3pPuv2v6/NvT100C3gH8Dtgvfe4MYBtwLrAL0FMjlouAr6T3dwVWAA+lj/8EuA94GfBi4FrgX9Pn/gz4N2D3NK5XAb+XPncTcGbFe7wIeAh4bxrH4cBjwCvS578IPAkclX6O3dJtf5c+/4Z0/8PTGP8R+GHF8QO4Pv3uan3Geek+uwBK3+cZ4Nhar28g3quAVel+BwMDwM1V8bw8vf/p9PuYnX5Pr0s/w/MxVbzujOHjAP81jUHp4z2BoYr/1uuBC4Gp6X+fB4CFdf4dPv9djvZvtPLfgm/53XzFYe24LCIGI+IJkhPyK9Ptfwp8NiJ+GhHbI+JLwHPAawAi4ur0dTsiaX74JXBkxXEHI+IfI2JbRAzVee/FkjaTnKheBSxKt78L+EREPBARTwPLgNPS5pytwN4kJ8jtEbE+Ip6qc/w3Aw9GxL+kcdwGfAM4tWKfb0XELenneLbq9e8CvhARt0XEc2kcr5U0r2KfSyLiiVE+IyQn/yeAlcDSiLihzuvrxps2o70NuDAifhcRd5Fcpe1E0iSS5PtXETGQfk8/Tj/DWH5Eklj+KH18KvCTiBgEjgBmRMTFEbElkv6UzwGnjXK8v06vDDdLeqyB97ecjNe2UcvHbyruP0PyyxLgpcB7JJ1b8fzU4eclvRv4AMkvWEiuDKZX7PtQA++9KiL+R43ts4BfVzz+Ncm/832BfwXmAldJmgZ8BfhQJM1d1V4KvDpNTsN2SY/RSJyzgNuGH0TE05IeJ/kV/2ADrx82PSK21Xmu8vWjxTsjvV+5f+V3NOL9SK6e7m8gthEiIiRdRdIf80PgnSTf8XB8s6rim0ySbOr5+4j4cLNxWPacOCwLDwEfi4iPVT8h6aUkvzSPJfk1ul3Sz0maY4a1U7J5kOQkNWx/kqavR9IT8EeBj6a//K8DNgKfr/GeDwE/iIjjR3mv0eIcEYekF5Fc7Qw0+PpGVL6+brzpFcc2kqT5i3Tz/nWO+RjwLEkfSHWfQSPxXgl8T9JykqbK4f6rh4BfRcT8Bo7RDJf3LoCbqqyWKWlH8vCt2R8YnwPOlvRqJV4k6b9LeglJG3sAmwAkvZe0w7dDrgTeL+kASS8G/jfJCKxtkl4v6ZD0RPoUSdPV9vR1j5C0uw/7NvBfJP2xpCnp7QhJf9BgHF8D3ivplZJ2TeP4aUQ82IHPWEvdeCNiO0lfz0WSdpd0EEm/004iYgfwBeATaWf25LQTfFeS/2Y7GPk9Vb/+9nS/lcDaiBi+wvgZ8JSkC9LO98mSDpZ0RJuf+xFgXtrEZjnxl221XEfSqTl8u6iZF0dEH0k/x+XAb0k6q89In7sH+AfgJyT/0x8C3NKZsIHkpPevJE0lvyL59TzcZDYTuIYkadwL/IAXmlI+RdIf8FtJl0XEfwJvJGmDHyRplvs4SSfxmNK+iL8l6Wd4mOQX/Gjt+W1pIN5zSJoEf0PS8fwvoxzur4ENwDqS/pWPA5Mi4hngY8Atab/Da+q8/krgOJLkORzfduAtJP1gvyK5slkJ7NHkR612dfr3cUm3jbqndczw6AczM7OG+IrDzMya4sRhZmZNceIwM7OmOHGYmVlTxuU8junTp8e8efOKDsPMrGusX7/+sYiY0ci+4zJxzJs3j76+vqLDMDPrGpLqVRPYiZuqzMysKU4cZmbWFCcOMzNrihOHmZk1xYnDzMya4sRhZmZNKWw4rqS5wJdJKpbuAK6IiE9V7SOSqqUnkiwUdEa6slmmVt8+wIq1GxncPMSsaT0sWbiARYfNzvptzcy6QpHzOLYB/zMibkvXaVgv6fq07PawNwHz09urgX9O/2Zm9e0DLLt2A0Nbk2UaBjYPsezaDQBOHmZmFNhUFREPD189pGsJ3EuyrGalk4EvR+JWYJqk/bKMa8Xajc8njWFDW7ezYu3GLN/WzKxrlKKPI13G8zDgp1VPzWbkOsn97Jxcho9xlqQ+SX2bNm1qOZbBzUNNbTczm2gKTxzp8p7fAM6PiKeqn67xkporT0XEFRHRGxG9M2Y0VG6lplnTeprabmY20RSaOCRNIUkaX42Ia2vs0g/MrXg8h2RZzMwsWbiAnimTR2zrmTKZJQsXZPm2td25Cj55MFw0Lfl756r8YzAzq1LkqCoBnwfujYhP1NltDXCOpKtIOsWfjIiHs4xruAM8r1FVdUdw3bkK/u082Jo2kT35UPIY4NDFmcRiZtaIwtYcl3Q08CNgA8lwXIC/AfYHiIjPpMnlcuAEkuG4742IMcve9vb2RjdUx60ewQXJ1c0lpxzCopsWJsmi2h5z4f135RilmU0EktZHRG8j+xZ2xRERN1O7D6NynwD+Mp+I8jfaCK5Fz/bXftGTdbabmeWk8M7xiWzUEVx7zKn9onrbzcxy4sTRptW3D3DU8hs5YOl3OGr5jay+faDh1446guvYC2FK1fNT0u1mZgVy4mjDcB/FwOYhghdmmTeaPEYdwXXoYnjLZUmfBkr+vuUyd4ybWeHG5dKxeVmxdiPHb/8BH5y6ill6jMGYzqXbFrNi7dSGRmGNOYLr0MVOFGZWOk4cbeh96noumbKS3bUFgDl6jOVTVrLsKYA3NHSMRYfNdg0sM+sqbqpqw7KpVz+fNIbtri0sm3p1QRGZmWXPiaMN+/JYU9tb4tnjZlYyThxtUJ2hsfW2N2149viTDwHxwuxxJw8zK5ATRzuyHjJ7w8UvlBwZtnUo2W5mVhAnjnZkPWS23ixxzx43swJ5VFW7shwyu8ecOvWqPHvczIrjK44SW3fguQzF1BHbhmIq6w48t6CIzMycOGoryUim8++ZzwVbz6R/x3R2hOjfMZ0Ltp7J+ffMLyQeMzNwU9XOSrQOxsDmIQY4mjVbjh6xXV7G1swK5CuOaiUZybT69oG6Nee9jK2ZFanopWO/IOlRSTVXJpJ0jKQnJf08vWVfGrYkI5lWrN1Yc3F1QTHL2JqZpYq+4vgiyep+o/lRRLwyvWX/s78k62DUW6sjwLWtzKxQhSaOiPgh8ESRMeykJOtg1GuOmu1mKjMrWNFXHI14raQ7JH1X0ivq7STpLEl9kvo2bdrU+ruVZB2MUdfqMDMrkJJlvQsMQJoHfDsiDq7x3O8BOyLiaUknAp+KiDHHovb29kZfX1/HY83b6tsH6q/VYWbWQZLWR0RvI/uW+oojIp6KiKfT+9cBUyRNLzis3Cw6bDa3nPgYv9r3Am559hQW3bTQBQ7NrHClnschaSbwSESEpCNJEt3jBYeVnxLNKTEzG1Zo4pB0JXAMMF1SP/ARYApARHwGOBX4c0nbgCHgtCi6bS1Po80pceIws4IUmjgi4vQxnr8cuDyncMqnJHNKzMwqlbqPY8Irw5ySktTtMrPycOIos6LnlHgFQjOrwYmjzIqeU1KSul1mVi6lHlVlZLtQ1Fjcx2JmNfiKw+orQx+LmZWOE4fVV3QfSzV31JuVghNH0cp8Mmy0jyWPz+COerPSKLxWVRa6plZV9cxwSH7RF1BUsWV5fYZPHpwmjSp7zIX311zOxcyaMG5qVY1742HUUl6fwR31ZqXhxFGk8XAyzOszuKPerDScOIrUiZNh0X0keZ3Qy9ZRbzaBOXEUqd2TYRk6jPM6oRc9GdLMnucJgEUaPundcHHStLPHnOSE2+jJsAzVc9v9DM2+lxOFWeGcOIrWzsmwJH0kq7cfxYrnLmPw2SFm7dbDku0LWJRrBGaWJzdVdbMSdBivvn2AZdduYGDzEAEMbB5i2bUbWH37QG4xmFm+Ck0ckr4g6VFJNQfiK3GZpPsk3Snp8FwCK7rDuVEl6DBesXYjQ1u3j9g2tHU7K9ZuzC0GM8tX0VccXwROGOX5NwHz09tZwD9nHlHeHc7tJKkSdBgPbh5qaruZdb+iVwD8oaR5o+xyMvDldLnYWyVNk7RfRDycWVB5djiPtab4navG7nQuuMN41rQeBmokiVnTemrsbWbjQdk7x2cDlXUm+tNt2SWOPDucx5p1PVpSSa2+fYAVazcyuHmIWdN6WLJwAYsOm935WOtYsnABy67dMKK5qmfKZJYsXJBbDGaWr6KbqsaiGttqFteSdJakPkl9mzZtav0dm+1wbqepabQk1UApjzJ0TC86bDaXnHIIs6f1IGD2tB4uOeWQXJOXmeWr7Fcc/cDcisdzgMFaO0bEFcAVkBQ5bPkdj72wdtG+Wh3OYzU1jWWPOXUK981p6MpntI7pPE/ciw6b7URhNoGU/YpjDfDudHTVa4AnM+3fgOY6nNst8DfaqKgGrnxy6ZjulhFmZpabQq84JF0JHANMl9QPfASYAhARnwGuA04E7gOeAd6bS2CNdji32x8y1qzrMa58Mu+YbveKyszGpaJHVZ0+xvMB/GVO4TRvtKamRtVLUg2U8si8Y7oMJU3MrHTK3sdRbs30h7RijCuf4X6FzEZVlaSkiZmVixNHO/Is8FdHph3TnbiiqlL08GEza58TR7vGccXWdQeey8HrP0yPtjy/bSimcteB53JEC8cbHj483LQ2PHwYcPIw6yJlH1VlBTr/nvlcsPVM+ndMZ0eI/h3TuWDrmZx/z/yWjue6Vmbjg684rK7BzUMMcDRrthw9YrtaHO7rulZm44OvOKyuesN6Wx3u2+njmVkxnDisriULF9AzZfKIbe0M912ycAGnTv0xN089jwd2fSc3Tz2PU6f+2HWtzLqMm6q6XSMVdFu06LDZzH7o28y9bQX7xCYe1QweOnwJRxw2WiX8UY43+RbePGUlu2x/FoA5eozlk1eyy+Q/BMbnAAOz8chXHN0s67VD7lzFERs+wkw2MUkwk00cseEjrR//houfTxrDdtn+bOMlWsysFJw4ulm7tbLyPr4nFJqNC04c3SzrE3Gnj1+CNdLNrH1OHN0s6xNxp49fgjXSzax9ThzdLOsTcaePX4I10s2sfR5V1c2yrpXV6PGbGdnVRomWlupcZTjqzGyicuLodkXXysppzY6W6lx5PRGzTLipyuprZLhv1iO7Ui3VucoqtrxWRfTqi1ZShSYOSSdI2ijpPklLazx/hqRNkn6e3s4sIs4Jq5ETb05DbFuqc5VFbFnPncn7fcxaUFjikDQZ+DTwJuAg4HRJB9XY9esR8cr0tjLXICe6Rk68OQ2xbanOVRax5XSFldv7mLWgyCuOI4H7IuKBiNgCXAWcXGA8Vq2RE29OQ2xbqpuVRWx5TWL0ZEkrsSITx2ygcnm5/nRbtbdJulPSNZLm1juYpLMk9Unq27RpU6djnZDWHXguQzF1xLahmMq6A899YUOzQ2xbbLdfdNhsLjnlEGZP60HA7Gk9XHLKIaOPqspi+G9ekxg9WdJKTBFRzBtLbwcWRsSZ6eM/Bo6MiHMr9tkbeDoinpN0NrA4It4w1rF7e3ujr68vq9AnjKOW38irnrqeD+6yill6nMHYm0u3LWb97x3PLUvH/M+ws+pRTpBcAXTTXI5an2HSFNj1JTD0284N+R0P35V1FUnrI6K3kX2LHI7bD1ReQcwBBit3iIjHKx5+Dvh4DnFZqtMLOY3abt/myTC3tcyr57b07AlbnoahJ5LtnRryW4L17M3qKTJxrAPmSzoAGABOA95ZuYOk/SLi4fThScC9+YY4sc2a1sNAjSTR8sJLGbXb576WeeXcmU8e/ELSGNahZFj4HB2zOgrr44iIbcA5wFqShLAqIu6WdLGkk9LdzpN0t6Q7gPOAM4qJdmLq9EJOWbXb15rjcfz2H/Cab/237OdAuBPbJqBCZ45HxHXAdVXbLqy4vwxYlndclhj+td6xJqBjL6zdbt/mCKzquRwnTbqZ5VNWsjtbkg1ZzhjfY04616LGdrNxyiVHbFSLDpvdueaejNrtq5vUPrjLKnbXlpE7dar5qFpGydCszJw4LF8ZtNsvWbhgRB/HLD1We8csmo/ciW0TkBOHdb3qJrVHNYOZ1JjLk1XzkTuxbYJx4rBxYUST2p2/c/ORWYZcHdfGHy8YZZYpX3HY+OTmo/q8uJW1yYnDrE25zVrvBC9uZR3gpiqzNgzPWh/YPETwwqz11bcPFB1abS7Xbh3gxGHWhhVrN3L89h9w89TzeGDXd3Lz1PM4fvsPRl+ZsEjtznT3qoSGE4dZW3qfup7lU1YyZ9JjTBLMmfQYy6espPep64sOrbZ2yr54VUJLOXGYtWHZ1Kt3mqW+u7awbOrVBUU0hnYWt3Izl6XcOW7WpMrO8Pt3qz1LfV/qzF6vJc9RTu3MdHdBR0s5cZg1obqE++COvZkzaeckoUZnqRcxyqnVocou6GgpN1VZ9yhBx2x1CfdLty3mmarldZuapd5NzT85rS9v5ecrDusOJZl/UF3Cfc2Oo2FrUpF3zqTHm29q6qbmHxd0tFShiUPSCcCngMnAyohYXvX8rsCXgVcBjwPviIgH847TSqDJZWezmpRXa1XENTuOZv3uLa7D3m3NP56RbxTYVCVpMvBp4E3AQcDpkg6q2u19wG8j4uXAJ/Ga4xNXE7/Ms5yU1/FVEY+9ECZNGblt0hQ3/1ipFdnHcSRwX0Q8EBFbgKuAk6v2ORn4Unr/GuBYScoxRiuLJuYf1FpKdmjr9o5Mylt02GwuOeUQZk/rQcDsaT1ccsoh7V3NVP+T9j9xK7kxm6oknQN8NSJ+2+H3ng1UXqP3A6+ut09EbJP0JLA37DzWUdJZwFkA+++/f4dDtcI1sdJedT/EWNub1dFVEW+4GLZXrVa4fUs2qxWadUgjVxwzgXWSVkk6oYO/+GsdJ1rYJ9kYcUVE9EZE74wZM9oOzkqmiVLps6b17Pz6UbYXqps6x6H1kW0lGBFnnTPmFUdEfFjS3wJvBN4LXC5pFfD5iLi/jffuB+ZWPJ4DDNbZp1/SLsAewBNtvKd1s9E6Zism0V3fM5MLp76Na7a8bsQuv3tuG6tvHyhX5dpu6hxvdWRbSUbEWec01McREQH8Jr1tA/YErpF0aRvvvQ6YL+kASVOB04A1VfusAd6T3j8VuDGNxewFVTWUdh96mOVTVnL6breO2G3z0NbyVa7tprkRrc456aa5KtaQMROHpPMkrQcuBW4BDomIPycZIvu2Vt84IrYB5wBrgXuBVRFxt6SLJZ2U7vZ5YG9J9wEfAJa2+n42jtU4Me2y/Vn+iit32rVTneQd002rFbbarNZtzXE2pkbmcUwHTomIX1dujIgdkt7czptHxHXAdVXbLqy4/yzw9nbewyaAOiegfaJ2vahOdZJ3TL0muLKt1Ndqs1o3NcdZQ8a84oiIC6uTRsVz93Y+JLMm1TkBParpNbeXspO8WhlLmLfarNZNzXHWENeqsu5X58T00OFLOjtZLyOrbx/gqOU3csDS73DU8huTPpgy9gu02qzWTc1x1hCNx77m3t7e6OvrKzoMy9O3PwDrvwixHTQZXnUGvPkTrFvzWebetoJ9YhOPagYPHb6EI076s6KjfV51tV1Ikts9k09DNUeeCy7anF+ANmFIWh8RvY3s6yKH1v3uXAV3fC1JGpD8veNrAByx4WvAEAhmsomZGz4C8/Ysza/derPcH5k8nZls2vkF7hewEnBTlXW/es06679YvuaeKvU66i/Z8nb3C1hpOXFY96s3rDO2195eomGg9Trq+37vePcLWGm5qcq6X73hnppcO3mUqLlnycIFNfs4lixcAIe+wYnCSslXHNb96g33fNUZpW/uyaTarlnGfMVh3W+0len2f025JtHV0NFqu9XKMomwLHFYRzhx2PhQPft6uBrr8InqlCsm3omqLMUFyxKHdYybqmz8yXrWdbeUCC/LJMKyxGEd48Rh40+WJ6oylgKppyzFBduJo1uS9ATjxGHjT5YnzG769dzEcruljKObkvQE48Rh40+WJ8yy/IpvRI3RZkPsyroDz801jHUHnstQTB0ZR0wdO45uStITjBOHjT9ZVmMty6/4Rhy6mHWHfJSBmM6OEP07pnPBlvfx7nUvzXUxq/Pvmc8FW8+kf0dFHFvP5Px75o/+wjyStJvCWlLIqCpJewFfB+YBDwKLI+K3NfbbDmxIH/6/iDipeh+znYw2PLddx144coQQlG5uSKXz75nPwHOXjdy4I1nMKq+5IoObhxjgaNZsOXrEdo21LkrW63h4tFfLirriWArcEBHzgRuov7LfUES8Mr05aVjjDl0M778rqST7/rs6dyLoshLh9Wph5bmYVb2yKmOui5L1Oh5uCmtZUfM4TgaOSe9/CbgJuKCgWMyaU2/FvhKaNa2HgRpJIs/FrEYtqzKaLK8cobv6q0qmqMSxb0Q8DBARD0vap85+u0nqA7YByyNidb0DSjoLOAtg//3373S8Zl2p5ZN2Bw03ia1Yu5HBzUPMmtbDkoULGmsqyzJJe0nblmWWOCR9H5hZ46kPNXGY/SNiUNLLgBslbYiI+2vtGBFXAFdAspBT0wGbjUNtnbTHuy7rryqTzBJHRBxX7zlJj0jaL73a2A94tM4xBtO/D0i6CTgMqJk4zKy2RZNvYdGuF8Nu/bDrHJh8IZBfU1v1KocDm4dYdm0y5qUjCazVOlhZN4WNY4UsHStpBfB4RCyXtBTYKyI+WLXPnsAzEfGcpOnAT4CTI+KesY7vpWPNUneugtV/ATu2VmycBD3TYOi3uZwsj1p+Y81+ltnTerhl6RvaO3j1yChIrhpKPGChrJpZOraoUVXLgeMl/RI4Pn2MpF5JK9N9/gDok3QH8O8kfRxjJg0zq/DdC6qSBsAOGHqCvGZjtzOya/XtAxy1/EYOWPodjlp+487zTzwyqhCFdI5HxOPAsTW29wFnpvd/DBySc2hm48vQE2PvM3yizegXeqsjuxpq4mplZNR4KfFe4OfwzHEzy3QI6pKFC+iZMnnEtkZGdq1Yu3HEaDCAoa3J5MXnNTuTf7zUvyr4czhxmI1nPXs1tl+GQ1BbXeWwoSauZicJ1mna6r9mWe2msLIquInOCzmZjWeveCv0fX7EpgBU8Xjb5N3YJeMhqK2scthQE1ezI6PqXFnN0uOdH+2VpYInL/qKw2y8unMV3PG1EZsCuGXHwSMKDi7deiartx9VTIyjaLiJq5nyMnWurAZjb6BGU1hZFVxs01ccZuNVjeYMAfP4DUdvGVn48Cc5Fj1sVCaTF2tM+nsmpnLptheSTZ51vFpW8ORFJw6zklp9+0B7J81RmmWqlfVk2UoT16gqmrZ2PNnP4I69uXTbYtbseKFyb551vFpW8ORFJw6zEurIbOs6tZiGm2UqdcXJslPS+ldrhr/jHU3W8SrLcN4Ci226j8OshBoaijqWGiOOtk3ejf/DaSO25V30sCxaGu01XobztslXHGYl1JF1NGo0Z+xy7IUcvf0ofpJ30cOy/Eqv0nRT2GjDYEvwefLixGFWQh1bR6NGc8Yich5uOp5W2htlGGzbfVJdxE1VZiXU6mzramPWempSS8cbR/WknumptVIE/HbKPiy7dgMDm4cIXuiT6poJhU1y4jAroVZnW1ca7mDv1Mms5eONo5X2Lt36Dp6JqSO2PRNT+egzp7bfJ9VF3FRlVlLtDkUdrYO9leO2fLxxtNLel54+kicmbeGDu6xilh5nMIaH89aeQFnWYc7tcuIwG6c60sHeieONo5X2Zk3rYc3mo1mz5egR2ydLbK+xttF4HebspiqzcareSavVk1nLxzt0cbKw0h5zASV/G11o6c5V8MmD4aJpyd+Ch73W63s6/dVzO9In1S18xWE2Ti1ZuGDEJEJo72TW1vFamaxWwtFYo5VB6X3pXhNmVFVRS8e+HbiIZJW/I9MFnGrtdwLwKWAysDIiljdyfC8da5bo9BDRXIecfvLgOn0jc5NihtZRzSwdW9QVx13AKcBn6+0gaTLwaZKlZfuBdZLWePlYs8Z1utZTx2tHjWYcjcYab4paOvZeAEmj7XYkcF9EPJDuexVwMuDEYTYRjKPRWFkocsJhmTvHZwOV/2r60201STpLUp+kvk2bNmUenJllrNHV/UrWgZ6HTs/RaVZmiUPS9yXdVeN2cqOHqLGtbodMRFwREb0R0TtjxozWgjaz8mhgstexAAAKzElEQVRkNNYELTrYkSKYbcisqSoijmvzEP3A3IrHc4DBNo9pZt1krNFYrRQdLGnBxWZ0eo5Os8rcVLUOmC/pAElTgdOANQXHZGatyKo5qdkO9FpXKNeeBRft0VXNXJ2eo9OsQhKHpLdK6gdeC3xH0tp0+yxJ1wFExDbgHGAtcC+wKiLuLiJeM2tDls1Jza69XesKZbgFvIuauTpVBLNVhczjyJrncZiVSJbzMaonCULSgV5vZvpF0xilq7RjcWUx4qn6mK///Rn8+y82dew9umEeh5lNFFnOx2h27e16Q3w7GFdHlv1t4JjfWD/QdMXkTilzH4eZjQfNNic169DFyRXCRZuTv6N1dNca4tvhuLIY8VT0KKpqThxmlq1G52PkYcQQX9hp1H8H4spixFPRo6iqOXGYWbbaqY6bVTzvvwsuehJOuaJzcaUjx+7f7V3cPPU8Tpp084in2xnxVPQoqmru4zCz7LVSHbdBbXVEdyquik76ScCcSY+xfMpK2Aprdhzd9oinTlc6bpcTh5l1rSw6oltSY5jv7trCp6b8E3+jq3no8CUccdgJLR9+tHLuRXDiMLOu1enlcVtWZySWBDPZxMwNH4F5e7Z1dZNrZeIxuI/DzLpWaTqNxxqJNVwGZZxw4jCzrlWaTuNGhvmOo3VEnDjMrGsVXXrjeTsN861hHK0j4j4OM+tarXQaZ7YA0vAIrXplUIqYt5IRJw4z62rNdBrnMgqr2TIorSqwPLwTh5lNGB0ZhdXICTvDeSvPx1B5VTNc2Xf4vTPmPg4zmzDaHoXVaIn4rJezHW0Bqxw4cZjZhNH2KKxGTth5LGebZcXhBhS1kNPbJd0taYekuvXfJT0oaYOkn0vyAhtm1pa2R2E1csLO42og64rDYyjqiuMu4BTghw3s+/qIeGWjC4yYmdWz6LDZXHLKIcye1oOA2dN6mlvTopETdh5XAwVXHC6kczwi7gWQNNauZmYd1VbpjmMvHHuobb3Fojp5NZDXyK06yj6qKoDvSQrgsxFxRb0dJZ0FnAWw//775xSemU0ojZywG0kunYqloNL0mSUOSd8HZtZ46kMR8a0GD3NURAxK2ge4XtIvIqJm81aaVK6AZM3xloI2MxvLWCfsgq8G8pBZ4oiI4zpwjMH076OSvgkcSWP9ImZmxSnwaiAPpR2OK+lFkl4yfB94I0mnupmZFaio4bhvldQPvBb4jqS16fZZkq5Ld9sXuFnSHcDPgO9ExP8tIl4zM3tBUaOqvgl8s8b2QeDE9P4DwB/mHJqZmY2htE1VZmZdJesyI1kfvwllH45rZlZ+WRcdLLioYTVfcZiZtSvrMiMFFzWs5sRhZtaurMuMFFzUsJoTh5lZu7IuOlhwUcNqThxmZu3KuuhgwUUNqzlxmJm169DF8JbLYI+5gJK/b7mscx3XWR+/SYoYf2Wdent7o6/Py3eYmTVK0vpGl6/wFYeZmTXFicPMzJrixGFmZk1x4jAzs6Y4cZiZWVOcOMzMrClOHGZm1pSiFnJaIekXku6U9E1J0+rsd4KkjZLuk7Q07zjNzEqjRGXVi7riuB44OCIOBf4DWFa9g6TJwKeBNwEHAadLOijXKM3MymC4rPqTDwHxQln1gpJHIYkjIr4XEdvSh7cCtSp1HQncFxEPRMQW4Crg5LxiNDMrDZdV38mfAN+tsX028FDF4/50W02SzpLUJ6lv06ZNHQ7RzKxAE6WsuqTvS7qrxu3kin0+BGwDvlrrEDW21S2sFRFXRERvRPTOmDGj/Q9gZlYWJSurntnSsRFx3GjPS3oP8Gbg2KhdabEfmFvxeA4w2LkIzcy6xLEXjlw6FiZeWXVJJwAXACdFxDN1dlsHzJd0gKSpwGnAmrxiNDMrjZKVVc/simMMlwO7AtdLArg1Is6WNAtYGREnRsQ2SecAa4HJwBci4u6C4jUzK9ahiwtLFNUKSRwR8fI62weBEyseXwdcl1dcZmY2tjKMqjIzsy7ixGFmZk1x4jAzs6Y4cZiZWVOcOMzMrClOHGZm1hQnDjMza4pqV/vobpI2Ab8uOo4apgOPFR1EExxvthxvthxvc14aEQ0V+huXiaOsJPVFRG/RcTTK8WbL8WbL8WbHTVVmZtYUJw4zM2uKE0e+rig6gCY53mw53mw53oy4j8PMzJriKw4zM2uKE4eZmTXFiSMjkt4u6W5JOyTVHWIn6UFJGyT9XFJfnjHWiKXRmE+QtFHSfZKW5hljVRx7Sbpe0i/Tv3vW2W97+v3+XFLuq0iO9X1J2lXS19PnfyppXt4xVsUzVrxnSNpU8Z2eWUScaSxfkPSopLvqPC9Jl6Wf5U5Jh+cdY1U8Y8V7jKQnK77bYtaGHUtE+JbBDfgDYAFwE9A7yn4PAtOLjrfRmElWY7wfeBkwFbgDOKigeC8Flqb3lwIfr7Pf0wV+p2N+X8BfAJ9J758GfL3k8Z4BXF5UjFWx/FfgcOCuOs+fCHwXEPAa4Kclj/cY4NtFf69j3XzFkZGIuDciNhYdRzMajPlI4L6IeCAitgBXASdnH11NJwNfSu9/CVhUUByjaeT7qvwc1wDHKl1TuQBl+u87poj4IfDEKLucDHw5ErcC0yTtl090O2sg3q7gxFG8AL4nab2ks4oOpgGzgYcqHven24qwb0Q8DJD+3afOfrtJ6pN0q6S8k0sj39fz+0TENuBJYO9cottZo/9935Y2/VwjaW4+obWkTP9eG/VaSXdI+q6kVxQdTC2FrDk+Xkj6PjCzxlMfiohvNXiYoyJiUNI+wPWSfpH+KslEB2Ku9Us4szHdo8XbxGH2T7/jlwE3StoQEfd3JsIxNfJ95fqdjqGRWP4NuDIinpN0NsnV0hsyj6w1ZfpuG3EbSc2opyWdCKwG5hcc006cONoQEcd14BiD6d9HJX2TpKkgs8TRgZj7gcpfmHOAwTaPWddo8Up6RNJ+EfFw2vzwaJ1jDH/HD0i6CTiMpB0/D418X8P79EvaBdiD4pozxow3Ih6vePg54OM5xNWqXP+9tisinqq4f52kf5I0PSJKVazRTVUFkvQiSS8Zvg+8Eag52qJE1gHzJR0gaSpJZ27uI5VSa4D3pPffA+x0xSRpT0m7pvenA0cB9+QWYWPfV+XnOBW4MdKe0gKMGW9VH8FJwL05xtesNcC709FVrwGeHG7eLCNJM4f7tyQdSXKOfnz0VxWg6N758XoD3krya+c54BFgbbp9FnBdev9lJKNW7gDuJmkuKnXM6eMTgf8g+dVeWMwk/QA3AL9M/+6Vbu8FVqb3XwdsSL/jDcD7Cohzp+8LuBg4Kb2/G3A1cB/wM+BlBf87GCveS9J/r3cA/w78foGxXgk8DGxN/+2+DzgbODt9XsCn08+ygVFGOJYk3nMqvttbgdcVGW+9m0uOmJlZU9xUZWZmTXHiMDOzpjhxmJlZU5w4zMysKU4cZmbWFCcOMzNrihOHmZk1xYnDLGOSjkgLAu6WVgu4W9LBRcdl1ipPADTLgaS/I5kh3gP0R8QlBYdk1jInDrMcpHWf1gHPkpSR2F5wSGYtc1OVWT72Al4MvITkysOsa/mKwywH6VrnVwEHAPtFxDkFh2TWMq/HYZYxSe8GtkXE1yRNBn4s6Q0RcWPRsZm1wlccZmbWFPdxmJlZU5w4zMysKU4cZmbWFCcOMzNrihOHmZk1xYnDzMya4sRhZmZN+f8UdRJgIzQJUgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Plot Linear Posterior Predictive fit\n", + "plt.scatter(x_1s[0], y_1s)\n", + "plt.scatter(np.tile(x_1s[0], posterior_predictive_draws), y_l)\n", + "\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')\n", + "plt.title('Linear Posterior Predictive Fit')" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0.5, 1.0, 'Quadratic Posterior Predictive Fit')" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEWCAYAAABmE+CbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8XXWZ7/HP0zRpU9QGSBHSFgrKdESKFlJkKI4eUSsIJVSpd8QjAuMBdMYppygHkUFb6JwBq+NLoHLUUYEKWOuAp6KISkewKZeWixwugm3DJQELAilt0+f8sdZud9J9WfuyLjv7+3698mqy9tprPXsnXc9ev8vzM3dHRERkTNoBiIhINighiIgIoIQgIiIhJQQREQGUEEREJKSEICIigBKCiIiElBCkImY2zczczMbW8ZhfNLNl9Tpe2szsATN7Z9pxRGFm3zWzS8Lv325mD1d5nG+b2f+qb3Qlz7e/mb1kZi1JnbMZKCGMAmZ2mpmtN7NXzOxpM/uWmU1MO65CzOydZrYxf5u7f83dT6/iWLeb2ZbwwjBgZjeZ2X41xudm9sZajuHub3b322s5Rl48p5nZUPgaXzSze83shHoceyR3/527T48Y0x0jnnuWu/9LvWMysyfMbDB8/bmvLnf/s7u/xt2Hwv1uN7OK/4ZkOCWEBmdmXwAuBRYAE4GjgGnAL8ysNeFYzMyS/ps6291fA/wN0AFcnvD5d6r1rqnE838fvsYO4DvAcjPbq97nz7ATw4t/7qsv7YBGKyWEBmZmrwO+Apzj7v/X3be5+xPAfOBA4KPhfjubBcKfh31KN7OFZvaYmf3VzB40s5PzHmsxs38NP4E/Drx/RAy3m9lXzWw18ApwkJl9ysweCo/3uJmdGe67B/BzoCv/056ZXWRmP8g75jFm9l9mttnMNpjZaeXeC3d/HrgRODQ8xkQz+76Z9ZvZk2Z2QS5Zmdkbzew3ZvZC+LquD7f/NjzcfWFsHwq3nxB+Mt8cxnVYXqxPmNn/NLN1wMtmNjbc9u7w8XFmdoWZ9YVfV5jZuPzfQ/j8p4H/U+Y17gCuAdrD97ng88vEO9PM7g5/N9cD4/MeG/l3MTW86+o3s+fM7Jtm9ibg28Dfhe/R5nDf/Kanh/LvYsL3ZMDMDg9/Pirv93ufVdG8ZnlNl2b2VeDtwDfDmL5Z6fEkoITQ2I4m+A99U/5Gd3+J4ML73ojHeYzgP9REggTzA9vV9PIZ4ARgJtANfLDA8z8BnAG8FngSeDZ8zuuATwGXm9nh7v4ycBzQV+zTnpntH8b+DWAS8Fbg3nIvwMw6gQ8A94SbvhG+noOAdwCnhrEA/AvwC2BPYEq4L+7+9+Hjbwljuz68iF0DnAnsDVwJrMxd1EMfIUiUHe6+fURoXyK4a3sr8BbgSOCCvMf3BfYCDgjfw1KvcSxwOvAS8Eih55eK18zagBXAf4TP+XH4nhU6VwvwnwS/z2nAZOA6d38IOIvwrsXdOwo8/drwPcmZAwy4+91mNhm4GbgkjOGfgRvNbFKp116Ku38J+B3h3aK7n13tsZqdEkJj6yT4jzbyIgTwFMEFtSx3/7G797n7Dne/nuBic2T48HzgCnffEH4KX1TgEN919wfcfXt4l3Kzuz/mgd8QXHzfHvE1fQz4pbtfGx7rOXcvlRCWhp9S7wtf8z+FF7MPAee7+1/Du6b/TZC4ALYRXEC73H2Lu99R4Lg5nwGudPe73H3I3b8HvEpwkd8ZQ/j+DBZ5PRe7+7Pu3k+QcD+R9/gO4Mvu/mqR5wMcFb7GpwkutCe7+wtFnl8q3qOAVoLf5zZ3vwFYU+ScRwJdwAJ3fznC+5TvR8BcM5sQ/vzRcBvAx4Fb3P2W8O/tVqAXOL7E8VaEdxObzWxFxBikCkoIjW0A6LTCbcf7Af1RDmJmp+Y1MWwmaHbpDB/uAjbk7f5kgUPkP46ZHWdmd5rZ8+Hxjs87XjlTCe5YojrX3TvcfbK7fyy86HYCbSNifZLgUy7AeYABf7BgRNB/L3H8A4Av5F2QNocxduXts6HwUyHcb2Qc+c/td/ctpV4gcGf4Gjvd/Sh3/2WJ55eKtwvY5MNLHBf6fRI+58kiHzZKcvdHgYeAE8OkMJddCeEA4JQR8R1D8PdaTE/4+jvcvafSeCQ6JYTG9nuCT3/z8jda0FZ/HPCbcNPLwIS8XfbN2/cA4GrgbGDvsAngfoILJgSfuqfmPXf/AnHsvMCETSk3Av8KvD483i15xytXb30D8IYy+5QzwK67gJz9gU0A7v60u3/G3bsImla+ZcVHFm0Avpp3Qepw9wnufm3ePqVeU1+BOPKbyWqtPz/y+aXifQqYbGaWt3+h32fuOPsX+bARJeZcs9FJwINhksgd9z9GxLeHuy+OcMxSVMe/DpQQGljYbPAV4Btm9j4zazWzaQRtwwPAD8Nd7wWON7O9zGxf4PN5h9mD4D9TP4CZfYqwYza0HDjXzKaY2Z7AwjJhtQHjwuNtN7PjGN6X8QywtxUfFvtD4N1mNj/sMNzbzN5a5pzDhEMRlwNfNbPXhknvn4AfhK/xFDObEu7+F4LXP5QX30F5h7saOMvM3maBPczs/Wb22ojhXAtcYGaTwn6OC3NxxKRUvL8HthP8Psea2Tx2NQ2O9AeCBLI4PMZ4M5sdPvYMMCXskyjmOoLf+z+w6+4Agtd+opnNsWDAwviwM3tKwaNEN/L3JlVQQmhw7n4Z8EWCT+R/Bf5EcDfw7rATF4JOxPuAJwja86/Pe/6DBO3rvyf4TzUDWJ13iquBVeHz72ZEB3aBeP4KnEtwQf4LQfvxyrzH/0hwkXw8bDLoGvH8PxM0MX0BeJ4gmb0lynsxwjkEd0aPA3cQXJSuCR+bBdxlZi+FsX3O3f8UPnYR8L0wtvnu3kvQLv/N8PU8CpxWQRyXELSRrwPWE7yHl5R8Rg1KxevuWwnuJk8LH/sQRX6fYVI9EXgj8GdgY7g/wG3AA8DTZjZQ5PlPEfxNHc3wv7cNBHcNXyT40LCBYMh0rdeirwMfNLO/mNnSGo/VtMy1YtqoEraHfwWYHV5cRUQiUUIYhczsE8A2d78u7VhEpHEoIYiICKA+BBERCTVU7ZPOzk6fNm1a2mGIiDSUtWvXDrh72YmqDZUQpk2bRm9vb9phiIg0FDMrNgFxGDUZiYgIoIQgIiIhJQQREQGUEEREJKSEICIigBKCiIiElBBERARQQhARkVBzJYR1y+HyQ+GijuDfdcvTjkhEJDMaaqZyTdYth5+dC9vCZWtf2BD8DHDY/PTiEhHJiOa5Q/jVxbuSQc62wWC7iIg0UUJ4YWNl20VEmkzzJISJRZZsLbZdRKTJNE9COPZCaG0fvq21PdguIiJNlBAOmw8nLoWJUwEL/j1xqTqURURCzZMQRESkJA07heJ3CeuWB6OQXtgY9DUce6HuKERk1GqeO4RKh53mEsgLGwDflUA0mU1ERqnmSQiVDjvVvAURaTLNkxAqHXaqeQsi0mSaJyFUOuxU8xZEpMk0T0KodNip5i2ISJNpnlFGEFz8o44Syu2nUUYi0iSaKyFUqpIEIiLS4JqnyagaWj9BRJpIancIZjYV+D6wL7ADuMrdvx7nOVfcs4klqx6mb/MgXR3tLJgznZ6ZkwvvrPUTRKTJpHmHsB34gru/CTgK+B9mdkhcJ1txzybOv2k9mzYP4sCmzYOcf9N6VtyzqfATNA9BRJpMagnB3Z9y97vD7/8KPAQU+bheuyWrHmZw29CwbYPbhliy6uHCT9A8BBFpMpnoQzCzacBM4K64ztG3ebCi7ZqHICLNJvWEYGavAW4EPu/uLxZ4/Awz6zWz3v7+/qrP09XRXtF2zUMQkWaTakIws1aCZPBDd7+p0D7ufpW7d7t796RJk6o+14I502lvbRm2rb21hQVzphd+gtZPEJEmk+YoIwO+Azzk7v8W9/lyo4kijzKCmuYhVDSiSUQkA9KcmDYb+ASw3szuDbd90d1vieuEPTMnJ3JRzo1oynVi50Y05WLQOgsikkWpJQR3vwOwtM4fp1IjmnpaVmt+g4hkUuqdyqNRyRFNmt8gIhmlWkYlVNsP0NXRzhEv3sp5Y5fTZQP0eSeXbZ/P2te9R/MbRCSzlBCKKNsPUMIVhzzCoWuX0W5bAZhiA1zauoz7D5kGj00Jl+UcQfMbRCRlajIqouKZzXlmPfaNnckgp922Muuxb2h+g4hklhJCERXPbM5XqllI8xtEJKPUZFREV0c7mwpc/IvObM43sUyzkNZZEJEM0h1CERXPbM6nZiERaUC6QyiiqpnNOVp+U0QakLl72jFE1t3d7b29vWmHISLSUMxsrbt3l9tPTUZx0fKbItJg1GQUBy2/KSINSAkhDqXKU+QlBFVEFZEsUUKIQ4TyFLXMhBYRiYP6EOIQYfnNWmZCi4jEQQkhDmXmIaxZeSXXv/IZHh/3Ue5oO5e5Y+7YuVukmdAiIjFQk1EcSsxDWLPySg5dewHtY3YVvlvcugy2wcodx0SbCS0iEgMlhLgUKU8x9e4luxW+m2BbOW/scm71d0SbCS0iEgM1GSVsH+8vuL3LnmPRvBnqUBaR1CghJOxZm1Rke6eSgYikSgkhYRsOX8Cgtw3bNuhtbDh8QUoRiYgElBASNmvumdx/xCU8zSR2uPE0k7j/iEuYNffMtEMTkSan4nYiIqOcituJiEhFlBDSpIqoIpIhmoeQFlVEFZGM0R1CWkpVRBURSYESQloiVEQVEUmSEkJaIlREFRFJkhJCWspURAXU6SwiiVKnclpKVEQF1OksIolTQkhTkYqoQORlOEVE6kVNRlkVd6ezmqNEZAQlhKyKs9M51xz1wgbAdzVHKSmINDUlhKyK0ulcLc2BEJEClBCy6rD5cOJSmDgVsODfE5fWp/9AcyBEpAB1KmdZqU7nWkycEjYXFdguIk1LdwjNKM7mKBFpWLpDaEbl5kDkWXHPJpasepi+zYN0dbSzYM50LfUpMkopIUhRK+7ZxPk3rWdw2xAAmzYPcv5N6wGUFERGoVSbjMzsGjN71szuTzOOphNx2OmSVQ/vTAY5g9uGWLLq4QSDFZGkpN2H8F3gfSnH0HwiDjvt2zxinzLbRaSxpZoQ3P23wPNpxtCUig47HT7yqKujveBuxbaLSGNL+w5BSliz8kqevuiN7PjyRJ6+6I2sWXllfQ5cdHipDWs2WjBnOu2tLcP2aG9tYcGc6fWJQ0QyJfMJwczOMLNeM+vt7+9PO5zErFl5JYeuvYB96WeMwb70c+jaC+qTFI69ELACD/iwZqOemZP5/qwnuXP853h83Ee5c/zn+P6sJ9WhLDJKZT4huPtV7t7t7t2TJk1KO5zETL17Ce22ddi2dtvK1LuX1H7ww+YDXvix/OakdcuZtf7Lw5LSrPVfVs0jkVEq8wmhWe3jhe+G9vGB+pxg4tQi2/Oak1TzSKSppD3s9Frg98B0M9toZp9OM54sedYK3w09a531OUGU2cqqeSTSVNIeZfQRd9/P3VvdfYq7fyfNeLJkw+ELGPS2YdsGvY0Nhy+ozwmiFM/Tus8iTUUzldO0bnnR8hGz5p7JGoK+hH18gGetkw1HLGDW3DPrd/5yxfOOvXD4Mp5QtOaRSlyIND5zL9K5mEHd3d3e29ubdhj1MXLNZAgutvUqcR1BpIt4iaSVf5z8EhcQDE9dNG+GkoJIBpjZWnfvLref7hDSkvKayZHrFEUowb1k1cO8Z+g3nNe2nC4boM87uWz7fJasalNCEGkgGmWUlpQ7bOtZp6j7xVtZ3LqMKWMGGGMwZcwAi1uX0f3irfUKV0QSoISQlpQ7bOtZp+j8th8zYcSciQm2lfPbflxVbCKSDiWEtKS8SE096xS9nsJzI4ptF5FsUkJIS5xrJkdQzzpFVuSupth2EckmdSqnKa41kyPIdfbWZahoBcNTRSS7lBCaWE/LanrGXQzjN8K4KdByIVBFgqpgSU4RyS4lhGY1ch5EbtU0GH4hjzAPYedzlABEGpr6EJpVlMJ1EZfaFJHRQQmhWUWZB1FJtdN1y+HyQ+GijuBfJQ2RhqOE0KyizIOIOnmu1jsJJRORTFBCaFZR5kFEnTxXy7oJ65az/afnDEsm2396zvCkoIQhkgglhGYVZR5E1MlzNZTheOXnFzJ2aMuwbWOHtvDKz8NzqB9DJDEaZdTMyo0MijqcdOKU8ILN7tvLGD/4dOntKRcBFGkmSghSWpThpDVMTOvbsTdTxuxe4qJvx95MgdSLAIo0EzUZSe1qKMOxrO3jvDJiZbhXvI1lbR8PftCqbSKJ0R1Ck6r7CmdVTkx76/vP4MKfbOfzfh1d9hx9vjdX8GGOef8ZwQ4qiyGSGCWEJhR5cZwEBOf7LB9adWzh5KSyGCKJUUJoQqUWx0ljhbOemZNLnnfF0GyWvLqUvi2DdI1vZ8HQdHoSjE+kWSghNKF6Lo4TtyzdzYiMdupUbkL1XBwnbvVc6lNESlNCaEL1XBwnbo10NyPS6MomBDM728z2TCIYSUbPzMksmjeDyR3tGDC5o51F82Zksglm5F3L3DF3cEfbuTw2/mMqYyFSZ1H6EPYF1pjZ3cA1wCp393jDkriV68jNigVzpu/sQ5g75g4Wty5jgm0NHiy2hoOIVKXsHYK7XwAcDHwHOA14xMy+ZmZviDk2kWF3M+eNXb4rGeRELaInImVF6kMI7wieDr+2A3sCN5jZZTHGJgIESWH1wncxZcxzhXfIahkLVWmVBhOlD+FcM1sLXAasBma4+z8ARwAfiDk+kV0XVoq0VNarjEU9L+Cq0ioNKEofQicwz92fzN/o7jvM7IR4whIJjVz7eaR6lbGIusZ0VKrSKg0oSh/ChSOTQd5jD9U/JJE8hS6sORUU0avqPLX0T6hKqzQgzVSWbCt6ATX4x/tjP4+/sJFjFt9WeRHAGtaIEEmLJqZJtiVV/rrI8fp8bzZtHsTZVTZjxT2byh8v6mpzIhmihCDZVujCigWfvus5cqfAeQYZx6XbhjdHRS6bUcMaESJpUZORpG/d8uLlrYeVv94AGDtHG9VzYlqBMtsL+09k5Y5jdts1ctmMKteIEEmLEoKkK8rontyF9fJDd2+Xr+PInZFltl8Zvx1e2bbbflksAihSD0oIkq5KhmfGOHKnUJnt1jHGyWNX84Ux19NlA/R5Z7Ca25zP1nw+kSxSH4Kkq5KLfLGOZBtTc19CoTLbx/E7vjb2aqaMGWCMwZQxAyxuXUZPy+qaziWSVbpDkHRVMjyz0PrKAD5Uc19C3+ZB5o65g/PGLt95NzDBttDO8NpJY4e2aHKZjFqp3iGY2fvM7GEze9TMFqYZi6SkkuGZuZE71rL7YzUWufvka/7A4tZlw+4G9uSlwjtrcpmMUqklBDNrAf4dOA44BPiImR2SVjySkgqHZ64Yms0O31H4WDVcqM9rvX63SqpmRXbW5DIZpdJsMjoSeNTdHwcws+uAk4AHU4xJ0hBxeGau47fb9mbKmIHdd6jhQj1h8OloO2pymYxiaTYZTQbyG483httECsp1/F62fT6veNvwB2u9UBdLJu17aXKZNI00E0KhG/Ld6hub2Rlm1mtmvf39/QmEJVmVmxC2cscxLNx2Oht3dLLDjY07Omu/UBfry3jzyTVELNJY0mwy2ghMzft5CtA3cid3vwq4CqC7u1tLdzaxro52NuUlhZVbg1nEkzvaWX3Yu2o7eIGZyhz8XrjvR/UriS2ScWneIawBDjazA82sDfgwsDLFeCTjFsyZTnvr7iOMXn51e7SCc+UcNj+ooHrR5uDfR35R35LY+eqxGI9WZJM6S+0Owd23m9nZwCqgBbjG3R9IKx7JvlzZ6a/87AHevuXXu+YMDHVyxU8+DHw2WmnqqOKaGV2PxXjqvaCPCCnPQ3D3W9z9b9z9De7+1TRjkcbQM3MyJ41ZvducgYvtKu69+ar6niyu0tv1WIyn3gv6iKDSFdKATt/6g93mDEywrZy+9Qf1PVFcaxrU485DK7JJDFS6QhpO15jnKtpetUIdzfmluctYcc8mlqx6ePfV1uqxmppWZJMYKCFIw9nSvi8TBp8qvL3eJ6tyTYNC1VPPv2k9AD2FajJVeudRj2OIjKAmI2k4E467mO0t44dt294yngnHZaf9vFD11J2rrdVjNTWtyCYx0B2CNJ7D5gd/uHlNOWMraMpJQrFV1XZur8dqalqRTepMCUEaU8YvhvmT6EZur5tSS48mKStxSM3UZCQSg0KT6NpbW1gwZ3p9TpCbh/DCBsB3zUNIenJaVuKQulBCEIlBz8zJLJo3g8kd7RhBeY1F82bUb+JcVuYhZCUOqQs1GYnEpGfm5PrOnM6XlXkIWYlD6kJ3CCKNKK5Z1I0ah9SFEoJII4prFnWjxiF1oYQg0oiyMg8hK3FIXagPQaRRZWTo7Yqh2Sx5dSl9WwbpGt/OgqHp9KQdlFRFCUGkQRWtlZRwDHf85Ftcz3V0jRug75WYSpFLItRkJNKAcrWSNm0exNlVK6kuCwVV4N6br+Jiuyr+UuSSCCUEkQZUslZSghIrRS6JUEIQaUBlayUlJLFS5JIIJQSRBlSsJlJdayVFsKV934q2N6QmWrtancoiaaixINwVhzxC19rL2I8B+ryTy7bP521jH+Ejr/4KLtoB1gJHnAYn/Ft8r4GwFPlPz2Hs0Jad27JWirwmTbZ2tbl72jFE1t3d7b29vWmHIVKbkRcZCCZzRR2/X+D52zFacGzkvt2fjj0prFl5JVPvXsI+PsCz1smGwxcwa+6ZsZ4zMZcfWmRluqnwj/cnH0+VzGytu3eX209NRiJJq7UgXIHnjy2UDADWfjfaMatsFllxzyZOXXMAR235Oge9+kOO2vJ1Tl1zQOKjnWLTZLWalBBEklbrRaaSi5EPld+nhhLWWRntFJuiNZl8VPYnKCGIJK3WgnCVFI6zlvL71HDHUu1opxX3bGL24ts4cOHNzF58W3bvKArVasoZhWs/KCGIJK3WgnAFL1JF/isfcVr549Vwx1LNaKesTKqLZFitpgJG2doPSggiSau1IFyh58+7MuhAzt0RWEv0DuUa7liqWRmu4ZqZDpsfdiAX7KUZVf0JGnYqkoZaC9MVev5h86sbUXTshYVHPUW4Y8nVK6qkplLkZqYahubGUudp4pQiI45Gz9oPSggio0TVF8HcRbbKi2+lK8N1dbSzqUBSGNbMVMP4/1yTVO4uJNcklYu1ajUkzkahhCAyCtR8EUywlPaCOdO54yff4vNcR5cFE+uu4MMcM+ezu3Yq1dFdJs5STVI1JYQaE2cjUEIQGQViuwjGoKdlNSe0Lts5u3mKDbC4ZRljW94ChBfXGjq6617nqcZZ5Y1EncoiKaj3sMusFLuL5FcXDyt1AQQ/54/WqaGju651nmqYo9GIlBBE4lJk9m8cwy5rvggmWcAtyqf/GobmVjPyqahaZ5VHkaHieUoIInEo8ckyjmGXNV0Ek/4UHOXTfw1Dc3tmTmbRvBlM7mjHgMkd7SyaN6O6prO4S1dk7A5Exe1E4lCiKNqBz1xKof91Bvxp8fujHb9Au/aKodnVjTJKuoBbrcX9khT3e5PQex+1uJ06lUXiUOKTZaRhl6UUGZLZc+JSehZWcUFNuoBbXKN18pNk+57BtsG/RD5+wWG7x15YsLz32HoNNc1Y8Tw1GYnEoUSzSM1t3PVu1661tlIWjGx6GXw++IrYDFOsX+eCx9/Ewm2ns3FHJzvc2Lijk4XbTmfF0Oz6xF3sPW7fM5V+BSUEkTiU6BSttY3bi3x6LLa9nDVvOIdBbxu2bdDbWPOGc6o6XllxtJsXSpL5yiTMYv061961gRu2Hs0xW5dy0Ks/5JitS7lh69H1K7NR6O9kTCtsfSmVfgU1GYnEoUyzSKWze/M9Qyf70l9ke+U+/+DBHLHtdM4bu5wue44+35vLts9n7YMHs3puVSGWVsOks6KiJMNCbfWhYsNzh4r0sdZtOO9h8+HPdwbrVvhQUINqbBtsfXn4frW+PxEpIYjEJabZv4u2nsKi1mVMsK07t73ibSzadgpfr+J4fZsH2cQxrNx6zLDtFtcchjjazYvVGcpXohR4sX6dFrOCSaFua1evWw73/WjXuhU+tHsyyEmgX0FNRiINpvd17ynYrt37uvdUdby6TuSKIo4+i4PfW36fEosFFevX+cjbptZvTkMh5Zq68iXQp5NKQjCzU8zsATPbYWZlh0KJyC4L5kzn1pZ3DGvXvrXlHVVfpOo6kSuKWteDKOSRX5Tfp9iaBhSfu3BJz4z6zWkoJOqn/oSK6KXVZHQ/MA+4MqXzizSsakpOx3q8Smv9xDHstNyFNcIFtVi/Ti39PWUVa+pq3wva9ki8flKqE9PM7Hbgn9090mwzTUwTyZisTDIrNsELgjuDrBakS+j9izoxLfN9CGZ2hpn1mllvf//uIytEJEVJ1PqJolgz1Lyrgxm/WUwGULREx4qh2amsOR1bk5GZ/RIKjoL7krv/NOpx3P0q4CoI7hDqFJ6I1ENWZto28loFI0ajxbbATwSxJQR3f3dcxxaRjMjSspIJLvITpzTXtsh8k5GIZFgcI4bSkpEy1GmubZHWsNOTzWwj8HfAzWa2Ko04RKRGNZSpLivJC3SGylAnPi8kj8pfi0j2JD16KekS4HlGVln9b387iRvXbhrWbNTe2lLT/IdRM8pIRJpQ0qOXInaO13vp00JVVm9cu4kPHDE5vslwJaiWkYhkT9KjlyJ0jscx+qdYB/Kv/9jP6oXvquqYtdAdgohkT9JrNEToHI9j6dM0O5ALUUIQkexJevRShM7xOC7eaXYgF6ImIxHJnsPms+aJvzD17iXs4wM8a51smLGAWXHOMygzj6HmpU8LWDBn+rBmKIi5sGAZSggikjkr7tnE+WsOYHDbrhUe2te0sGjqpkQ6VwuJ4+JdtLBgy2q4PPlZ10oIIpI5ac7WLabeVWbzjzvsGCOH3ObmREDsSUF9CCKSOVnrbM3paVnN6nHn8qfxH2P1uHODT/L1lmLBQCUEEcmcrHW2AsnNZk6xYKASgohkTuKruEWR1Cf3pIfc5lEfgohkTjXt9SNLQNSjfX+YpD65H3th4bIdo3gJTREZLSpdQjPbFuDIAAAF7ElEQVSiSpauTGQNgaRKfae4toMSgohUL8URMfkSGZWU5Cf3lNZ2UB+CiFQvI0toJjIqKc5S3xmhOwQRqV5GltCMYxZxQaNkVbZidIcgItVLcURMvkyOSmpASggiUr2MLKHZM3Myi+bNSGUNgdFETUYiUr0UR8SMVMmoJClMCUFEajPK29WbiZqMREQEUEIQEZGQEoKIiABKCCIigXXL4fJD4aKO4N96VzFtAOpUFhHJSAmOtOkOQUQkIyU40qaEICKSkRIcaVNCEBHJSAmOtCkhiIhkpARH2pQQRESaoLR1FBplJCICKsGB7hBERCSkhCAiIoASgoiIhJQQREQEUEIQEZGQEoKIiABKCCIiElJCEBERAMzd044hMjPrB55MO44COoGBtIOogOKNl+KNl+Kt3AHuPqncTg2VELLKzHrdvTvtOKJSvPFSvPFSvPFRk5GIiABKCCIiElJCqI+r0g6gQoo3Xoo3Xoo3JupDEBERQHcIIiISUkIQERFACaEqZnaKmT1gZjvMrOhwMjN7wszWm9m9ZtabZIwj4oga7/vM7GEze9TMFiYZ44g49jKzW83skfDfPYvsNxS+t/ea2coU4iz5fpnZODO7Pnz8LjOblnSMI+IpF+9pZtaf956enkacYSzXmNmzZnZ/kcfNzJaGr2WdmR2edIwFYioX8zvN7IW89zd763O6u74q/ALeBEwHbge6S+z3BNDZCPECLcBjwEFAG3AfcEhK8V4GLAy/XwhcWmS/l1J8T8u+X8BngW+H338YuD7j8Z4GfDOtGEfE8vfA4cD9RR4/Hvg5YMBRwF0NEPM7gf9MO85SX7pDqIK7P+TuD6cdR1QR4z0SeNTdH3f3rcB1wEnxR1fQScD3wu+/B/SkFEcpUd6v/NdxA3CsmVmCMebL0u+3LHf/LfB8iV1OAr7vgTuBDjPbL5noCosQc+YpIcTLgV+Y2VozOyPtYMqYDGzI+3ljuC0Nr3f3pwDCf/cpst94M+s1szvNLOmkEeX92rmPu28HXgD2TiS63UX9/X4gbIK5wcymJhNaVbL091qJvzOz+8zs52b25rSDGWls2gFklZn9Eti3wENfcvefRjzMbHfvM7N9gFvN7I/hp4i6q0O8hT65xjYmuVS8FRxm//D9PQi4zczWu/tj9YmwrCjvV6LvaRlRYvkZcK27v2pmZxHc3bwr9siqk6X3Nqq7CWoKvWRmxwMrgINTjmkYJYQi3P3ddThGX/jvs2b2E4Lb9lgSQh3i3QjkfyKcAvTVeMyiSsVrZs+Y2X7u/lTYDPBskWPk3t/Hzex2YCZBO3kSorxfuX02mtlYYCLpNSmUjdfdn8v78Wrg0gTiqlaif6/14O4v5n1/i5l9y8w63T3twnc7qckoJma2h5m9Nvc98F6g4OiDjFgDHGxmB5pZG0EnaOIjd0IrgU+G338S2O0Ox8z2NLNx4fedwGzgwcQijPZ+5b+ODwK3edi7mIKy8Y5og58LPJRgfJVaCZwajjY6Cngh18yYVWa2b64PycyOJLj+Plf6WQlLu1e7Eb+Akwk+obwKPAOsCrd3AbeE3x9EMJLjPuABgqabzMYb/nw88P8IPmWnGe/ewK+AR8J/9wq3dwPLwu+PBtaH7+964NMpxLnb+wVcDMwNvx8P/Bh4FPgDcFDKf7fl4l0U/q3eB/wa+NsUY70WeArYFv7tfho4CzgrfNyAfw9fy3pKjPbLUMxn572/dwJHpx3zyC+VrhAREUBNRiIiElJCEBERQAlBRERCSggiIgIoIYiISEgJQUREACUEEREJKSGI1MDMZoXF4MaHs9MfMLND045LpBqamCZSIzO7hGBWcjuw0d0XpRySSFWUEERqFNYGWgNsIShHMJRySCJVUZORSO32Al4DvJbgTkGkIekOQaRG4XrO1wEHAvu5+9kphyRSFa2HIFIDMzsV2O7uPzKzFuC/zOxd7n5b2rGJVEp3CCIiAqgPQUREQkoIIiICKCGIiEhICUFERAAlBBERCSkhiIgIoIQgIiKh/w+1qxVAB3sNdgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Plot Quadratic Posterior Predictive fit\n", + "plt.scatter(x_1s[0], y_1s)\n", + "plt.scatter(np.tile(x_1s[0], posterior_predictive_draws), y_p)\n", + "\n", + "plt.xlabel('x')\n", + "plt.ylabel('y')\n", + "plt.title('Quadratic Posterior Predictive Fit')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}