statistical modeling with TensorFlow

data.csv

modeling.ipynb

tf_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### デモデータの生成"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "n_data     = 50\n",
    "a_true     = 100.0\n",
    "b_true     = 1.0\n",
    "sigma_true = 2.0\n",
    "d_X = np.arange(n_data)\n",
    "d_Y = np.random.normal(a_true + b_true*d_X, sigma_true, n_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 可視化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x25a02e07630>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAFiZJREFUeJzt3X2MXFd5x/Hfj8TgLYU14I1xvE6dIocmsWiIJmkkKpQiAhtKY/qCEkRFRCO5FUGlqBElRWr6okhUtLiiLZGWxkqCICG0pVgUuZiUyq2UkK7JCxsMjsEJXtu1jVIvATbBTp7+MXed2e283bt3Xu6Z70dazdwzd2fOVdbPnDz3nOc4IgQASNeLBt0BAEBvEegBIHEEegBIHIEeABJHoAeAxBHoASBxBHoASByBHgASR6AHgMSdPegOSNLatWtj06ZNg+4GAFTK3r17fxARE53OG4pAv2nTJs3MzAy6GwBQKbaf7OY8UjcAkDgCPQAkjkAPAIkj0ANA4gj0AJC4oZh1AwCjZt/Ree2aPabDJxe0Yc2Ypras04Xrx3vyWYzoAaDP9h2d1/Seg5pfOKX146s1v3BK03sOat/R+Z58HoEeAPps1+wxjY+t0vjYKr3IPvN81+yxnnwegR4A+uzwyQW9bPXSzPnLVp+twycXevJ5BHoA6LMNa8b09DOnl7Q9/cxpbVgz1pPPI9ADQJ9NbVmn+YVTml84pecjzjyf2rKuJ59HoAeAPrtw/bi2vfF8jY+t0tH5ZzQ+tkrb3nh+z2bdML0SAAbgwvXjPQvsyzGiB4DEEegBIHGkbgCgS/1czVomRvQA0IV+r2YtU8dAb3uH7eO2Zxva/tT2YdsPZz9va3jtZtsHbH/H9lt71XEA6Kd+r2YtUzepmzsk/Z2ku5a1b4+Iv2pssH2RpOskXSzpXElftX1BRDxXQl8BYGAOn1zQ+vHVS9q6Wc06DOmejiP6iNgj6aku32+rpHsi4tmIOCjpgKTLV9A/ABgKRVazDku6ZyU5+vfbfjRL7bwia9sg6VDDOXNZGwBUWpHVrMOS7ika6G+T9BpJl0g6Kumvs3Y3OTeavYHtbbZnbM+cOHGiYDcAoD+KrGbtd/GyVgpNr4yIM19Htj8l6UvZ4ZykjQ2nTko60uI9piVNS1KtVmv6ZQAAwyTvatYNa8Y0v3BK42OrzrT1snhZK4VG9LbXNxz+uqTFGTk7JV1n+yW2z5e0WdKDK+siAFRTv4uXtdJxRG/7bklXSlpre07SLZKutH2J6mmZJyT9riRFxGO275X0LUmnJd3IjBsAo2ox3dM46+bayyb7PuvGEYPPmtRqtZiZmRl0NwCgkEFNobS9NyJqnc5jZSwArMCwTKFsh0APACswLFMo2yHQA8AKDMsUynYI9ACwAv3e/7UIAj0ArMCwTKFsh0APACvQ7/1fi2DjEQBYoX7u/1oEI3oASByBHgASR+oGAJYZhs1CysSIHgAaVGGla16M6AGMrGYj98aVrpLOPO6aPVbZUT2BHkDSWqVhFkfu42Orlozcf/TsKf3Cq1++5D2GbaVrXqRuACSrXRqmVY2a+YXTQ7/SNS8CPYBktSs41qpGzctXnz30K13zItADSFa7gmOtatRcfO7wr3TNixw9gGS127N1ass6Te85KKke/J9+5rTmF06d2QGqyoF9OUb0AJLVruBYFWrUlIURPYBkddqzNbWReysEegBJG5Vg3g6pGwBIHIEeABJHoAeAxJGjB5CE1CpOlokRPYDKS7HiZJkI9AAqr12pA3QR6G3vsH3c9myT126yHbbXZse2/QnbB2w/avvSXnQaABq1K3WA7kb0d0iaWt5oe6OkqyR9v6H5akmbs59tkm5beRcBoL1WdWuqXHGyTB0DfUTskfRUk5e2S/qQpGho2yrprqh7QNIa2+tL6SkAtNCu1AEKzrqxfY2kwxHxiO3GlzZIOtRwPJe1HS3cQwBJK2O2TKdSB6Mud6C3/TOSPiLpLc1ebtIWTdpke5vq6R2dd955ebsBIAGtdnkqUlyMUgetFZl18xpJ50t6xPYTkiYlfcP2q1UfwW9sOHdS0pFmbxIR0xFRi4jaxMREgW4AqDpmy/RH7kAfEd+MiHMiYlNEbFI9uF8aEf8jaaek92Szb66QNB8RpG0ANMVsmf7omLqxfbekKyWttT0n6ZaIuL3F6V+W9DZJByT9RNJ7S+ongArIm29vtzEIK13L44imKfS+qtVqMTMzM+huAFiBxnx7445N7fLtrX7nzRdO6Kv7TuR6r1Fke29E1Dqdx8pYAKUokm9vtcvT/mM/JndfIoqaASjF4ZMLWj++eklbN/n2ZrNlbv+vJ1q+Fymd/BjRAyhFmatTW73Xi88yxcsKINADKEWZq1NbvZclUjoFEOgBlKJVvr1IWqXVez37XDAdswBy9ABK02p1apG8erP3ajcdE60xogeQ276j89q+e79u+vwj2r57f9sceZmbglC8rBgCPYCeBu4yyxyUmR4aJaRugBGXt7BYY+CWdOZx1+yxpucXnXbZCsXL8mNED4y4vCPuvPVp2BRk8Aj0wIjrdeAmrz54pG6AEZe3sNjUlnWa3nNQkpbUobn2ssmm78+mIINHUTNgxBUpLCaJMgRDoNuiZozogRHXasTd7qbrB6+6gMBeIQR6ALkLi6FauBkLoClmy6SDQA+gKWbLpINAD6ApVqGmgxw9gJZYhZoGRvQAkDgCPQAkjkAPAIkjRw+MCDbVHl2M6IERUObmH6geAj0wAsrc/APVQ6AHRkDeUsRIC4EeGAGUMxhtHQO97R22j9uebWj7C9uP2n7Y9ldsn5u12/YnbB/IXr+0l50H0B3KGYy2bkb0d0iaWtb2sYh4XURcIulLkv4ka79a0ubsZ5uk20rqJ4AVoJzBaOs4vTIi9tjetKzthw2HL5W0uHvJVkl3RX03kwdsr7G9PiKOltRfYOSUNS2Scgajq3CO3vattg9JerdeGNFvkHSo4bS5rA1AAUyLRBkKB/qI+EhEbJT0GUnvz5rd7NRmv297m+0Z2zMnTpwo2g0gaUyLRBnKmHXzWUm/mT2fk7Sx4bVJSUea/VJETEdELSJqExMTJXQDSA/TIlGGQiUQbG+OiMezw2skfTt7vlPS+23fI+mXJM2TnweK27BmTPMLp87s1yp1nhZJqQMs1830yrsl3S/ptbbnbN8g6aO2Z20/Kuktkj6Qnf5lSd+TdEDSpyS9rzfdBkZD3mmR5PTRTDezbt7VpPn2FueGpBtX2imgKno9el6cFtn4GddeNtnyMxpz+pLOPO6aPcaofoRRvRIoaHH0PD62asnoud389CJfDHmmRR4+uaD146uXtJHTB4EeKCjv6LnTF0MZ/3dQJKeP9FHrBigo74yYdlMly8qtU+oAzRDogYLyFgpr98VQ1nx5Sh2gGVI3QEFTW9Zpes9BSfWA/fQzpzW/cErXXjbZ9Px2aZUyc+uUOsByjOiBgvKOntulVSgjjF5iRA+sQJ7Rc6epknn+7wDIg0AP9FGrL4a88+WBPAj0GDnDWiIgb259WK8Dw4ccPUZKKiUCUrkO9AeBHiMllbK/qVwH+oNAj5GSStnfVK4D/UGOHiOlXyUCep0/p9QB8mBEj8rbd3Re23fv102ff0Tbd+9vm6fuR4mAfuTPKXWAPAj0qLS8QbVoiYA8Xyb9yJ9T6gB5kLpBpRWpv15kGmOecsT9KhVMqQN0ixE9Kq0fNyXzjtApZ4BhQ6BHpfUjqOb9MiF/jmFDoEel9SOo5v0yIX+OYUOOHpXWjxoxecsRL/aLwI5hQaBH5fU6qFJwDFVHoAe6wAgdVUagR0+VuUK016tNqQaJVHEzFj1T5grRXq82pRokUkagR8+UuUK016tNqQaJlBHo0TNlLmbq9cIoqkEiZR1z9LZ3SHq7pOMRsSVr+5ikX5P0U0nflfTeiDiZvXazpBskPSfp9yPi33rUdwyRZvntMiss9rpaI9UgkbJuRvR3SJpa1rZb0paIeJ2k/ZJuliTbF0m6TtLF2e980vZZpfUWQ6lVfvuCdS9tuZgpT5EwqdjCqGGragkMSsdAHxF7JD21rO0rEbG4VPABSYsrR7ZKuicino2Ig5IOSLq8xP5iCLXKb+8/9uOmK0Qltbzx2So4511t2q+qlkAVlDG98nckfS57vkH1wL9oLmtDwtpVa2w2/3z77v1NK05++v4ntXDq+ZZVIvPMZe9HVUugKlZ0M9b2RySdlvSZxaYmp0WL391me8b2zIkTJ1bSDQxY3lowrW58PnToZGkzX7i5CrygcKC3fb3qN2nfHRGLwXxO0saG0yYlHWn2+xExHRG1iKhNTEwU7QaGQN78dqsvBsulBWdKBQMvKBTobU9J+iNJ10TETxpe2inpOtsvsX2+pM2SHlx5NzHM8ua3W30xvH7jeGnBmZurwAu6mV55t6QrJa21PSfpFtVn2bxE0m7bkvRARPxeRDxm+15J31I9pXNjRDzXq85jeOTJb7cqEiYpd5XIvJ9BDh6jyC9kXQanVqvFzMzMoLuBIUC9GaB7tvdGRK3TeRQ1w1Bh5gtQPkogAEDiCPQAkDgCPQAkjkAPAIkj0ANA4gj0AJA4Aj0AJI5ADwCJI9ADQOJYGTsCKCsAjDZG9InLu9MSgPQQ6BPXapu/Ipt5AKgmUjeJa7fNXxGkgYDqYUSfuDJ3WiINBFQTgT5xZe60RBoIqCYCfeLybvPXDhtuA9VEjn4E5N3Mo1UefsOaMc0vnNL42Koz57LhNjD8GNFjiXZ5eDbcBqqJQI8l2uXhy0wDAegfUjdYotN0TPZ0BaqHET2WKHM6JoDhwIh+hDW76Tq1ZZ2m9xyUVB/JP/3Mac0vnNK1l00OuLcAimJEP6Ja3XSVRB4eSAwj+hHVeNNV0pnHXbPH9MGrLiCwAwkh0CckTx2asmvgABheBPoBK6tI2GIqZnxs1ZJUTKu0C4ufgNHRMUdve4ft47ZnG9reafsx28/bri07/2bbB2x/x/Zbe9HpVJRZJCxvHRoWPwGjo5ubsXdImlrWNivpNyTtaWy0fZGk6yRdnP3OJ22ftfJupqnMImF569Cw+AkYHR1TNxGxx/amZW37JMn28tO3SronIp6VdND2AUmXS7q/jM6mpsw8eZFUDIufgNFQ9vTKDZIONRzPZW1ooszFSaRiALRSdqD/f0N8SdH0RHub7RnbMydOnCi5G9VQZnAmFQOglbJn3cxJ2thwPCnpSLMTI2Ja0rQk1Wq1pl8GqVsMzo2zbq69bLJjcG41U4dUDIBmyg70OyV91vbHJZ0rabOkB0v+jKQUqRWfZxolAHQM9LbvlnSlpLW25yTdIukpSX8raULSv9p+OCLeGhGP2b5X0rcknZZ0Y0Q817PeD9CgNslut6KVQA+gmW5m3byrxUtfaHH+rZJuXUmnht0gR9WsaAWQF0XNChjkJtmUEQaQFyUQCujHqLpVaogywgDyYkRfQK9H1e1KIzCNEkBejOgL6PWoutMNV6ZRAsiDEX0BvR5V561bAwDtMKIvqMj8926nY1JCGECZGNH3Qd5yxNStAVAmAn0f5J2OyQ1XAGUiddMHRaZjcsMVQFkY0fcBi5wADBKBvg/IuQMYJAJ9H5BzBzBI5Oj7hJw7gEFhRA8AiSPQA0DiCPQAkDgCPQAkLsmbsWVt8zeo7QIBoEzJBfp22/xJ6jpwswk3gFQkl7ppVVfm0/c/mauw2CC3CwSAMiUX6FvVcn/o0MlcgZua8ABSkVygb1VXxnKuwE19GgCpSC7Qt6or8/qN47kCN/VpAKQiuZuxi3VlGm+6Lu7lmmef11bv0+lGLDN1AAwbR8Sg+6BarRYzMzM9/5xeB+HGmTqNXybM1AHQC7b3RkSt03nJjejb6XVhscaZOpLOPO6aPUagBzAwHXP0tnfYPm57tqHtlbZ32348e3xF1m7bn7B9wPajti/tZeeHDTN1AAyjbm7G3iFpalnbhyXdFxGbJd2XHUvS1ZI2Zz/bJN1WTjd7a9/ReW3fvV83ff4Rbd+9v+Xc+k6YqQNgGHUM9BGxR9JTy5q3Sroze36npHc0tN8VdQ9IWmN7fVmd7YXFvHq3C6naYaYOgGFUdHrluog4KknZ4zlZ+wZJhxrOm8vahlaZK2DZSQrAMCr7ZqybtDWd1mN7m+rpHZ133nkld6N7h08uaP346iVtK8mrs5MUgGFTdER/bDElkz0ez9rnJG1sOG9S0pFmbxAR0xFRi4jaxMREwW6sHHl1AKkrGuh3Sro+e369pC82tL8nm31zhaT5xRTPsCKvDiB13UyvvFvS/ZJea3vO9g2SPirpKtuPS7oqO5akL0v6nqQDkj4l6X096XWJyKsDSN1IrYwFgJQkvzKWmjIA0J1KVq8sc+47AKSukoGe3Z8AoHuVDPTUlAGA7lUy0DP3HQC6V8lAz9x3AOheJQM9c98BoHuVnV5JTRkA6E4lR/QAgO4R6AEgcQR6AEgcgR4AEkegB4DEDUX1StsnJD1Z8NfXSvpBid2pklG9dq57tHDdrf1cRHTcuWkoAv1K2J7ppkxnikb12rnu0cJ1rxypGwBIHIEeABKXQqCfHnQHBmhUr53rHi1c9wpVPkcPAGgvhRE9AKCNSgd621O2v2P7gO0PD7o/vWJ7h+3jtmcb2l5pe7ftx7PHVwyyj71ge6Ptr9neZ/sx2x/I2pO+dturbT9o+5Hsuv8saz/f9tez6/6c7RcPuq+9YPss2w/Z/lJ2nPx1237C9jdtP2x7Jmsr7e+8soHe9lmS/l7S1ZIukvQu2xcNtlc9c4ekqWVtH5Z0X0RslnRfdpya05L+MCIulHSFpBuz/8apX/uzkt4UEb8o6RJJU7avkPSXkrZn1/2/km4YYB976QOS9jUcj8p1/0pEXNIwpbK0v/PKBnpJl0s6EBHfi4ifSrpH0tYB96knImKPpKeWNW+VdGf2/E5J7+hrp/ogIo5GxDey50+r/o9/gxK/9qj7UXa4KvsJSW+S9I9Ze3LXLUm2JyX9qqR/yI6tEbjuFkr7O69yoN8g6VDD8VzWNirWRcRRqR4QJZ0z4P70lO1Nkl4v6esagWvP0hcPSzouabek70o6GRGLe2im+vf+N5I+JOn57PhVGo3rDklfsb3X9rasrbS/88puPCLJTdqYQpQg2z8r6Z8k/UFE/LA+yEtbRDwn6RLbayR9QdKFzU7rb696y/bbJR2PiL22r1xsbnJqUtedeUNEHLF9jqTdtr9d5ptXeUQ/J2ljw/GkpCMD6ssgHLO9XpKyx+MD7k9P2F6lepD/TET8c9Y8EtcuSRFxUtJ/qH6PYo3txcFZin/vb5B0je0nVE/Fvkn1EX7q162IOJI9Hlf9i/1ylfh3XuVA/9+SNmd35F8s6TpJOwfcp37aKen67Pn1kr44wL70RJafvV3Svoj4eMNLSV+77YlsJC/bY5LerPr9ia9J+q3stOSuOyJujojJiNik+r/nf4+Idyvx67b9UtsvW3wu6S2SZlXi33mlF0zZfpvq3/hnSdoREbcOuEs9YftuSVeqXs3umKRbJP2LpHslnSfp+5LeGRHLb9hWmu1flvSfkr6pF3K2f6x6nj7Za7f9OtVvvp2l+mDs3oj4c9s/r/pI95WSHpL02xHx7OB62jtZ6uamiHh76tedXd8XssOzJX02Im61/SqV9Hde6UAPAOisyqkbAEAXCPQAkDgCPQAkjkAPAIkj0ANA4gj0AJA4Aj0AJI5ADwCJ+z+Nd7jzWatBEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(d_X, d_Y, alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### モデルの定義"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.reset_default_graph()\n",
    "\n",
    "## データの型定義\n",
    "X = tf.placeholder(tf.float32, shape=[None])\n",
    "Y = tf.placeholder(tf.float32, shape=[None])\n",
    "\n",
    "## 確率変数の型定義および初期化\n",
    "# 制約のないraw変数\n",
    "a = tf.Variable(100.0, dtype=tf.float32)\n",
    "b = tf.Variable(1.0,   dtype=tf.float32)\n",
    "sigma_raw = tf.Variable(0.2, dtype=tf.float32)\n",
    "\n",
    "# raw変数から推定に使うparametersを作る\n",
    "sigma = tf.exp(sigma_raw)\n",
    "mu = a + b*X\n",
    "\n",
    "# 計算しておくと後で便利なparameters\n",
    "mu_upper = mu + 1.96*sigma\n",
    "mu_lower = mu - 1.96*sigma\n",
    "\n",
    "## Costの定義\n",
    "cost_func = - tf.reduce_sum(tf.distributions.Normal(loc=mu, scale=sigma)\n",
    "                            .log_prob(Y))                            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 最適化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EPOCH: 100, cost: 114.066\n",
      "EPOCH: 200, cost: 114.062\n",
      "EPOCH: 300, cost: 114.062\n",
      "EPOCH: 400, cost: 114.062\n",
      "EPOCH: 500, cost: 114.062\n",
      "EPOCH: 600, cost: 114.062\n",
      "EPOCH: 700, cost: 114.062\n",
      "EPOCH: 800, cost: 114.062\n",
      "EPOCH: 900, cost: 114.062\n",
      "EPOCH: 1000, cost: 114.062\n",
      "a: 99.644905, b: 1.015577, sigma: 2.368601\n"
     ]
    }
   ],
   "source": [
    "n_epochs = 1000\n",
    "n_epochs_report = n_epochs // 10\n",
    "\n",
    "train = tf.train.AdamOptimizer(0.03).minimize(cost_func)\n",
    "sess = tf.Session()\n",
    "sess.run(tf.global_variables_initializer())\n",
    "\n",
    "for epoch in range(n_epochs):\n",
    "    cost_val, _ = sess.run([cost_func, train], feed_dict={X: d_X, Y: d_Y})\n",
    "    if (epoch+1) % n_epochs_report == 0:\n",
    "        print('EPOCH: %d, cost: %.3f' % (epoch+1, cost_val))\n",
    "\n",
    "a_est, b_est, sigma_est = sess.run([a, b, sigma])\n",
    "mu_est, mu_upper_est, mu_lower_est = sess.run([mu, mu_upper, mu_lower], feed_dict={X: d_X})\n",
    "print('a: %f, b: %f, sigma: %f' % (a_est, b_est, sigma_est))\n",
    "\n",
    "sess.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 結果の可視化"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x25a04022278>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzt3Xl02+d95/v3QxIgwQ3gLi6iSEmkJIpabFOy7Di2G2dxHcdO2qa1m7huuri5x8m0Z6bTNMnMpDP35J7O6b3tdCb39FRzmxsriS27WRo1ixsnk1pxEtmWHWvjLpGiuIgUF4AECBI/AM/88QMkkAJIrFzA7+scHRE/wsDvl1Bf/Ph8nuf7KK01QgghslfOep+AEEKIzJJCL4QQWU4KvRBCZDkp9EIIkeWk0AshRJaTQi+EEFlOCr0QQmQ5KfRCCJHlpNALIUSWy1vvEwCorKzUTU1N630aQgixqbz11luTWuuq1Z63IQp9U1MTZ8+eXe/TEEKITUUpdTWe58nQjRBCZDkp9EIIkeWk0AshRJaTQi+EEFlOCr0QQmQ5KfRCCJHlpNALIUSWk0IvhBBZbkMsmBJCiK2ga8zFyxfHGXF6qXfYeLi9hn219oy/r9zRCyHEGugac3H89AA3XAYluQW4vAbHTw/QNebK+HtLoRdCiDXwz78cx/BYmBm34HIq7DYLdpuFly+OZ/y9pdALIUQG+XzQ1we/7PWCLw+b7db3SgryGHF6M34OMkYvhBAZEAjA6Ch0d4PWUO+wMe83AMvN58wt+Kl32GK/SJrIHb0QQqSR1nDjBrz2Gly8CCUlUFEB92yvwWMYeAwDrTUur4HLa/Bwe03Gz0kKvRBCpInLBW+8AW++CXl5UF0NltANfJPDzmOtzRRZLEwtLGC3WXjm/uY1mXUjQzdCCJGi+Xno74fhYSgqgpoYN+lNDju1NrOw33vv2p2fFHohhEiSzwdXr8Lly7fu4JVa77O6nRR6IYRYxfKFTu/fV0Optt8MWsvLITd3vc8yNhmjF0KIFYQXOrm8BrX2Aq5PGXzxWwO8fMZ1M2jdyEUe4ij0SqkvK6UmlFIXI479hVJqRCn1TujPIxHf+6xSql8p1aOU+kCmTlwIIdbCyxfHsdssFORYGBtVzE1bKMm30Osdvxm0bnTxDN18BfgScGLZ8b/RWv/fkQeUUm3AE8B+oA74kVKqVWsdSMO5CiHEmhua9GLxFzAxB1YrFBdDUOcxOR97odOg08WZkXEm571UFto4Vl9DkyPzs2tiWfWOXmt9GpiO8/UeB05qrRe11gNAP3A0hfMTQoh14fNBby/4nDamZv0UFZmFHmDe8FNZGH2h06DTxaneATw+gwpbAR6fwaneAQadme9pE0sqY/SfUkqdDw3tlIWO1QPXIp4zHDomhBCbQiAA167Bq6/C4CA8uKsGf4650CmoNW6f+fWx+uhzKM+MjFNksVBstZCjFMVWC0UWC2dGMt/TJpZkC/3fAbuAw8AY8P+EjkebWKSjvYBS6hml1Fml1NkbN24keRpCCJEeWsPEBPz0p+aK1tJSM2jdWRFa6GS1MOVdoMhq4bHW5phDMZPzXgotS0fFCy23hnp8gQALfn/GrydSUtMrtdY3P5qUUv8T+G7o4TCwPeKpDcBojNc4DhwH6OjoiPphIIQQa8HpNHvSTE+D3W62LYjU5LDHPcZeWWjD4zMott5KaucNPxW2Al4bGuVrF7o5WruN99zfls5LWFFSd/RKqdqIhx8BwjNyTgFPKKXylVLNQAvwRmqnKIQQmTE/D+fPw89/DouL5orWgoLUXvNYvdnTxu27NdQzMufmjZHr/M3rv6Qwz8Lh6ur0XECcVr2jV0q9ADwIVCqlhoEvAA8qpQ5jDssMAn8EoLW+pJR6CegE/MCzMuNGCLHR+Hzm+Pvly2bAms4VreGeNmdGxhlyzTI862ZkzkO5rYBPHTnEsW315Kzx8tlVC73W+skoh/9hhed/EfhiKiclhBCZsLx1cKqLnWJNoyy3mcM3Z0cnsOTm8GR7Kx9q2Ul+Xi6Li+m7nnhJCwQhRNYLtw7u7ISFBSgrM3vTpCI8jbLIYrk5jfKfeq5Qmm/lJwPDLPj9vHdnI7+1vxVHQX56LiRJUuiFEFlttaA1WZHTKLXWTHm9nBufZMEf4K7aap46uJftpWl6sxRJoRdCZCWPx9zCb3R05dbByZqc91JhK+CGx8u58RtMLyziyLeyp6KMz913JL1vliIp9EKIrOLzwcAAXLmS/qA1Un5uLj8dGmXcM48tL4+jdTVU2Aoozrem/81SJIVeCJEVAgEYGYGeHnNMvrIScjLQn9e1uMg/dvbxwytDKKC13EFbVTm+QBCPYfDendtXfY21JoVeCLGphVe0dnaac+HTEbRG4wsE+F7fIN/q6mchEOB9O7dzT/02Lk3O3Jx181Bzw7o2L4tFCr0QYtNyOqGrC2ZmzKC1tDT97xHUmteGRvn6xR4m5710hILWhlDQeqCmKv1vmmZS6IUQm06mg9awSzemeO5cF5dnXOx0lPKpIwc5UF2ZmTfLICn0QohNY62C1pE5N189382bo+NU2Ar49NFD3N+49ita00UKvRBiwwsEYHjY7A+fqaB10Oni1asjvDU2weicB2tuLh9r38MHW5vJ3+h7Ba5CCr0QYsOKDFoXFsxNuNOxonV52wIjEODv377ENdccAa3ZXlpCfWkRd9ZWbfoiD1LohRAbVKpBa7SCDixpW+Be9PE/377EoHOWhUCAuuIiDtZUUppvxe0zODMyviFn0SRKCr0QYkNJR9AarQ/Nqd4B8nNybrYtmPDMc258kpmFRaw5OdzfWMe24qKbrxG5WUi6aA1uNzgcaX3ZVUmhF0JsCIuLt4LW/PzUgtbIPjTAzb87b0zTWuHgtaFRRt0eCvPyOFJXw4TbQ7F16YrWlfaFTYbbbX6I1dVBa2vaXjYuUuiFEOtqedBaVZV60BruQxMpR8GU18srV2bJzcnhQHUFLeUOFvwBHPlWPIYBmHfy84Yfj2HwUHNDaieCmS24XOZd/L33rv3dPEihF0Ksk0wErWGR2/n5g0H6pp103ZgmoDW1JUXsqyynrCD/ZkF/rLUZYMmYfqqrXP1+s2NmQQHcdVfmpoLGQwq9EGLNzcyYQavLZYas6V7Reqy+hu/0XGFkzk3vlBOv39yz9ffvaKOmqDBmQU9H8BoMmtenFLS1QUNDapubpIMUeiHEmlketGZq61S3z8+VGRfDcx5KrBYe3FHPh1qb01rQl9Pa/ODy+WDnTmhqMrOGjUAKvRAi49IZtK5keHaOr57v5uzYBJWFNv746GHua6zL+IrWyKC1pcX8ENtIpNALITImE0FrNK6FRV7s7OWVK9coyM3l4wf28khLU8YXO4WD1rIyOHRofYLWeEihF0KkXSaD1kiLgQDf7R3g292XWQwE+MCuRj7a1oI9w2MmGylojYcUeiFEWmU6aAWzdfDpoRGev9DDlHeBupIitpcWU2y1MuNdyFihjwxa9++H+vr1D1rjIYVeCJEWHo85RDM2ltmg9cLEJM+d62LAOcv20mIaSiupLymm0JJ3cwXsYxHBazosD1qbm83umZuFFHohRErCQevAQGZbBw/PznHifDdvhYLWP7n7MMOzbuYN/20rYNPZo2ajB63xWLXQK6W+DDwKTGit25d970+BvwKqtNaTSikF/C3wCDAP/K7W+u30n7YQYr35/eYerZlsHQzgXFjkxUu9/GjgVtD6wZYmrLm5fOnNc7etgE1Xj5rNErTGI547+q8AXwJORB5USm0H3gcMRRz+VaAl9Odu4O9CfwshssRa7dG66A/wz31X+Hb3ZYxAkId37eA32nYvGX+PXAEblmqPGr/fHIfPz98cQWs8Vv2/R2t9WinVFOVbfwP8GfCdiGOPAye01ho4o5RyKKVqtdZj6ThZIcT6Wqug9dWrIzx/sYdp7wJH62p46uBe6kqKb3vusfoaTvUOAKn3qFm+onWzBK3xSOpzWCn1GDCitT6nln7U1QPXIh4Ph45JoRdiE/N4oKcHrl+H4uL0BK3R+sXPLhqcOG8GrbvK7PzJ3YfZX1UR8zWaHHYea21OqUdNZNC6a5e5onUzBa3xSLjQK6UKgc8D74/27SjHdIzXeQZ4BqCxsTHR0xBCrIHlQWu6NuFe3i/+utvDF197k2nvIlWhoPVd2+Nb0drksCcdvGZD0BqPZO7odwHNQPhuvgF4Wyl1FPMOfnvEcxuA0WgvorU+DhwH6OjoiPphIIRYH+GgtacnM0FruF98Xo7i7esTDMzMkpujOFhdyWfv68C6Ritay8s3f9Aaj4QLvdb6AnDzFzel1CDQEZp1cwr4lFLqJGYI65LxeSE2D61hfNwch89k0Dru9jDlXaBnaoZAULO73MHeijLchsHonPu2IZ10TZUMr2i12aCjw2zJsNmD1njEM73yBeBBoFIpNQx8QWv9DzGe/n3MqZX9mNMrP5Gm8xRCZFg4aHU6zTvcdOzRurxAB7Tm1avDvDk6gdfvp76kiIPVlZSE9mjNUyrqFoCpLoCKDFrb282hmmwJWuMRz6ybJ1f5flPE1xp4NvXTEkKsFbfbnAsfDlrTuUdrZIE+Pz55M2htLC2m3FZAXWhFq9tn4DGMJXu6QuoLoLZC0BoPWRkrxBaVzqA11h6tZ0bGyVE5fPV8F29fv0F1oY1/e+wO7m2o5apr9rbZMt/tG6TUsrQshRdAxfMbQyS3G+bnzWmSu3dDYWHy17fZSaEXYovx+2+1Dob0BK2x9mj92dAo3+zso8CSx1MH9/LI7qabQWu02TKxFkAlMqSzsACzs2a+cPgw2NO/x8imI4VeiC0ik0Hr8j1ae6Zm6J6cIag1j7Q08dF9LZTkrz5mEmsBVDxDOuGgtbDQXNG6VYLWeEihF2ILSDVoXU14j9Zrs3P0TTtZ8AeoLCzgDw7v50j9trhfJ9YCqJWGdAIB87q2atAaDyn0QmSxdASt8XAt+uibdjLmnqfUauWeploebWlKKkCNe0jH56dE2ZiaMsfgd+zYmkFrPKTQC5GFFhfN/VkHB83mXOks8JGhqCUnhyHXHN1TM0uC1sjWKImGqNEsH9JxzftxeQ0+dlcDDxzZ2kFrPKTQC5FFEg1aEy3C4WmUuUox6JxlwDlLXo7i0ZYmPn5gL5ZlYybxTLuMR3hI52dXxxmZ9lLnsPHMQw0caZWkNR5S6IXIAsGg2To4kaA1mSL82rUxRuc8DMy4CGpNS7mDHfYSbBbLbUUeVp52mUihNwwo9Nl5fJedffskaE2UFHohNrnpabPAz86aUwnjDVoTKcIBrfnXwWG+1zeALxCkoaSYAzUVlFitBLWOudFHtGmXiWwMIkFrekihF2KTWh60Jto6ON4ifG78Bs+d6+Kqa45yWwF7yh1st5fc/P5KG30kuzGI1maB9/vNFa0StKZGCr0Qm8ziIly+DFevpha0rlaEfzE8ytfO93DdM0+RJY+nDuzhUE0l/9w3iNtnxLXRRzIbg4RbBzc0yIrWdFFme5r11dHRoc+ePbvepyHEhrY8aC0rS21Fa+QYfWQRfnBHPS9fvsrrI+Pk5eTQVllGbUkxC34/j7U2AyQc4Mbz/HDr4IoK2LtXVrTGQyn1lta6Y9XnSaEXYmMLB62dnWZzrnSuaI0swo6CfDw+g9NXR1gMBNhhL+XQtkryQ4Pibp9BkdXCE/tb0/PmIYZhLugqLDS38KuslKA1XvEWehm6EWIDWx60pvsut8lhZ7u9lJ8MXuPkxV5mFha5p6EWrTUNpcVLdnhKJESNR2TQeuAA1NZK0JopUuiF2IBSDVrj9cvrNzhxvosh1xx7Khz86T13sreynJOXepMKUeMhQevak0IvxAaysGCuaE01aF3NVdcsJ8518874DWqKCvl3x+7knoZtN1e0JhOixkOC1vUhhV6IDSATrYOjmfYucPJSLz8ZuEahxcLvHtrHw7t23LbYKVZzsWR3eYoMWqV18NqTQi/EOgoGb7UOTnfQGsnr93Oq5wrf6blCIBjkg63N/Pq+3ZSsMGYSrblYoiKD1iNHJGhdL1LohVgnmQ5awVzR+pOBa7xwqRfnwiL3NGzj4wf2sq24KP1vFvm+AbPA5+ZK0LoRSKEXYo253dDTYwatpaUZDlrPdTE0O8eeijL+7N672FNRlpk3Cwnv0WoYErRuJFLohVgjy4PWbfHvx5GQQecsJ853cW58km1FhfzpPXeyrcjG66MTvHJlKOlWwauZmzP3aJWgdeORQi9Ehvn9cO2aGbQqldmg9YWLPfxkcJgii4VPHGrjA7t3MDI7l5ZWwbFEBq133CFB60YkhV6IDFnLoPU7PVc41XOFgNZ8qLWZX9/XsqQbZTpaBS9nGGbOUFwMR4+ahV6C1o1JCr0QGRAZtDocmQta/9fANU6GgtZ3ba/lt9v3sq146ZhJqq2Cb3vfiKD14EEJWjcDKfRCpNHcnBm0TkxASUlmglat9c0Vrddm3eypKOMz995Fa4ygNdlWwbe/rwStm9WqhV4p9WXgUWBCa90eOvZ/Ao8DQWAC+F2t9agyl9X9LfAIMB86/namTl6IjWKtVrQOOmd57lwX5yduBa3H6rct2aN1uXSscp2bA6/3VtBqS70TglhDq3avVErdD7iBExGFvlRrPRv6+t8AbVrrTyqlHgE+jVno7wb+Vmt992onId0rxWbl98PQEPT1mQGrw5GZoHXKu8DJcNBqtfDRthY+sGsHljjfLNkNusNBa2Wl2To43t2rxNpIW/dKrfVppVTTsmOzEQ+LgPCnxeOYHwgaOKOUciilarXWY3GfuRCbwFoFrd2TU3ztQi89UzMA3L+jnt87vH/JMEw8El3lKkFrdkn6R1Mp9UXgdwAX8Cuhw/XAtYinDYeO3VbolVLPAM8ANDY2JnsaQqy56WmzN/zcXAaD1mCQlzr7+E7PFYyguUfr7nI7Gpicn6fYmpk5jMuD1rq6zPyGItZW0oVea/154PNKqc8CnwK+AET7zI86NqS1Pg4cB3PoJtnzEGKtrFXQ+nYoaB2edVNWkM+d26qoCAWnbp+R8rTI6O9rtg42DGhpMYNWS2K/NIgNLB2/bD4PfA+z0A8D2yO+1wCMpuE9hFg34aB1cBAKCjIXtA44XTx3rosLE1NsKy6kvaqCPRUOciNuqdO9+QfcClq3bzdn00jQmn2SKvRKqRatdV/o4WNAd+jrU8CnlFInMcNYl4zPi81qedBaXZ2ZceqpeS/PX+zl1atm0Pp7h9t4/64dfLOrPzQt8lahX21aZCKhq9drzvOvqoI775SgNZvFM73yBeBBoFIpNYx55/6IUmoP5vTKq8AnQ0//PuaMm37M6ZWfyMA5C5FRa7ai1fDz7Z7L/HPvFYIaHtuzk1/fu5uiUNCa6LTIyM2+V2p1IEHr1hPPrJsnoxz+hxjP1cCzqZ6UEOtlasos8PEGrYlOWxx0uvj58HUuTExy1TnHYiDAfdvr+NiBPVQXLV3RmujmH6u1OpCgdeuSlbFCkFzQGu8ddNjAjJOvnOtm0DmL2zAoK8inraqcj+zdybxhcPJS720fGIlMi4zZ6sDjZWbGHIravVuC1q1ICr3Y0iJXtCYatCbSLOzKjIu/+vnbTMx7KbZauLehlvqSIjyGn5f7r+ILBlPuLhmt1YHT48embGzbJkHrViaFXmxJy4PWqqrEx6njaRYWGbTm5eRwuKaSXeUOckNvVmjJ4+yNKdqqKlLuLhk5pm8lD+e8n2CewSff20B7S2LXJrKLFHqxpYSD1s5OM5RMJWhdqVlYZNCqNTy+ZxcajREI3izy4eeDotCy9CSSmUbZ5LDzyM5mTl8ZZyrgZe8OGx850kBbnTSI3+qk0IstY3nQmmrnxWizYtw+H7a8XJ79wU9wLfp4d2Mdv91uBq3hMf3I53sMg9ZyO/OGP6XukuGgtTLPzuces1NbK0GruEUKvch6mVrRGjkr5oZnHiMYpH/ayemhUdoqy/ncffvYXe6I+vzIWTRA0t0lwyta/X5zRWtjowSt4nZS6EXWSiVojVeTw05Qw4nz5orWuuIiPnPvXRypq4naOjjWLJpEplGGhfdobWyUoFWsTAq9yDqGYe7RGitoTbZl73KT816ev9jDq1dHKLVa+YM79vO+nY3kJTFmksg0Sq8XOkdd9M2P47N6abLYyK2oYZ9NxuJFdFLoRdYIBuH6dXMc3u83x+GXB62Jzn2PZt4w+Hb3Zb7bO4AGPrxnF7+2bxdFGR4z8fnMYZqpgIt3FgaoLrewraAAl9fg+OkBnrm/mX21UuzF7aTQi6wQb9CaykbZgWCQVwau8eKlXmYXfdzfWMeT7bevaE23cNCal2euaH3x0jjVDgt2m3nu4b9fvjguhV5EJYVebGrhoHV83GxXsFrQmsxG2Vprzo5N8NXzXYzMeWirLOfpZUFrJsQKWkd/7qXWvvQaSgryGHGmt6ulyB5S6MWm5PWaQevQENzwubjgGmfy2upj7olulH1lxmwdfPHG6kFrpFRzgNlZ8xp37ICdO5cGrfUOGy6vcfNOHmBuwU+9Q9JYEZ3MtBWbimHA5cvw6qswNgYei4sfjw7g8RlLxtwHna6o//2x+ho8hoHbZxDUGrfPwGMYHKtfOiVnct7Lf3/jHf79j15jyDXHR/bu5J7ttbwxOs6LnX0xXx9u5QDxnlMkr9f87aS4GO67D/bvv302zcPtNbi8Bi6veQ3hrx9uz1CjfLHprbo5+FqQzcHFaoJBs7B3dy8NWk9e6r3tDt3tMyiyWnhif2vU11rpbnt50PpoazN3bavilYFrFFksS+a5xwpwkzmncNBaXAxtbVBevnJLhq4xFy9fHGfE6aXeYePh9hoZn9+C0rY5uBDrSetbe7S63bcHrcmMuUebyugPBvnRlSFe7OwLBa31/PaBPVQV2jh5qTehADeRc1oetMa7onVfrV0Ku4ibFHqxYUWuaC0tjR60Jjrmvlw4aD1xvovROQ/7q8p5+lAbu8puFdFEP0ziOSdZ0SrWkhR6seFEBq2rrWhNdBemSP3TTk6c7+LSjWnqSor483d10FFbfVvQmuiHyWrnNDtrrtptbLw9aBUiE6TQiw3DMG61Ds7Li691cKK7MAFMeOZ54WIPp4dGKbVa+cM79vPeFVa0JvphEuucavLtXL9u/mZy112yR6tYOxLGinW3PGgtKzO3u0s3j2Hwra7LfK9vAAV8sLWZj+yNb0VrKtMlw0FrSQns27d60CpEvCSMFRue1rdWtEYLWtPFHwzyypUhXrrUx6zPxwM76nmy3Qxa45VIL5qwZINWIdJNCr3IuGhTAeuL7PT0wI0bsYPW5RK9q9Za8+boOF+90B0zaE31PaK/r1ngAwFobYXt2yVoFetLhm5ERnWNuTh+egC7zUJJQR5Ot5+RGwbHyppprbZTUhLf60Q2I4tnLnv/tJPnznXROTlNfUkRTx3cFzVoTeU9ookMWnftMsNkITJFhm7EhvDyxXGzyOdbcDlhesqCCkCfd5y7SuK/U463GdmEZ57nL/bw06FRSvOt/OGd7byveTu5cYyZpNLwzOs1i3x1NXR0EPcHmBBrQQq9SJtoQzQjM16Kcwu4et0MXW02sKk8Jr2JNeBabS67GbT2872+QRTwa3t38ZG9uyhMYMwkmcVXkUHr3XdDRUX81yTEWlm10Culvgw8CkxordtDx/4K+BDgAy4Dn9BaO0Pf+yzw+0AA+Dda63/J0LmLDSRyiKbWbvZI/x+vDOCdy2Eq4Ke8xHJzJo3bl9h+qBB7Lnu5rYDv9w/y0qVe3D7jZtBaWWhLeLw9kfnyfr9Z4PPy4NAh2LZNglaxccVzR/8V4EvAiYhjrwCf1Vr7lVL/Ffgs8BmlVBvwBLAfqAN+pJRq1VoH0nvaYqMJD9HYbRYWF8EzY2HeBSo3QNBi4A1AYc7SOeixCnG048vnsnt8BkOzc4x75pmcX6C9qoKnD+1jZyhoTWaDkXjmy4eD1mDQXNEqQavYDFa9B9Fanwamlx37odbaH3p4Bgj/S3gcOKm1XtRaDwD9wNE0nq/YoEacXmx5eUxMmIuefD6otOdBruax1maKrBamvAsUWS081toMELXD48+vjUQ9Dtx8nQHnLG+OjXPpxjQFuXl87r4O/uKBu28WeVg63p6jFMVWC0UWC2dGxmNeQ3ih0/JzDX8wzM6as4Tq6uCBB8xVrVLkxWaQjjH63wNeDH1dj1n4w4ZDx0QWMwwo1Da6+w2K8y0UFZkLgsJDNNHmoMdqFPaD/iGaHKVRA9H3NDVw3e3hrbEJ7PlW/ujOdh6KEbQmM94O0efLz8+bfXckaBWbVUqFXin1ecAPfD18KMrTos7fVEo9AzwD0NjYmMppiHUSXtHa1QVNuTUM5g4QyAVNHh7fym0CYhXi6YUF2izlS47n5SheHx7jn7ovoxT8+r7dfHjPzhWD1lSbnYEErSJ7JF3olVJPY4a0D+lbk/GHge0RT2sARqP991rr48BxMOfRJ3seYu2FV7R2doLHY7YsOFRmx+6Iv+dMzHC1oIB5w0+x1UJQay7PuLg4MYURDPLgjgZ+u72VijiKdSrNzvx+cxzeaoXDh82mahK0is0sqUKvlHoY+AzwgNZ6PuJbp4DnlVJ/jRnGtgBvpHyWYsOYnSXmitZE2gTEKsS/uruRN0cnuO720DvtxO0zcBTk84lD+7ivMf5RwGSanUUGra2t5qKnPJmALLJAPNMrXwAeBCqVUsPAFzBn2eQDr4RWGp7RWn9Sa31JKfUS0Ik5pPOszLjJDl6vuYXf0BAUFq7cOjgesQqxEdRcm3Uz4Jyl0JLHfdvr+PCeZprLEt+IO5EPnvCK1vAerbKiVWQTaYEgVrS8dbDDkZnOi+Oeeb5+oYefXRvFkZ/Pb7W38lBTQ1wrWlMRGbTu2SNBq9hcpAWCSEkwCKOjZuvgQMBsrZuR1sE+g2909fP9/kFyFPzGvt18eM8ubJbM/mj6fOYwTWkpHDtmXp8Q2UoKvVgiWtCaibniRjDIv1y+yj929uHxGTzY1MCT++MLWlMRGbTecYcErWJrkEIvbnK5oLc3sdbBidJa8/rIdb56oZvr7nkOVJsrWpvoI7rFAAAUCklEQVQTbAWcqGDQnCopQavYiuRHXeD1Qn8/DA+bTcdSDVpj6Z2a4blzXXRPzdBQWszn7zvCHduqVmwdnA4ulxm0NjVJ0Cq2Jin0W5hhwOCgWeSt1vj2aE2GGbR287NrYzjy8/nkXQd4zxoHrUeOSNAqti4p9FvQ8qC1oiIzQavbZ/DNdQpanU4JWoUIk0K/hUQGrW63WQDXImj9laYGnmjfc1vLg3STFa1CRCeFfotwucwVrZOTYLdnZhxea82Zket87Xw31z3zHKyu5OlD+2hylKb/zSJI0CrEyuSfQ5YLB61nul1cdI4zG0h+0+uV9E7N8JVzXfRMzbC9tJj/8O4jHK5Zm6B1cVFWtAqxEin0WSoctF6+DNcXXLx2Y4Aia/ybcMTrutsMWn8+PIajIJ//464D/IqsaBViQ5FCn2UCAbN1cOSK1h93j1NkTW7T61jmfD6+2dXPD/oGyc3J4TfbWnhsz05sGR4zWVw07+IlaBUiflLos0Q4aL10yVzRGhm0JrsJRzRGIMDzF3v4Qf9VjGCQZkcpTx3Yw6FtGVhdFUGCViGSJ4U+C7hc5h381FT0oDUdm3CEg9b//51OprwLVNoKOLytCmtuLq8OjWIvyE/rmH9YMGgWeIC9e6GhQYJWIRIl/2Q2sfn5Wytai4piz6RJZRMOgJ6pGZ4710nPlJPSfCtHamtoLls6kyaVYaBYJGgVIj2k0G9CPh9cvWoGrRaLGUiuNLklmU04IHrQ2jU5fdtvAskOA8UiQasQ6SWFfhOJFrTGu6I1kU04wkHr9/sGzf/WXsIdtdXsLrcz6V1IeRgoFglahcgMKfSbgNbmQqfOTvNuN2OtgwMBXr58lW909uMxDGqKC2mrLKfcZu7jeqp3gI7aKs6O3QCSGwaKRoJWITJLCv0GtzxozVTr4F8MX+drF7oZ98xzuKaKqiIblpyc26ZkDs16khoGikaCViHWhvyz2qDiDVpT1T05zYnzXfRMOWm0l/Af3n2UO7ZV8aU3z2HPty55bngsPpFhoFgkaBVi7Uih32ASDVqTdd3t4WsXuvnF8HVKrVY6aqsptlromZqhrMCalimZ0Xg8ZkO1bdvMvjTFxaleiRBiNVLoN4jIoDUYzNwerXM+H9/o7Ofl/kHycnL4wM5GPIZBaX4+hZa8m+0R0j0WL0GrEOtHCv06W8ug9QeXr/KNzj68hp/3NG/nif2t/MuVIXIzOBYfDlrz8809WqurJWgVYq1JoV9HmQhaB52uJcX57rpqRt3zfO1CNxMeL4drqvidQ3vZYTcXPK3UHiGVsXgJWoXYOOSf3jqIFrQOOl2cuXyrQCfTRnjQ6eJU7wBFFrNL5TXXHC/3DzLnM9hhL+E/vvsoh7dVLflvMjEWHxm07tpl3s0LIdbPqr9EK6W+rJSaUEpdjDj2UaXUJaVUUCnVsez5n1VK9SulepRSH8jESW9WPh/09cGrr8LEhHkHX1x8q0B7fMaSNsKDTldCr39mZJwiiwWN5szwdc6MXMcXCNJRV81fve/dtxV5MNsjeAwDt88gqDVun4HHMDhWn/g0H48HxsfB4YB3vxva2qTIC7ERxHNH/xXgS8CJiGMXgV8D/j7yiUqpNuAJYD9QB/xIKdWqtQ6k5Ww3qUDg1h6tWt++R2u4QKfaRnhszs119zxXZlzkKMX+qnJayh24Fn3kxpi6k2x7hEiRQes995g5gxBi41i10GutTyulmpYd6wKi7R70OHBSa70IDCil+oGjwC/ScbKbTbxBa6pthH2BAD/oH+T1kXGMYJCdjlL2V1Vgs+Th9hmrDsMkOxa/PGitqcnMVFAhRGrSPUZfD5yJeDwcOrblJBK0JjtOrrXmZ9fG+PqFbibmveytKKPYaqG6qJD8vNybwzCptCeIRoJWITaXdP/zjHY/p6M+UalngGcAGhsb03wa6yeZFa3JtBHumpzmuXNd9E07abKX8p/uP8ChmqrbZt0k254gFqfTzBqam80/MgYvxMaX7kI/DGyPeNwAjEZ7otb6OHAcoKOjI+qHwWbi893ao9VqTWxFayLj5KNz5orW10euU16Qz7MdB3mgqeHmGHw62hNEIytahdi80l3oTwHPK6X+GjOMbQHeSPN7bCirBa3xWq1Azy36eKmzj3+5fBVLTg4P72rEZrHQOTnNxLw3qemY8VhcNO/iHQ4JWoXYrFYt9EqpF4AHgUql1DDwBWAa+B9AFfA9pdQ7WusPaK0vKaVeAjoBP/Bsts640Rpu3ICXz7h4fWQct/ZSXWzjWG56C244aP1GVz8Lhp+HdjZyb8M2/vXqCDkqZ8l0zMdam9P23pFB6513StAqxGamtF7/UZOOjg599uzZ9T6NuDmd5h38hWsuXh0foLTAsmRsPR0Fd3nQese2Kn7n4D4a7SWcvNR7W3jr9hkUWS08sb81pfeNDFpbWyVoFWIjU0q9pbXuWO158k84AeGgdWQECguh1ztOaUHq89+X67xhtg6+FbQe5FBN5c3vpzodMxYJWoXITlLo47A8aK2qMocx0l1wzaC1i9dHxim3FfCpI4fYYS/mjdEJfjo0crM1QrrbFkjQKkR2k0K/guVBa2Xl0s6L6Sq4s4s+/jEctObm8GR7Kx9q2cmY272kd026WwiHV7Ta7RK0CpHNpNBHEQ5aOzvB6zV7p0cbp05m/nskXyDA9/sH+WYoaH3vzkZ+a38rjgJzzCRWa4RUWwjLilYhthYp9MuEg9aZGbN3S0lJ7Ocm2ycmqDU/uzbK1y/0cGPey1211Tx1cC/bS5e+WbpbCAeDMD1t/laybx/U10vQKsRWIP/MQzwes7Pk6Ki5ojXe3vCJFtwfDwzxwsVeZhYWcRTk80d3tvP+XTuiPjedY/EStAqxdW35Qu/zwcAAXLmS+IrWRIzOufn7ty5w8cY0Bbm5dNRWU1Voo3tqhtYKR9QPi1SHhmBp0Lpnj/khJoTYWrZsoQ8EzGmSPT3Rg9Z0cS0u8o+dffzw8hBKQWu5g/bqCvJCb6aUijkdM5UWwhK0CiHCtlyh19rc9KOz0yyGZWWZGaf2BQJ8r2+Qb3X1sxAI8N7m7XgMg9riInIifmVYbTpmokNDErQKIZbbUoXe6YSuLrMQ2u1m2JpuQa15bWiUr1/sYXJZ0BptRWuq2/bdfN/QilalJGgVQiy1JUrB8qA1ntbBybh0Y4rnznVxecZFs6OUTx05yIHqWyta0zHmvpzWt/Zo3blTglYhxO2yutCvZdD61fPdvDE6ToWtgE8fPcT9jfVLhmggPdv2RXK7zQ+xujpoaZGgVQgRXVYW+kDA3Pijt3ftglZrbi4fa9/DB1ubyV+hT3E6+sUvLMDsrAStQoj4ZFWhjwxaFxZir2hN1WIgwPcjgtb372zkN9tasBdkdszE7zcXPBUUmK2DM/UbihAiu2RVoR8chEuXzDvctQhaj9TV8PEDe2kozWwXsMigta3NbB2czOYmQoitKasK/eKi2T64INQ1YPn+qanswnRxYooT582gdaejlE8fOUR7dUUaz/52ErQKIdIhqwp9pEGnK2rnx0Q3BRmedfO1C928OTqOLS+PfZVl7K+qoNia2f/pJGgVQqRL1hb6WJ0fw6tQV7vbdy0s8lJnHz+8MoQlJ4dmRymt5Q5K8q3MG/60b90XtrBg3sU7HHDvvebfQgiRiqwt9Ct1flzpbr+2pJjv9Q3wra7LLIaCVktuDoGgTvtOUpEig9a77pKgVQiRPllb6Ffq/Bjtbl9rzYuX+rgy42LSu0BHbTVPHdxHQ2kxX3rzXEa27gMJWoUQmZe1hX6lVajf7RtcUrgnPPO8Mz6Jc2GRXWV2Pn308JKgNd1b98GtoNXnuxW0Wq1Jv5wQQsSUtYV+pVWo4cId1Jrz45OMuj0U5OVytL6Gf3/PXbetaE2mdcFKGYAErUKItZS1hR5ir0LdX1XG8bcuMeb2kJuTw56KMqoKbXxk787binz4dRJpXRArA3h4RzN27JSVwaFDErQKIdZGVhf65Rb9Ab7bN8C3uy+zGPCzs8zOtuJC6kqKV51jn0jrguUZQGGeBe88/GJknM8/ZpegVQixprZEoQ9qzemrIzx/sYcp7wJH6mp46uBe6ksys6I1PONHa3NzcaVge20eHu3NWOdMIYSIZdVCr5T6MvAoMKG1bg8dKwdeBJqAQeA3tdYzSikF/C3wCDAP/K7W+u3MnHp8LkxM8ty5Lgacs+wqs/PHdx9mf1XsFa3pWE1babMx7TYoyLFQVmYO0bh9fhpsqfedF0KIRMXT0/ErwMPLjv058GOtdQvw49BjgF8FWkJ/ngH+Lj2nmbjh2Tn+r9fe5C9efZ05n8Gf3H2Yv3zoXasW+VO9A3h8xpKx9UGnK+73dbuhxVZDINegrNqgrFzj9hm4vAYPt8vtvBBi7a16R6+1Pq2Ualp2+HHgwdDXzwH/CnwmdPyE1loDZ5RSDqVUrdZ6LF0nvBrX4iIn3+rlRwPXKMjN5akDe3mkpQlrHJPTV1tNu5LwitayMvjN99t5t7eZly+OM+L0Uu+w8VtHGthXm95VtEIIEY9kx+hrwsVbaz2mlKoOHa8HrkU8bzh07LZCr5R6BvOun8bGxiRP4xavL8AL565w8vxl/MEgH9jVyEfbWrAn0AVspdW0sUTu0Rq5otXhsEthF0JsCOkOY6PNJdHRnqi1Pg4cB+jo6Ij6nER8/fWrPPfLXu6oruETdyYXtCayMGr5itb6elnRKoTYmJIt9OPhIRmlVC0wETo+DGyPeF4DMJrKCcbr48d24Ag6KA+UYy9J7jXiWRgVuaJ11y5oapIVrUKIjS3ZDfZOAU+Hvn4a+E7E8d9RpmOAa63G5wssubTXlKf0GuGFUUVWC1PeBYqsliUdKt1ucwerigq4/35obZUiL4TY+OKZXvkCZvBaqZQaBr4A/CXwklLq94Eh4KOhp38fc2plP+b0yk9k4JwzKtrCqHDQWl4uK1qFEJtPPLNunozxrYeiPFcDz6Z6UonqGnPx8sVxOq96sQVtPLAr+Z2kIkW2Du7ogKoqWdEqhNh8kh262TC6xlwcPz2Ay2tQWViAx0h87vtywSBMTcHsLLS3m8M00rZACLFZbfoWCC9fHMdus2C3WZj0QJHFgj8nuU1BJGgVQmSjTV/oR5xeau2pbwridsP8vDlNcvduc5NxIYTIBpu+0Nc7bLi8BnZbcpuCLCyYQzTSOlgIka02/Rj9w+01uLxmL5mg1ngMA49hcKx+5b4yhmFOlfT7zRWtd98tRV4IkZ02/R39vlo7z9xv9pXpnPBSZLHxwK7Ym4IEAuB0msFqe7u5y5OsaBVCZLNNX+jBLPb7au10d8PoKNij1HitzQJvGOYY/I4dErQKIbaGrCj0qwnv0drQIEGrEGLryepCLytahRAiSwu9YZidJQsLZUWrEEJkXaF3uczAVYJWIYQwZVWhLy42C7wErUIIcUtWFfqGhtWfI4QQW82mXzAlhBBiZVLohRAiy0mhF0KILCeFXgghspwUeiGEyHJS6IUQIstJoRdCiCwnhV4IIbKc0lqv9zmglLoBXE3DS1UCk2l4nc1ErnlrkGveGhK95h1a66rVnrQhCn26KKXOaq071vs81pJc89Yg17w1ZOqaZehGCCGynBR6IYTIctlW6I+v9wmsA7nmrUGueWvIyDVn1Ri9EEKI22XbHb0QQohlsqbQK6UeVkr1KKX6lVJ/vt7nkwlKqS8rpSaUUhcjjpUrpV5RSvWF/i5bz3NMN6XUdqXUT5RSXUqpS0qpPw4dz9rrVkoVKKXeUEqdC13zfw4db1ZKvR665heVUlm1vY5SKlcp9Uul1HdDj7P9egeVUheUUu8opc6GjmXk5zorCr1SKhf4f4FfBdqAJ5VSbet7VhnxFeDhZcf+HPix1roF+HHocTbxA/9Oa70POAY8G/r/NpuvexF4j9b6EHAYeFgpdQz4r8DfhK55Bvj9dTzHTPhjoCvicbZfL8CvaK0PR0ypzMjPdVYUeuAo0K+1vqK19gEngcfX+ZzSTmt9Gphedvhx4LnQ188BH17Tk8owrfWY1vrt0NdzmIWgniy+bm1yhx5aQn808B7gG6HjWXXNSqkG4IPA/xd6rMji611BRn6us6XQ1wPXIh4Ph45tBTVa6zEwiyJQvc7nkzFKqSbgDuB1svy6Q8MY7wATwCvAZcCptfaHnpJtP+P/DfgzIBh6XEF2Xy+YH94/VEq9pZR6JnQsIz/X2bJnrIpyTKYTZRGlVDHwTeBPtNaz5g1f9tJaB4DDSikH8G1gX7Snre1ZZYZS6lFgQmv9llLqwfDhKE/NiuuN8C6t9ahSqhp4RSnVnak3ypY7+mFge8TjBmB0nc5lrY0rpWoBQn9PrPP5pJ1SyoJZ5L+utf5W6HDWXzeA1toJ/CtmPuFQSoVvzrLpZ/xdwGNKqUHMYdf3YN7hZ+v1AqC1Hg39PYH5YX6UDP1cZ0uhfxNoCaX0VuAJ4NQ6n9NaOQU8Hfr6aeA763guaRcaq/0HoEtr/dcR38ra61ZKVYXu5FFK2YD3YmYTPwF+I/S0rLlmrfVntdYNWusmzH+7/0tr/TGy9HoBlFJFSqmS8NfA+4GLZOjnOmsWTCmlHsG8C8gFvqy1/uI6n1LaKaVeAB7E7HA3DnwB+CfgJaARGAI+qrVeHthuWkqp+4CfAhe4NX77Ocxx+qy8bqXUQcwgLhfzZuwlrfV/UUrtxLzjLQd+CXxca724fmeafqGhmz/VWj+azdcburZvhx7mAc9rrb+olKogAz/XWVPohRBCRJctQzdCCCFikEIvhBBZTgq9EEJkOSn0QgiR5aTQCyFElpNCL4QQWU4KvRBCZDkp9EIIkeX+NyBb1bk9fpVEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.fill_between(d_X, mu_lower_est, mu_upper_est, color='b', alpha=0.2)\n",
    "plt.plot(d_X, mu_est)\n",
    "plt.scatter(d_X, d_Y, alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}