count_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Count regression\n",
    "\n",
    "This notebook demonstrates the use of tfp for regression on count data. In good statistical tradition, we do not use a testset. We compare again the R-Code and use a standard dataset.\n",
    "\n",
    "* Poisson Regression\n",
    "* Zero Infated Poisson Regression\n",
    "* NB work in progress\n",
    "\n",
    "#### Links:\n",
    "Using plain TF\n",
    "https://github.com/gokceneraslan/neuralnet_countmodels/blob/master/Count%20models%20with%20neuralnets.ipynb\n",
    "\n",
    "\n",
    "."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TFB Version 0.7.0-dev20190511\n",
      "TF  Version 2.0.0-alpha0\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "\n",
    "import tensorflow as tf\n",
    "import tensorflow_probability as tfp\n",
    "\n",
    "sns.reset_defaults()\n",
    "#sns.set_style('whitegrid')\n",
    "sns.set_context(context='talk',font_scale=0.7)\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "tfd = tfp.distributions\n",
    "print(\"TFB Version\", tfp.__version__)\n",
    "print(\"TF  Version\",tf.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading some count data\n",
    "\n",
    "We use some traditional count data from: https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta see also: https://stats.idre.ucla.edu/r/dae/negative-binomial-regression/ The data have been processed:\n",
    "\n",
    "\n",
    "```\n",
    "dat <- read.dta(\"https://stats.idre.ucla.edu/stat/stata/dae/nb_data.dta\")\n",
    "dat <- within(dat, {\n",
    "  prog <- factor(prog, levels = 1:3, labels = c(\"General\", \"Academic\", \"Vocational\"))\n",
    "  id <- factor(id)\n",
    "})\n",
    "summary(dat)\n",
    "m.p <- glm(daysabs ~ math + prog, data = dat,family = poisson)\n",
    "\n",
    "df = data.frame(model.matrix(m.p))\n",
    "df$daysabs = dat$daysabs\n",
    "df$X.Intercept. = NULL\n",
    "write.table(df, \n",
    "            file='~/Documents/workspace/dl_book_playground/nb_data.csv', \n",
    "            col.names =FALSE,row.names = FALSE)\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = np.loadtxt('nb_data.csv')\n",
    "X = d[...,0:3]\n",
    "y = d[...,4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((314,), (314, 3), array([[63.,  1.,  0.],\n",
       "        [27.,  1.,  0.],\n",
       "        [20.,  1.,  0.],\n",
       "        [16.,  1.,  0.],\n",
       "        [ 2.,  1.,  0.]]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape, X.shape, X[0:5,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Definition of the Model\n",
    "\n",
    "$$\n",
    "    Y \\thicksim \\tt{Pois}(exp(w^{T} \\cdot x + b))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Definition of the custom parametrized distribution\n",
    "def my_poiss(t):\n",
    "    return tfd.Poisson(rate=t) #Equivalent to log_rate=t when using linear activation \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: id=0, shape=(4,), dtype=float64, numpy=array([0., 1., 2., 3.])>,\n",
       " <tf.Tensor: id=6, shape=(4,), dtype=float64, numpy=array([0., 1., 2., 3.])>)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_poiss(np.array((0.0,1.0,2.0,3.0))).mean(), my_poiss(np.array((0.0,1.0,2.0,3.0))).stddev()**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Definition of the model\n",
    "#Using the functional API for clarity\n",
    "inputs = tf.keras.layers.Input(shape=(3,)) #we need to estimate 3 (2 cat + 1 numerical) coeficient (excluding intercept) \n",
    "rate = tf.keras.layers.Dense(1, \n",
    "                             activation=tf.exp, #We use the exponential \n",
    "                             bias_initializer='zeros', \n",
    "                             kernel_initializer='zeros')(inputs)\n",
    "\n",
    "#p_y = tfp.layers.DistributionLambda(my_poiss)(rate) #<------- Here we hand-in the the custom distribution\n",
    "p_y = tfp.layers.DistributionLambda(tfd.Poisson)(rate) \n",
    "\n",
    "#Glueing input and output together. Note that output is a tf.distirbution \n",
    "model = tf.keras.models.Model(inputs=inputs, outputs=p_y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def negloglik(y_true, y_hat):\n",
    "  return -y_hat.log_prob(y_true)\n",
    "\n",
    "model.compile(optimizer=tf.optimizers.Adam(learning_rate=0.01), loss=negloglik)\n",
    "hist = model.fit(X, y, epochs=600, verbose=False, validation_split=0.0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.422600928036271"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEbCAYAAAA1T5h7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XecVOW9x/HPbyuwwC7SpCmIiF2JJfZuvCaxa9QYE5NobnrxGq8lRo3XEomJsSTGxFhijz2KIghWQCyIdBaQXhfY3nef+8eZmT0zO30XZnf2+3695gV7zjMzz9mZPb/ze9ox5xwiIiIAOZmugIiIdB0KCiIiEqKgICIiIQoKIiISoqAgIiIhCgoiIhKioCAiIiEKCiIiEqKgICIiIQoKIiISoqAgIiIhCgoiIhKioCAiIiEKCiIiEpKX6QpEY2bNeAGrMtN1ERHpRvoDrc65tM/t1hXvp2BmrYAVFxdnuioiIt1GRUUFgHPOpd0K1CUzBaCyuLi4uLy8PNP1EBHpNkpKSqioqOhQC4v6FEREJERBQUREQhQUREQkREFBRERCFBRERCREQUFEREK66pDUtE1esJHtNY0cOnoAew7pl+nqiIh0K1mXKfz+jcVc88I8Zq7YlumqiIh0O0kFBTP7kpnNMrNKM1tqZufFKDfOzGaaWW2g/LjOrW5ivfJyAWhoatnZby0i0u0lmyk8BrwIlABXAI+Z2aAo5Z4B/gPsArwKPN0ZlUxFr3zvkBqaW3f2W4uIdHvJBoXdgWedc63OuXeArcAYfwEz2xsYC0x0ztUDdwLjzGx8Z1Y4kcJAplCvTEFEJGXJBoX7gG+aWZ6ZnRzYtiCizL7AEudcE4BzrhFYEtgexszK4z2AtFfCU6YgIpK+ZIPC63jNRvXAJODnzrnaiDJFtF/quhLo26EapkiZgohI+hIOSTWzXfD6CS4DXgH2AyaZ2XLn3Dxf0Rogcgxof6A68jWdcyUJ3jPtbCGYKSgoiIikLplMYSxQ6Zx70TnX4pz7HJgFnBBRbiEw3szyAMysANgrsH2nCWYKaj4SEUldMkFhKdDPzM4wz/7AcYA/S8A5txhYAVxlZoXAVcAy59ySzq50PMoURETSlzAoOOcqgAuBW/D6CF4D7nTOvW1m15nZA77iFwFnA+XAWcDFnV/l+HrlK1MQEUlXUstcOOcmA5OjbL8t4uelwBGdU7X0FOYpUxARSVfWLXNRmB8cfaRMQUQkVdkXFPI0T0FEJF1ZFxRCfQpqPhIRSVnWBgX1KYiIpC7rgkKoo1nNRyIiKcu6oJCfawA0tygoiIikKuuCQo55QaHVZbgiIiLdUNYFhbxgptCqTEFEJFVZFxRCmYJigohIyrIuKOTmeEGhxan9SEQkVdkXFAKZQos6FUREUpZ1QSEnkCkAtCowiIikJOuCQp4vKKgJSUQkNVkXFPyZgpqQRERSk3VBIdinAAoKIiKpShgUzKw64uHM7LwYZd82s3pf2dc7v8rx5ar5SEQkbQlvsuOc6xv8v5kdiHd/5jfjPOUy59zTnVC3tOSYOppFRNKVavPRJcDLzrmqHVGZzhCc0QxqPhIRSVXSQcHMDO+ey08kKHqvmW0xsymBzCLaa5XHewDFyR9CuBz1KYiIpC2VTOE4oA9R7tXsczUwBtgNmAJMMrO+ccp3OvUpiIikL5WgcAnwrHOuKVYB59xs51y1c67OOXcnUA0cHqVcSbwHUJHykQRo9JGISPoSdjQDmFkBcD5wRoqvv9OXpcvxhTktiicikppkM4Wv4l29z4hVwMxKzOxUMys0swIz+xUwAPioE+qZNDUfiYikL9mgcAnwpHPhZ1kzu8Q3FyEfuB3YCmwEvg6cvrNHKoUFBaUKIiIpSar5yDl3QYztTxAYjeSc2wIc2nlVS094n0IGKyIi0g1l3zIXWvtIRCRtWRcUwpbOVp+CiEhKsi4oaEiqiEj6si8o+DKFZgUFEZGUZHVQUPORiEhqsi8oqPlIRCRtWRcUdI9mEZH0ZV1QgLYmJM1oFhFJTVYHBXU0i4ikJjuDQqBfQc1HIiKpyc6gEGw+UlAQEUlJVgaFYF+zhqSKiKQmK4NCW6aQ4YqIiHQzWRoUvMNq1tLZIiIpydKg4P2r5iMRkdRkZ1AwNR+JiKQjYVAws+qIhzOz82KUHWJmb5hZrZnNN7MjOr/KiQVnNWtIqohIahIGBedc3+ADOAqoA96MUfwBYAUwEJgIPG9mhZ1V2WRpRrOISHpSbT66BHg52n2Xzawf3n2Zb3bO1TnnHgWqgBM6XMsUaUaziEh6krpHM4CZGXAx8OMYRcYB5c65Tb5t84B9gckRr1We4O2Kk61XNJrRLCKSnlQyheOAPkSc4H2KgMqIbZVA3zTq1SHKFERE0pN0poDXdPSsc64pxv4aoF/Etv5AdWRB51xJvDcKZBJpZwv5gTGpLZqnICKSkqQyBTMrAM4HnohTrBQYYGZDfdv2BxamX7305Od6mUJTizIFEZFUJNt89FWgApgRq0Cg8/lV4AYz62Vml+JlCm93tJKpCmYKjc3KFEREUpFsULgEeNK58DGeZnaJmb3u2/QjvA7nbcA1wHnOuYZOqWkKgkGhSbPXRERSklSfgnPughjbn8DXpBQYeXRa51QtfW3NRwoKIiKpyMplLtoyBfUpiIikIjuDQp6aj0RE0pGdQSFHzUciIunIzqCg5iMRkbRkZ1BQ85GISFqyMigUaEiqiEhasjIoaEaziEh6sjIo5ClTEBFJS1YGBS1zISKSnqwMCgW5WjpbRCQdWRkUtPaRiEh6sjIo5Kn5SEQkLVkZFAq0IJ6ISFqyMigEm4/UpyAikpqsDgpNaj4SEUlJVgaFvEDz0fqKeiLuCyQiInEke49mM7PfmNl6M6sys5kxyq00s1ozqw48Hujc6iYnuMwFwKR5GzNRBRGRbimpO68BPwOOBw4DNgAT4pQ9yTk3q6MV64i9h/UP/X/l1poM1kREpHtJGBTMLBe4FjjKObcusPmTjrypmZUnKFLckdcfM6iIg0YWM3dtBXWNLR15KRGRHiWZ5qNRQC/gYjPbZGYLzewbccq/FCj3opnt3jnVTN3gfoUA1DUpKIiIJCuZ5qPhQAmwe+BxEDDZzOY655ZElP0m8ClQANyCFyC+5CJ6e51zJfHeMJBJdChb6JWfCygoiIikIplMoS7w763OuXrn3IfAm8CpkQWdczMCZSqBK4FxwJhOq20KegeCQr2aj0REkpZMUFgKNKXx2i7wyIjeBcoURERSlTAoOOdqgBeAa8yswMwOwcsSpvrLmdluZnakmeWbWREwEfgi8NjpQpmCgoKISNKSnbz2E2AksA14Gvhv59xiM7vONxehH/AgUA6sBMYCZ0X2J+ws6lMQEUldUvMUnHNbgTOjbL/N9/8FwAGdV7WOaWs+0lIXIiLJysplLgB65XmHpo5mEZHkZW1QUEeziEjqsjYoqE9BRCR1WRsUNE9BRCR12RsU1HwkIpKy7A0KgUyhudXptpwiIknK2qAQ7FMAZQsiIsnK2qAQbD4C9SuIiCQre4OCMgURkZQpKIiISEj2BgVf85HuviYikpysDQqFeW2HpkxBRCQ5WRsUzEzLZ4uIpChrgwL4JrA1ap6CiEgysjsoaP0jEZGUJBUUzPMbM1tvZlVmNjNGuXFmNtPMas1slpmN69zqpqZXvnd4CgoiIslJNlP4GXA8cBhQDPw0RrlngP8AuwCv4t2lLWOCzUe1Dc2ZrIaISLeRMCiYWS5wLfAD59w651yrc+6TKOX2xrsF50TnXD1wJzDOzMZ3dqWTNay4NwArttRkqgoiIt1KMpnCKKAXcLGZbTKzhWb2jSjl9gWWOOeaAJxzjcCSwPYwZlYe74GXjXTYQSO9l5m7trwzXk5EJOslc4/m4UAJsHvgcRAw2czmOueW+MoVAZURz60E+nZGRdOx7/D+ACzbXI1zDjPLVFVERLqFZIJCXeDfWwPNQh+a2ZvAqXiZQFAN0C/iuf2B6sgXdM6VxHvDzsoWinsXAN7y2Q3NrWErp4qISHvJNB8tBZqSKLcQGG9meQBmVgDsFdieEUWFbUGgVktdiIgklDAoOOdqgBeAa8yswMwOwcsSpkaUWwysAK4ys0LgKmBZRBPTTlVU0JYI1WgEkohIQskOSf0JMBLYhjfM9L+dc4vN7Doze8BX7iLgbKAcOAu4uDMrm6qiQl9QaFRQEBFJJJk+BZxzW4Ezo2y/LeLnpcARnVO1jvM3HylTEBFJLKuXuSjIzSEvxxtxVNOgPgURkUSyOiiYGX2Cs5rVfCQiklBWBwWAvoF+hWplCiIiCWV9UOgTCArKFEREEsv6oFAUyhQUFEREEsn6oNC/lxcUquoVFEREEukBQSEfgMq6ZCZli4j0bNkfFHp7mUKlMgURkYSyPygEMoUKZQoiIgllf1DoreYjEZFk9ZygUK+gICKSSPYHhcDoo8o69SmIiCSS/UHB13zknMtwbUREurasDwqD+xYC0NjSqs5mEZEEsj4o7FrcK/T/DRX1GayJiEjXl/VBYWBRAQW53mFuVFAQEYkrqaBgZm+bWb2ZVQcer3ek3M5kZqFsYeGGygzXRkSka0slU7jMOdc38Di9E8rtNHsMLgLg8VmraG1VZ7OISCwZaT4ys/J4D6C4M9/vW1/eHfD6FMqqGzrzpUVEskoqQeFeM9tiZlPM7MBOKLfTHDl2YOj/a7bXZbAmIiJdW7JB4WpgDLAbMAWYZGZ90y3nnCuJ9wAq0jqaGIoK89ilqACAtdtrO/OlRUSySlJBwTk32zlX7Zyrc87dCVQDh6dbLhNGDugNwFplCiIiMaXbp9DayeV2uKH9vRFImys1LFVEJJaEQcHMSszsVDMrNLMCM/sVMAD4KJ1ymTIoMLO5rLoxwzUREem6kskU8oHbga3ARuDrwOnOuSozu8Q3FyFmuc6vduoG9/OCwhaNPhIRiSkvUQHn3Bbg0Bj7ngCeSFSuKxjc1+to1pBUEZHYsn6Zi6Bg89GWKgUFEZFYek5QCDQfVdU3U9/UkuHaiIh0TT0nKAQyBYCtNepsFhGJpgcFhYLQ/8vUhCQiElWPCQp9C/MozPMOV53NIiLR9ZigYGa+uQoKCiIi0fSYoABtnc33TV+W4Zq02VxVzz/eW8HmKs20FpHM61FBYfxQb22+NdvqWLGlOsO18Vz2z4/4v9cW8Z1/domJ3yLSw/WooHDL2fuTn2sATF20KcO18QTvBrdId4UTkS6gRwWFwrxcTttvVwBKN3WNTEFEpCvpUUEBYESJt4T2Rq2WKiLSTo8LCrsWe0tov1daxt/eWY5zumeziEhQzwsKgfsqANz++mKmL9mcwdqIiHQtPS4ojB5UFPbz9x75mIsfnEVLqzIGEZEeFxT2GdafW8/ZP9S3ADBzxVaN/hERIcmgYGZvm1m9mVUHHq/HKDfEzN4ws1ozm29mR3RudTvHJV/enSuOHRO2beqiTTS1dJm7h4qIZEQqmcJlzrm+gcfpMco8AKwABgITgefNrDBG2Ywa5ssUAO6eWsptkxZlqDYiIl1DpzUfmVk/vFtw3uycq3POPQpUASd01nt0phERQQHg4Q9W7vyKiIh0IakEhXvNbIuZTTGzA6PsHweUO+f8U4XnAftGFjSz8ngPoDi1w0jdPsP6M3pgnx39NiIi3UqyQeFqYAywGzAFmGRmfSPKFAGRvbWVQGS5LiE3x3jlZ8cwJmI0kohIT5ZUUHDOzXbOVQeahe4EqoHDI4rVAP0itvUPlI18vZJ4D6AijWNJWf9e+fzguD3Ctl360Idc/OAsmtXpLCI9ULp9CtHOmKXAADMb6tu2P7AwzffYKXYpKgj7+b3SMmau2Mr7y8oyVCMRkcxJGBTMrMTMTjWzQjMrMLNfAQOAsLWenXNVwKvADWbWy8wuxcsU3t4B9e40w4p7Rd3equUvRKQHSiZTyAduB7YCG/FGGJ3unKsys0si5iz8CK/DeRtwDXCec65L3+bsgBHFnH/IyHbbzSwDtRERyay8RAWcc1uAQ2PsewJ4wvfzJuC0TqvdTmBm/OGCg1i3vY6ZK7aGttc2tGSwViIimdHjlrmI5VtH7B72c01jc4ZqIiKSOQoKAQP65If9XNugoCAiPY+CQkBJn/BRSNUNzcxfV0Fjs4amikjPoaAQMLhf+BJN90xbxtfvfZ//ff5zLastIj2GgkLA4H6FXP/VfUI/BzOEF+es48CbJvO3d5bv8DroLnAikmkKCj5XHLcH504Y0W57TWMLt7++eIe/vzIS6cpaWh2vfr6e1VtrM10V2YEUFCL07ZVwlO4O06ygIF3Yv2au5KdPzuG4idMzXZUu6/3SMi5/9CO+KKvJdFXSlrkzYBc1dnDm1u9TpiBd2b8/WZvpKnR533roQwC+KKvhrf85IbOVSZOCQoQJu5XslPf589RS1pWHp+HNLQoK0nXpzoTJW75FmULW2GdYf0aU9GZdeV27fc65Tln+YvmWav40dWm77c2t+qOTrqtJFy09gvoUIuTn5nDbuQdE3dfQSXMWttc0Rt2u5iPpynrinJ115XWcff8H/GvWqkxXZadRUIjiuHGD+OmJe7bbXt1Js5yr6qO/TpOCgnRhjT2w+ejq5+by2ZpybnhpfqarstMoKERhZlx12nhG7RJ+H+fHZq7i1tcW0tTSyscrt1FR15TW61fFCC5NPfBKrKeqqG2iojb1749zjvqmzCzW6O9TaO0hFzBLNlZlugo7nYJCHPdd/KWw+y3c81Ypf3/vC7790GzOf2Aml/xjVsqvWdvYTHWMTKGzMpGuTHe0g7rGFk68622O/8P0lD/zK5+dy5dumcLyLe1uaNhh9U0tvPr5erZWR1/t3n/R0tQJ/V/OuS4fXGJl9dlMQSGOg0aV8N7VJ7bbHlxie/66yFtSx3f/9GUccNObPPfJmqj7sz0o/GHyEva7cTJvzN+Y6aqEWVlWw9XPzWX+up1yF1gWrK9gW00j5bVNzP5ia+In+Lw4Zx21jS3csQMmU97x+mJ++uQcLvnHh1H3+zuaO9rp3NrqOPevMzj1T+90auazbHM1900rDWVhK8tquOvNJWyurE/r9TqrH7E7UVBIIC83h0N2H9AprzVx8hJaWh2fri6Puj9WBpENPl29nfumL6OhuZXnuth492899CHPfryWc/8yI265itomrn3hc6Yv2dyh9/MPKKhtTO+EmOyghNZWx29fns8ht0xJWO9HZqwEYHGMJhN/n0JHO52Xbalmzupylm+p4d2lW+KWbW5p5aU561izLfFM6lP/9A5/eHMp1780D4ArHvuYe6ct43uPfpTgmTvGe6VbOO+vM5izentG3j8dKQUFMzvSzFrN7JoY+1eaWa2ZVQceD3RONTPr8mPGxNxX39TCuvK6ThmZkc2ZwpMfrg79v64puePcWWtBrd3uDT9O1JF626RFPDV7Dd99OPUTTG1jM++VbqGhuYVyX19UWVXyNyb0B4KcJEdGL91cxWMzV7G1ppGX5qxL+r0S6eicBX92kOhTfnPhJn75zGcce+d0ymI0bYVeK/Bir36+AYDSzV4zmz+rd86l9d1K1PQZ7TUvfWg2n6zazjkJLji6kqSDgpnlAH8i4t7MUZzknOsbePywQ7XrIk4/YBizrz+Zhb87jTMOGh62b9rizRx9xzSueOzjDr9PrA7oZK3ZVstLc9Z1yUlGGyva0veyquhDcv1++K9POPmud7pUoPx41ba0n3vlM3O59KHZ3PraIspr245/QwrNGrVhN35KLir428Q3lKfXhNL+vTueKfiz4twEc39Wbm2bCPbBsrIOvW9LoNnqa/e8T0NzallaoqakqYviZ2Krt9ayemstN72ygGWbO79PqLOkkin8APgQWLSD6tKlDenXiz4Feew7rH/Y9p88+SkA7yRIgZPx4LvLE14J+a0sq2HSvA2hzrpT/vgOv3zmMx58d0WH69LZ/FeGGxOcCOubWnhjwUZWlNXw9OzVccsCTJq3gZ888WnY766xuZUfPf4Jd0eZJJiuyOvAhz/4gnP+8kFYwIvljQVeP8pjM1dR7ht1tLky+c+7ztfUFCtTaG113PNWKZPmbWj3nM1V6QeFORFNnh298PBnS7UJ+hSamtt+8zUduE1uRW0TyzZ7zVYLN1Ty7tL4ASbyyj9RUEh0YXjcxOmcef/7PDJjJefc/0Fylc6ApIKCmQ0EfgncmETxl8xsk5m9aGa7RytgZuXxHkBx8oewc11w6EgG9W2794L/exO8mvKn+c457p++LPRHGs+abXX8NBBkknHCH97mx098yvOfem30wS/tEzthok2izsFPV2/nplcWhE7U9b6rsoq6prjPr/FlB8lkCj9+4lNem7eBW19ru175z9z1vD5/I3dPLW13lZs232f9yart3PyfhcxZXc5vUhzDvt0XFKKNbonVtFHjO8HHurh+p3QLf5yylB8/8Sl1jS1hfRbrK+rTbjp5+bPwpqey6sYOjSTzB8aaBJ9xra+5sS7JTuncQNTsld92iqtubA4LZttr42eskZ3pndEhHjzujrYK7EjJZgq3Anc756L3kLb5JjAaGAesxgsQHV8XogsZ1LeQaVcdH3Xf0XdM4+IHZzH2uknsed0kVpbVMHnBRiZOXsKPn0juZD9rxbaUZzYH20+DEl3RdHRY6Px1FRz8uze57sV5Mcuc+5cZPDJjJde94JWpi+hQXRVn+WX/iWzJxqqkr0pLN7d1kFbWt5101kdZsqTtvZL/42z1nUzP+2tbG3Gqw0Mr6tpORpHvP/uLbRx261v8a+bKds/znzwtRvPRFl8fxaKNlWEnssbmVh7/cDWH/N/UpC5SghqaW9p9x77xt5n88PHkL2Ai+U/IiYJCve/7EO/E7P9eR8uk6ptawi4yEg3siAxAPWUkUsKgYGYTgMOAvycq65yb4Zyrd85VAlfiBYd2vbTOuZJ4D2DnjA1MU/9e+ey9a79227fXNoWGqza3Oq554XMe/mBlyq9/2yTvitc5x9JNiU+KkWE33h/Ov2atYv+b0h8W2tjcytfvfZ/6ptawzuOgmcu3cu0Ln4d+fq+0LFCn8GPwn8AjPRoYBQPw+vyN/CTJgJqXk8Pa7bXUNDTTOz83tD3YkRzNujj7IsUK1YlOapG218S+Sv7+ox9RVt3ADS8vCNte3dDM715dGPo51qVWL99xz1tbwZSFm8L23/DSfLbVNCa8SPHPH/hwxbaoo6SmLtrUbluy/BM/42WD22oa+XhV28idyIsLP38m1dTi2vUZ1DW2hGVmWyKaamsbm3n4gy9C6541RPwddfakwdVba7vk0jbJZArHA+OBdWa2EbgQuD6JkUWOxAMLuq3/Pn6PhGVmrdjGh1+k3jn50PtfAPD0R2v4yp/e5bcvt2+e8Hf0RZ4f6uNc0dzw0nzqm1r54eOf8EVZDTe9soCVvrXf/zhlKde+8HnMSUWPzPgi5ms3tbRy8d9n8dTstnkYwTQ+8g80VkdbRW0T/3g//D3eXBj75OP/Q12xpZpjfj+dCx+cGXaiiRcU/BlLQW78P4dYrS6pDisN9i9A+IkMYk+W+sPkJcz2fZf8TRvvl5aFhmv6r6pvfGUBr8XJCFpbHcs2V9Ha2r5Jyd/ct3STF8D3HNJ+Wfl0T5T+jCZeUD3t7ndZsL5t5NCKsmqOnzidm15Z0K7s2xFDbu+btizs54bmFqp8GeRf317Oe6VtfYFPz17Dzf9ZyNF3TKO11bXLFDo7KBw3cTp/fquU1lbXpdaVSiYoPAjsCRwceLwC/Bn4X38hM9stMGQ138yKgInAF4FH1jn74BH8/ORxaT/f39YZy7WBphf/STaoLqx92cJS58irj8/WlEe9W9ZFD87kkRkr+c7DswHvRH3PW6U8NXsNf5yyNOwPKOilOevDfvafTD5c0T4ABtP4YH2Dxx1rJMySTdEziK/86Z2wDCLI30xUGTihzl9XGfo/ELcjeJVv7HtjS2vcK7fWmG396bcPJ5tl+E9e0HaCmr54M9966MPQOv6pNIftd+NkTvnju9w1ZUm7wLa5soG122sDr+ntK+md3+41UhkY4eefc1Ado/PYORcWPAAmzdvIqq21PDJjJYs2hE8e/cXTn4X9fO+0ZWEZal1ja7us5LGZbf1v/nkcn67e3i67TbX56PbXE4/JueetUs5/YAZH3TEtrWVPdoSEZybnXK1zbmPwAdQB1c65CjO7zpcx9MMLIOXASmAscJbL0hsPmxlXnroXc2/8CqW3ns6PTxib0vOP2XNQ3P1X/Xtu3P3+E9H22saYozJKN1Vx9v0fcNzE6e1OeJsCI1+CV8szlreNxrhv+rKoi4D1KcgN+9l/NRWtSaiyvpljfj8tdEU8pJ+3bEh1jJPX4o3RZ4kv3VTNja8saHcMlXXRX2eTLxBURgluwcxl8oLwZrQz7n2ft2I0i8Rq5nDOOxl//5GPeCDiXt61gc7Ngrzof2rxgkK8K9Pg7/3ut0oB7zOsbmimrin5E1fwNe6fvrxdv8uZ973PSXe9w9w15aHvWp/C9ivtb62O3lnb1OI1L071ZXn1TS18sKyM+qYW1mxvCwqxfgeJOpVTHZ5a19TSLhPbFlix+I35G0NNnQBrttd2OFP42zvJjQL8dHU5ZdUN/OHNJXznn7PbNfntbCnPaHbOXeacuyPw/9uCcxGccwuccwc454qcc4Odc2c757IyS/Ar7p1Pfm4OZx48PNRUksxzvn3k6LBtt50Tvlx3tFm/znnpfkNz+KiSOavLueyR2VHfa9ritqufz9bEn1UZeeX10mfr25XpHREUgp11m6vqufk/C9uVh/Dmm0F9C4C2ZhLnHNtqGkMZh/8PM5opC72T+G2TFnHcndNjLk3hvx/Giogbnjz54Wr2++1k7p++LKxJBmDhhkq+/2j7oYWNza1hI2YiPfDOCt5avDm0/ERzSysbK+o58vZpnHHv+zEzkJrGlrBsy/8d2ve3bzBxcvTlLIKZl3/5hv1vnMwnac6l+FFEH0NlfTONza387Kk51AYuOIoiPnuAD5aXtWt6un/6MsZd/zrXvTiPyx/7OJRx3vjyAi75x4fc9MqC0AUJwCtz1/Pxyvb1jhXwg/xIDi9mAAAUxklEQVTfq2QGT1TVN7W7QNhe00h9Uws/fPyTsO3/+/w8LnggfMJZZH9Gc0srd725JDQyq6NDav41axXvLN3Cr5+Lf0G4o2mZi06y9679+ej6UxJmDDOvPYmpVx7P4H6FYduH9i+M8QzP2u21vDF/I6f88V0uenBWu6adyHHkjc2tNDa3hnU8nvfXmTFfv6KuqV2qDrRLaYsKwq8Wg800L89pH0CiCR53dX0T68rrOOzWqXzplin8/b0V1De18M6S+PM9/ufZuazdXsuD765g9bZafvnMZ1HL+a98319WFjak8roX59Hc6pg4eUnM97l/+rKwE3lkp2SkewJX7AAbKuo4+vfTOOL2t6ioa2LxxqqYQaGl1YWaJWoamsPKtTrvKv7RGSvb3clrXXkdDc0t7YZNJppAFUusPp7V22pDmULkBQHAnW8s4flP14WuouubWtr9XoOZ6DMfe82gT3/Uvjn0/AfafzejNV/6PTJjJZc/+jGbq+rZlMTM8CufndvuO76ttpFNUebNNDa3tvvd+vtn6ptaeGr2au6dtoxfPP0Zby3aFLPPKVXxLj52Bt15rRPtUlTAT0/ak8/XVvB+lNR25IDeDCv2luN2vj74Hx4/llG79In72sf8vu1m6XNWl/OjBMMBL3hgBnPXJj+I67oX5rE5yh/WzBVb+a/9dw39HJkNBZtUkl1rPxQUGpr59b/nUhZofrht0mKOGjso4evUNLZw5TOJr6QiO5d/8fRnjBzQm49XJrcGzcTJS5g4eQk/P3kcV566V0oLqr0xf2PYlXAiNQ3N5Ofm8LV73ou6/8YonaoVdU28/Nn6hGPtO8MLn3oBtaggjwNGFDMvIju76t9zefajNXx5j124N6JzF7zJhdH6g6LZWt1ARV0TewzuG7XZL9LURZtoeLaFXyTZvxc8liP22IVZK7zl7+MNRPCbvGAjFXVNFPfO56dPzgkbfRUtu+yuFBQ6WZ+CPB6//Mtsrqzn7Ps/YHC/QnYfWMSGijruOO/AULkh/Xrx928fSl6uceL4IVE7guNJNCs4lYAA3lWQf5nwoJc/W8d/7b8rj3zwBZ+vq2h3Ejr7/g/4/jFjkp77MLiv9x5LN1WzdFP41enX730/7nODf8izozQ1RIoWXM5/YGbMq7mjxg5kxvL2K5be81YpPzlxLN/4W+wsK1KsZjS/3BwLZQUfrdzGASNLWJnid+Dq5z5PXChF/XvlhXXS+/UpzOVvlx7Cko1VTFm0KWxI8uyVsT+Xv7y9POr2SL95aR6Pz1pNbo7xxi+O5aMYAXxgUQFbfXcvfK+0jAsOHQVASZ98fnfW/vz8qTlx3+viw3dj1optOEfCWe/9CvNobGmlobmVGcvKOH784JSH4547YQQVdU28tThxJmfmZZDJNkd3NgWFHWRI/17MuPbkuGVO3Xdo6P+7DezDxYfvxlNJLOuwowRT668dMIxDRw/g5v8s5K1Fm1mwvoKb4pzoHno/+a6jQf0K0q7fieOHMCvKCKdkxUvvn7j8y6woq+Hku95pt+/oO6Z12v2JH/rOoew+sA9jBvXla/e8x+KNVfzjvS/47tGxF13sDH+79BBGDejDnDXbGTu4L4s3VDJ/fWW7vqvLj92DP06JfpIsKshjeElvhpf0Zr8R/aPOU+mIx2d5r9fS6vjave/HHKb51QOGtbs95sMfeN/BYcW9GTmgd7SnheQYfHnMwNDPsYJPUFVDM/sN78+C9ZWs2lbLK1H62vx+cNweYUvN7DW0L3+88GCe+2RtUkHBOXj+k7Wcf8hIcjIQGNSn0IXcfu4BPPzdw0IdeoV5OYyLMjY8HfsN75+wTHPgyvWXp4zjwsNG0TdwhfS1e+Jfwfv5A100wdFH8dx+7gH8+rTx7bZ/49BRnDNhBCftPYQLDhnJzWful3S9EjEzxg7uy90XHtxuX5lvhM3NZ+7Xofc9eFQJew7pR26OceFh3tXtx6u2h9bQStfE8w+M2tF51wUH8flNX+G0/XZl3+H9ueTLu3PEHgO57OgxXPWV9r9j/6S/SP6XH9KvF0//4IgO1TlyJJtfZEB49HuHh/6fl2sMieiTC/apDS/uxYRRJdxy1n48eOkhTPnVce1ee1hxb3Yt7sUuRcldoBw2egDDS7xAU7qpmt+/Ef9eFt87egx///ahoZ+DTabH7zU47vMKfSPUrn7+8w53XKdLQaGLOXH8EObffBpzb/wK7119IucdMjLl14j2Zbr+q/uERv7Ec8o+Q9lzSF/6FOTxy1Nit9PmxbiCueiwUTz6vcN563+O58JASh90wvjB9O/VPjn1n4gPGlnMxYfvxonjh7QrN6CogD9deDD/vOwwJl5wEN85anS7VWvjOfvgxGXPnjAi5sluzg2n8p2jRvOdo0YzemD8PqBoigpyKenT9hnsNbT9rPho9hhUlLDMBYeOYsktp4dtO3XfoZx3yEj692o/vwC8wQ2RM/P3HBr7IiRyvZ4j9hjIw989jMPH7BK2ffckfzfDS3q3G3ARzbkTRnDU2LYr+4NHlTA6xu+kuE8+ZsalR47mK/vtyrgov+Ngs91TV7T/nPNyjPu+OYG+hXlcfswYrj19b24+c39GBILC85+uDVu7Kppe+TkU++Z0DAh85oP7FTLp58dGfc6fLzqYz377ldDPE3YrIVMrBKn5qAsyM+9L1Tufy44azepttbz2+QZO3mcIr32+Ieokml+fNp77py/jxPFDuPGMfSnMz+W+aaVMW7yZuy+cwAEji3nqiiP416xV/GvWqlBTymVHjWa3Xfrw70/WMrR/IX+68KDQl/Gyo0bzXmlZ1BVgrzl9b6Yv2czC9ZXUNLbQ2NzK2QcP56S9h4Se//vzD+SsCcP55t+9iVXXnr4PQ/sXkp9rYc0xB4ws5tenjefBd1dw8eG7AdC/d/hX888Xtb+CBy+gnDthBN99xFvR/dZz9mdQ30Lqm1oYXtKb+6cvY49Bfbnuq3szc8XW0DDbf//wSC6IMuIFwrOqw0fvwuyV2zh3wggG+K4sJ15wEB+u2Mrwkt5c+Wz0ju8/XXgQC9ZVcsDIYl75bD2/PGWvsHbi8VGWSjlmz0Hces7+zP5iG79+7nOOHTeIB751CPvdODlU5thxg8KG7v7Xft5AgIK8HH531n488sFK7rl4AvuPiL+upJnx928fyhWPfczijVX893F7cMJeg7nz/AO5+rnPOXrPgfQtzGPyAq/9vDLKPclPHD+Ew0bvwr3TSinMzaFfr3yOHDswav/QSXsP4ZwJI3ivdAv//mQt3z5yd75x6CgmzdtAXVML17/Yfl7MmQcN5/bzDiA/N4cXfnwUn67azun7D2PckH6ccV/7ob4Hjypp9xpXnrpXWJPYKft6Fxx7DC5i1C69WbOtjl75Ofznp8dQVOg1kZ2+/7CwzyoYFILOOng4L8doRirMy2X0oD6h7/mBI9s+h32H96dfr7yw+RK983M5dtxgeuXn0Lcwj+qGZn5ywp5RX3tnsK44t8zMyouLi4vLyxOtv9czXfTgTD5auZ3Xfn4M5bVNjB3cl8H9CmlsbiU/1xJeYdQ3tZCfm8P22sawFV9jaWj2TvpLN1Xzh8lLKN1cxas/O5Zdo3RMR3LOcefkJQwr7hWam7GyrIa8XOOVud4f1Y8DfwDOuVDdK+ubOPCmNwG45Mu7cWvEPA6/1lbHc5+uZa+h/aKeFPx1+fNbpRQV5HHFcXvwr1mruOGl+fzma/tw+bHhy5Y8PXs1pZurue6r+7C5qp5BfQvJj7EMRkNzCxvK63nqo9WhCUsXHjqK359/YNTyftMWb6K+qZWxg/sye+U2Ttt3KEP6t/+9LtpQyZ+nljK4XyHXf20fVm6toaK2iSWbqvjGoaPChh6nyls5lVD7tbfmVjVjBhVRkJfDyXe9zfItNbz+i2PZZ1jiZkiAj1duw8ybRfy3d5dzzel7s9/wtpNjTUMzRRGT4X765KftFt675ez9ufSIqIsts6Wqgf6987j2+Xm8MGcdZx08nN+fd2C734VzjrXb6yjMy+H5T9dx3pdGhH7HNQ3NVNQ10Ss/N25z0vx1FaFAt8+w/jx5+Zc5buL0qMuSrLjtq+TkGEs3VVFW1cARewwM6xtYsaWa+6cv5yv7DeWtRZu46PDd+NJu3t0dF2+sZN32Ok7eJ34zbCwlJSVUVFRUBNaQS4uCQjfU2NxKVX0TA5M4oe8Ira1uh3eAtbY69rhuEuC1Jydqj03X9prGsAygoxqaWyjMS/8E3RXVNjaztbox4bDpjtpe08gjM1Zy+Jhd+M/c9VQ1NHPXBQclFfCaWlpjBu3OUt/UQl1jC/1755ObY6zdXssngcylrqmFXz3zGfsM68evT9t7h9YjHgUFyWozl29lzfZaLjhkZMbaV0W6k84ICupTkC7ryLEDOZKBiQuKSKfR6CMREQlRUBARkRAFBRERCVFQEBGREAUFEREJUVAQEZGQrjpPoRWw4uL40/RFRKRNRUUFgHPOpX3B31WDQjNeFhP9hr3xBSNJajcU6Jp0LF2TjqXryZbjgI4dS3+g1TmX9hy0LhkUOsLMygE6MqOvq9CxdE06lq4nW44DMn8s6lMQEZEQBQUREQlRUBARkRAFBRERCVFQEBGREAUFEREJUVAQEZGQrJunICIi6VOmICIiIQoKIiISoqAgIiIhCgoiIhKSVUHBzIaY2RtmVmtm883siEzXKRYzu9nMFppZq5ld5Nuea2b3mlmFma03s8sjnvf9wPYKM7vPzDL+GZpZoZk9bGZrA/Wabmb7+fb/xszKAo8bIp77dTNbYWbVZva0mfXZ+UcQzsz+amYbzKzSzOaZ2dd9+7rVsQCY2ZGB79k1gZ+743fsbTOrD/xuq83sdd++bvWZmOc3gd9xlZnNDGzvGp+Lcy5rHsALwF+A3sB3gHVAYabrFaOu3wJOBWYBF/m2/xz4CBgIHAxsAw4K7DsY2AocFNj/CfCzLnAsRcANwEggF7gSKA3sOxP4AtgN2B1YBZwR2LcrUA6cAvQDXgXu6gLHszfQK/D/QwN1HNBNjyUn8B37ELimG3/H3vb/nfi2d8fP5OfAFGBE4PM5pCt9Lhn95XTyL7of0AgM9W1bDJyW6bolqHfYlz3wB/wN388PALcH/n8H8BffvouAmZk+hijHVAC0Br68TwNX+/ZdAzwV+P8PgUm+fUcAGzJd/4hjORSoA/btjscSqNefgUd8QaHbfcfiBIVu9ZngXTRtAMZE2dclPpeMp4WdaBxQ7pzb5Ns2D++PuTvZF5jv+9l/DPH2dSVHApudc1tJ/Xh2NbMBO6WWcZjZX8ysDu/KbSqwiG52LGY2EPglcGPEru76HbvXzLaY2RQzOzCwrVt9JsAooBdwsZltCjQhfyOwr0t8LtkUFIpof6e2SqBvBurSEZHH4T+GePu6BDMrAR4ErgtsSvp4nHM1QAtd4Jiccz8O1ONUYKrzLs+627HcCtztnCuP2N4dv2NXA2PwmommAJPMrC/d7zMZDpTgNXXtDnwXeNDMxtNFPpdsCgo1eE1Ifv2B6gzUpSMij8N/DPH2ZZyZ9QJeBl51zv0zsDnp4zGzIrz0uksck3OuxTk3FTjFzE6jGx2LmU0ADgP+HmV3t/uOOedmO+eqnXN1zrk7A3U6nG70mQTUBf691TlX75z7EHgT7+KjS3wu2RQUSoEBZjbUt21/YGGG6pOuhXj1DvIfQ7x9GWVmuXjtu2uAq3y7Uj2ejc657TuwqunIA8bSvY7leGA8sM7MNgIXAteb2QN00+9YhNbAv93pMwFYCjTF2Nc1PpdMdyB1cifOC8B9eG12l9K1Rx/lB+r5bqCuvfCC9M/xRorsgjfSYCvhIxDKgAPxOnE/oguMDAnU7Z/AG0B+xPYzgRV4af9ueCNF/KNDtgMn4aXCr5Dh0SGBenwr8G8ecB5QDxzQnY4F6BOoU/DxDPB/eDeF71bfMbzmllOBQrxBDL/C66zt150+E9/xPI03SrIAOCRQx727yueS0V/ODvhlDwUmA7XAAuCITNcpTl0fAVzE4wS89PY+oCLwxb884nmXB7ZXAPcDuV3gWHYP1L8OL6UNPo4N7L8h8AXfCtwQ8dwzAn/INYETV58MH0sRMA1vKGMF3tC/c3z7u82xRPm+BUcfdavvGDAY+DjwndoGvAUc3F0/k8BJ/ZXA8ZQSGHHUVT4XrZIqIiIh2dSnICIiHaSgICIiIQoKIiISoqAgIiIhCgoiIhKioCAiIiEKCiKdzMxGm5kzs0GZrotIqhQUJKsFbs7S4Ls5S/CRm+m6iXRFCgrSE1zvnOsb8WjJdKVEuiIFBemxzGxl4LaIHwVui/i2mY317R9iZs+Y2ebArUbvDqwEG9w/zsxeCazxv83MXoh4i6+ZWWng9onPBlbpDN6+9KHA8yoDa+qfsHOOWiQ+BQXp6a4ALgOG4N1I5wUzs8C+pwL/jsVbpvlo4DYILcM8Fe/GJ6PxFl+7N+K1z8C7c9s4vAXOfhTY/h1gArAX3gJ1ZwCrO/WoRNKkoCA9wS1mVu57zPDt+5tzboFzrg74Nd7drA4ysxF4q2te6Zyrcs6tx7uD2WWB530db/nm651zNc65Rufc9Ij3vdk5V+Gc24y3ANohge2NeKt27gPkOOeWO+dW7IDjFkmZgoL0BDc450p8j6N8+1YF/+Ocq8ZbaXMEMBJods6t85VdgXfPjj54K8Mud/FXlNzg+7//JimP4y01fj9QZmaPm9mu6R6cSGdSUJCebvfgfwK3dxyIdx+OtUBeIGMIGg1sd87V4gWTPXxNTUlzzjU75+5wzk3Aa1rqD9yZ/iGIdB4FBenpfmBm+wY6kO/A61f4PJAhTAf+YGb9zGwYcDPwaOB5r+HdKOl3ZlZkZgVmdmIyb2hmJ5nZwWaWh3fvj3q8eweLZJyCgvQEt0aZp7BbYN9DeCf6LXh3VzvPORe81eM38e6+tgLvJi8fAtdDqKnpFOBLeJ3EG4GfJFmfoXid2OV4ty/NA67p2CGKdA7dZEd6LDNbCVzlnHsu03UR6SqUKYiISIiCgoiIhKj5SEREQpQpiIhIiIKCiIiEKCiIiEiIgoKIiIQoKIiISIiCgoiIhPw/cXJ8tyNZy7sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(hist.history['loss'])\n",
    "#plt.plot(hist.history['val_loss'])\n",
    "plt.xlabel('Epochs')\n",
    "np.mean(hist.history['loss'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<tf.Variable 'dense/kernel:0' shape=(3, 1) dtype=float32, numpy=\n",
       " array([[-0.00728608],\n",
       "        [-0.43613392],\n",
       "        [-1.288468  ]], dtype=float32)>,\n",
       " <tf.Variable 'dense/bias:0' shape=(1,) dtype=float32, numpy=array([2.660357], dtype=float32)>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For comparison the results of the r-code?\n",
    "\n",
    "See also: https://stats.idre.ucla.edu/r/dae/negative-binomial-regression/\n",
    "\n",
    "* null deviance: only the intercept?\n",
    "* residual deviace all included\n",
    "\n",
    "```{r}\n",
    "> summary(m1 <- glm(daysabs ~ math + prog, data = dat,family = poisson))\n",
    "\n",
    "Call:\n",
    "glm(formula = daysabs ~ math + prog, family = poisson, data = dat)\n",
    "\n",
    "Deviance Residuals: \n",
    "    Min       1Q   Median       3Q      Max  \n",
    "-4.2597  -2.2038  -0.9193   0.6511   7.4233  \n",
    "\n",
    "Coefficients:\n",
    "                Estimate Std. Error z value Pr(>|z|)    \n",
    "(Intercept)     2.651974   0.060736  43.664  < 2e-16 ***\n",
    "math           -0.006808   0.000931  -7.313 2.62e-13 ***\n",
    "progAcademic   -0.439897   0.056672  -7.762 8.35e-15 ***\n",
    "progVocational -1.281364   0.077886 -16.452  < 2e-16 ***\n",
    "---\n",
    "Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1\n",
    "\n",
    "(Dispersion parameter for poisson family taken to be 1)\n",
    "\n",
    "    Null deviance: 2217.7  on 313  degrees of freedom\n",
    "Residual deviance: 1774.0  on 310  degrees of freedom\n",
    "AIC: 2665.3\n",
    "\n",
    "Number of Fisher Scoring iterations: 5\n",
    "\n",
    "> logLik(m1)/314\n",
    "'log Lik.' -4.231346 (df=4)\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Fitting a glm using the convenions functions in TFP "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor: id=37559, shape=(4,), dtype=float64, numpy=array([ 2.65197378, -0.00680838, -0.43989746, -1.28136411])>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# First build the model matrix, appending an intercept term\n",
    "Xm = np.ones((np.shape(X)[0],np.shape(X)[1]+1))\n",
    "Xm[:,1:4] = X\n",
    "Xm[1:4,]\n",
    "# Doing the fit \n",
    "fit = tfp.glm.fit(\n",
    "    model_matrix=Xm,\n",
    "    response=y,\n",
    "    model=tfp.glm.Poisson()\n",
    ")\n",
    "fit[0]#2.65197378, -0.00680838, -0.43989746, -1.28136411"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other then Poisson\n",
    "\n",
    "Other than poisson cases.\n",
    "\n",
    "### Using zero inflated poisson\n",
    "\n",
    "We define a network with two outputs. One for the poisson mean, one for the prob that there is a zero. \n",
    "\n",
    "First, we define a function which produces the zero inflated poisson distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Definition of the custom parametrized distribution\n",
    "def my_zero(t):\n",
    "    mean = tf.exp(t[:,0]) #Needed so that Poisson has batch_shape=(5,) \n",
    "    #print(mean)\n",
    "    s = tf.math.sigmoid(t[...,1:])\n",
    "    probs = tf.concat([s, 1-s], axis=1)\n",
    "    #print(probs)\n",
    "    return tfd.Mixture(\n",
    "          cat=tfd.Categorical(probs=probs),#batch_shape=(5,) event_shape=()\n",
    "          components=[\n",
    "            tfd.Poisson(rate=mean),   #batch_shape=(5,) event_shape=()\n",
    "            #tfd.Poisson(rate=0.00001*mean), #Hack to produce 0\n",
    "            tfd.Deterministic(loc=tf.zeros_like(mean))\n",
    "        ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Logging before flag parsing goes to stderr.\n",
      "W0611 11:15:48.761339 140333529929472 deprecation.py:323] From /usr/local/lib/python3.5/dist-packages/tensorflow_probability/python/distributions/mixture.py:154: Categorical.event_size (from tensorflow_probability.python.distributions.categorical) is deprecated and will be removed after 2019-05-19.\n",
      "Instructions for updating:\n",
      "The `event_size` property is deprecated.  Use `num_categories` instead.  They have the same value, but `event_size` is misnamed.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.     1.   ]\n",
      " [ 1.     1.   ]\n",
      " [ 1.     1.   ]\n",
      " [ 1.     1.   ]\n",
      " [ 0.999 10.   ]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<tf.Tensor: id=37615, shape=(5,), dtype=float32, numpy=\n",
       " array([1.9872233, 1.9872233, 1.9872233, 1.9872233, 2.7154417],\n",
       "       dtype=float32)>,\n",
       " <tf.Tensor: id=37617, shape=(), dtype=float32, numpy=2.7182817>)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Testing the distribution\n",
    "t = np.ones((5,2), dtype=np.float32)\n",
    "t[4,0]=0.999\n",
    "t[4,1]=10\n",
    "print(t)\n",
    "my_zero(t).mean(), tfd.Poisson(tf.exp(1.0)).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"model_1\"\n",
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_2 (InputLayer)            [(None, 3)]          0                                            \n",
      "__________________________________________________________________________________________________\n",
      "dense_1 (Dense)                 (None, 1)            4           input_2[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "dense_2 (Dense)                 (None, 1)            4           input_2[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "concat (TensorFlowOpLayer)      [(None, 2)]          0           dense_1[0][0]                    \n",
      "                                                                 dense_2[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "distribution_lambda_1 (Distribu ((None,), (None,))   0           concat[0][0]                     \n",
      "==================================================================================================\n",
      "Total params: 8\n",
      "Trainable params: 8\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "## Definition of the model\n",
    "inputs = tf.keras.layers.Input(shape=(3,))\n",
    "mean_logit = tf.keras.layers.Dense(1,bias_initializer='zeros', kernel_initializer='zeros')(inputs)\n",
    "p_zero_logit = tf.keras.layers.Dense(1,bias_initializer='zeros', kernel_initializer='zeros')(inputs)\n",
    "\n",
    "# Combining the path\n",
    "tt = tf.concat([mean_logit, p_zero_logit], axis=1)\n",
    "p_y_nb = tfp.layers.DistributionLambda(my_zero)(tt)\n",
    "\n",
    "#Glueing input and output together. Note that output is a tf.distirbution \n",
    "model = tf.keras.models.Model(inputs=inputs, outputs=p_y_nb)\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def negloglik2(y_true, y_hat):\n",
    "    # return -tf.reduce_mean(y_hat.log_prob(y_true))\n",
    "    return -tf.reduce_mean(y_hat.log_prob(tf.reshape(y_true,(-1,))))\n",
    "\n",
    "model.compile(optimizer=tf.optimizers.Adam(), loss=negloglik2)\n",
    "hist = model.fit(X, y, epochs=1200, verbose=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.8213007666882435"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEbCAYAAADajfNFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XeYlNXZx/HvvbvsLkvZFZBeBUQRBEQI2BB7N2piiya2GBOjrxoTFVuMoiaaiJGo0RhrLEk09ooCinSkSO+dZWnb2L5z3j9mdpnZNjPbZnbm97muuXzmOWdm77OLc885z3nOMeccIiIiCZEOQEREooMSgoiIAEoIIiLio4QgIiKAEoKIiPgoIYiICKCEICIiPkoIIiICKCGIiIiPEoKIiABKCCIi4qOEICIigBKCiIj4KCGIiAgASZEOoCZmVoY3WeVGOhYRkRakPeBxztXrs93C2Q/BzMYC3wITnHOP1lC+EegMeHynXnPO3RB2UGYewNLT08N9qYhI3MrJyQFwzrl6jf6EnEXMLAF4ApgXpOpJzrnZ9QnGT256enp6dnZ2A99GRCR+ZGRkkJOTU++RlXC6FdcDcwB9bRcRiUEhJQQz6wjcAowBJgWp/q6ZGTATuMU5t6mG9wv21V9JR0SkmYU6zjQRmOScC/ZBfjnQFxgIbOZAchARkSgXtIdgZiOAUcCNweo652b6DovM7DYgB+gHrK9SLyPIz8xGvQQRkWYVypDROGAQsM33ZT8dKDOzvkFmEDnfQ0REWoBQEsJzwJt+z58E1gCP+Vcys95AD2A+kAz8Adjge4iISJQLmhCccwVAQcVzMysE8p1zOWY2Aejt6ym0w5s8DvHV/xY434Vzo0MDeTyOD7/fQU5hKace3oWu6anN9aNFRFq8sO9mc85d5Xf8sN/xMmBo44RVP2Zw21uLKPM4endIU0IQEQlDTK1lZGakt24FQHZBSYSjERFpWWIqIQCkp3kTQm5haYQjERFpWWIvIfh6CDlKCCIiYVFCEBERIAYTQoYSgohIvcRcQjhwUVkJQUQkHDGbENRDEBEJT+wlhLRkALZlF0Y4EhGRliXmEkLFNYSt+wrJyi2KcDQiIi1HzCWEYwd0qjxekZkXwUhERFqWmEsIXdNT6ZHRGoDMHA0biYiEKuYSAkA33xpGCzbtY6eGjUREQhKTCWHcoQcD8O/5WxnzyJcs317vPadFROJGTCaE0f06VB47B899vS6C0YiItAwxmRCG9gzcfbN1cmKEIhERaTliMiGkJScxuFv7yucpSUoIIiLBxGRCAPj9eUdUHu/dr70RRESCidmEMLpfB+4+63DAO9tIRETqFrMJAWBE7wzAu4xFQUlZhKMREYluMZ0QendMqzz+dGlmBCMREYl+YSUEMxtrZh4zuzNIvUvNzJnZpQ0Lr2EObptC61beC8rfrt0TyVBERKJeyAnBzBKAJ4B5Qeq1Ae4BljUstIYzM352TF8A3v5uKx8t2RHZgEREolg4PYTrgTnAiiD17gVeAHbXVsHMsut6AOm1vTZcpx3RpfL4xte/a6y3FRGJOSElBDPrCNwC3B+k3qHAmcDkhofWOEb0ygh47pyLUCQiItEt1B7CRGCScy47SL0ngTucc3VuV+acy6jrAeSEGFdQZsYx/TtWPt+VX9xYby0iElOCJgQzGwGMAp4PUu98oMw592kjxdZoXrx6FEkJBsCMNbWOZImIxLVQegjjgEHANjPLBC4B7jazZ6vUGw+cYGaZvnrHAM+a2d2NGnE9pCQlcuIg7wqot/17Mc45yj1Ow0ciIn5CSQjPAQOA4b7H+/iGhqrUuxdv4qioN99X58nGCrYhxg3qXHn85YosTvjTVC7++yw8HiUFEREIISE45wqcc5kVD6AQyHfO5ZjZhIqegnMur0q9EiDHOZfftE0IzdlDu1UevzBjA9uyC5m3cR/rdkVFeCIiEZcU7gucc1f5HT9cR70T6xdS0+jQJpmxh3Rk1vo9zFp/4Ca1VTvzGNilXQQjExGJDjG9dEVVg7pW/+DPzNEWmyIiEGcJ4ZJRvaqd25WnaagiIhBnCeGwru04dXCXgHNZSggiIkCcJQQz4/mfHh1wbmeuhoxERCDOEkKFCWcdVnmsHoKIiFdcJoTrT+jPMz85ClAPQUSkQlwmBICu6akA5BWVkV+s3dREROI2IXTPaF15vCO7MIKRiIhEh7hNCAe3TSE5ydv8dbv24/E4Xp+zmSVbgy3oKiISm8K+UzlWJCQYw3qmM2/jPuZs2ENmTiG//2A5ABsfPTvC0YmINL+4TQgAR/ftwLyN+3jx2410aJMc6XBERCIqboeMAE4/omvl8d79JZXHWhZbROJRXCeE4b0y6NMxrdr54jJPBKIREYmsuE4IAMcP7FTtXFFpeQQiERGJrLhPCL84oX+1c4VKCCISh+I+IfTqkMZzV44MuKhcVKohIxGJP3GfEABOO6IrU39zYuXzwhL1EEQk/igh+LROTqw8/m7zPm55cyFLt+VEMCIRkealhOCTnJTAwe1SALjn3aW8u2g7506eEeGoRESaT1gJwczGmpnHzO6spfx9M8sysxwzm2NmYxsnzObRLiXwPj3djiAi8STkhGBmCcATwLw6qt0FdHfOpQN/AN4xM2tYiM2nx0Gtg1cSEYlR4fQQrgfmACtqq+CcW+acK/NLAl2B6nd+RalHLzqy2rmHPlwecBeziEisCmktIzPrCNwCjAEmBan7IXAqkAz8zTm3v4Y6wZYUTQ8lrsbWI6N6D+EfMzaweW8Bz1XZelNEJNaE2kOYCExyzgVdG9o5dw7QDvgRMLcBsUXExzcfX+3c58t3RiASEZHmFTQhmNkIYBTwfKhv6pwrcc69DfzGzA6voTyjrgcQsfmeg7u356vfjIvUjxcRiZhQhozGAYOAbb5LA+lAmZn1dc7dEOS1yUA/6rjuEI36dmxDYoJR7tE0IxGJH6EMGT0HDACG+x7vA08Cd/hXMrM+ZnaumbU2sxQzuwnoASxo5JibXEKCcdFRPQLOLdi0L0LRiIg0j6AJwTlX4JzLrHgAhUC+cy7HzCaY2bN+1ScAO4FM4GLgXOdcixyAH9oj8Lr2Rc/MjFAkIiLNI+wd05xzV/kdP+x3vAloUTei1eXInhmRDkFEpFlp6YpaHNkznQtHBA4bfb9VaxuJSOxSQqiFmfHni4cFnPNf20jbbIpIrFFCqIOZkZxU/Vf07sJtjJr4JdNWZUUgKhGRpqGEEMRTl40IeF5UWs4tby1id34xV71Y17JOIiItixJCEEf2DJxt9M532yIUiYhI01JCCKJbemseOO+IyucT/vd9BKMREWk6SgghuKDKTWoiIrFICSEEVTfOERGJRUoIITAzVj54Bml++y6LiMQaJYQQpbZK5Kpj+kY6DBGRJqOEEIYbTuwf6RBERJqMEkIY2qe24qjegWscebREtojECCWEMP3nhmMCnn++PFPLWIhITFBCCFNigvHS1aMqn9/w2nd8/H1mBCMSEWkcSgj1MO7QgwOeP/f1ughFIiLSeJQQ6sHMaOt3b0JKknc66uY9BXy6NFPXFUSkRVJCqKept59Yebxwi3d7zXGPT+WG1xbwv4Va70hEWh4lhHo6uF0KVx/bF4DSckd+cRkV15b/u2Br5AITEaknJYQG+MkP+lQeD7n/swhGIiLScEoIDTCgc1tSathAxywCwYiINFBYCcHMxpqZx8zurKGss5n928x2mtleM3vPzHo2XqjRafpvx1c7l6CMICItUMgJwcwSgCeA2rYJawPMAIYAXYD1wIsNDTDadU1P5YoxvQPOlZR7IhSNiEj9hdNDuB6YA6yoqdA5t8E591fn3C7nXCnwNDCmprpmll3XA0iv6XXR6rxhgfslfLdpH09PW0tuUWmEIhIRCV9ICcHMOgK3APeH8d7HAMvqE1RLc1i3dgHPyzyOP326it+/HxfNF5EYEWoPYSIwyTmXHUplM+sD/BG4p6Zy51xGXQ8gJ8S4okL71FbcfPLAaue1/7KItCRBE4KZjQBGAc+H8oa+3sQnwETn3JSGhddy3HbqoXx88/G1lj/+2SqO/9NXbNlb0IxRiYiELpQewjhgELDNzDKBS4C7zezZqhXNrA3wEfCuc+6pRo20BRjcvX2tZZOnrmXL3kLue29pM0YkIhK6UBLCc8AAYLjv8T7wJHCHfyUzSwbeAZY55yY0cpwtxg+Hd688Tq7hHoXM3OLmDEdEJGRBE4JzrsA5l1nxAAqBfOdcjplN8OspjAVOAy41s3y/R+9a3zwGPXzh0MqtNj0eV22vhMKSsghEJSISXFLwKoGcc1f5HT/sdzwdiPs7stKSk7jwqB68NHMjZR7Hsu25pLduVVleUFIewehERGoXdkKQ4Hp3SKs8PuepGQFlhUoIIhKltJZRE8hIS2bMIR1qLCsoVUIQkeikhNBEHjhvSI3ny7V5johEKSWEJjKoazsuHdUr0mGIiIRMCaEJ3X/uETUOHZ0x6WuKNHQkIlFGCaEJtU5OZOIFQ6udX5mZxweLt0cgIhGR2ikhNLFeB6XVeL5QPQQRiTJKCE0sOSmBJb8/jbvPOjzgvEcXl0UkyighNIP2qa247vh+Aec+XZbJIx+voLhMPQURiQ66Ma2ZWJVtNWev38vs9Xs5qE0yN4zrH6GoREQOUA+hGf3qxOof/Is2h7TFhIhIk1NCaEa3nzao2v7LFverP4lItFBCaEYJCcb4QZ0DzikhiEi0UEJoZqP7Bd6otnpnPou3eIeNcgpK+WrlTkrKPJEITUTinBJCM2uX2oqZd55U+XxtVj7n/+1btu4r4PJ/zOaal+bzly9WRzBCEYlXSggR0D2jNQ9XuYP5vUXbWbY9F4Bnp6+LRFgiEueUECLkhyO6Bzx/7LNVEYpERMRLCSFC0pKT6NI+JdJhiIhUUkKIoJtOGhjpEEREKoWVEMxsrJl5zOzOWsr/bmbrzMyZ2ZjGCTF2XT66N+MOPbjGsvzismaORkTiXcgJwcwSgCeAeXVUWwhcA2xtYFxxISHBePma0Yzsc1C1siVbdQeziDSvcHoI1wNzgBW1VXDOPeucmw5oxbYwTLpkeLVzHt2KICLNLKSEYGYdgVuA+xvjh5pZdl0PIL0xfk5L0atD9T0Tvly5MwKRiEg8C7WHMBGY5JzTOEYT+cvFwwKev/jtRgpL1NESkeYTNCGY2QhgFPB8Y/1Q51xGXQ8gp7F+Vktx4VE9efXa0QHn/jZ1bYSiEZF4FEoPYRwwCNhmZpnAJcDdZvZsk0YWh44feDBrJ55JuxTvNhXPTF9HuXZWE5FmEkpCeA4YAAz3Pd4HngTuqFrRzJLNLBUwINnMUq3qzjBSp6TEBO47dzAA5R7Hpj37IxyRiMSLoAnBOVfgnMuseACFQL5zLsfMJlTpKXzuK+8NTPcd92mKwGPZRUf1pH2qt5cwacoaMnOKcE49BRFpWhaNHzRmlp2enp6enR2/17CveWkeX63Mqnw+rFcG7914bAQjEpFol5GRQU5OTo7vWmzYtHRFlPrRyJ4BzxdvyWbm2t0RikZE4oESQpQ6c0hXHrkwcInsdxdti1A0IhIPlBCilJlx2eje9OrQuvLcnvwSvli+k7Jy3cYsIo1P1xCi3OqdeZz2xNcB5yZeMITu6a3JLSrl/OE9IhSZiESbhl5DUEJoAR7/bBWT/W5S69I+hZ25xQAM7taeiRcMYUTv6gvkiUh80UXlOHDzyYH7JlQkA4DlO3K54OmZzR2SiMQgJYQWIDkpgWUPnB7pMEQkxikhtBBtUpJISdKfS0Sajj5hWpAzh3SttayoVCujikjDKCG0IPeeM5ifjq15JZCduUUArMzMZXd+cY11RETqooTQgnRsm8Ifzh/C6z//QbWyzJwilmzN5oxJ3zB64pQIRCciLV1SpAOQ8B3Tv1O1c5c8N5vRfTsA4HHgnEMLzYpIONRDiCFzN+6tPC7R3cwiEiYlhBbq7V+O5YoxvTnl8C41lheVKiGISHg0ZNRCjezTgZF9OpBdUMKxj37F/ir7LxeXlQOtIhOciLRI6iG0cBlpySy877Rq55dty41ANCLSkikhxIDkpAR+OLx7wLmrX5pHWbmHfftLIhSViLQ0Sggx4uEqeycAXPnCXEZNnMLqnXkRiEhEWholhBiRlpzEf24YG3Bu1vo9lHkcpz3xNdkF6imISN3CSghmNtbMPGZ2Zy3lA81slpkVmNlsMxtYUz1pGkf3qX0J7IVbtJS4iNQt5IRgZgnAE8C8Oqq9BXwAdAA+BN5sUHQSFjPj5pMG1Fi2K8+7nMX27EJ+9s+5fLB4e3OGJiItQDg9hOuBOcCKmgrN7DCgP/CYc64I+BMw0MwGNThKCdltpw1i1UNnMKpvYG9hR7Z3raPf/XcJ01fv4qY3FkYiPBGJYiElBDPrCNwC3F9HtcHAKudcKYBzrgRY5Ttf9f2y63oA6WG3RCqlJCXynxuOCTj3xJTVzF6/h++35UQoKhGJdqH2ECYCk5xzdQ1EtwGqTn7PBdrWJzBpuBl3jOe4AQfWPbr0udnkFZVGMCIRiWZBE4KZjQBGAc8HqbofaFflXHsgv2pF51xGXQ9AX2MbQc+D0njtuh8w9pCOlec8fltoZ+UVRSAqEYlWofQQxgGDgG1mlglcAtxtZs9WqbccGGRmSQBmlgwc6jsvEfTA+UfUeH7iRzVeDhKROBVKQngOGAAM9z3eB54E7vCv5JxbCawHbjezFOB2YK1zblWjRixhO7RLO24+ufoM4O827wNgVWYe1708j/l+q6WKSPwJmhCccwXOucyKB1AI5DvncsxsQpWewqXAD4Fs4HzgsiaJWsJ226mHVju3ZW8hHyzezumTvmbKiizueXdpBCITkWhhzrngtZqZmWWnp6enZ2frZqrGVFxWzpNT1rAjp4j/LdxWY50Hzz+CK8f2bd7ARKRRZGRkkJOTk+O7Fhs2LV0RR1KSEvndGYfxqxP711rn3veW8f1WXdMXiUdKCHFoYJd2XDCiR63l7y3axpz1e7j+lfmszdLCeCLxQkNGceyjJTu48fXv6qzTu0MaX/9ufDNFJCINoSEjqbezj+xGrw6t66yzeW9BM0UjIpGmhBDn3rx+LP0PblNnnQWb9jVTNCISSUoIca5HRms+veUEzh3WvdY6P352ZjNGJCKRooQgtEpM4KnLRvDJ/x1fY7nHQX5xWTNHJSLNTQlBKh3erT2TLx9Bp7Yp1cqG3P8Zp/5lOgUl3sSwYNNeZqzZ3dwhikgTUkKQAOcc2Z3595xSY9marHzmb9xHTkEpFz0ziytemMOG3fubOUIRaSpKCFKjkbVsx1lQUsbW7AMzj1buqLriuYi0VEmRDkCi08DObWucXfTHT1fR86ADU1WTEvWdQiRW6P9mqdEdZxzGOUd249ELh3JY1wPbXGzYvZ9v/K4dzFq3JxLhiUgT0J3KEpLb/r2Id76reUG8H43syY3jB9AjozXJSfqOIRIpDb1TWQlBQrIrr5ib31jIrPW19wh6ZLTmy9+MI7VVYjNGJiIVtHSFNIuD26XwxvVjmHv3yYzqW/MF523ZhczZsJdo/JIhIsGphyD18sKMDTz4Yc27o3ZLT+WT/zuejLTkZo5KJL5pyEgiZsGmfVz0TO3LWgzp0Z4eGa05Y0hX3l+0nYcvHEq39LoX0xOR+lNCkIhasGkvu/NL+MWrC4LWHT/oYF68enTlc+ccP39lAcVl5bx89WgSEqwpQxWJebqGIBE1sk8HTj+iK3MnnBy07tLtgTexLduey5QVO/lmzW7mbNjbVCGKSIiUEKRRdG6fyoZHzuIP5x9Ra52qvdGi0vLK49JyT5PFJiKhCSkhmNkzZrbDzHLN7HszO6eWekea2bdmlmdmC81sROOGK9HMzLhyTB+G9UyvsdzjywfOORZtyeZHz86qLEswDReJRFqoPYQngX7OufbA1cBrZhYw99DMWgHvAi8DGcCjwLtmpqkmccTMeOsXY7n9tEMZe0jHgLK9+0vIKSjlyhfm8sO/fRtQdtWLc3nu63XNGaqIVBH2RWUzOxr4BhjpnFvud34I8KVzrovfuVXAr51zX4T5M3RROUbMXLubP362isVbQvtbbnjkLKwJegvlHkeC0STvLRItmu2ispk9bWaFwDxgCrCiahXfo+r7VxtUNrPsuh5AzWMO0uIcM6AT7914LEfXsnpqVd9vywn5vbdlF/LF8p14PHV/qdlfXMYJf5rKJc/N1k1zInUIOSE4534FtAVOBaa46v9nrQIKzexXZtbKzH4C9AfSGi1aabFeumZ08ErAeZO/Zeu+guAVgZMen8bPX5nP299trbPeu4u2sS27kLkb9rJnfwkA//hmPde/Mr9ywx8RCXOWkXOu3Dk3BTjFzE6vUlYCXAD8BMgEzsPbk6i2IppzLqOuBxD610RpEdqmJLHsgdP5zamHMuaQDozu16HWuif9eXrlDKT84jJ+/OzMGu+KLi7zzkx6Y+7mOn92uV8PouJrzEMfreDz5Tt5dvr6cJsiErPqO+00Ce+3/wDOue+cc8c65zoClwOHAsHvWJK40CYliZtOHsib14/llTp6DCVlHn76z7kAPDNtLfM27uOFGRsqh4byi8sCEkRZkCEj/3HM8ip1d2QXhtkKkdgVNCGYWVszu8L33yQzuwgYj/fCctW6Q80sxcza4Z1ltMg5t7Txw5aWLrVVIhPOOqzW8rkb9jJgwsfMXn/ghrX8kjJenrmRc/76DS/M2FB5vqw8yHUBvwvJpeWewB5DPWIXiVWh9BAccA2wFdgDTAAuc859b2Y/MbNP/OpeA2T56nYFrmrccCWWXH9CtU5mgDKPC9i17a9T1nD/+8vYuCfwGsPyHbl17u3s30MoKfdQUnbgJrhg15idc7oQLY3i3YXbeGbauqj+9xR0C03n3H7gpFrK/gX8y+/5rcCtjRadxLypt5/It2t386ORPUltlUhBiXdG0O78kmp1/+HXK6hq/OPTeOKSYZxzZHda1bGtZ0lZYEIA7y5wvTukkVhlLaWi0nLOeWoGXdun8uq1ozVlVeptf3EZt7y1CIChPdI5bmCnCEdUMy1dIRHVr1MbrhjTp3JTnbTkJN7+5TH8/tzBYb/XrW8tZuDdn/C3qWtrrVNa7qG4/MCSGW9/t5Xxj0/jtn8vqlb3oyU7WJuVz4y1u9mZWxx2PCIVcotKK4+3hDiLLhKUECTq9OnYhquO7ceaiWfSPjVoJ7aaxz5bxfmTZ3Dxs7N4f/H2gKmlJWUeikurr5v03qLt1c4VlR1IHE5XG6QB/HuldfVgIy16I5O41yoxgcX3n8Y1x/YjOSmBW04ZGPJrF2/NYe7Gvdz8xkIe/nhl5fl73l1KSYgL6fkP9ZaWKSFI/RX5fQmJ5n3HozcyEbxLTdx37mBWP3Qmt5xyKF//dnxl2fTfnhj2+63MzKt1GY0Lnv42YLjJ/+JfSXk5zjnW78qnrEpCieaLhBId9vv1UpMTo/dalBKCtCi9O6ax/A+ns/7hs+jTsQ3v3Xhs2O/x2Gerajy/cHM2j322ih05hdz59hLe/u7APZVFpR7+M38rJ/15esA9EDe9sZB+d33MN2t2hd+YCCsqLQ8Y25amU1hyYPgxMSF6P3ajNzKRWqQlJ1XurjasVwYbHz2bjY+ezZlDuob0+h05RXWWj33kK96ct4VFfj2JZdtz+N3bSwB4edYmwNsz+GCx99rDlS/MDbsd9ZFXVMqHS7aT18AP8nKP4+Q/T+foB6ewb3/1GV2Nrazcw+Sv1tSZOHfmFrEtxBsFg61fFW0K/BJC1Zsjo4kSgsSMx348jMmXj+DFq0c1+nvf8fb3Ac9Lyz0B48L1VVxWzqdLd7A3xA/lW99axK9fX8hv/7OkQT93d34x27ILKSn38PHSHQ16r6rW7crn1VkbAzZAen3uZh7/fHWtiXN/cRk/ePhLjn30q6C/i398s55hf/icWev2NGbYTcp/YkOZJ3o3g1JCkJjRNiWJc47szvhBnfnwpuPo2Kb2rTj+7+TQL1DX5Pb/LOateYFrKDnn+H5rDne9syRggb6CkrLKD8e1Wfn8e/4WfvTMTKatyuKxT1dxw2vfccnfZ7Fka3aN3x7fXbiNMQ9/yTdrdjFlRRYAny7LbFD8TXnZ4+Q/T+fe95bx5JdrKs9953eDYU3W7coPue5DH60gr6iMq19qnl5ZY/CfZRTNPYTw5/SJtABDeqQz446TmPjxcgZ2bsemPQWcN7w7c9bvoV+nNpx2RFeSEow/f7G6Xu//3qLt1aaq5haVcfnzs8krLuONuVtY//BZlJR7GP/4NAzjvV8fyyl/mV5Z/6oX51Uer8nK57zJ33LbqYdyc5VkVXFDU13DUmXlHpLCmM7ov2WpVVu1Prgvlu/k/cXbue+cwRzcLqXGOh8u2c4dZ3iXJ0nwu+mvpMxTbaaN/2dkqEPsjdFDay7FfgmhNNhSKxGkhCAxq3VyIg/9cGjAueG9DuwbctPJA8nMLeJfc7zf9K89rl/AGknhGvbA5wHPD5nwMScd1rnyprZnpgXfEe4vX6yuTAjbswtDuvth0ZZsLn9+NucP78EjFw4N/gIImHpbn3ssfv7KfADKPR6e/snIGusk+X2y+2+RWlBSRnJSYO+toQnK31crd/LB4h3ce85gOvj1Er3LkAQmp+bi30OoOkstmighSFx78PwhAHTPaM2N4wcwul8HfvFq4y3Q+9XKrMrjNVl5Ib1mV14xKa0SOObRr0Kq/8Hi7RSUlPPG3M3cduqhdGqbjJmRX1zG7HV7OG5gJ0rLPbRKTKi8I7ymm/OqKiv3cNnzs1mwaR+z7jqZLu1Tq9Wpa0OjJL8P3kS/hJBfXEZGWmBCyC8+MMZeXBbaB2aCeT/kl+/IpVeHNNqntgLgmpe8ySqvqJR//Mx7Pck5x0//OZcNu/fz6S0nkFdUytUvzuPsod24qYHDh6HwT8DBVueNJCUEiWsJCcbECw58qz79iK6sfuhMAFbvzCO1VQIXPj2T3KIynrtyJC9+u5FZ6+t3MfPbtaG9btn2HPbUsJZTbXIKD8w4evLL1XyweAeTLx/ByzM3MmVFFqcN7sK3a3eTmGDcddbhTFuVxQWcZGiCAAAP9klEQVQjelS+puo38k179vP6nM10apvCvI3e8fyL/z6L6X73gFTYsjdwVpD/7B//taH8h4H8Z9x88v0OEhMs4AOzsDS0TYsSzJi+ehdXvTiPQV3a8dmtJwSUV1xvAe/v6Js1uwE45pEvyS3y/oyVmXnNkxDq6CF8unQHC7dkc/tpgyJ+F7MSgkgVFePbQ3p4d3Kd/tvx7C0oof/BbTmqz0F8uHg7o/t1ZPLUNeQVlVV+0FRISUrgzCFdebeG5TBC4X9toS6z1u3hB/06BHz7fm22d/jL/3rD58t3Vh7f9Y53ttRnyw6cm/C/7+nXqQ1j+3cE4NS/fF3tbu5Ne2pff+e+95Zy4/gBdGmfGrDcR5LvBqx1u/J5Y+6WyvMVU2bX7crnl//6DoBbTzm0srygpJxt2YV0bJNc2aMB+HRpJlP9elxlHlc5DLdqZx5b9xWwbV/N01b9eyAVyaBCRe+pqpqudQD8ffo6MnOLuO+cwSEveBiQEKr0EG54zfs76J7emp8d0zek92sqSggiQRzUJpmDfGPRndqmcNWx/QAqx86LSsuZtiqLL5ZncUT39lwyqhdtUpK49rhDOHfyjCaL67LnZzfqe2189Gw+XZpZ69IezjnMjP/M3xJw/pVZm3hl1iYeOO8Ijh1wYBXPihuwHv4ocPv1VZn5jOzTgdWZB4bQvl17IKnOWreHu/+3lGP6d+T1n49hbVY+f5++jv8sqL5V6pwNB/bLOO6PU2ttX25h7b2OXXnF7M4v5s15W/j1+AF0z2jNpCmreXraOl742dEM6NyWg9K8ySm7oIRHPvEuhXL8wE6cdFiXWt+3wpcrdjLZ7w54/4Tg36NamZkb9L2amhKCSAOltkrkjCHdOGNIt4DzQ3ums+i+U0lLTiK/uIzUVgkMvu+zCEUZ3LlPzajzmkC/uz6u8/X3v78s4PniLdk89OFyvvT7Vg/eLU8vG90rYOho7sYDH+wfLvHeFzFz3R5Kyz0BM7Pqq647stdm5Vfu0Pf16l3MuOMkJk3xTpm98oW5mMGRPdJ579fH8V+/pPTUV2sZ2LkdvToc2Dbe43Es3LKPAZ3b4Zwjt7CMa1+eH/Dz9uQXc8akr1mZmcejfpMAomF5dYvGdVjMLDs9PT09O7vmNWdEWqq+d34EwDlHduPxHw9j+updPPjhcob3yuCy0b158MPlrPR9c26XmsSgLu2YH2Refkt09tBufPR98BvihvfKCLhjPFw9D2pNr4PSwrruM6J3Bgs3V/+Zi+47leF/+CLg3MHtUvjo5uP465drOH7gwTz/9Xrmb9rHoC7tWLUztEkEFa4Y07varLhwZWRkkJOTk+Pbmz5sSggizWjJ1mw+WZrJDSf0Jz2tVUivKS338PaCrQzvnUHHNim0SUkkMcF4ZeYm8opK+etX3uGIKbeNo01KIq/P2cyrszeRXXDgW/HxAzvxyjWjWbUzjzMmVdv9tlncfPJA/up3s5pUt+GRsxrUU1BCEJFardiRy/pd+zn58M6VF2inrsxiTVZe5bLgndqmsDu/aTcAmnHHeDbvLeDy5+c06c+JBamtEvjXdWMY2eegsF+rhCAi9ebxuIAbtXKLSpm6Mov/e/PADnKvXjuavftLGNYzg+4ZrUlMMHbkFLJmp3c3uZ+O7cPNbyxk8dYczh7ajTvPPIyT/zKd1KQEume05pwju/Hrk7xTOzfs3s8Hi7fz3wVb2bw3cObSNcf248dH92TbvkK+XLmzcmbS4vtO4/g/fUVuURlXH9uX1Tvzgk7h/UG/DgEXnFuavh3TmHLbuLDuPodmSghm9gzwQ6ANsAm4yzn3YQ31jgKeBgYDmb56b4cdlBKCSMQ455i1bg9tUpLolp5K5xpuSKsqv7iMTXv2M7hb+5CHPIpKy1m9M48h3dMpKfcETDHNKSjlgQ+Wce6w7ow/rDPOOfbsL6FTW+8yGVv3FTBt1S7ueXcpw3tl8PAFQznrr96hsOTEBFZPPJNyj2PKip0s2LSP575eX/nefTumccmo3hzTvyO3vrWI9bv3A/DYj47kqa/Wcv0Jh3DFmD6szcqv84L23Akn85cvVvPNmt3VVmm98Kge9DoojWmrsli89cCF+gfPP4J73/NefH/t2h+wJiuPBz5YTlWf3nI8h3VtH9Lv0V9zJYTDgI3OuSIzOxqYAvRzzu2rUm8p8CrwGHA88DHQxzm3u+p7Bvl5SggiElReUSntfHcoL9i0j79NXcttpx5aeQ9JhQWb9rEyM5fLR/cOSFhrs/J59JOVXHtcv8r7MPxt2rOfV2dt4sqxfejTsQ2b9xRw4TMzOXtoVx7w3eUOMG1VFl+v3s34ww6mtNwTMB31ljcX8u6i7XRLT2XWXSdX+xlrdubx0syNXHf8IfTr1KZyem99NPuQkS8hfAOMdM4tr1KWBxzpnNvge74ZuMg5F9qdNgfeRwlBRKJSuB/Y2QUl/HfBVs4c2o0eGa2bMLKGJ4SQ70Mws6eBq4FU4ENgRQ3VJgOXm9kfgXG+c8uqVjKzYJ/06UHKRUQiItxv7xlpyVx3/CFNFE3jCvmKhXPuV0Bb4FRgiqu5a/EJ8HOgCO9w0c3OudrveRcRkagR1iVs51y5c24KcIqZne5fZmYdgA+AW4EUYBQw2cyq3WnhnMuo6wHUfrukiIg0ifourZcE9K9yrj+Q65z7ny9xLAFmAyc2ID4REWkmQROCmbU1syt8/00ys4uA8XgvLPtbDbQzs3PNawhwAvB91fcUEZHoE0oPwQHXAFuBPcAE4DLn3Pdm9hMz+wTAOZcDXAI8COQCHwF/cs5Na4rARUSkcelOZRGRGBGrS1d4AEtP1+xTEZFQ5eTkADjnXL2uD0drQijDO5xVnx0jKrJIPM1UUpvjg9ocHxrS5vaAxzlXr71uojIhNETFTW/17TK1RGpzfFCb40Mk2xzZHZ1FRCRqKCGIiAighCAiIj5KCCIiAighiIiIjxKCiIgASggiIuITc/chiIhI/aiHICIigBKCiIj4KCGIiAighCAiIj4xlRDMrLOZfWpmBWa21MzGRDqmhjKzFDN70cy2mlmOmU01syP8yu8xs92+x71VXnuOma03s3wze9PM0pq/BQ1jZmPNzGNmd/qeJ5rZU77fxXYzu65K/Wt953PMbLKZtZh/476dBu/xxZ9nZrN852O5zUeZ2WwzyzWz1b4dGWOqzWb2gJkt9/07vtTvfL3baGY/8H3GFZjZZ2bWuVGCdc7FzAN4B3gaaA38DNgGpEQ6rga2qQ1wL9ATSARuA9b4ys4DNgC9gT7AJuBcX1lXIBs4BWgHfAj8OdLtCbPtCXj35Z4D3Ok7dzMwD+gIDAf2AsN8ZcPx7uo3zFe+ALgp0u0Io703A18APXxtHxkHbV4K3OFr7zhgP9ApltoMXAGc6vu3fGmVv3fYbQRSfZ9tV/g+654F3m6UWCP9y2rEX3o7oATo4nduJXB6pGNr5HYmAx7fP5I3gd/5ld0JvOE7vgH42K9sDLAj0vGH2dYbgCeBl/wSwmzgYr86zwKP+I4fBZ72K7sUmBXpdoTY1kRgB9CvhrKYbLMv3jz/NgObgVGx2GZgWpWEUK82AmcAy/3KugPFQNuGxhi13ax6GAhkO+d2+p37HhgcoXiaylggyzm3B2/blvqV+be3prKuZnZQs0TZQGbWEbgFuL9KUbhtbil//154v/ldZmY7fUMMF/vKYrXNAJOBy80sycxO9p1bRmy3uUJ92xhQ5pzbjnczsYENDaheu+pEqTZU32EtF2gbgViahJllAM8BE3ynqrbZv71tgKyKAufcfjMr95Xva/poG2wiMMk5l21m/ueDtbm2smjXHcjAO/TXB+9QwWdmtpjYbTPAJ8ArwANAOXCJc67AzGK5zRXq28Ym+6yLpR7CfrzDRv7aA/kRiKXRmVkq8B7woXPun77TVdvs396AMt//YIm0gN+HmY3AO2zwfA3FIbeZlvX3L/T9d6Jzrsg5Nwf4HO/Yc0y22cw6AB8AtwIpeP/mk81sKDHa5irq28Ym+6yLpYSwBjjIzLr4nRsCLI9QPI3GzBLxXi/YAtzuV7Qcbxsr+Le3prJM51xL6B2MAwYB28wsE7gEuNvMniX8NreUv/9qoLSWslhtc38g1zn3P+dcuXNuCd5x9ROJ3Tb7q28bA8rMrBvehLCmwRFF+kJLI1+0eQfvmGQqcCUxMMvI165/Ap8CraqcPw9Yj3eWUW+8M478ZxntA07C25V8nxYyywhI88Vf8XgLeAjv5uM345111AHvsMoeAmdm7AaOxHvRfR5RPPukhna/iXeWXDIw0vf3OyxW2+z7e2YD5wKG90MuC29CiJk2A618n0lf+z6XUvF+Ga9XG32v3w5cjneW0TNollGNv/guwGdAAd4LU2MiHVMjtKkP4PAOKeT7PY73ld/r+4e0B7i3ymvP9SWJ/b4P1bRIt6eev4OXODDLKBFv0s/BOyvnuip1r/OdzwH+BiRGOv4w2tkRb+LOx/tt7+I4aPPpwCK8s402AbfHWpt9/35dlceJDWkj3lmDy3yfdZ/jN7uyIQ+tdioiIkBsXUMQEZEGUEIQERFACUFERHyUEEREBFBCEBERHyUEEREBlBBEGp2Z9TUzZ2adIh2LSDiUECSmmdk0Myv2bRLk/0iMdGwi0UYJQeLB3c65tlUe5ZEOSiTaKCFI3DKzjb4tK+f5tqycZmb9/co7m9lbZpbl28J0km/V2YrygWb2vpntMrO9ZvZOlR9xtpmt8W2B+G/firMV26K+4Htdrm/vgxObp9UitVNCkHj3c+AqoDOwAnjHDmzA8Ibvv/2B0cCxwMNQuZz4FLwblfTFuwjfU1Xe+1zgaLwblwwDfuk7/zNgBHAo3gXezsW7U5hIRCkhSDx40Myy/R4z/cr+7pxb5pwrBH6LdzeqYWbWA+9Ksbc55/Kcd1eq+/EmD4Bz8G5lerdzbr9zrsQ5N7XKz33AOZfjnMvCu2jdSN/5Erwr0B4OJDjn1jnn1jdBu0XCooQg8eBe51yG3+MYv7JNFQfOuXy8q8b2AHoCZc65bX511+PdcyMN7yq061zdq0Pu8Dv239TkNbxLmv8N2G1mr5lZ1/o2TqSxKCFIvOtTcWBmbfEuQb0N2Aok+XoKFfoC+5xzBXgTySF+w0shc86VOecedc6NwDuc1B74U/2bINI4lBAk3l1vZoN9F4sfxXsdYYmvZzAVeNzM2vl2pXoAeNn3uo/wbnzyBzNrY2bJZjY+lB9oZieZ2XAzS8K7nn0R3v2ERSJKCUHiwcQa7kPo7St7Ae+H/C5gKHCRc87jK7scSMI7VDQf7+5Wd0Pl8NIpwFF4LwhnAjeGGE8XvBess/Fui5oE3NmwJoo0nDbIkbhlZhvx7tD130jHIhIN1EMQERFACUFERHw0ZCQiIoB6CCIi4qOEICIigBKCiIj4KCGIiAighCAiIj5KCCIiAsD/A9ies/MwJF/lAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(hist.history['loss'][200:])\n",
    "plt.xlabel('Epochs')\n",
    "np.mean(hist.history['loss'][1000:]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing with R\n",
    "\n",
    "```\n",
    "m1.zi = zeroinfl(daysabs ~ math + prog | math + prog, data = dat)\n",
    "summary(m1.zi)\n",
    "logLik(m1.zi)/314\n",
    "Call:\n",
    "zeroinfl(formula = daysabs ~ math + prog | math + prog, data = dat)\n",
    "\n",
    "Pearson residuals:\n",
    "    Min      1Q  Median      3Q     Max \n",
    "-2.3410 -1.1702 -0.6026  0.5864  7.3907 \n",
    "\n",
    "Count model coefficients (poisson with log link):\n",
    "                 Estimate Std. Error z value Pr(>|z|)    \n",
    "(Intercept)     2.5879178  0.0614173  42.137  < 2e-16 ***\n",
    "math           -0.0052060  0.0009291  -5.603 2.10e-08 ***\n",
    "progAcademic   -0.3040946  0.0567438  -5.359 8.36e-08 ***\n",
    "progVocational -0.9397286  0.0794916 -11.822  < 2e-16 ***\n",
    "\n",
    "Zero-inflation model coefficients (binomial with logit link):\n",
    "                 Estimate Std. Error z value Pr(>|z|)  \n",
    "(Intercept)    -1.920e+01  1.706e+03  -0.011    0.991  \n",
    "math            1.264e-02  6.926e-03   1.825    0.068 .\n",
    "progAcademic    1.674e+01  1.706e+03   0.010    0.992  \n",
    "progVocational  1.766e+01  1.706e+03   0.010    0.992  \n",
    "---\n",
    "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 \n",
    "\n",
    "Number of iterations in BFGS optimization: 13 \n",
    "Log-likelihood: -1192 on 8 Df\n",
    "> logLik(m1.zi)/314\n",
    "'log Lik.' -3.797503 (df=8)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([[-0.00513406],\n",
       "        [-0.2881807 ],\n",
       "        [-0.9255564 ]], dtype=float32),\n",
       " array([2.5703218], dtype=float32),\n",
       " array([[-0.00844655],\n",
       "        [-0.31162283],\n",
       "        [-1.2921721 ]], dtype=float32),\n",
       " array([2.6434429], dtype=float32)]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.get_weights()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using negative binomial (Under construction)\n",
    "\n",
    "* Problema in TF only has the default parametrization of the NB. Which is not suited. Therefore we stack to distributions (using the gamma-poisson case). Much time has been wasted to construct a single NB distribution from scratch. \n",
    "\n",
    "* Question. It seems that in NB there is a globa parameters $\\Theta$ it's not modelled depending on data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_gamma(t):\n",
    "    mean = t[...,:1] #First column mu\n",
    "    #disp = 1e-3 + tf.math.softplus(0.05 * t[...,1:]) #Second column sigma\n",
    "    disp = t[...,1:]\n",
    "    disp=1.03\n",
    "    eps = 1e-10\n",
    "    # translate these into gamma parameters\n",
    "    gamma_shape = 1. / (disp+eps)\n",
    "    gamma_scale = (mean / (gamma_shape+eps))+eps\n",
    "    return tfp.distributions.Gamma(gamma_shape, gamma_scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEXdJREFUeJzt3XusZWV9xvHvw4BMOpczQqvWiFxsNR2nohgTIAGMRq1StIkRkaQtRFQMWGwyCSPVNsSSTMFUDFgNWC/UkCmNF0RU1HqJqaKkqcgwUkgBbyimI3MOMwOtML/+sdfQ7eGdmTNn7XP27DPfT7IzZ71r7f3+1lz2M++73r12qgpJkmY7ZNwFSJIOTAaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmg4ddwFzkeQxBmE2M+5aJGlCrAZ2VdW83+czCd8HkWQXkKmpqXGXIkkTYXp6GqCqat4zRRMxggBmpqamprZt2zbuOiRpIqxZs4bp6elesy5eg5AkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpaVI+KNfLMRtufuLn+zeePsZKJGlyOIKQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpKZ9BkSSS5NsSbIryVlD7RcnuSvJw0nuTPK6Wc/74yT3JtmeZFOS31qIE5AkLYy5jCDuAS4Cvjer/XHgDcAUcD7wiSTHAiR5BvBJ4K3A7wIrgfeOqGZJ0iLYZ0BU1Ser6ivAo7Pa31dVd1TVrqr6FnAncEK3+0+Ab1fVV6vqYeBvgbNHXLskaQGN5G6uSVYD64AtXdNaYPPQIXcAz0jy1Kp6qPH8bfvoYmoUdUqS5q73ReokAa4BbqyqH3bNK4CZ3cdU1Q4GU1Ir+/YnSVocoxhBbASeBbxiqG0HsGr3RpIVwDJge+sFqmrN3jroRhiOIiRpEfUaQSRZD5wOnFFVjwzt2sJgymm3dcAvWtNLkqQD01yWuR6WZHl37GFJlic5JMm5wDuAP2q88X8WODnJy5KsBP4KuH7UxUuSFs5cRhDXAo8ApwDXdT+fCvw1gyWsd3Wfddie5BKAqvoF8GfAPwIPds95z+jLlyQtlH1eg6iqc4BzGruO3cfzbgJumldVkqSx81YbkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpaSR3c10Kjtlw8xM/37/x9DFWIkkHBkcQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNcwqIJJcm2ZJkV5KzhtqXJbkqyXSSB5KcN+t5b+7ap5NcncRAkqQJMdc37HuAi4DvzWq/ADgROA54DXB5kuMBkrwQuBx4dbf/pO54SdIEmFNAVNUnq+orwKOzdp0NXFFVW6vq+8ANwO4RxlnAP1fV7VW1FbiiO16SNAH6fqPcWmDz0PYdwCuH9n1p1r61rRdJsm0f/UzNt0BJ0vz0vSawApgZ2p4BVs5hnyTpANd3BLEDWDW0vRrYPod9v6Gq1uytk26E4ShCkhZR3xHEFmDd0Pa6rm1f+yRJB7i5LnM9LMny7vjDkizvlqxeD6xPckS3eulMYFP3tE3AG5O8IMmRwPrueEnSBJjrCOJa4BHgFOC67udTgQ8CtwH3MbggfXFV3Q7QrWraANwC3Mtgiew/jLJ4SdLCmdM1iKo6BzhnD7sv7B6t530E+Mh8Clsox2y4+Ymf7994+hgrkaQDm59sliQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1NT3+yAm2vB9mSRJv8kRhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ19Q6IJCckuTXJTJK7k7y+a1+W5Kok00keSHJe/3IlSYtlFPdiug74J+Bk4BTgC0m+CZwNnAgcBxwFfC3JbVV1+wj6lCQtsFFMMR0N3FBVu6rqm8BW4FgGAXFFVW2tqu8DNwBnjaA/SdIiGEVAXA2cneTQJC/v2u4E1gKbh467o2t7kiTb9vYApkZQpyRpP4wiIL4IvAV4FPgC8BdVtRNYAcwMHTcDrBxBf5KkRdDrGkSSI4CbgHOAzwHPZ3AN4r+AHcCqocNXA9tbr1NVa/bRj6MISVpkfUcQzwFmquozVfV4Vf0AuBV4KbAFWDd07LquTZI0AfoGxN3AqiRnZGAdcCqD6w3XA+uTHJHkeOBMYFPP/iRJi6TXFFNVTSd5I/B3DALhV8DlVfWNJN8CngvcB+wELnaJqyRNjt6fg6iqW4BbGu2PAxd2D0nShPFWG5KkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkpoMCElSkwEhSWoyICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU29AyID707yQJKHk3yna1+W5Kok092+8/qXK0laLIeO4DXeAZwGvAT4OfCirv0C4ETgOOAo4GtJbquq20fQpyRpgfUKiCTLgHcBJ1fVz7rmf+9+PRu4oqq2AluT3ACcBUxUQByz4eYnfr5/4+ljrESSFlffKaajgOXAm5I8mGRLkjO7fWuBzUPH3tG1PUmSbXt7AFM965Qk7ae+U0zPBNYAR3eP44FbktwOrABmho6dAVb27E+StEj6BsQj3a+XVdWjwHeTfBl4BbADWDV07Gpge+tFqmrN3jpxFCFJi6/vFNPdwK/3sG8LsG5oe13XJkmaAL0Coqp2AJ8GNiR5SpIXMxg9fBW4Hlif5IgkxwNnApv6FixJWhyjWOZ6AfAx4FcMlrm+raruSnIP8FzgPmAncLFLXCVpcvQOiG4Z62sb7Y8DF3aPiTK8tFWSDlbeakOS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpKZRfKPcQW/4C4bu33j6GCuRpNFxBCFJajIgJElNBoQkqcmAkCQ1eZF6PwxfjJ7LMV6wljTJRjaCSHJSkl1JNnTby5JclWQ6yQNJzhtVX5KkhTeSEUSSQ4D3A7cNNV8AnAgcBxwFfC3JbVV1+yj6lCQtrFFNMb0V+C4wNdR2NnBFVW0Ftia5ATgLMCAkaQL0DogkRwLvZDBauHJo11pg89D2HcAr9/Aa2/bRzdQ+9kuSRmwU1yAuA66sqtlv8iuAmaHtGWDlCPqTJC2CXiOIJC8CXsLgesNsO4BVQ9urge2t16mqNfvoZxuOIiRpUfWdYjoNeB7wsyQweBN/LMkxwBZgHfDD7th1XZskaQL0DYhrgE1D2x8A7gGuAP4cWJ/kXxmsYjoTeFnP/iRJi6RXQFTVTmDn7u0kjwDbq2o6yQeB5wL3dcdc7BJXSZocI/0kdVWdM/Tz48CF3UOSNGG8F5Mkqcl7MY2B92uSNAkMiEUylxv9SdKBxCkmSVKTASFJajIgJElNBoQkqcmAkCQ1GRCSpCaXuS4gl7ZKmmSOICRJTQaEJKnJgJAkNRkQkqQmA0KS1GRASJKaDAhJUpMBIUlqMiAkSU0GhCSpyYCQJDUZEJKkJgNCktTUKyCSHJ7kY0l+mmQ6ydeTPH9o/7uT/Hf3eE//ciVJi6XvCOJQ4F7gROAI4CbgswBJXgu8GTgBeDFwXpIzevYnSVokvb4Poqp2AO/dvZ3kauB9SY4EzgY+VFU/7vZ9qGu7afbrJNm2j66m+tQpSdp/o74GcRLwy6raCqwFNg/tu6NrkyRNgJF9o1ySNcA1wCVd0wpgZuiQGWBl67lVtWYfr70NRxGStKhGEhBJlgM3Ap+vqo92zTuAVUOHrQa2j6K/pWRPX0t6/8bTF7kSSfpNvaeYkiwDNgE/AdYP7doCrBvaXte1SZImwCiuQVwLLAfOraoaar8eeHuSZyd5NnB+1yZJmgC9ppiSHA2cCzwKPJRk965XV9XnkhwP/EfXdmVVPWkFkyTpwNR3meuPgOxl/3sZWgYrSZoc3mpDktQ0smWuGq3h1U2uaJI0Do4gJElNBoQkqcmAkCQ1GRCSpCYDQpLUZEBIkppc5jqBvMGfpMVgQEyAPQXCXI8zOCTNh1NMkqQmRxAHAT+VLWk+DAg9iVNVksApJknSHhgQkqQmA0KS1OQ1iIOMF6wlzZUBIWDun7WQdPBwikmS1OQI4iC2v6OGPU1P7W+7pMlgQGhRLEZYGEjSaDnFJElqWvARRJKnAdcBpwL3AudV1a0L3a8W1kJc1F7sEYAjDmnvFmME8WEGwXAkcAXwqSSHL0K/kqQeUlUL9+LJKmArcFRVPdi13QVcVFW3DB23bR8vNQUwNTU1rzpmHn1sXs/Twli9/P8HrnP5s9nT8cPte3ut2cfN5bX255g+xx8s/H3ZP6P4/Zqengaoqpr3QGChA+IE4EtV9bShtn8Bvl1V7x9qm0tAFDAzjzJ2p8r0PJ57IFpq5wNL75yW2vnA0jung+F8VgO7qmreibzQUb6CJ7+pzwArhxuqas1CFbA7fBayj8W01M4Hlt45LbXzgaV3Tp7P3Cz0NYgdwKpZbauB7QvcrySpp4UOiHuApyZ5+lDbOmDLAvcrSeppQQOiqh4GPg+8J8nyJH/KYATxjYXsV5LU32IsJ3g7g89B/Aq4D3h9Vf3PIvQrSephwQOiW976qoXuR5I0Wt5qQ5LUZEBIkpoW9INykqTJ5QhCktRkQEiSmgwISVKTASFJalrSAZHkaUm+lGRnks1JThx3TfOV5PAkH0vy0yTTSb6e5PnjrmtUkpyUZFeSDeOupY8MvDvJA0keTvKdcdfUV5ITktyaZCbJ3UleP+6a9keSS5Ns6f5+nTXUvizJVd2/pweSnDfOOudqL+dzcZK7ur93dyZ5Xd++lnRAsLS+rOhQBudyInAEcBPw2bFWNCJJDgHeD9w27lpG4B3AacBLGNyC+cLxljMS1wGfAdYAbwGuS/Lb4y1pv9wDXAR8b1b7BQz+PR0HvAa4PMnxi1zbfOzpfB4H3sDg7935wCeSHNunoyW7zHWuX1Y0qZI8BXgU+J2q2jruevpIcj7wBwz+Yt9VVRvHXNK8JFkG/BQ4uaruG3c9o5LkYeAFu88pyY8Z3DJnogI9yTeAD1fVpm77VuDvq+qGbvvDwENV9a7xVTl3s8+nsf/fGJzfp+bbx1IeQfw+sG13OHTuANaOqZ5ROwn45RIIhyOBdwJ/M+5aRuAoYDnwpiQPdtMAZ467qBG4Gjg7yaFJXt613TnOgkZkLbB5aHvJvD8kWc0I7py9lL/7b05fVjSJkqwBrgEuGXctI3AZcGVVbUsy7lr6eiaDaZiju8fxwC1Jbq+q/xxrZf18kcE006UMpjHeWFU7x1vSSMx+j1gq7w9h8P5wY1X9sM9rLeURxJL8sqIky4Ebgc9X1UfHXU8fSV7EYK7+2nHXMiKPdL9eVlWPVtV3gS8DrxhjTb0k2X296y+Bwxn8eV2d5A/HWthozH6PmPj3h85G4FnA2/q+0FIeQTzxZUVD00zrgI+MsaZeujnuTcBPgPVjLmcUTgOeB/ysGz1MAY8lOaaqzh9rZfNzN/DrcRcxYs8BZqrqM932D7q5+5cymJKZZFsYvCfs/l/2xH+ZWZL1wOnAKVX1yL6O35clO4JYol9WdC2DOe5za2msLrgG+D3ghd3jc8AHgIvHWdR8VdUO4NPAhiRPSfJiBqOHr463sl7uBlYlOaNbwrsOOJUJCockh3Uj70OAw7r3g0OA64H1SY7oVi+dyeA/YAe0PZ1PknMZrKL7o6p6aCSdVdWSfQBPB24BdjK4qHbiuGvqcS5HA8VgGmP70OOUcdc2wnP8OLBh3HX0PIcjGQTddgaj2DPHXdMIzulVwPeBh4EfAevHXdN+1v/x7t/O8OOlwDIGF+CngZ8D54271p7ncx/wv7PeHy7p09eSXeYqSepnyU4xSZL6MSAkSU0GhCSpyYCQJDUZEJKkJgNCktRkQEiSmgwISVKTASFJavo/7HdLcJRH1kwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "my_gamma(np.ones((1,2))).sample(1000)\n",
    "plt.hist(my_poiss(my_gamma(np.ones((1,2))).sample(1000)).mean().numpy().reshape((1000,1)),100);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Definition of the model\n",
    "#Using the functional API for clarity\n",
    "inputs   = tf.keras.layers.Input(shape=(3,))\n",
    "t        = tf.keras.layers.Dense(2, \n",
    "                             activation=tf.exp, \n",
    "                             bias_initializer='zeros', \n",
    "                             #bias_initializer=tf.keras.initializers.Constant(value=2.61),\n",
    "                             kernel_initializer='zeros')(inputs)\n",
    "gamma    = tfp.layers.DistributionLambda(my_gamma)(t)  \n",
    "p_y      = tfp.layers.DistributionLambda(my_poiss)(gamma) \n",
    "#Glueing input and output together. Note that output is a tf.distirbution \n",
    "model_nb = tf.keras.models.Model(inputs=inputs, outputs=p_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_nb.compile(optimizer=tf.optimizers.Adam(learning_rate=0.001), loss=negloglik)\n",
    "#model_nb.predict(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist = model_nb.fit(X, y, epochs=400, verbose=False);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 226,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.912700449239676"
      ]
     },
     "execution_count": 226,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEbCAYAAAA8pDgBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXe8FNX5/z/PltsbvaM0K9UGitj92WssWPLVxBKNGmNi/JoYo8YYE2O+1iRGjRVbNIkaY1RQUSwoIqCIKNIEAWm3cOu28/tjdnbPnD1Ttt3dvfd5v168uDszO3N2duc85+kkhADDMAzTO/EVegAMwzBM4WAhwDAM04thIcAwDNOLYSHAMAzTi2EhwDAM04thIcAwDNOLYSHAMAzTi2EhwDAM04thIcAwDNOLYSHAMAzTi2EhwDAM04thIcAwDNOL8SQEiOgmIlpGRDEimqnZX09Em4joldwPkWEYhskXXjWBFQCuBPChzf6bAHyVkxExDMMw3UbAy0FCiFkAQETXqfuIaDyAAwDcB+CMTAZBRBEYAqklk/czDMP0UuoAxIQQnuZyHRm/UeJuAD8FMMbuACJqcjmHHwDq6+vrczAehmGYXkFzczOQpW83KyEQ9w9sFkLMIyJbIeCF+vp6NDW5yQqGYRjGpKGhAc3NzVlZUDIWAkRUDcMXcITbsUKIBpdzNQFgLYBhGKabyUYTGAdgFIAFRAQAlQAqiOgTIcTEXAyOYRiGyS+ehAARBWHY7X0AgkRUAWAZgJHSYWcCOCn+P8MwDFMCeHUoPACgA8AMAI/F/z5ACLHJ/AegGUBICLElP0NlGIZhcg0JIQo9BhBRU319fT07hhmGYbwTdww3u/ldneCyEQzDML2YXOQJFJSFaxuxcnMrhvepxAFj+xd6OAzDMCVFyWsC//h4Pa75xyd44oOvCz0UhmGYkqPkhUDQRwCAcDRW4JEwDMOUHqUvBPzGR4jECu/gZhiGKTVKXggE4kKANQGGYZj0KXkhEPQb5qBIlDUBhmGYdCl5IRDwmeYg1gQYhmHSpfSFgN90DLMmwDAMky4lLwQS5iDWBBiGYdKm5IVAwhzEmgDDMEzalLwQCPo5T4BhGCZTSl4IBDhPgGEYJmNKXwj4OESUYRgmU0peCJgZwyE2BzEMw6RNjxECERYCDMMwaVPyQiDAGcMMwzAZU/JCIBEdxHkCDMMwaVPyQoDzBBiGYTKn9IVAImNYoBj6JTMMw5QSrkKAiG4iomVEFCOimdL2i4hoFRG1ENHXRPTz/A5Vj+kYBjhXgGEYJl28aAIrAFwJ4ENl+2wAewkh6gDsB+AcIjomx+NzxcwTANgkxDAMky6ujeaFELMAgIiuU7avkV4SAAFglO4cRNTkcpl6t3HYIWsC4VgMlfBneiqGYZheR1Y+ASI6m4h2ANgAoArAP3IyqjQwfQIAawIMwzDp4qoJOCGEeBLAk0Q0HsCpAFptjmtwOk9cU8hIGzCjgwBOGGMYhkmXnEQHCSGWAmgD8MtcnC8dgpImEGbHMMMwTFrkMkQ0AGBMDs/n7aKyTyDCmgDDMEw6eAkRDRJRRfzYIBFVEJGPiP6HiAaRwWQAlwF4M98DVpE1Ae4uxjAMkx5eNIEHAHQAmAHgsfjfBwHYF8AnAHYAeA7AXwHcl59h2hOUfALcZ5hhGCY9vISIng/gfM2uuQCuyO1w0oejgxiGYTKn5MtGqHkCDMMwjHdKXghwxjDDMEzmlLwQ8FuEAGsCDMMw6VDyQoCIpJ4CrAkwDMOkQ8kLAUDuKcCaAMMwTDr0DCFgagLsE2AYhkmLHiEEEs3mOTqIYRgmLXqEEDAjhMJsDmIYhkmLHiEETE2AzUEMwzDp0UOEQLzPMAsBhmGYtOgRQiDAPgGGYZiM6BlCwMfRQQzDMJnQM4RA3BwULRJNQAjBTmqGYUqCHiEE/GayWJFkDF/8+EJM/e3r2LKjq9BDYRiGcaRHCIFg3BwULRJz0Oxl32J7Wwj3vbWy0ENhGIZxpEcIAbOIXLHVDoqJ4hoPwzCMSo8QAsXmEzAhkPtBDMMwBaRnCIEi8gnEpDEQywCGYYqcHiIEiidZLCqZgFgGMAxT7LgKASK6iYiWEVGMiGZK2/+XiJYT0Q4i+oyITsrvUO0xfQLRItAEoqwJMAxTQnjRBFYAuBLAh8r2KIDTAdQDuATAo0Q0KrfD80YxVRGVncHEUoBhmCLHVQgIIWYJIWYD6FS23y6E+FQIERNCzAPwGYC98jROR4pWEyjgOBiGYbwQyMVJiKgOwHgAy2z2N7mcoj6b6xdT2QiLMsJSgGGYIidrxzAZNo/7AbwghPg8+yGlTzJEtPBCwOoYZinAMExxkwtN4HcAhgM40u4AIUSD0wnimkLG2kAxlY3gBDGGYUqJrIQAEV0N4DgAM4QQHbkZUvokQ0SLwDHM0UEMw5QQrkKAiIIA/DBMR0EiqgAQAnAegCsATBdCNOZ1lC6Y5qBi0AQ4T4BhmFLCi0/gAQAdAGYAeCz+90EAfgVgCIDlRNQa//eLvI3UgUCxRgexFGAYpshx1QSEEOcDOF+zqyA5AToCiR7DxWAOSv7NjmGGYYqdHlU2wkkTePWzTTjsj3Px/spteR2LxRzEMoBhmCKnRwgBM1nMySfwg8cXYtWWNpz1wPy8joWTxRiGKSV6hBAwy0Y0toXwt3dWY9329oKNxRIiyqoAwzBFTo8QAqYmsGJzK25+aRmOu3tewcbCmgDDMKVEjxACpk/ApKUzUqCRcHQQwzClRY8UAoUkxmUjGIYpIXqEEPD7i+djFEGqAsMwjGeKZ/bMgmLSBNgcxDBMKcFCwIVQJL0EtBiXjWAYpoToGULAn95067W8xN/eWY3xN7yK2cu+zejcrAkwDFPs9AghYJaStkMo5Z1bu7xFD9380jKEojFc9NhHnscSY6cAwzAlRI8QAkEXc1BIqSnU5lEIZIJcNoLlAcMwxU6PEAJ+FyHQGe5GIRCThQBLAYZhipseIQTcfAJd4ajltVdzUCbEWBNgGKaE6BlCwMUnkKoJRG2OzB6L5SkHmoAQAu2hwmVAMwzTs+khQsBFE4h0nyZgNQdlf77Ln1qEvW6eja82t2Z/MoZhGIUeIQTS9Ql0hvOnCVjNQdlLgf98shGd4Rhu+c+yrM/FMAyj0iOEQEBTNkJekXcqmkA+exHnyyfg46QDhmHyQM8QAhpNQG41qa78I3lsQykLHzU/IRuIhQDDMHmgRwgBnTlIXu2r5qDu0wRyd50iqpHHMEwPwtPUQkQ3EdEyIooR0Uxp+2lE9AERdRHRffkbpjNBzQwZjuRHE3hl6Sa8snSj7X751GwOYhim2PG6vlwB4EoAHyrbtwO4DcCDuRxUuug0gXAsORu/+9VWyz5TE3jondW45/UV2nPqTDnrtrfjklkLccmsj21bWMZymCwmuGk9wzB5JuDlICHELAAgouuU7W/Et08AMDjno/OIzicQiRoT6IamDjy9YJ11X0zgm6YO/PolI+LmwHH9MWVkH8sxuiJzX0sT/zdNHRjRtyrlGLlsRLbWoIilGB1LAYZhco8nIZAtRNTkckh9NufXZQybjuH1jR0p+6Ixgca2UOL15h1dlv0doSiOuevtlPfJk7JdbkIuy0aYggxgcxDDMPmhR7gbdRnD4fgEuqMznNg2cXh9fJ+zT+C5j9djzbZUc49s6rHLTcilY1g2aRVR3xyGYXoQ3aIJCCEanPbHNYWMtQF9dJAxge6IN50fUFuO6jLj46qmHvndQgist7H3WzUBvfzMZcYwawIMw+SbbhEC+SaoMwdFrJpAbUUgYTYKR4VllS5PsLf853M8+M5q7XWi0srcbk7OZZ6AHMXEMoBhmHzgSQgQURCAH4b5KEhEFQBCMBbRwfh5/PHtESFEt1Y8KwukrspXbW3FH2d/gaZ2UwgEE6Gk0VjMskqXJ1g7AfCvReuxoakz8dquO5nFHJRlTlpYdgxzs0qGYfKAV03gAQDnxf+eAeAxAIcC2BnAw9JxFwK4CcCNuRmeNyqD/pRtVz692PK6riKQMBuFoyIlV+Dixz6yRP+oXPXMEstru4QzebOTT+AfC9fjuYXrcfsZkzCsoVJ7jDxG9gkwDJMPvIaIng/gfJvdj+RoLBnjJXyytiL5UaMxkXAcA0ZV0dfS6CMM2E/wXn0CP33WECo/e3YJnrxomvYYWdCwT4BhmHzQI6KDAKAi6PxRasuDiV7EkVjMEiHUEUq/qqjstJWJpekTWKuJQtJdw6VlAsMwTEb0mKmlptxZqampCCR6EYejwiIEMukvYKsJpBki6iQowhbHMGsCDMPknh4jBKrKnIVAreQTMMxB2QkBW59AmiGiTodYzUFeR8YwDOOdHiMEql00gdqKYKLvQDgas/gEMmk8H7OZ4dPXBOz3WUJEOTqIYZg80CPyBACgpjw1QkimtiKQKPWgagIPzNOHhTqhagItnWF0hKKWKqJe0gScBIUsqNgaxDBMPugxQsDNHFSnJIu5lY5wQ44CikRjOOqOt7G1tQvHThiS2O5JE/B4jRy2JigIkWgMi9Y1YcKwelRoQnoZhikMPcYc5OYYHtpQKWkCVnNQJsgT9MotbdjY3IlwVOCtL7cktmdrDpJrB0VzIAVeWboJlz/5Mba2drkfnGP+8OoXOP2+9/HDJz7u9mszDGNPj9EE+tWUOe7fdXAt/rvUDBHNgSYgTcpLv2lO/F1dFkhkKXurHWR/kBwiaueDSIdLZi0EYAiwv5y7d9bnS4e/vr0KAPDG8s3del2GYZzpMZrAjw4fh6H1Fbb7ywP+RIjovBVbcdsrX2R1PbmO0KeSEJCrlnrJE/DqGLYrU5EJn21oydm5GIYpbXqMEOhfU453rz0Mk0ekFiw1yzL4pYyrjnD6CWIysiKxbGNyUpXDTbMNEZVrB+XCHJQ4Vx57LDMMU1r0GCEAGAlVVWVWp+P+o/vh8Qv2A6BvPpMpsibQ0pFc/XutHWTipC3ImkAuzEGJ82Zb2Y5hmB5Dj/EJmNRVBBN/+wh44sKp8MXNQHbdwHQsvekodISi2PeWOdr9siZgp1V4mbedjolYNAH3c3mFNQGGYUx6lCYAAHWVAenvYEIAAEgki3mhpjyAftX2zmZZE2i3qT3kxSfQ3BHGuQ9+gDc1DtNcO4YT52UhwDBMnB4nBGolTaChMmjZl44mAAA+H9m2kZRX03YF6OxkgCoc3vlqK773yIKU42SzTS5X79FcqhUMw5Q0PU4IyOag+irrSj4Tn4CuaxmQXE0LIRzMQamTbVN7CGc9MN/TteVchlw6hlkTYBjGpOcJAckcVJ+hJjCgtlx6j/4WmRN8KBpz7DJ255wvcfKf3sX2thAA4P9mf4n5q7Z7GofqGI7FBH701CKc99CHWWkG7BNgGMakxzmGnc1BzjJvcF0Fzp46EsdOGJzY1r+mTFtlNBITWLahBRubO2zPFxPAnXNWAADunPMlfn3SeEuLSjciSojox1834sUlGwAAS9Y3Ya+RfTyfy3pejg5iGMagxwmBOqmDWEOVIgRczEE1FQH86PBxlm0j+lZhjabxS3NHGMfePc/xfDs6k8LDPEc6zedlx3A0JrB4XVPitT+LinKsCDAMY9IDzUHJiV+tJ+SmCejMJMP76Pv/frlph+tYvm1JrvpDEcNvkM78K6/YhQCWSZm+IclUtHhdE65/fim27Oj+mkAMw5Q2noQAEd1ERMuIKEZEM6XtfiK6h4iaiWgDEV2Yv6F6Q3YMq5E9dpE+JnohUKU9dtXWNtexmH4AAOgIG5O2UwKZqSWEIjE0tYesjuGYsGQmhyJJITDz/vfx+Py1uPDR1AgjhmEYJ7yag1YAuBLAzcr2ywBMAzAawAgAbxDRAiHEktwNMT3khvJqS0a7SB8TnRA4YvdB+MOrqXWGnHoD69gc1wqcrEFdkRgqgn6c+7cP8NGa7RazTVQINLaHpGOTEUmdcQGzZH2yhhFj5dP1zRhQW47BDvWlGKY34kkTEELMEkLMBqB6Nc8G8AchxDYhxGIAfwcwM+UE3YhsDlIX/m7JYjohsOvgWtxx5qSsx7WppRNdkaijJtAZjkIIgQ9Xb0+x28diwuJjkDUBmXR8DrnmqQ+/xn8+2Viw69uxeF0TTrj3HUy79fVCD4Vhio5sfQJ7AFgqvf40vs0CETU5/QNQn+U4EtRKfoADx/a37HMLEf3+gTtrtx84doCnaw+UQktVhLCah3S0h6IWW79MOBqzZCZ3SUJA1n7S1VByxeJ1Tfj5Pz/FZU9+jMa2EDrDUfz25c8x94vCl45+cfGGxN+FFJIMU4xkKwSqAch1iVsA1GR5zqzw+QgvXXEg/nbePthn576WfU5C4EeHj8P3p4/S7nPzJQBGbsF71x6GedccmqhaqtIZjjlqAu2hKNq69IlnLZ3WMFVZCOzUL+m3WB33VXT3ZLdyc2vi75bOMO56fQXuf3sVzn94QcEnXvmeZ9tMiGF6GtkKgTYAtdLrOgCt6kFCiAanfwByasweP6weh+8+KGW7U4joMeMH25qL7ITAhGFJBaYzHEXA78OIvlW44EC9MOkIRR19Ah2hqG3T+6Z2qxZhZw5q6gihKxLFCfe+g7MfmJ+oOdSdE3FMAO+t3JZ4bafddBfyZy/0WBim2MhWCCwDMF56PT6+rSjxOcTWlwXsb4WdEBg/rC7xt2yvt1McOsLOQqA9FEFbSC8EmqVy1YBVCMjzWlN7GG99sQVLv2nBeyu3YdVWQyZ352K8MxxFe5e7/6K7kD96ocfCMMWG1xDRIBFVxI8PElEFEfkAPAngaiLqS0STAJwB4On8DTc7nGrmlDk4je3MSA1V+iqjPpvjO8NRCIdMgfawvSagDl02B8kVRpvaw5ZaRqYfIZe1h9zoVD5HoSdeqzmIhQDDyHjVBB4A0AFgBoDH4n8fBOBPABYAWA3gFQD/W8jwUDfkCeCXx+1u2VfuoAnYaRB2gkMNTTXpCEUds3U7QlG02vgEVCyagDTJNXeE0RVO7vu2xUggUyOfHpy3Cuu2586JLJ+9MxxDm+TELrQJRpZ/hRZIDFNseA0RPV8IQcq/uUKIqBDiciFEvRBiiBDiwXwPOBtkp+DkEQ04c58RiddO5iA7TaA86EuUmZi5b/JcTuYgp5RhJ5+ASiianGRlTaCxPYStbcnM4Yse+wjPLVyfYg76zX8+x0l/etfTtbxgyVuIRNEeKiZNIPl3FwsBhrHQ42oHOTFlZLL/8B5D6yxmAichYGfeKQ/4cfGM0Th4l/7Yc2jSSWynOXSEnfME2l32y8irfVkTaGoPY1ur1Yl89bNLcPT4wVBxC1lNh05pPJ2hqEXgZtvPOXskxzALAYax0KuEQF1FEAt/eQR8RKgqC1hWiE4+ATvKAz4E/D7svZM1FFWVGXUVAbR0RuI+AXs6QhFEo97GIZtYZMHR1BHG1tbUGkJehUumdIatmoCMrN14CbcNR2PY1hrKWXavXDS10Kap7ubD1dvxxAdrcfX/2xUj+upLoDC9mx5XQM6NfjXl6BNvGymHDnqZnFTs/AiqT2BIvZE3YPgE7Cfj1q6oxZbuRMjiGE5ub24PpWgCxjHdKATC1olWzn0wTWtvfrEZ33v4Q6zakhJRjEseX4hpt76OD1d767vghnzPS0kTmPvFZqz41r1QoRNn/PV9vLB4Ay567KMcjYrpafQqTUBFNqPYOXOdsDMhyeagyqAf9fGS1m4hone/vsLRLCXTZQkRlX0CYVQEUzUBuz7IuaJTE5FkImsCphD43sNGsbt1jQsx5ycHW45/Pd5v+cYXP8PLV87Iemzy/SmV6KCFa7fj/Pg9Wn3rsRn9PmWWe6h6y/ROep0mIFNV5s/q/eUB/ftlpaJPVRCVQeO4jnhtICe8rlTtooNaOsPYrCkpvWBNblbVJk3tIby5fHNiUpVX/2pJa7kpjzqZfbU5VRMw2daWXmnsWExg9rJvsb7RGvUkC8xS0QQ+kLSgXPpuGEalVwuBq47YBcP7VOJ/9t8po/eXB901gbrKpBDodAkRTQe7PAG7GkVXPr04NxeOc/p97+N7jyzAn99cCcDq/N28w1pnUNYM0mltme7k989F3+Cixz7Cgb9/07Jd1lJKJTpoUG3SH7Kx2Xs3OoZJl14tBAbWVWDeNYfi1yeNdz9Yg71PIPl3dXkAlWVJTSBXTd7lkEy7RLDDdxuYk2vpWBFfwf957lcArBNtc7s1u1nOgu4IR/H5xhZ4Id06P88v+ka73aIJlIg5qEYqCvhNk30L03S46d+f5eQ8pcIXm3ZgjYe+H72dXi0EgMx8ASb25iCrT6AiYQ6K5cwmvb6xAy99sgHhaEzr9PURcN9398ao/tWO58m2ppD52WRzkFriokPxERxzl3NbzmzHAlg/lyygSsUcJH+nG3IkBB5+d02KgPbK5xtbEsUJS4Fvmjpw1J1v45Db51pyVoqBpz78Gn987YuCF1Y06fVCIBvsNAFZCFQE/UmfQCjqWQhUu/grVm9tw+VPLsJf5q7UmpgG1VUg6Pe5hr7qzCPRmEiZyO0w74FFE1Dem45TWjUXma8j0ZjrQ1Mp3bMOGxNQPoVANCZw/9srcxLVFMmDEADgWLbEjnXb23HMXfNw6O1zUwR6sbJA+g422ZjTOsNRLP2muVsn465IFD//56e4542v8Opnm7rtuk6wEMgCeyGQ/LuyzI/KsuREGfYwCZUHfBg7UF+RWxUOd7++QmtnN8tZu0UbdYVj2NzSiWc/WpdYMX3/kQWY/OvXLOGJ6xvbccMLS/GFEmWS1ASSk0NTihDwvhJTJ+ktO7qwuaUT+/32dVw662PH91ZKPpqWjuQ15bHlMzrowXmr8NuXl+OMv76f9bnksNYNTbnzCaTjkzGRhVqp9LGWhZ2dtn/Ogx/g+HvewVMfruuuYVl+359+UxydAFkIZIGdOYgs5iCfJToo5GLnLg/4MPuqg20dyLKtGDBWjDqfwPA+hhBwa6nZGYni5D+9i5899wluetEoAPvWl1sgBHD7a8m2mj94fCEefX8tjrrzbcukWhGfeOUEsWw0AdVmv76xHQ++sxrb20J45bNNjvkOVWXJe9PSmRxDd2kC//xY75O44qlFOOL/3kprRR+Rfie5nHgzEQIxSyh1zoaSV+TcGbshL1zbCAC47dXl+R9QHPl7bczQNJdrWAhkgX10UPLvCtknEIoiEnOehC4/dCxG9quyTSqrKU9N7dBNjMP6eNMEOsNRbIiry898ZL8i+mxD0pkrh3zKn81EnWgXfd3kOAYZdaW+Zlu7RV3foamt9NaXW3Dyn97F8k3JMe6QhIDFJxCN4YtNO3Do7XPx2PtrPI9rU3MnrnluiaOpZ+32VJt5W1cE/16yAV9tbsX3H1ng+XqyYG8P586mnUk12SIxXaeFsPlbh1vHwVwSlp5/tUdIoWAhkAVeMoYrg35LdJCdOaiqzI8BteU4Z5oRrmr34NVUBFO26R7sYQ1GiYAyG23FRM3ulbEbg9w7oTzgw/OLvsHKLfZOw6/TqFaqCpCvt7VZHL6yY3Nraxc+39iC8x76EIvXNWH+quQEbZqDWrsilryJrkgMlz/5MVZvbcOvXvAeLfPTZxfj7x+tdzT16O6lXMJj+aYdiHg0R8krdlmTavVYYNCOSAad1eQFSa6i2/KNvHBw034yqRaQKfJYiiX/o1dnDGeLl2SxiqA/sXpv7YrYhj0+ceFUTBhWn+huZqcJ1Go0AfPQ6jJ/ouxEQhNwcQx3KsXdvDjJWjutmsCt//3c9T1eUTWBtdvbscugZPO6xvYQRsbbae57yxxbQWWag656xpofEYrEMnr43Jy96uQeiwn4fJRSxykUjdl2sJORJwtTy/rL3JW47dXluOesKTh+4lCvQ7c9r1fkt3gVYoVG/l24ad8BX/ethWUh3MTmoNLHU9mIMn+iEXxLR9g2Tj3o91kmB7vJrW+1vpENYG1yYzqGnfokAKkVPuWVnjwE2RYsm1qCfl9GE4sdOnOQfG3T6RyKxBzNFGZP5tnLvrVsD0ViGYUFO3WlA4CtSr0m83vessO6vctB85LRaQK/f2U5hAAuf3JRyvGxmEgRbrrvJZOVfCn2aLZoLy5jzrUm0B6K2C6m5N93I5uDSh+7H488X1QG/aiLm3CcslXV/sd2msCgunLbczRUJU1FXqODZE2AyD56pkoyychN7yMxYdthLRPUe7SxqcNiZmlqD2Hh2u34+OtGx/O02IS4hqMx234PTrjZjVUzlnkfdZqAF3SagBM/+fti7HXzbMxfleztrPsuM6kmKyzmoOw0ASGE554ZWV1H+lsnDOXPlEufwMK1jZh002v4xb+WavfLQpgdwz0YNVmsVmPHV1F/iPZCwL68sikE+lWXJfwQrtFB0gTrJ0I4Yi1BYSLb5W95OdlGOhyJJcxD08f2c7yWF9SVZjgaQ5ckqOat2Irv/OV9zLx/vuN5THPQzv2s5ZNDkZjrql6H22oxrEyOjW3G9VUhkIkmEIrGXM0wzy/eAAD4xb8+TY5J8550fQKrt7ZZBHO2msBVzyzGlJtne84az5SYi+CShbFdv5BMuPYfnyAcFXjqw6+1++X7XyyJiywE8oAlWazMj7pKd9eLapccGl/JqzgJgZHxevFyjoFrnoAU2unzke1KVRYC67YnQx3D0VhiZdcnBxqBOnFFYsJisnpu4XpP5zGd1+qUFcpUE3Cx46vjPugPb2LBmu0aTcBbuKzq7G/XNOZ57bNNuO5fn1pX1tLbdBN2Oqa7fy1aj0Nvn4vf/Cfp88nWJ/D84g0IRWK4/nn9SjkbVm5pxQ0vLMXj76+x+gQ090Fe/ORSE3DTcrLVpPIBO4bzgE8xB3nRBNSV5m9PmYDLn1qEPYfW4ckPkquKukr7c31nr+HYfUgdpo/tn9hW5neLDrLW+rdOZsmHp9ImgzkUjaE1ngzm5K/wiro6ikRFRlmq7fGH0Zz0yvw+hKIxQxPI4KGXBbsQIsWvIGtQJj9+ejHGD6uzbPNawE6drHX34OLHFwKwVsO12u81mkAak5AueipX0UGmkDN/f/IiI1MO/+Nbib+vP36P5LU0Y5a1y1z6BMKKGa884LP83lTBbAYQFBLWBDLgphO+JVKiAAAgAElEQVT3xIuXT7c/wBId5EN1md919RlUVpoj+lbhhcump1Q4dYr2qSoL4H/23xljBjhrAi9ePj0xYcuZtX6yCgF5NVVp85A2d4QTx+VCE1A1kWhMoDMDtdl8GM0JwBRiXRmbg+zHaLetKxJLmIXkbV5QJy6nhLtlkmlFfpfO3PCXuSsx8/73PcWo6+5SrjKuYzGBWEzg+Hvewf63vp51uKR6v+QItrBGCORLE5DHse8tc3Ds3fMs2pOqSWWSt5FrshYCRDSRiN4loh1EtIiIpuRiYMXMeQfsjInDG2z3y4ut8oAfROSqDditRvzKhFUWsP/B6iLddEJgwrD6hOP4lpeTqr4vRROAZZ8O+eHtV5OZEJCddCmaQCyWkSZgPmzmytVcLYejMcu99moekc11nSGdrT11WzgaQ5ey3asdWF1xO5XekIWamybw2rJvMX/Vdtwx+0vXMei+80zyDHRE4vWpvtrcisb2MB5+d3VW51M/q3y/ohrtR85wz0YT6AxHcfp97+GXzxu+GPl30NoVwfJNOyzlIexqYxWSrIQAEQUBPA/gUQANAH4H4Hkiyl24SAki29nNEE03v4CdA1d9EJ3MO6rAkK8vQ0TYe6c+qe/3EUKyY1jaZ2cLlhPHMtUErI5HJd5epFd7yMScrKKKEOgIRy3RW14nZVnAqj2UAX3iTyiS6tD1qgmoWeBeBaG8sHRy4qohrTp0v8hc2bSjMWERWJ+sz66Ojio05aQ6vU9ANoOmPiMtnWFPrT2fWbAOC9Y0Ytb8rxGKxByTL4FUraTkhQCAXQFUCyHuF0JEhRDPAOgEcHD2Qytd5AfdLC1RW56ZJqCqqkFHTSB1n5356IYT9kjZ5iO9JnD7q19YykbYkalPQDal6K6fSUih+bCZk7A5tg9Xb8cqKbvZqxCQJwp1Qn58/lpc+kRqcbtwNJYyATW1hzz1B0jVBOyFgNVfYb2+7Xs8rH51+RS5yhNYvmkHPpbKiSzNspiamokv3y+dH0P+DnUa9GG3z8WRd7xtCbmV+XpbOy5/8mPM+TyZh7K1tUtrFpS/H3VRUAwZ2NkKAULqgsEHYE/LQURNTv8A1Gc5jqLCqgkYK1B3TUD/VchF0QBjUrfLFdBpAnYaBhGlVCT1+6D4BATWbG3DvW9+5Th2wHCG13lwgOuQJ2LdpKyrF+SG+bCZz5idgPIaty/Pmau3teH7jyzAfz7ZCAC2kS6RmEiZiK98ejGm/+4NfOmyylRDhFUhIJvQ5GMtpjWHz+YSOQwAWj+WkybQ3B52LJKnJlBd9NhHib+3ZekTUO+zvHDQTbSyn0lnljc1Jbv6Uhc8ugAvfbIR81ZsTWyzE+5hyyIn1TFcaLIVAl8A6CCiHxJRkIjOATAGQJXL+3o0ciy4aY7J1CdQr0QDlQV8mHXBVByx+0DceeZk13M41Q5S/QV+soaIxoRzD2CZmvKAbUE9N6ydv1IfCtnJ55Wm9jAueGRBwizQt1ovOL0KAVkTuPKpRXhj+WZc9qRzaWsgNX/A5PH31zq+T9UgOpQicvJkYmdSkFfHFcp3k2tNIBoT2OeW2Tjgd2/YCgI3LSKbuv6qmcXa0lTjE/BYXtxO5q3QPBd2zm3ZRKR+VyWvCQghQgBOAXAOgE0ATgQwB8A3ynENTv8AFEdh7Rwhh2j2rym3/G+HXYSCOlGX+X0YN6gWD563Lw5QkrO05iCHPIEUIeAny4P61pdbcKG0WnOipjxgq824Ia/+dQX2MimatmxjC15fvjnxup+NJvDx2kZPtl95PmxJQyjZOVLtQm4BY1Jau81akK+1K5pyjO4awnKM8crvoxSzoE5rVNH6BOLnFELgy293JFayTe2hxPVkE4nlvS7+hGxaf6q/G7mlqU74yEJg1dY2XPDIAryxPHXc6UTv2C1WZMuAeg8yyeDONVlHBwkhPhZCTBdC9ANwNoBdACzMemRFxrETBgMATtt7uOuxI/pWYc5PDsb7Pz8sMdHO3HeE43u81rORJ27VoaWTI2kJASJPTW+A1JVlTUXANTHNDos5KD4RyCGp6fQjsMPOHHTFU4tw5B1vuwqaTFdsdqvfCod7demsjy0CDLCWwVDLe8gTiy46KOinFC3RS0SMLpTWvNbdr3+F/3fH24nosjZJSNnF/LtpAm5OVSeczEH6PAG5FEkYry/fjO8/Yix4hI15zQ2735D8udR7UPKaAAAQ0QQiKieiWhjRQYuFELlPByww/3fGZDxx4VTccoq3pvRjB9ZgSH0y63fSiAZcesgYjHbp+euGvNpW6w3pVndOeQXqPqcQURXVvDW0odK1RIUdS9Y34frnl2Ld9vaEEKpyaa+ZLm5O628anZ21mab4293PcofkKN1KulEyNahmO0vRP2lOCSWEgA9+dcHgyRyUus2cxO6YY4SY/u0dI7RTbuJjl1Pilm2sVrR1YltrF/489yusbzTKlKtahCyU9D4B+2vJE/Wczzdj5RZvJlE7IWDRBDTVZgtNLpLFvg9gM4D1AAYDOD8H5yw6KoJ+TB/b37Z8tBf+9+jd8MbVh1i2jR5QjbP2G+n4PtlUZNUErE+p7sF2qiKqM994VcnVktY/O2pXlLtkJ9txzXOf4PH5a3HRYx8lJk0nc0kmuAkBN+GXaZKU3cRgmlW2tnZh5v3v45F4nLzddeRJ1oji0lfJtJqDjHOV+X0pvxVzweAUoqjVBKIxzJqf6s+QNRXzNyqEwMK12xOd5txWvekIgUffX4vbXvkCB/7+TbR1RVLMbrI5SCd8nEJu1e9AzkR2YoeNOUjWBNR7sGR9Eybc+Cr2vnl2wfo358IcdJUQoj7+77tCCO9tpBi88dNDcOupExyPkSdE+WFWzUH66CD7r1gVENGY8BwCKLe5LAv4sOfQesfwVS8s37QjkVylW03arTC9oLblVHGL389UCNhNsh3hKDpCUdz80jLMX7UdN/7bKMpn1xRdzuz2+aw2cNkcJLTmIJ/WHPTL5z/FuOtexjkPzresVk10msCX37bil5poKLmlqPmZX/50E77zl/dxVrzYn9s91Jn9lm1owf8+9wnWbLX6SOSQ0heXbEhNFuuSHcPOGcMqmSbEySXWZRZ93ZgYgyqQ2rui2NEZwba2ULc2t5HhshEF4KHz90F9ZRB/OG2ip+PlyU/2HXjRBJxMNKoNPxJNDWm0Q25zaYaaqualP5w2EUfsPsjT+Uz++tYqAHpz0CUHj0nrXDJBl8YhbglpduYgr8k+6gN+31srsecNr+CFePVPk/U2ZqkdXaomYOMYlvME4ol/wUCqTyASi+HvC9YjJoB3v9qGOcusPghALwRM84uKrKmYY3vwHeO7XLaxBZ3hqOvkesxd8xJ9f03O+Ov7eOajdTj3bx9YtsuO/m2a+Hw7c5mJkzlIpw17MdvYaX3PLlyPW+JF+NSxyOPI1JyaLSwECsBhuw3C4l8didP3cXYWm9iZRtRJX7eScLL9qkIgHI1lJATMXAa/jywTx6C6Cpw8JbMOWLpn7pQpwzIWBG7No9wS0uw0JK++AjUnA0j9jCfe+w5+/4q+6bla48nWJyCPTdIE1AXDZxtaLOe47MmP8fKnGy3HkCY+yE7oyeMz79UeQ5LF85asa/KUbXyZknRnTqyqcJR/p53h1KQ8GZ05qN3h+9aN087UYxmrwzEPJcx9ihCIm8CCfsqo2VEuYCFQINL5wr2aQbyE/cmoK/doTHie1GTziikQiMiyEq0uD7i2t7TjU00GaXW5H9cesxsG1jqH2+pwayHY1hVFNCbw4pINWkegna9EZ0bRUa1pC6ryyfpmLF6nt6bKpgYjOkj2CThHB5VpzEGLvk69zrMfrbO81q0fdCvoznBUqwnI2tyCNds9mRrbPJYICVmEQNRx8aLTBOwmdSGEtiJsU4d7MpuXpEY1Z6EjXoeqO1tcqrAQKAG8RsroZICTWNBrAt7MG7JjuKpcP76a8szDRneV+gqbmBNpJmF1fh85mt/aQhH8a9E3+NFTi1IcgdGYsF0Be60F5EUIOCFPWn4liku+H7LZwsknYHLy5KHYf7SRb7Jd6XSlW6g0azq2tXSELY5h87qy3X3Lji5PtnYnH9a5D36QcErLi5UOFyGg++7scj1C0ZhW4HvpAmZqAk4LnxRNIJLUBAoFC4ESwKyNfthuAx2P0z3oThqHmk0c1ZQ5sEMOEa0u009w1eX+jITAsRMG4/fKhO33UcKRnUljk4CPcPo+I3DwLgO0+1u7Injts02WbY1tIdw1Z4Vj1rTXTmE6c1A6yEKAlHwOa+0lkfK3Lk/ApLYiiJMmGya7RiXjVffTkc0+Js0dYcukak7QspYUjgnb7GkZp7LO73y1NeGUljPLO8POixfdvhYbJ25nOKY1B3kpvW36bfo7VNNVz22agzJdLOUCbipTAkwZ2QcLrjsCfapcSk9ontpJw+uxy6AafPlt6kSmrljCMeF5gpXNQXaaSk15wLXRvcqgunL8+Zy9U1Zv1WX+hEDLpPKiOQnarTTbuiIWP0dHKIrvPbIAi9c14b9LN2rfAzg7GGWc8gK8IE/00ZiwvJZDCzvCUdz44meYNrpfonFKWcBne8+qyvzoE3eypggBzfE6TaCpI2zZbk66liSpiLPd3sRL1nlrVwQh6b53Rtw0geS+ls4w1m1vdwjnjOrNQWloAv1ry7HBJspLvQfmPWJzEOPKgNpy1xaHOidwwO/Df688SHu8uvowJhdvE6zsp7AzdRg+gfQmP7MInfpR5Al6jNQ+0yumELATSm1dUYtZa1tbV8I+v3yTfVkJryUtvKwkdWiLuClmO9Uk9ch7a3DJrIWJVpxVZQFbTaAi6E/kUOzoiuCHTyzEQ/EEMK/Nd5rb7cxB1oQtLwsMeZx2E/vabW3Wz5+GT+C8hz7EcXe/Y9vjuCMU1WosXr4/M0HNrpDi6q1t+OvbqyzbTEGdbXh1NrAQ6AXYTQC6PAGvjmH5nHaaQNDvS/vHbRahIyKLaaBKEgJ3zZyCaaP7pmVHDSQ0Af17WrsiFmH23kp9CWGVFs3KWIdd/L8bOlObrjqpDnM1Whn02/4GKsv8lj4QL3+6Cb9+aRmiMeFZCDQq5bHNyb5LMVnpOnypyFnwdtE2a7e1W36nbuYgc/W9Zmub1iEu0xGOakunePEJmNqZnWnngkcXpGxL+gRYE2AKgG5C7PCYtSlPKk5OT6eJZEh9Rcq2CslPIU8I8jVG9a/G0xfvjxMmeg8/dTMHtYcilmzwj9ZsT/ztVPzPS+ggAAysS/2sXtA53b0KAZMKByFQVebXZlM3toe0PgEds+avxUZJyIUS5iBruQQvmoAs+O20rDXb2hQh4E0T+PeSDbbHWM+VKlB0ZjA7yvw+S49jE7mPhYlpynPLY8knLAR6MboVS4fHED0hgMkjGlBV5scFB46yPc4pCXJk39SK43LxMdlOWqOZDNPJsDTPFbRZpbV2WSeSNVuTSVFO98SrEPjNyeMxaUQDjtozveQ5nSYQjQnPDmkAqCzz2dqcK4J+1FcGUyb8zS1dns+/ROkKljAHyVnNUeFpNf3lt6244qlFAKzO29lXHYSpo/oCMBq6WPIEPPoE3vpyi+v1O8J6c5BXjQ8wnqsLDhyFEyYlFyl2GqjpE2BzEFMQdPZ6WRM4ebL9SjsqBJ69ZH98eN0RGOSwynXSBHS2fdlEZTU5pU6G6QgBcw60C99r64pYVpdrpFLObQ41XeyiTFR2G1yLFy6bjnOn7eTpeBO78Ntr/vGJ53M4moPi+xqUvhVbWruyqJwaNwdJv6WFXzfi6meXeHr/v5dsQCgSs5iDRvWvxm6DjbDhxvaQxdS04ttWx3aZkahAa1fENgdDptPGHOQ1AABIapvy57erOcbmIKag6DQBs37LWfuNwJ0zp6TsN0Mdj95zMIJ+n8Vhq2PnftWYMKweI/pWpuwb1S+1oupuQ5L5AfLqSXcdL5UwTRKagM2KTBUCm3d4Wwmr9WImDq/H9cfvgfMP2Nl6/fhDnm79I53wSxcnIWD6c/ooJqHNLZ2eorB05w1rfAJydI2XZL/OSDShZVWV+RHw+xJhyTs6I5aVf1ckhvsVh6tMJCawYM12T0LNCBHV1xp6ZelGrNrS6lpCwnyu5M9vV5I6oQkU0BzEIaK9GJ0QMBt+28X+z//F4QhFYujn0iTHxOcjvHDZdMSEwNjr/mvZ16CEvB4/cQguO3Rs4rXV75A6ecr2Yx8B5x8wCusb2/HastRSzOa57CKsuiIxT9m/teUBS2aoGjdfXxnEBQeOSqm0aY413eqomeQX1FYELGaqijJnTQBIjWgxNAF3k1Pf6jJsUQRmWOMTkDHLWDhNyp2haMInYC4AauNhyS2d4bQa0DS2h1IK0NnREYqmlGgHgDeWb8YbyzejvjKID35xuOM5TG12gCTs7MbbydFBTCHRCYHmjjCIgFP2GqZ9T21F0LMAMPH5SDv5qqv7e8/ey7Lyle3YOqEkm5qCfh9+dcIeuOHEPVOOM85F8ffoxxiKxDxFRtUpZhPVHGSWbahVqpaa6n66fZir0sw0/uEhY3DchCHWcwT9tklYFXEhowrZzS1diHoIF67TVGfVhYjKBPyknWhl2kPRhJZl5qSY935HZ8TTd2Umwc1bsRV3zP7S9XgA+OmzS2zbRALG8/Hwu2scz2E+V9cctWtim130UrJ2EJuDmAJQbvPD2390P+w5tB4A8NMjd8nb9d1KKdhFByX2+6xCAACCNpOdaTrSFUUDjInLbXW59059UsJhVcewuXiWJ3sfJTWRepeEP5V0NIEDxvTDNUfvliLcKx00AfPzqELWq09A1ztbVzZCJhoTruaPjnAUW+J2frNiqClYVXOQHbKfQ1cmwu7evqpkjqvYFfkzMU2OA+sq8I9LD3A81hQCnCzGFAQB/UO+i1S354rDx+HIPbxHtDxz8TQcMKYfXrhsuuuxNRUBxzIBbmGofosQcDb3mNhdLhx11wR+d+oEjRCwagLmPZUzquUx1aRp40/HJ5AQhMo9cAoRNc1B6v1tbAu5+gRuO21iimYEAKGIgBDC1rzWFYlZBPyQ+gr8+IhxlmOOuWseXl1qTMZmeK0pcJo7wp5qXFVoJvlRUme/nW26/GVb118OuHDLmDcFZRmbg5hCYFf8bMwA68Nh16Rdx9TR/fDkRdMwaUSD67E15QHHfgOygNCt2mTHcMDv7PhNYBOtFI4K12Jw5QF/ik1f9QmY86ZsDpK1k3Sc2YDeF2KHOfmnaAKehID1OuGo3kFq8tfv7o0z9hlh+ZymDTwcjWFHV0RbDhwwomZkwVgZ9KeYsADgi2+NTG3TkWxey2vZEN0CY5wUkbazJjDBGF9mDYRM5PvvKgSKIDqIHcO9GDub7ZgB1tDNq4/aFYvXNbkWsEuX6vIAfnvqBPSpLsMJE1MnAbkvrqs5yCUZzMRuCg550ATKAr6Ulbk5USVICIHkCtlNO1HxUVKYpKMJmKtJ9R5Ultn7BCptzEGhSMx2sn34e/vi4HFGIT7ZJzCwthxbdnQhHI3huLvn2Y6zMxKz+IN8PrJtTm+ctyJ+rfRMaWpfZQAYN6gmETgwrE8lyvy+FDOgrsNZOshCwK0wnBk0VEhzEAuBXozdikeN3+9fU45XfqyvP5QNNeUB1FcGbdtruoaIUqom4GReUt8jE47GPGgCPsfonrKADzecaGSKWjSBNMsEVwb9idyEdHwC5uSvrj4NTcA+WQxIFTZdkdRqmtPH9sNJk4fh0F2TiwFZ4Ayqq8BnG1rQ2B7Guu36DmmAIWBkwRjwkeN9NTUBnRPaCZ1Ql02d9ZVBVJf7EWrXV/YEgFOnDMM/F32T1nXLpO/ba09yNgcxBWGGTVnlTJq2ZEI6E5yuPpHfYg6ilG067HYL4dx8HDDqGtnF+Y/sW4VF1x+ZcKhX20Q56VAFXKX03nSigxLmIK1PQHc8Jd6jNQcpdvcnLpyGM5RueLLJyPzdeAnHlH0CPiLHnhmDFJ+AExVBH4bWV+CXx+2OZk0jGFkI1FUGtd+nqQlMHdXXNkrOiXQ0AZOSjg4ior2IaD4RtRDRl0T0nVwMjMk/k0c04MmLpuJVZZXfHW3uTpg01NVMImsqro5hX7LwnBNONnm3iqBlfp/tQ11bEbCM0e+jxOSuhkOq5xg7sAY/OGh04nVlmdWm7NWNELTxixjmoNRxVzhUgg25+ARM5DBS04Hrpf6UHB3k95GlZpTKwLry+Hjdp6vhfarw3s8Px4UzRmvLP8t+gIPHDdCaoTqkEtxOZio73ITApYeMwc79rCVTSj066DEA/wLQAOAiAI8RUf8cnJfpBg4Y0x+7DEqaf4ZqirrlmrKAD/eclZqNrCKn6uvMQbIQKPcwQbjh1tow4PdZVtnypFSviZIxTULqKk811xBZBYW8Oo0JgcEei8+ZZgi1WZCdY1hefatCwGx76EZUyoTtm0b4q/x5/T6Cz0d44sKp2mMHxPNSdAK+vjJo+W3Ihd5UIdCvugyVZX7M+cnBeOmKAzGyX5VeCITMsE1n4WSHW3TQT4/cJSWqqtSTxXYC8HchREwI8RaAbQAsFcWIqMnpH4D6HIyDyRD54RrZL7WoW67x2gvZTROQ7f9eV2xOtYzcms2r1+xXnTSb6ZyWphBQ/RTqWH1EFpu9bA6KxoChDaklNwDgT2fvZXntFB2k0yYsPSEUc4xaNG+YzRguPWQMAj7CfqP6pggfJ2Qt0BRQ08f2x+gBqRE7Tvkkr111EN7/+WGJ152SSa9PtfU7MaOXxg6swfhhxpQjC3LzfpiO4qDf50n7UJHvf8BHloC08oAPAb8v5XeYaS/uXJCLK98L4GwiChCRmU/9WQ7OyxSAGeP0foJc4tXaJGsCulBJ+UHyKgScru0l9lyuQipPMnWVqROVacNWNQHVwUmwCopKaeKJCYF+mnaFZQFfivZhjk1nDtJFvMi5DCmagGTSGVhbjke/v1/K+wEjkuyDXxyOJy6cmnLdGeP6475z99K+Tw6blRcFqh/C7yPLeZ+9ZH/rZygPoLYiiB8cNBo+Au46a3Ji33XH7YHJUqjy/x69W8o4Ki3akPU3FMzQHCTXySIiywSfMBEqUrnUo4P+C8MkdBOAKIAzhRDt8gFCCMegcdYG8s+AePieHbd9ZyKWrG/CxZJtOl94VXxlTUDnwJNNHBUeHXDyxPnkRVPx/sptuOeNr7TH7ja4NqWrmDx59ZU0AZ05yHzg1cnxjjMn48R73028NjQBG3NQTGjLdPSpCqb6GjTRQUTGa60QkCZ+NURUdgc8dP6+GOvQzc0cnyrspozsY1tixOIYlt6m9hyoCPgsmuq+O/fF0PqKRMMcc9V97TG74bLDxlo0smENlXj+suloag9hzbZ2i0BInj95r9UIpaCPLGbGgI+w57B6LHGpRqqGWJcHfInIM7NQn+qbKllzEBH1BfBvAFcBKAewL4B7iUgf88cUjGcunoYz9hmO/145Q7v/jH1H4JZTJnRLlIJXx7Occap7j2Xi9BhpdOpew7D/6H44Z+pIHDCmP/5n/51Tjjlx0lDMGNcft58+KWWffH/kJDpHc5ByTycOb8B71yZNGERWQSGHa8aEPlmvT1VZymrSnBBlc0RdRRBEpDV1yULArmS1MX5v35f62xlQW277ewpaQkSTf6vdx3TfqyUb2ywHQmSbR9BQVaYVAIA1q7gqmFrvSdYEhvWpxK9talOp75ORzWR9413c1O+ulKuIjgHQIoT4V/z1J0Q0H8AhAD7N8txMDhk9oAa3nZY6qRUCr+Ygt+AUqybgTQiUB/x46uJpidc6W+xZ+43E/mP6oVkTXSJP6HJHLl35BNMcpMtdkFfeRNbEJnniiQqBU6YMwz1vfIUqyazTUBVMES6mIJEnIdNMpYvYsQgBByHqlnuhXt9kQE2qoBozoBo3nrhnoo8xYF0Vp2gCGg1QDhXNNpLNSRMI+H2W/X4fWQSQV2TNzPzNqI76dHNJckm24udLALVEdAIZjAdwEFgAMA547V3rht/iE0j+lP9yzl44ZYq3+G6dGm6upHX75IdVFgI6c5BZKlu3mpXP7VP6KcufRQiB0QNqMOcnB+HNqw9JbB9SX5m6mtTkCdSUG2PQhb9WW4SA/eRml2imonZt619TbtFK+lQF8fpPD8GMcQOsjmHpY6hhqXohkLscVzkcVxWEZX6rPyLos++fcfPJ41FfGcRjGt+JfN4+dkLAozkzH2R1ZSFEM4AzAdwMoAXAfwDcJoSYm/3QmJ5KrtIQ/DbRQcdMGII7zpyse0sKOnOFuXLTOevkCUhe/escw6fvPRxH7D4QFx6Y6meRr6uGiMoTp1m6YezAWgyqq8BZ+43EkPoK/PyY3VLMNLroINPJ3N7l7BiurwzihhP2wEGaBEKvmoDqt+lfU66t9Gr8LYeIJrerjmGdL2jmviNStmWKvNJXhUBt3JRm4vcR+laXabXH707bCYt/daT2/g3vk3QU940HE6gRciXdY1gI8aoQYrIQolYIsZMQ4vZcDIzxhhlv/9002xYWkqH1+nDDdLETAumgm+DKbKJsAOCkycMwYVg9jhk/GINqnUNERw+owYPn7YsDx6WmzVgb4lg1gYbKpIYxdVQ/y/tuPXUC3rv2MAysq0j1CWiqiJpayK6Da6GiVjT93vRRuPGE1AbpXqtqjlaqcvarKVMmfr0fQJ5Tz5++s+UcuhDNo/YcjEe/vx/mXXOop3E5ccyEwYm/1dBatelRIJ5h/fY1h+L1nx6cci4709ROUoJanyo7TaBw5iCuHVTinDBpKKaN7of+mjDCYuO20yZi1vy1uNPjKt2NXAgBM4RPLiJmagLyQ23+WRH0499XHAjAWnfea40Y+boy1jwBH1798UFobA9hj6F1tu9VNRVzIim3aALGuK4/fg/4yHDWPjDPsMfr4u91Ga5eNYEBSrmRmvIA2iQNxBI/rySLmVx5+G8FcJMAABAbSURBVDh809iBF5dsAKD/XokIB9uUPEmXvXfqi2cunoaBdRW4c4618Yw5YZs0xF8PTjOhcmTfZO6NGe6b6hMoYU2AKTwDasu7pdRDtpyxzwi8ePmBGD3APtxQ5s4zJ8NHxgSmw29jR08XeUKqKQ9YHloTnR9jiDQZqIlJ6eAjSjGP7Dq4FtNG93N4l/1EIk8o5iQ6oLYcd86cgiP3SK58db4KnanDaxVU9TdIZO0gptrXk5/DOt7T9xmeMv58MnV0P4zqX50iVM3v9AcHj8bA2nJPkUE6dpISME3fkfrdeayOnRdYE2CKlpOnDMORewyyzRiVJ+Z0G7jLyHboJy6cqhWoOhE7cXgDfnbUrqgM+jG8T+aZ1kZ0UPIKXlfeaqkMvU9AjXhxFpw6TSCdJit9qoJolKKq7MxBsvlDtbrJWlV3ZtKq5j9z5f/zY3bHtUfvlvFCa5BU9sPuHE45PPmGhQBT1DiVDMikbIQO2RTkReDIXHbo2IyvazJuYI3ls3iNy1cjeso0QkC9L34XwZmNOQgwEssuf3IRLj1kDABl9S9N6PLY1cQpeQzpNuHJBvW+y+agbDTt3YfUYa+RDWjpjGD/uHantt5UTWndCQsBpmTJlTlIxi4E0HOacxrcceYkvPXFFlx5xC6Yv3JbYnum0ThJh3aqY9hEDsEs1wkBzco7HU1gysg+eFdKhJPHUmYRAslrq59X9ml0owxINQel2Q/aDr+P8I9LD4AQSaEmJ0IetecgbWe17oJ9AkzJIq8SdRPaiZOGAgD21DhX7bBr55iPyeiUKcNx58wpqCkPwG8TMumE30faEMxyB02gS0oa0yXY6Sb8bJyWlvFJJiA5WU69plUIdJ8UkLUWovQ7mTlBRJbfq9y74jcnT8i6r3E2sCbAlCzyY6Mzbdx66gQcvvtAHJRGUTy1ho5Jvicj62Tu/VoBPyVW9+b75BX3UXtaezjLfgSdM1s1ezz8vX2zmqDk81nMQZKwTamoqdQ+6i5kB3h9ZTCvpig5g9tW++wmWAgwJYvc/lBXt726PICTJqfXGcruwc+/EMhs9Rv0+dCJZOljwPgM8645FNvbQolOZyZ7jeyD0/YejoqgD7sNdteQZozNXWsQWQhYO6+pmoC9gMgncnFANTw018hCIFemzExhIcCULKFI0r6d7zjrfK/WAoopIpP3ySvoEX2rMEIT6kpE2sJ49ufP3X2VtTXZJ+DoGO5GTcBacyl3piAdcoXcQod3s0+AKVnCUlSP116u6fKzo3ZFQ1UQ931377yc3yRTk0vAJgSzmDhu4hDUVgTws6N2TWyTo7DU1X6hfALyvazNs9D30oKzu2BNgClZLEIgTxPgZYeOxaUHj8l7qGKmtWMy9SV0J/eeNQXhqLAIalkTUCuHyt9ld66Sg0rSYD7p0PR3KBQsBJiSRbbbOuUTZEt3xKrLmgClEY9qMQcVqSZARChTauPI35caMy/f7+4NEZWil3qRJlCcvxqG8cDBuwzAWfuNxI+PGFfQZJtcIE/mAt5rCMgaRLGag3TImoDcRtRk/9H9UBH0JZLOugOLOSiNctW7aYrzuXHzSUYJisN3G5j2e3MNawJMyeLzEW49tXc3sfNb4vBLRwjI0UGdmlXxrAunoj0UsTSQyTeyOcguX0TmvnP3xqz5a/Gbk8enfa1zp+2EfUf1xej+3upo5RMWAkyv59gJg/Hyp5tw6l7phZPmi/TMQfqM3GyoqwigpTOC0/ce7n5whsh5Amb/XRm/j7pVAADWMF2zGY8TR48fjKPHD3Y9TgcReQrR7Q5YCDC9nttPn4Qz9hnhWrUznwjJApSWOcimSmc2vHTFDLz5xWaclkchIAssnSZQCCzVZHPYvazY6T2flGFsqCoL4JBdC2+bzQTZmZmrSJqR/apw3gE75+RcdshjVR3DhcLaljP/JayLhdIxIjJMD2ag5NhOx9GYy2SuQqEzBxWCgD89c1BPgTUBhikC+lSX4ckLp6I9FMXYgd6FQLHmBqRDXZGYXmR56sUx3FPIahlBRK3KP0FE38nV4BimN3HA2P44Yo9B7gdK9K0u3dDYXx63O8YMqC6aCK+QpJHUsibgDSFEIr6JiCYCmA/gtWwHxTCMN35x7G74dH0TjsowSqWQXDhjNC6cMbrQw0ggm6WqepEmkEs97BwALwghduTwnAzDODCkvhJzf3ZooYfRI+i0VPZkIZAWZLj6zwLwQ5v9TS6nqHfZzzAMk1d2leL2c9VVrBTIlSZwEIAqAK/m6HwMwzDdyuQRDbjjzEnoX1Oe0r+5J5OrT3oOgL8LIcK6nUKIBqc3xzUF1gYYhikop0zJX4JcsZK1ECCiMgCnATgh++EwDMMw3UkuMk2OBdAM4L0cnIthGIbpRnIhBM4B8KQQwnvBE4ZhGKYoyNocJIQ4PRcDYRiGYbqf0i88wjAMw2QMCwGGYZheDBWDKZ+IYgCovp6jRBmGYbzS3NwMAEIIkfGCvliEQASGVtKS4SlM6dGcmxHlDB5XevC40oPHlT7FOrZMx1UHICaEyNi/WxRCIFvMshRuSWndDY8rPXhc6cHjSp9iHVshx8U+AYZhmF4MCwGGYZheDAsBhmGYXgwLAYZhmF4MCwGGYZheDAsBhmGYXgwLAYZhmF5Mj8gTYBiGYTKDNQGGYZheDAsBhmGYXgwLAYZhmF4MCwGGYZheTMkLASIaSESvEFE7ES0lomkFGsdcIuokotb4v/9K+35JRFvj/67P8zhuIqJlRBQjopnSdj8R3UNEzUS0gYguVN53QXx7MxHdS0Q5/W04jOtGIgpL961Ved/xRLQqvu9pIqrK8bjKiehhIlof/+xvEtGe0n7b7y6fY3MaFxGdT0QR+Z4R0UjpvVPjz0I7Eb1KRANzNa74+f9CRBuJqIWIPiWi46V9BblfTuMq9P2SrrN//Pd/bfx1QZ/JBEKIkv4H4J8A/gygEsB5AL4BUF6AccwFMFOz/UQAqwGMBLATgLUATsjjOM4FcCSA+fJ4APwIwAIA/QBMBrAdwKT4vskAtgGYFN+/EMAV3TSuGwHcZ/OewQCaABwBoBbASwD+mONxVQO4HsBwAH4APwGwwu27y/fYXMZ1PoBXbN5XEX8Gzo0/E/cB+EeO79luACrif+8Tvw99Cnm/XMZV0PsVv44v/tv/AMC18W0FfSYTY8vHSbvrX/zHFAIwSNq2HMBRBRjLXOiFwNMArpFeXwvgqe4eT/wHeIb0+j4At8b//h2AP0v7ZgJ4v5vGdSPshcAlAF6WXk8DsDHP960MQCz+4Nl+d909NmVcTpPa0QCWSa+HAugCUJOnce0DoAPAHkV2v+RxFfx+xT//XQAekYRAUTyTpW4OGgegSQjxrbTtUxhffCG4h4i2ENFsIpoY37YHgKXSMYUan9M4Cj3GmUS0nYgWEdGp0nbduAYTUZ88jmV/AJuFENtsru90z/I5NnlcADCdiLbFTWyXSMdZxiWE2ACjWdO4XA6GiP5MRB0wVrJzAHyuXhsFuF824wIKeL+IqB+AHwO4QdlVFM9kqQuBaqR2I2sBUFOAsVwDYBQMVXg2gJeJqAapYyzU+JzGUcgx/h2GGj8QxsrxYSLaRzcuIUQbgGi+xkZEDQDuB/AL3fXhcM/yOTbNuN4CMAHAAADfA/ArIjrFZszquHOCEOKH8XMeCWCOMJarBb9fNuMq9P26BcCdQogmZXtRPJOlLgTaYJiEZOoAtGqOzStCiA+FEK1CiA4hxG3xMeynGWNBxucyjoKNUQixTAixSQgREUK8CuApGLbllHERUTUM+3jOx0ZEFQBeAPCSEOIh3fXhcM/yNTbduIQQq4UQa4QQMSHEBwDuBmBOat32TAghokKIOQCOIKKjNNfu9vulG1ch7xcRTQGwL4AHNLuL4pksdSGwAkAfIhokbRsPYFmBxiMTi/+/DMaYTAo1PqdxFMsYgeR9A/Tj2iSEaMzlBYnID8OevQ7A1S7Xd7pnOR2bw7hUbO8ZEQ2BMYGsyNW4NAQAjFGvjW6+Xw7jUunO+3UwgF0BfENEmwCcCeA6IrpPvTYK9UzmyzHTXf9gRAfdC8PL/10UIDoIQAMM9bMchgPvKgAbYUjyEwGsgmEmGgkjeiKf0UHB+L14O34/KmAI+x/BiEzoCyPiYBuskQhbAUyE4XhcgNxHB9mN60QYTbZ9AA6DEdExNf6ewQAa49trALyIHEeUxK/zEIBXAASV7bbfXXeMzWFcRwMYEP97LwDrAZwWf10BYAOAs2FEu/wFOYx2iX/Wc+P/BwB8B0AnDHNLwe6Xy7gKeb+q4p/d/PcMgN/Ef/MFfSYTY8zHSbvzH4BBAF4F0A7gMwDTCjCGAQA+gqGubQfwOoDJ0v7r41/wNgDX53ksjwAQyr9DYKje9wJohiGgLlTed2F8ezOAPwHwd9O4no5PDjsALAFwuvK+E+KTSVv8AarK8bh2io+lI/79mf9muH13+Ryb07gA/BHA5vjrLwD8UHnvtPiz0A7gNUjRczkYVzWAN2AI62YYoYunePmt5/l+2Y6rkPfL5jkwo4MK+kya/7iKKMMwTC+m1H0CDMMwTBawEGAYhunFsBBgGIbpxbAQYBiG6cWwEGAYhunFsBBgGIbpxbAQYJgsIaKdiUgQUf9Cj4Vh0oWFANOjIKO5T5fSQKQ1XoKBYRgFFgJMT+Q6IUSN8i9a6EExTDHCQoDpNRDRmnj7wwVEtCOuNYyR9g8komeIaHO8reOd8Sqe5v5xRPRivGfEdiL6p3KJ44hoRbwd4N/jlTLNVpF/i7+vJV7T/pDu+dQM4wwLAaa3cRGMTlMDYTQc+ScRUXzfU/H/x8AoAz4dwG+BROnjOTAafewMoxjYPcq5T4DR0WocjIJgl8a3nwdgCoBdYBQOOwHA1zn9VAyTISwEmJ7IzUTUJP17T9r3VyHEZ0KIDgA/g9GtaRIRDYNR4fInQogdwugwdQMMgQEAx8MoQXydEKJNCBESQrypXPcmIUSzEGIzjCqZe8e3h2BUt9wdgE8IsVIIsSoPn5th0oaFANMTuV4I0SD9O0Dat9b8QwjRCqPa5TAYDd0jQohvpGNXwehXUQWjqudK4VxxcaP0t9wUZBaMstB/ArCViGYR0eBMPxzD5BIWAkxvYyfzj3j7z34welCsBxCIawQmOwNoFEK0wxAeoyXTkWeE0TXtd0KIKTBMRXUAbsv8IzBM7mAhwPQ2LiaiPeIO39/B8At8EtcA3gRwOxHVxjtM3QTg0fj7/gOjMc6viaiaiMqI6FAvFySiw4hoMhEFYNSs74TRX5dhCg4LAaYncosmT2BkfN/fYEzsW2B0nfqOEMJsN3g2jK5Uq2A0CfoAwHVAwnR0BIzOVF8D2ATgMo/jGQTD6dwEo1VkAMC12X1EhskN3FSG6TUQ0RoAVwshniv0WBimWGBNgGEYphfDQoBhGKYXw+YghmGYXgxrAgzDML0YFgIMwzC9GBYCDMMwvRgWAgzDML0YFgIMwzC9GBYCDMMwvZj/D+MeOPrWVOFfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(hist.history['loss'])\n",
    "plt.xlabel('Epochs')\n",
    "np.mean(hist.history['loss'][0:]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 227,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<tf.Variable 'dense_31/kernel:0' shape=(3, 2) dtype=float32, numpy=\n",
       " array([[ 0.00079479,  0.        ],\n",
       "        [-0.37138182,  0.        ],\n",
       "        [ 0.5108788 ,  0.        ]], dtype=float32)>,\n",
       " <tf.Variable 'dense_31/bias:0' shape=(2,) dtype=float32, numpy=array([-1.7670381,  0.       ], dtype=float32)>]"
      ]
     },
     "execution_count": 227,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_nb.weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "W = model_nb.weights[0].numpy()\n",
    "b = model_nb.weights[1].numpy()\n",
    "W,b\n",
    "mean_disp = np.matmul(X,W) + b\n",
    "np.mean(np.exp(mean_disp[:,1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}