TensorFlow Distribution as Cost Functions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TensorFlow Distributions as Cost Functions\n",
    "\n",
    "Sometimes you want to fit a function were your standard mean square error or mean absolute error doesn't not do a great job. Sometimes you want to find parameters to probability distribution. Here I show how you could do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/stephen/miniconda3/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n",
      "/home/stephen/miniconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n",
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/stephen/miniconda3/lib/python3.6/site-packages/tensorflow/contrib/learn/python/learn/datasets/base.py:198: retry (from tensorflow.contrib.learn.python.learn.datasets.base) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use the retry module or similar alternatives.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as stats\n",
    "\n",
    "import statsmodels.api as sm\n",
    "\n",
    "import keras\n",
    "import keras.backend as K\n",
    "import keras.optimizers as kopt\n",
    "from keras.models import Model\n",
    "from keras.layers import Input, Dense, concatenate\n",
    "\n",
    "import tensorflow.contrib.distributions as tfdist\n",
    "\n",
    "from sklearn.metrics import r2_score, mean_squared_error\n",
    "\n",
    "from rosey.keras_utils import keras_r2\n",
    "\n",
    "import matplotlib.pyplot as graph\n",
    "import seaborn as sns\n",
    "\n",
    "graph.style.use('fivethirtyeight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rsq = 0.558\n",
      "MSE = 7.780\n",
      "Mean loglikelihood = -2.4497965129009436\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAD1CAYAAABgOJMiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJztnX18VNWd/z9nZjKTZMAAiQQhEU2MrlarXVjEBPGpdLVq0V2CWKtARenaraxuXz+tdOmvdtH6q6srVl0QBR9aMdhVFMQtPv0QEC32RystxRB8SKAEEyDAkMxkZs7vj5k7uXPvOfdh7sPMnTnv18uXZJ7uufNwvvf7PZ/z+RJKKQQCgUAgcAtfvgcgEAgEgtJCBB6BQCAQuIoIPAKBQCBwFRF4BAKBQOAqIvAIBAKBwFUC+Tx4X1+fkNQJBAJBkVNVVUXkf4uMRyAQCASuIgKPQCAQCFzFc4Gnvb0930NwnFI4R0CcZzFRCucIiPO0C88FHoFAIBB4GxF4BAKBQOAqIvAIBAKBwFVE4BEIBAKBq+R1H49AIChcVndEcO9HR9EVSaAu7MeiCcPR2hjOekxgywYEX1qO83q7QatrEZsxD/HmaXkascAriMAjEAhUrO6I4PbNfehPpPZ4d0YSuH1zHwBkgk9gywaEVjwIEosCAEhvN0IrHgQAEXwEmohSm0AgUHHvR0czQUeiP0Fx70dHM38HX1qeCToSJBZF8KXlroxR4F1E4BEIBCq6Ignd20nvAeZjeLcLBBIi8AgEAhV1Yb/u7bR6NPMxvNsFAgkReAQCgYpFE4ajwp/l64gKP8GiCcMzf8dmzAMNhrIeQ4MhxGbMc2WMAu8ixAUCgUCFJCDQUrVJAoLgS8tBhKpNYAIReAQCAZPWxrBKPq0k3jwN8eZpaG9vR1NTk0sjE3gdW0pthJCnCSEHCCE7ZLeNIoRsIIS0p/8/0o5jCQQCgcDb2LXGsxLA5Yrb7gbwFqW0CcBb6b8FAoFAUOLYEngopRsBHFTcPB3AM+l/PwPgGjuOJRAIBAJv46SqrZZS+lcASP9faCwFAoFAAEIp1X+UkRci5BQAaymlZ6f/PkwpHSG7/xClNGudp6+vL3PwUmmwJBBYZf0BPx7/vAzdUYLaEMVt4wdxxWj2hs9cGbljK8a+8zLK+g5isGoU9l1yLQ6dPdnWYwiKG7nYpKqqKkub76SqrZsQchKl9K+EkJMAaG5nNqqIKQX1TCmcIyDOMxdWd0Rwf8eQh9r+KMH9HeU4aUyVrgLNKIEtGxB6/fmMHU6w7yDGv/48xtSO4UqlxWdZXDh9nk6W2l4FMDv979kA1jh4LIGgJNDzUFvdEcE5bfsxcsVenNO2H6s7IqaPUegebIEtG1B553UIz74ElXdeh8CWDfkeksAktmQ8hJAXAFwMoIYQ0gXgJwB+DqCNEHIzgC8AtNpxLIGglNHyUDPiKG2EQvZgE47YxYEtgYdSej3nrsvseH2BwAsY6V9jlbqwH52M4FMX9mtmQ2bGQatHg/R2M2/PN1rZWCkEHqn/Eek9AFo92rNOEcKrTSCwASnb6IwkQDGUbeRS6tJCy0PNiKO0EQrZg62QszGnkbI9X283CCh86WzPi6VGEXgEAhsw0r/GDlobw1jSUoX6sB8EQH3YjyUtKWGBEUdpI8SbpyE694dIVteCgiBZXYvo3B8WxJV1KTtiF/ramxmEV5tAYAN2ZRtG4HmoLZowPGuNB1A7ShtF8mCzghOlx9iMeVlrPEDhZGNW0SujFVO2JwKPQGADWmsvbmHEUdot7BI6KMl2xPb2OoccI6KJQl57M4sIPAKBDWhmG0l7y21aGHGUdgO7hA4s7MjGCg0jooliyvZE4BEIbEAr27BiyuGGUs4J3Co9FovKy0gZrZiyPRF4BAKbsDvbcKpc5QZulB6d2NOTr0BmtIxmNNsr9AsWoWoTCAoUt5RyTmCkdbZV7FZ55VOubKeE3S1pvxVExiMQwNwVoltXk06Uq9y6ojcrdMhlXHarvPK5OdXOMpqT62t2IQKPoOQxU9Jys/xld7nKbbsZo6XHXMdlt8or33Jlu0QTbkr7c0WU2gQlj5mSlpvlL7vLVYW6ATHXcdntsFAsm1Pt2kjsJCLwCEoeM1eIbm8U5bkU5EK+r+h55Douux0WCtkqyAxurK9ZRZTaBCWPmZKW2xtF7VTKFeoGRCvjsnNPT7HIlQtpIzEPEXgEJY8Zqxk7bGnyJXUt1A2IhTSuYtmcmusFiyTyOK+3G7S61rHAKwKPoOQxc4Vo9Woyn3tzCvWKvlDHVWq4KT4RgUcggLkrRCvlr3xLXQt1AyJrXIEtGxB8/lGQyBEAAB12AmI3/EAEJIdwU04uxAWCvGNHu2avUEhSV14LaTs2IFr9TANbNiC0/OfwRY6AACAAfMeOILR0MSq//y1P9qApdNwUn4jAI8grXthlbSeFInXV2qVvVTJux2cafGk5SEIdjDMByKMN0AqZvRXVpm63ggg8JUqhZBletoXJhUKRumqVVaxmZXZ8pnpX2YWw/6jY+NEpMxHxBbNui/iC+NEpM20/lgg8JUghZRmFVHpymsCWDZj92Hdx5K0b8NkHC3B992bLe3NyRauswsu+RoaIoYsVOz5TOkw/EOd7/1GxsalhKuafMQ+fh2qQBPB5qAbzz5iHTQ1TbT+WEBeUIPle4JZTCA3U3ECpGKrr78GzHU8hOnUk4o25LdxaEQBo7Z1hScaDPuBojOJgNPVZ8dR4gS0b8OkHSzG2vwedoRosbJiJVbUtAICRwVTgMjReyr5ZOVY3kPvIfaVqJJKz/qkoBQ6LJgzH7QNTMp8XkMrGlziQjYuMpwQppCyjUEpPTmO3XY3VrFVrlz7LMWFYgGBQEQyU5TMpuNb198AHYHy0B0t3Lces7s0oI8CxODU8XhLRLsu5tc9HuRYW7DtYtOtL2Z87dTQbF4GnBCmUBW7AfluYQsVuxZDVdRQ9u5nWxjA+njkGh+aOS/0/xk5B5BcrrOAaTsbw88/aMDxIEEtmP1drvFzfNMCyNY4ZCtXfzimkz/3DKf34eOYYx36HotRWgtix+95OCqVds5PYbVdjR9ZqZpe+kZIoL4iO6+/Foah+4JLDczNwK+BIFKq/ndcRGU8J4vUso1AUeWaw24DS7azVSElUy93Z7HjtMADl7VMyQ7E4VhcaIuOxmUJvOSvh1SzDDcsZJ5ql2W0L43bWasQqSMtzbVGt+fGayciUv7sVwW24aN2jlu1f7PKRy1dL7UJFBB4byacPV6nglCJvaGJIlcOka3s7/arsNKDMhwOx3sWKMrhGqmqw8NTr8Pius1DXdRTfPq0cv+2K2T5e1u+u4a2Vtti/KM9pMAdVm9sN+LyACDw2Ukgy5WLFCUXetlfWYtKrS+BLxJj3u9X+2CxWs1ar2bn8+bWhcvzMF0FrOriygsGvdw8wS7p62YDe/azfXd1AD3PMrLUZvdeXXzC0t7ejqanJ8HsE5LeldqEi1nhspJBkysWK3WsbqzsiaFi/EpWcoCNRbIvJVuXYyufvj/qynm9Udadl3WPkfoD9++oM1TDHrVybMfL6VhECBTUi8NhIIcmUixW79/3c+9FR7tWxnGJbTLYqx9Z7vtGLMD25shE5M+v3tbBhJo77s+1fWGszbsilhUBBjQg8NlIqmyHzid2KvK5Ignt1LFEIzdLsxmp2rvd8oxdhetmAkWyB9btbM3YKPvzW7bqqODeykWJpqW0nIvDYiNdlyl5B2uS2bOoIAMCtGw/nLKuuC/uxsEFtjpiE+5sV3cRqdq73fKMXYXrZgJFsgfe7m3jNVTj+0IuIzr8HABBaep9KVu1GNmKHNLzYEOICm/GqTNlr2KUglPypAGDxnjbUR3vQVV6DPVfMwcRrrtJ87rZX1uK0N1ZibH8P9lXUYPfl+s8xi1MyXKtybL3nG1Xd6cmVjcqZeb87PUWZ3W23eYKNYmmpbRci8Ag8iV0KwqEJcipOq20xrO7a9spaTFqzBJXJlCihrr8Ho9YswYdATsGHNWFd373FVhmuPIjNrh6NU6feiLmxiTmp2pSBpTaUxM/OH5X1fK2LsKHzPQu3nXULFn/6IsJ9ParganX/k56iLN48Db72HSh79zUgmQR8PgxOuTyn91dspzAOodSADaxD9PX1mT54LnJGr1EK5whYO8+RK/YyDYwJgENzx1kalxEOf28G6vrVooSuihqM+K+Xsm7TO0/lhAWksoe/bluAYX1fqh6frK7F8Yde5L4eK0sCoGlBY1Vabeaz5J2vE2Xp8OxLQBjfFAqCyDPvqDIiQNuaR+s8z2nbz7QVqg/78fHMMRbOwn3snoOqqqqy6q5ijUfgSXhrDBRwxUZnLCPoaN2uBS97q+wzvhdFgicPDj7/KPfK3+3+TG42/9Nbw7FT1Sa2UxhHBB6BJ2EtXku40dhuXwVbCce7nYXkOce6SgaAzhC75bDWwjdvIiWRI8zHk95u17vAujFBSz5tpLdble/I13DsVLWJ7RTGEYFH4EnkSiYWTrfP3n35HBxXKOGO+4LYffkcQ8+XZxk8HjpzlmkZrukJ0+dz/Urd6Qk6O+tLlV8p2CpFO1VtYjuFcRwPPISQzwghHxNCthNCtjl9PEHpIMmq2XlPKvNxysF64jVX4cPpt6OrItUmuKuiBh9Ovz1LWCBlNJM2VajGwMoy5FT4Cc771hWmZbimJ8xk0lAgsMPpWcLpCZqZ9QGg6bUx+ftn5x4bsZ3COG6p2i6hlJovfgsEBuD1igGQtWYB2KsumnjNVcA1V+E4gBEAJsruy15AJ2jZsxGXvdWG8EAvaPVoTBn9j3hB1mJYYlb3Zvz8szaM6+8F3Z4SBmgJCSTkJqepIw5BAdBgOUhsQPU8Wl2rK41mSZKTy3+BbQeiOSn47DA41RJDmCmf2e0aXkjbKQrZKV/IqQWehzVxKnHLrFUKALN7D+DSUDUWNswEACzdtRzhtPSa9HZj6eHloEBWf/tZ3Zux7JPlGd84o/JpljJLDgEQ9QcRClLmfhW9QMDKICoTMTSsX4nV51yS03tqZYLWky2bbbpXjHtsCl3a7bicmhDyKYBDSF14LaWULpPuk8up29vbHR2HoLhZf8CPxz8vQ3dUEs8SzOrenNkU2hmqwY8bWnHnP/6t5ms82FGGI4lUvlAVAP61IYYrRhtb6xi5Yyvq1z0H/+CQ4WjEF8RxXxAnxo+pHh+HD3PO/F4m+Ox5/3acHO1VPS5WNQp/+sEDWedYG6L4ZWwjvvn736Cs7yC33CiRBPDqJbdkHj9YNQr7LrkWh86erHte5/37LczXTwI4+e+fw2t/p86klGO9bfyg4fdRzsgdWzH2nZezxjyl/2Lsj6pXCcaEknjt7waYn0OiLIjOK280dL7FwNW/K9d8j9xALsdWyqndyHhaKKX7CCGjAWwghPyFUrpR+SCjmvFS2ONSCucI2HueTU3A7enk4Zy2/WjZszEryxgf7cHST56C78tRzKvb1R0R/Kz9MAZl12F9ceDfd4dw0pgRhq4SK5/4MXyD2S7X4WQss8lUSQBJLPtkOQiATQ1TUR89yHxcWd8hbPeNxf0dQ1ewF3+xBd/Y9QyCnNdW0hmqwR2Vl+HrS26A9Iya9H960OpaZgbRGapBd9SHpqamrM9ydUcka6z7owT3d5TjpDHm1jsCWzYg9PrzmWwr2HcQ419/Hpc0hphlSmksaGrCYO0YEFn5bHDGPNQ0TzN0vlp45bfZvWkv+3bpPdLB6fN0XFxAKd2X/v8BAC8DmOT0MQWlzaIJw3Hfp22ZoCNRmYhx92fc+9HRrKAjEUvCsDouFwluZSKGZw78Bh/PHKOpsFKKERbvaeMGNCURXxALG2bmrFKLzZincnqWXpMlTLAqz5aEDKGli5nS8Ps/a2M+Tz6WePM0HH/oRUSeeUclKCgFCl3a7WjgIYSECSHDpX8D+AaAHU4eU1C4SCovSWm2/oAzP4LWxjDqB9QlK4AfHLQmZaMTNi9w9PjDqombNSYthZVyDPVRba1OMv3f56EazD9jHlal7YById48DR9+63Z8UV6T9Zprxk5hKtGsyLODzzyM0NLFGSk0i3H9vUK2rEOhS7udLrXVAniZECId69eU0jccPqbABG4pX1iLnfftDuKkMRFHjmd2gVlLGecjKYsevfeHaTgJoCYRyfybNZlKY9JSWNV1ZW807QzVYDwj+FAAkaoTcfu4Vjw7eqgcZXXSmXjNVVh9ziVZ35UlnPeC917qBb7Alg0oe3uN7noVrR6NJS1VBavYKgTy0RrdDI4GHkrpHgDnOnkMQe7YoXwxGrhY5ZeBJNFUmllxZjbrOrxownDc9t5hZrlNGrb0/pz68TuYuvE51biyA0cq6CknUZXUWTGmF2qbce/kc4bez9rhaIVaubewYWbWGhaQKn9974x5eOLumbiwI4J3bJ505Eq0wJYNCD429NmMvPCq1EIbY6yAscAXfGm5ftCRKfHyMYmuP+DHP/y//QU5mSspJGm3EiGnLmGsOjybCVxmyy96dvZ6mN2fIY33rq19OBhLnY8PqXKVnOn7NmHSu8vh40iepf8q77wOPkbGRQDEiQ8+mhxqpZAe07ZX1uKy9SsxeyClwlvYMDPTskF5Bbt6TCqbkav2FjbMxOaGqZnHGwnokaoaLDz1Ojw+4gLDEynrs6lf9xwGa8cg3jwt56ttrTUyipTIwa62ELmwuiOC+3YHMZBMfWcLTaLsJYQ7dQHi1jladXg248Zr1rmXN3HrOTPbCev96Xh/AbPEJY1LygA71s3iLqAmAQQv/hWAIVfm67u3ILn8F5k9PEAqg5l/xjxsbpiqeo9ydXhm7fmRjrOqtsXQa3A/m2En4Phjr3KfxxqL/MIAA/3wMTzlKIDo/IWGAo6TpeNicp/WQ7hTCxzDqvLFTBbDWuws91Fu+cWNlsR6sN4H3qI+6T2Q5b+m1U5bfp+UYQZfWp4VdICUFPvh9mfRGUmoLH9ytWdhbQYNJ2NYvKctazxacD+bY0cMW+mwXLRJ9DioP/s9pwAGL51uOOg46bIt3KftQwSeEsaq8sVM4GJNlPecFuNOlG60JNaD9f50lbMDilLyvLBhJgagfh+iJJBxM8i8ZiTBncxr4scwq3szcwKVvOoOzR2Hj2eOMXRlzzuOPKDqTaS8z4AAhtsJMP3U4nHQinC2N938hYjNvsPQazrtsl3oEmUvIQJPCWPV1NBs4FJOlACy5NXyidVO88ZcaW0M49unlUM6RT8Bftv8HUOS51W1LZh35nx86R+WcUY+WDYcN//NrVk2OUBq4tKazFfu/C9M37fJlgmUdxx5FqY3kcZmzGOWaAF1YOOZi/KzpqM5778xmpEoZf2rOyKGTFAXTRiOcl/2mReSRNlLCHFBiWNF+WJFsqm3UGu3eWMurO6I4Ne7BzKqtgQF7vBNwmlXBpiqtrqulGOCfMH/jtNvQlttMygIRoYIjsYolLN2JJ7EX06dgDN7X2equgJIYumu5ZgPYHXHlcz316gCkKX2SyKV8XS8vwA/bZyJC2dcqfm+vFDbjKsDw1DDsAGSBzYtgYhZubsRjMi4WYKY915ahxvbn4JvUFvI0toYxl/3d+PJvZWeULUVMiLwCCyRa+C696OjGEhmT7NKRV2+zRt5pZu5sYn4+CG1K/OK4DZM2rU84ygwPtqDpbtSpadVtS04GKUI+oCwD5DPjwejFMN2fKgpJZbWYc7enK1yA8wpAFmSb6nsMT7ag2XtTyHePRLxRvb7HtiyAdOeXYrq+DEkkV0yOe4PwifLSLW6e5qVuxvBiIyb9Zn+pKMNgUH2OJXv3xWjE7i9pbiEBPlAlNoElmCVLYzghYVas2OcuvE5lY2NfOEeSFnwDCg12gDqBvS7htRHe7LWLFZ3RLDw4ZcQWHq/qfbNkp0Mra5VBbvAIP95UoCr6+8BwZDcPONmcHp2lqUlEIk3T8PglMtBfb5UKdLnw+CUyy1daFzfvQV/3bYAsXe/g473F+D7h99XlY5Zn52WYETgDCLjKQLy1XfDygbUXHe3mx2f3vui9RizYzSycK90zF5bfR6u6t2uu3ESAHoDwwCkJs/VHRG899I6PPqX5Qiodhupx6M8zxXBbbiIUepinYf03HffXIrxigDnQyroNF7wiKobLK+cBh9B8JmHUbbpDZBkeuzJJMo2vYFk09k5BR8pKJanxzc+2oP//POTiJ5flZW9sT5TrgtE9WjV+3bLOD94KmMrG55LDZHxeJjVHRE0/Gofbtl42DEJqRZWVEROL9QakdbqPeYbdWx/Nd7tegv3s7o3Y+mu5Rgf7YEPqcnxtn1vYny0x1DgGR7vx6zuzagL+3HvR0fxkw61EWrWeIYNZ55ny56NmLRmCfeY8vOQP5eXGdRHe5jS+Pi5k9n7xJLJlDWOiSxND62ynhyWIOanjTMRL1MLRjZOvVH1/bhvd5D522LJw0MrHrTUqbWYEYHHo0gTgrTLXo6dElItrJTLWhvDuOe0mGNtgo0ERb3H/LaLPanzbmcp8Y6nXZyBlMuAMlAwvdsAsN7BciRw36dtWDRhOLo0AkHWC0F9nlrO1sp1FvlzeXuT9lXUqKTxgS0bUhkNZ2jc23Msbxnd98VScl4440rEv6tuMT43NpFr86TEaOATpBClNo/CmjTluLFWYrVc5uRCrZGgqPcYs4GVpcT7cOqNePf4uSBRoJ7R5I0Fz0wUAOoHetHaGMa9Hx3llogkSOQoc7y8gEUBROf+MKs8JH8uyx+OBkOouWk+rjgx+xisidgIuarazKjkWIKYeKNayNK1gt3ThvX5F8KGZy8hMh6PohdY3NjUVsjW60Y2++k9JpcNg8o+MBOvuQqv/d1AyoLI4KTaGarhZhfSRLpownD8tHEmIj5+uwXpscrx8l+7VjX5yp+7qrYF88+Yh89DNUjKMgPWOoaRCVd52WRF1ebEvi8zn38hbHj2EiLweBStyc+tyd/qBlQnMRIUF00YjjJGaiGt4dgdWJmTo+IxEV8QD505C59+c47mRNraGMaFM67EPV+9NbNJlfdY5XksbJiJ44qAxZuklc9dVduCs6cswTOL1nI3eAa2bEj1ktCABkMYvHS6qryV62J8vHkaonPV5TIri/tmbJ60Ap+Rzamlhii1eRTWngUAGBUieOB89yb/QrVeN7q5laUHe/aTfkyujVjaIMvihdpmbD/rFty5cxXqo704XlWD4IRmBP6wNVOa88+Yh8XpyTI6OqSpkmptDKP1jhkAZiCqoahSnsfmhqn48G8qmZtgJSSF1pzeA2g14WCd2VOUVL+z0jdV7jJtrIeqMeze98X6/G8Zd5x57rwNzwAsuawXK8KdugAxeo75klHbdex8f5Y8t2FA23HY7Lm3t7dju29sTm7STjEk/e0GfD4gmcwEBADMzZ1aGYT0WfKcq6nPh+gtP/L8ZGv2O1sILuu54LQ7tch4PEy+sg07GsgVAmbaXcv3aFxWXo23Tp2JVbUtWecO8LMjnoLurq19mu9ZLgFebz+JqjVCOjuRrsZpWYir0NILHNy1nSQ1HHSKaT+MEB2wEYFHYBqrDeTsxMokpdXuWr6GppyoTx7ItsKRAkh/AsxgfB74Qe5gjGJ1B7v9dy4B3oh9jpbijMSiAO8+I4IBix5suTYAzGf2r4UTnnTFgBAXCExTKHY3gS0bEHrqgexNe089YHjxlrcRtIyk7pOsgHqeXarZwwZIBRCtPUFSIJvVvRkd7y9A7N0b0PH+AjzyyQpM+8Vs5sJzLht0jewnyfVq28hkGZsxT7UZM15mXF2Wy34Yp/vwWKEQXNYLERF4XKRY1C0F05fk2UdB4vGsm0g8juCvHtV9quQ8rSTsB246vQK/3j2QmcjG9rP3vZwc7cGs7s2ax5GC8aIJw7nOBSnvM/Vu91wCPDeopBvVAfoBhA47IefJ8oXaZtzadHNacp2y07m16Wa8UNus+1yt8WsFS6f78FjBCbVdMSBKbS6RawnBKayUJoy4ADvN6o4I5vSr2yQDqU6YevA24I4q9+O3XbGs+3gbNQmApbuWI0iAtfVTcDCqfr26sB8jd2zF7PfWYk5vt2pjqOpv2VoKrxToI6m23KzPjVfa6QxVZ8p01zOcoTPPD4YQu+EHAPgtKVjlTZx4CoDU+9o5ugXPjs7uOfSOwTKs2dLU6o4It1xqNQPXOk8zxJun4YXa5qHfW5cfizjlVSsUarmRhch4XKKQLDWsliYKYf/O9lfXW3q+VjahvG9hA3+jZjgZw5K9q/HA+VWZPR/yctof37wFJ7+2Ml0ONIZ0dc/aRwKk+gLxPjdWaSeStu2RsoCMM7TidSmQcYhWboSNN09LZezf/xZCSxerPMlG7tiaef9YGA0CZkpT0veYh5UMnOe9Jp2nGdwoBRZyuZGFCDwuUUjqFjtKE7m0XbaLwJYNuO+Py/gml+ETuM+TSp2ffrCAWSarC/tVE5a0Y5+n/Q/39WS6lV6vKKcNjx6FL2HuyruzvDojOJB3QGWh/Nyk0o681DX/jHmZrqdSAAj8YSsz+wr8gT2xZibiY0eYWdrYd14GYL0Ma6Y0pWUbZTUD510oSudpBjdKgYVcbmQhSm0uYbe6xUpaXSjigFwJvrQcPo7JZZT4Qb/zA9XtylJnXX8PlsmUaUD2ZHXrxsNZgWZVbQsW72nj2ucDKfPQdxlGoGaI+IK459SZ2JyeMOQdUHlIn9vQd+Is+JofYT5PCgBmL4T0vNcCfQcRXrEXI0MEZQQYlB27jKS6rMrLg9d3b+GW8oxuBJXO+5FPVmD+vrfhRxIJ+LB07KWonHeHpYsh3vtQ1nfQ9KZXN35vXvtNi4zHJexUt1hNqwtGHJAjvEmBAnj/mjtUe1Yq77wOoaWLVRNnZTKGn3/WxiwXsuZ6VsmNImX/D8CYY7QGcfgy2UlXJKFrBCtRF/arvhOsp8kDq1lvMb3MvDNUA4pUN1VCgFFBAoL0/0nqdum7+t5L6xB4mt9CwGhzwbqwH498sgK37XsTASRTGRuSuG3fm7hh0zLN8erBex8Gq0aZfi03fm9e+02LwOMSdqpbrKbVhWzuaQTepBCpOhETrxlqSZ1dp2czrr8PdelxAAAgAElEQVQ3q1yotW6wqrYFv6m7KCsoEQBlm97AtlfWwkeAXv+wnM4p4gtizpnfy2RfdWG/oavVoC/1efKClJ+AGVj1LoSUCkyptw9v7FLrByDVZTVc5sOhueMQLvMhpnDP0Wo1beaiatGE4Zi/721mybDs3de445XQUpny3p99l1yr+7qscTr9e/Pab1qU2lzELi8pq2m13R5kbhNjqLJoMITArFshF1cbseZXBjG9dYN/6FN3CyWxKBrWr0TrqcdxQlIt0U6AgAwbDsJYHwFSWZN8HUaaMO796ChXsZV5bnqovM8+SZFyxlbA8xaTRARKBSb1+0EDgSz5OgXQ4x+GO06/KTN2Ca3WElqtps1sTm5tDMPP6b4KhlecHD2VKe/9OXTiKWB7e/Nx4/fmtd+0CDwexI620YVq7mkErUlTDrPtsgyKVBD71y0HsfKTft21lCUtVQi/xZ406wZ6sHhPG0I0rrqPVlSi/7FXUXnndcwxfRGqweoxLSAUmbbUUx97DnN6D6CzvBr3pO15WAzS1GSTy3eCdyHEXFhPJJAMnwBaVZH1nk/qOkfzuGZbTZu+qEr7zDFv10BLZaq51tTervm6PNz4vXnpNy1KbR7Ea2m1E7DkvnKMbs6t3X0WntqlH3Tqw360NoY121vzruT9/akyUWzGPNUakVSmkjKTnXUf46J1j2bWP04e6MGyT5bjeo2Nql2RRNZ3Qi7n/njTD0xvVOYKDyJHVe+53nfRTKvp2Ix5ptcqBi++mi0Lv/hqztmlz6WAVKaliMh4PIjX0up8EHxpue6+mS9CNWBdSM/q3ozFe9pQH+1BZ6gGP22ciQtnXAkgFTwCTz+YtUaRRMrFIAEffIzSz97yapy6Yi/qwufg78+ah5/tfBY1iWMAkOmLI02srCvxykQMz+58HPd/1oa7Txlqoy2N76EzZ6G1cQaA1P6m+2RdQof1fQkqKyExN0UiO3uk4eEgEfUmXF43TwD4tw8OojvqU30XWd/VC2dciXj3SPhkx9w49UbMTWdPBNniDq2LqtjsOwCk13SSScDnw+DFV2du5yE81PKLCDwexUtpdT7Qu3JVLohLSLY20sQ9PtqDJz55CrR7JOKNqR3o7zXdjJ90tGUyHKlsEEBS1bY64gvi7lOvyyyUDySBShrLPObE+DEs27UcH/5NJYCr+FfiSEnAV+x8AgHQzPPHR3vwH396EvEtVWhtnobZnatVUnP5RmXVusbynwOEZNZuSG83aCAA6veDyPYfsRSY8p49N1SNRHLWPzFLd3qtpofMUFPHk95DilSmqXdRFZt9h26gUT2Hs05Y6h5qbiFKbYKiZG9FNfN2CvWmSjmLGftwQvGhifvej47i2dEtaLzgEXSGalQ/IIKULJq1eRMAftqhfv3KZAxTNz6XGp+GegwAymRBRyIwODQ+rRISb+2G5XdHK8KaCkzlzv5g38EsSbQZWIICKeg4tTlZeKjlF5HxCGwl335R0vFbTpmZlbkAqeyDF3AktBRXQPYiN++xPiQRvPhXmb/lpTte+Y/0HsC2V9bigkgEIc5jtJDGp1VCMrN+QY4dReSxV7n3G1mcN0quKk2rfXvs7lgqMI4IPALbyHeDOPnxpeAiXwtZ2MBXh0lwFVfh4ai88zpEew+gM1SNtdXncV8jCR9i796AzlAN1lafhzn7N+q6GUSqatCwfiVCNLed5tLaRGzGPISW/zy7TOb3IzZj3lDHUQN8EarGxW37uRcOdi7O56LIKzTTXYE5ROAR5Iwyu4kMJnNqEKe8ct049UbMjU00nTUpSzaraluyAk2ZctVawc1nVOCnB2bi0b9kZ0oJnx++6PHMgvv4aA/+ad+bzDo1RWqtR+9xcuJlISw89To8vP2XOo9kI8nCMxBFXpX+m7mu4fdnrfEAQJQEUBkfQMe6Weh6qwbbrpiTtTEXsHdxPhe3czszLoH7iDUeDYxad+STfPX4Ye0wPxhjz+paJROWC/CkNUvQsmejaTsgreP4ADx+4QiMCrKLXfVhP/6jeRQunHEl7vnqrWmTTYJjVSeCVIZV6yC8H47y1fV+YBTAluFNeHzEBegMmd2amJYOXzo9e28Tp0dR8KXlQCwK6vOBAql1jXl3I3rzXUhW1yIJgi8Dw0ApxYmJY/Ah1W110qtLENiyIev3cEd9K1cSbZZc3M6FHNrbiIyHQ77LRkbIZ7nBqI8YoF0yYcqH0909pWyFlTWx6vt1YfZmRiAleZaer3V13doYRusdMwDMwPH0/WT2JYbOMxcIgCkHd+C67s1Y2KBel1Kq5OS30+pa1boGd0I+dmSoT1EymQkS8s2SI1fsxe73F+DE+LGs51YmYji2ahlun/iVzPv22IgLcLQpiSV7VyPc14NBDVWbEcyqNIUc2ts4nvEQQi4nhOwihOwmhNzt9PHswgs24/ns8WPUnkevZMKbKJUL9/Lj8XqlrAhu0x1PLlfXXENNxd88kxa98OwDMoF2/hnzMi0N4vDxWz9U1zI3znJNPpV/M74ndWE/VzBR2dej+j08O7oF50x5FJFn3sGffvCAqyUu0VLa2zgaeAghfgCPAbgCwFkArieEnOXkMe3CCzbj+Sw38LKYUSFiaFIPbNmArzx6F3jTsrLsJD8eL+BO3fgceMmVvMQm7yW0s+5jzH7su5qlSt4kt3PiN9FVMdT35omxX2c6E/RMuCgt2+UjTfiralNS7YpLfgU/5xmqNR29sXKOqfyeLJowHF3l7HJfZ4gtT8/X70HIob2N06W2SQB2U0r3AAAhZBWA6QD+7PBxLWOHH5rT5LPHD29B+IHz9TuRKkuESo4rNncqsyatgPufLSNw23uHVf1gHphcpTsOVqlSKulJayNIJkGra7Fx6o24+sjfov/8G7Je8/2q01WuAnOuOBcVTU1Y+PBLeHj7L5lXe8pAm6Cpyf5kTgYSevJ++Np3qDZOsnzsMNBvyImgtTGMbVfMQc2rS1CZGCr3HfcHcf/p6s22QH5/D0IO7V2cDjzjAHTK/u4CcL7Dx7SFXJQ2bmPn7mvemtbW7ih+2xVTBSMrtj0812hp3eLDqTdic2wiCOd1tQKumXHpKaNUAVK2NjK365zMTns5ciVdhZ9gSUsVRv7xv1H5xI/xn70HcMxfjnBiICv48FwU7mGs+QDpslkyibK31wAAM/iwNntmfU8AJGvVrtUTr7kKHwJoWL8SdQNDMvS2arUMvdB+DwLvQCg1tkCc04sT0grg7yml89J/3whgEqX0BwDQ19eXOXh7jq6vTrL+gB+Pf16G7ihBbYjitvGDuGJ04ZTaAGDkjq0Y+87LKOs7iMGqUdh3ybU4dPZk069z9e/KsT/KEwgPlanKfRT3nBaz9D6c9++3cBfMt//4SQDa5zVyx1bUr3sO/kGZ5LksiM4rbzR17nrj+MqjdyHYd1B1f6xqFMJfW8LwEEg9mwCZ78u3D2xWjTVKAjjiK0d14pju/iJp8+nJnM2nlPiwfeFS5nPl3983/3Afph76c9ZrUAA9Ey5C1xXfyXoe/7uQfZ7/OCaOu08b1HmcoFRpamrK/Luqqirr6+t0xtMFoF72dx2AfawHygepRXt7u+HHWqWpCbhde7+hI5g6x6YmxK69MdOOtyb9n1m6N+3l3JM93Q0kCZ7cW4nbW8bkcJQUtLqWvZHR58OZX34GAAi9/nzmCj3YdxDjX3kKY4/0pK7um5owWDsGRFZOGpwxDzXN07LOPUv5Fh4OkNSOfEkFxxvH3ooaTNpUiWjfIeb4y/oOoS4cYJZi68MBfDxz6L2pvPNH8A0qLHhoHMcD5TjpQnbAkCNlUIPv3sB+AE0yvyurOyK4v2Mog2059Bdmw7Tq37+HSWP+KSs77I4e1h0XQPDh0XI0NZ0CwN3fZT4R52kPTgee3wFoIoScCmAvgFkAvu3wMQVpzKzZ8Na0WFhdUGaVCAGAJJMIrXgQNBhS3weg7O01SDadnSkladX3Ves3sjUOaS1ncMrlKNv0RtaxjvuCuPuUmak9RKFqpotBAgQd62apshVW6cmoak+PBHyZjalZcPrO3PVBH6bv26Rv1UOTmc9dKq+ODBLmniyla/ePG2YCYK/9CARaOKpqo5TGAfwzgP8BsBNAG6X0T04eU5CCtcHz1o2HMYKzGZbVN4U3WVldUJYUSZSov34kFh3ab6K8DzAsFdfrPkpiUQT+sDVLGdVVUYNbZV5uCxtmqlRqkjOBDylngqW7lmNW92aMChGmgk+rf48Wo0Ik5bSQZunYSw33nVndEcE3vtiEpbuWY3y0Bz7wP8uEYgroT1CAQPVdkFy7pdcbH+3B0k+Wu7ZhWVBcOL6Ph1L6OqX0dEppI6V0sdPHE6TgOf4CbDcA1v6W755RYbnhHM9ZId48DaDa7YlZkN5uQw4NRiTlpPdAVkO5U89/JGutxci+mnB6s2s44MP13VtU5xo/dzIzYGh5vREAe749Fo9fOCLzeTz4tXnYOfGbGdcB6vNh8NLpzHYA9350lOmyrYQiFdCUHIrSrO/CqCDBfYzXq0zEbN0zxnIK8YJ7iMA8wrmgSNErh7HcAFi7xyfX8st1eu7AenLlwapRzMV7Gj4BiBxhXqUTxusoWd0RwWXl1Th5QLucpcxGWOVGuUotxllnqY/2YMqejQhteEp1rrQsxFxbmX5oOxZwxiVllOrP438hgv+leU5A6rPnlfIyQdDnw3N1l2FBwxzm8ZXHDv+2l/l6du0ZY6kqv7/pMChFRhpfiO4hgtwQXm02UyhXaEbKYUbWauSbLeW9UVjuAYGnH8TCh1/KnHt81TJNZ4V9l1zL3n3+nR9g8NLpmhsutRwa7v3oKO45VV0mUx1HITtnlRuDPmRKXrzyWGeoBvd/1sY8V9b+GQAY19+LN0f+P+zZmmpR3fH+Aszq3myLRLku7OeOlVbXIvLMu4iseBuJOf9iOKPlujfYZFHDytBjSWTtxwIKzz1EkBsi8NgIa13FqMGl3bAmUSVW1mpYayiBwSju3Lkqc+6VfV8ynytdJR86ezJ393ls9h2Izl+oueNfebUtlfU61s3C4j1tWDlmaqZM9qV/GHoCwzR3ubPKjY9NGZEpebHWfCK+IH7cMBPj+tkZAQ8aHo6L1j2KkweG1kyWfbIcr53we82reSMXNosmDMdPG9VjjZdlB1sz9kFOW9SYEawUknuIIDdEqc1GtPzd3C4NSMe7a2sfU6Fk9craiFqL29tGdpWspU6T7qu88zpdhwZlWW98tAdz9m/MavwmdbTUgmdWeX33FnS/2YbKZAxx+OBDMqNqAwD4CJBUv8+8fUKxJFCuNEdNpDqRHle0IJAwalzb2hgGZlyJe1714c6dq1Af7cXxqhoEZt3KDLZGvpssRwSzjde0MKOqLCT3EEFuiMBjI4Xo79bPOPSoIMEDk9lXtkYl2Dz3AHmJh+m2nMNVshGHBlYGFpa5XFsJtMFnHkbZ22swPv13AMkst4FlnywHSaqFEjxnaQAI9nOUexprJmYubFgu29nNEszjpEUNyykk6EPWGg8g3BKKBRF4LCKfqH0k5bGlJF9XaLzWBeEyHzfoGG0FwQoGSuuXVbUtqA758HDnaktXyUautrUysHoLLbgDWzag7O01XDUbIcjyNZOgPh9oxTDmGs/eihokktDNBpWYvbCx2hraTXhWR6zbhLDA+4jAYwHlRM0KOvm8QjM7UZm5olYGg0hVDW4f14pVo4fkyBV+gvO+dQWON86wchqZ42lNmrwMDNW1uuU1LYIvLedmLfXRXu59SFLEvvMDZqb2o/QGVWU2GPEF4dfIBs0Y13qxNTSv7CcCTfEhxAUW4GUUfgLDvV7MYFYxx8u0eLcbCVTyMZzZdQ6e+f7TiDzzDrBkNS6ccaXmQrVyT8/IHVt1ztg4rMXviC+IO+pbLYk7tEpfJLWjhnkfrR6NF2qb8S9n3YLPQ9WZbqbRuT/Epoapqj1Cn4dqcM9Xb82Yk7L2PrEEI7wLm3z2ahII9BAZjwV4E3WSAofmqp1/rZBLR1SzDtt6V9R6Y9BaqGZdgdevew6DtWPwQm0ztr+6XnchXItMJ81Vy1DZ14POUHXKzmbEBahQuGzfdvh9LP70RYT7enRLUNxMCvz1GxoMYePUG1PtGUZcgMcuuAAAcEP3Zjy2ahn29PWgs7war40a2kRKCHDtqeXsTGXpYoSWLsbs6lqcOvVGzI1N1C09FVNraDPWTwJvIAKPBdzs2ZOLYs5s6wK9QGVFtce6AvcPxtC/ahneG3cIj/5lqOw0rO9LxJ9ml4W01i3izdNwTpe6/XV/guLpXf2gSFm/3CcrcUklKF/7DgT+sFX1uhun3ohJa5agUscFAMhuRz3z069gUOb8Pqt7Mx6XHffkgR7ctu/NTPA6eaAH9eseTW04ZfjUSWOdsvZR7PyufsOzfPZqshMvtKAXmMfRtgh6yNsiGKWQ3GGVPwpgqAeLlR8F6xxHrtjLLOoQ2JtdaU0wvDEAqdKatrtBNzNDSIJwzTiT6fbOEqy+MvGyEP71K7fg8REX4LbD76ezph5uu4GO9xewF/WRncFI57m3ogZrRp6Hq3q3c1sTZD+HgFaPxk2j/zHr2Lzj6o2DRYL4cOe5t+HxERdwgwCzB08wlFOXTiPf8/b2dmz3jbU9OJ3Ttp/jAj4kjXczKBbS/OMkdp+n220RihorzdDM4lZ2pVUu442BACqHYyC190Wr0yiQcoDm2bsoy0Jam1Z7G5JZmYxk4AkgKwBwj8X5u66/B3Oiqf1Ai/e0aQaP1HMoSG83lh7KPrZZN2ot/DSJ+/64DL1nJLGqtiXr/ZZng4NTLmdmcWYxkumuP+DPasNgV2ait+4oMiJvIsQFFuFZytiNmYVlp+A5WCuzIGlS0nOIpkhNyEnO11BZFtKSTLNMMTOSZ9lteq7QLKTXYTkX6D3H7HHpsBNUIgm91+9PUGx/db3Kwqhs0xuIzZiHyDPv4PhDL+asZjMiOnn88zJucLKCnkBGKygKChcReDyCGXsTLXiKqVzHwCu9dUUS3EBBASSRClo+pDZkKl9Hae8CaLcY4GUU9dHejMv2rO7NqEwMMN2i9aiP9qjdqolPQ9cGnBztQcf7KS+2yvgAoiS7wMAbx+CUy0E5fXaUY5K4c+cqx1RsRtSR3VF2gdDq5mneBdc36oLcMpwdxxU4iwg8HsJqdsUy9gyteNB08JGPoV5jUuIFigR8qi8eQartgCQ7jjMW0HmS6YUNM/kZRfVo/EfzKLx2wu+x7JPlODF+LJMBUQDJ8AkpQ1KdLEN6/VW1LWi84BFUXfZrPPdvr2HlonXoLOdnM1L/mhMTxxAgFMlhJ2T84gYvnY5k+IRMwCEAfMeOoGzTGxi8+GrEy4yNCdAoIdqgYjOSbdeG2OHXaimYdbHz7dPK8evdA5oWO8JWp7ARgaeEcGJvh9akxAsUPlYnTQB+UBxP7wlilYWkBnLJ6trM3hfJi43ZtE1mqzN143MqhwECAOUVKUPSuT/kZi4UwP1NM5nZJs8JW8ross4vmQBCFZnSV2z2HUB5hXp9Kd2kLv7dH+JY1YlIpl9PjtwlosJPsJcT/PZWVHPOyjhGsu3bxg86VgpWXuz8tivG3D9n93EFziHEBUUIT+XjxN4OlsDiG3VB3PvRUdwaOQu3nXULFu1+ESOPfZlRmvEW6Y1IfSUHA2WZRVrEv//TNtQN9KoW0/nn3p1q2NY8DfFVyzCM4ajd4x+Gp2pacIjhgNAVSWSOLW8LbTQD0fpM4s3TgOZpOI5sGXmkqgb/XXUeFu9pw7M7n8Dxqho8O+o8zN6/UeWE8KNTZuIJ5hHMoWcmesXoBE4aU+WK0EarjGbFHkngHiLwFBlaKp/ZNu/tkJBPSsrjPzbiAjx9fjMmnRjAxv2DmazCqnkoa8/RmrFTcFnrlaZMTQmA0IoHsfVAFM+Oa8WjR9U2NnecfpPmOkdnOvjI1XPdm+ejevCYehzh7Ctxo/tt5JZBgS0bcKNMLTis70vMObYRK8ZMxVW927Pk5JsbpjLH7QRGna6twlNXGnEfFxQGotRWZGipfJzuqaJ1/D1Hk1g2NdXX5sXaFtzz1VvRN6xasz+OFq2NYbx2wu/x2QepxfvPPlig2cuGde4SJBbFaW+sxLOj1TY288+YhzVjp3BLN7xSYzjA+Wkp6mq5fCaskmllIoarD25H4wWPIHjxr9B4wSOa4/YyhaDwFFhDZDxFhpb01emeKnrHz74inoH29nNz3qQW2LIBF617NDMB1/X3YNzLD2Kwrz21dqIg3jwNWw9EcdHLDzI3aI7tT5XGlJkLADypoR7k7eUKvcWW85Jj2bfn8plwZeUDvdyNvBLFYD/j5v45gTOIwFNASHX883q7M9YrgLlJSW+jqVM9VaSxR3sPDPmkySZwu1VGTKEEgLK31yDZdDbzHOfGJuJdTnO6fRXsxfn6sF93QpMCqjSp37rxMKaXDcOoQXXw2VtRjVNX7M2aLM1+Jrzy3BehlJBg2dQRltteFDqs9/zej46KAOQRROCxmVyvKJnmkE89AFAKkkgM3aZjbW/WGNQOWN0/5a4B0vGVPmsjL7wKkGU8ZvrHcBflkQpKrOd1RRLM5nQRXxC7L5+DiiME0/dtyogEusprsOeKOQDYHUHlyCf1Wd2bUTl4XPWYKAng7nRLBCOTvvy7lGVsGh4OGgiAxIdau2Vk5Rqva0eH3Fy+305lWcUUSEsNEXhsxMoPgXkFH1f3jJTkz7wJOR9lCK3un5sbpmLRhOEq+xzS242TX1sJbFgNEjkKOmw4SH/EcJDVdI1mBKXVHRH4CJgKtIfOnIXF11yF115Zi0nvLs/IriXjzsG+dl3rGfmkvnhPG8qhzjqP+MqzskCtSV8ZyLKMTSNHQP1+JIedABw7qsowea9rtUOu1vf7vByeY/U7mWsgLYZyo9cRgcckWl9aK1eUZiTNeo91Q10kz054e/dPjvZmVEbBx9TByZdIAOkOneSYulOnPMgqs6H4uZOZnUEBtSJMmvykj0a+jlPhJ/ig73mE516Ki5JJ5p4a+XF4AVE+efOk1NUJtcrNSFM+lh0QSSRAQxUITXyC+e6zXjdXvz/pO896rvT9/u+vsZ9rR5bFI5dAKrKkwkCo2kwgfWk7I4msconUaMzKFaUZSTPrsUornG2vrDXVNM4MSgcEbl8a2Thz3StEeg8wHRfKNr2BxFl/q7adYSjCtBr2fdD3PM7c9joII+hkxqD8m7HpVj5581wUkvBhVvdm7vPkGAlkpPeAqWZ/uajB5N95Hlrfb6tZlhZmGx0CwtutUBCBxwR6X9pcfggSTFltIADqz34ua2JlTcyT1ixBy56NzACZK1L30e6VSzXNP1njzHWvEK0ezXVc8P9le+oxvpRnGk+WzZvkZu7fnAo6OYxLGUhXBLdhz1a+LxuQ8qRbumt5JvjoNeWT4AUyWj3aVDDJxe+PF7R5YzV6n/x2s511JXIJpHYHwlzHXuqIwGMCvS+tlf0FcjuYzCR6812Izrs7fRt/vwtzX4fCHVnvqk7PPFR+5cu7Ak8ZZrLHqbWPhocUvLhCgnSWQpJpQ5mBfubjWJPfrO7NWPrJct2gw7XRkQXSba+sxaRXl+DkgSFfNkopEoxXl9a+9CZ9+XdGyw7IbDAx6/enNyHrfb/1fhN6VQQtcgmkVi4OlVgZe6kj1nhMoFcjt7qwL8lqlU2Y9KS2Wu0C5PAmEaaiLr2bX2qz7CPIrJF0ciTJVNG4TXlugEwaPmw4cDwCX3JoTDQQAA1VpsQGskV8qZGc5nuA1KJ7aOli+Np3ZO3lYSn97vu0TeXdlnUuAL4I1eCklhaUbXpD1VBNyuZWd0Rw2fqVqtcqR0LTtVpvh31rYxh3be3DwRhVCSL2VdSg5qb5mffUSWkx7zsPZNvTtLfzzwPg/yasrgGZXc+0U/Xp5PpVsSMCjwmMfGndsg2Rw1N4SWsK0sTFu6rjlbIa1q9E5+TUqrH898WSJBtxQFDuV+l5+Tmc/N5aXfl0bMY83YZymXFDvZeHNfnVD/RyX4MCeHzs1/Hg1+bh45ljkGw6myvzvvejo5g9YK7JGwUy/nBaPDC5KvN9kwQRvA63Ti2a877zZlpyaP0mnFwD4o0FsEf16fbYiwkReExQqDumeROztKYAQNM+hZcx1XEm1Owr8F4gRweEQ2dPRs21N+o+Lt48Db72HVwVmxLWXh7l5Ee3cjzSkAo6d5/5XSxJv19aGzy7IgluBtjjH4bqxDFVPdvHGB8LM983p66+nf7Ou9VZV45dF4f5GHuxIAKPSfKR0eghTWChJ+8fWu9IE07G8PPP2rjmmQA/Y9LqmrmqtgVrxk7JqRldLgT+sNWUCEBPRccK1sf9Qcw/fR42NUzFEoOTa13Yz92UesfpN+HZnY/nND4Jo983J6++zX7nzeyTyceGZ7vw8tjzjRAXFAnx5mlAkr2qMK6/V3PiiJ87mdkNc221elugn8BSB1QWRpRBWt1MmbfrqOiUYg7q86EiEcMzB36DnXUfGz6vRROGY83YKUxz0VW1LZqKNDuxc9HcCmYX3O3qrJsPvDz2fCMyniLCqMW+ElY2QQBcfXA7FshuM1vbN4LRtQmuU0GoAjTanzV+o47bmUyRIayQ36/FUClqKhoV5qJA7uthcoxYCRXK1XcuJb9CrCIYxctjzyci4ykicm17oOd27OTVnNENfbxzi865E9H5C3Ul5zzs6Mra2hhmyoaBVEly/hk3o3fYiTmNz2i78kK5+hYL7gIjiIyniMi17YFWpuR0Yy29iSrraj88HDQYAjl2VHVuuTpu29WVVWuj5eC55yN093XIZXeHVmCUzlm5psJzp3YDseAuMIIIPEVGLm0PWAvtvEzJLoPFwJYN+MqqJxDtO8Rto6DaXxQ5kspy5t9jW2uHXMuTSngBlAC4+7TBXIaWer5OYCw077FCKfvz6zYAABAWSURBVPkJChtRahMoFtr55SBpkmvZsxG731+AjnWzcNn/mY1tr6zNPEbPAUF6TGjFgwj2HYQPNNNGQWklY0cZTA+7urI6tbjPC4DS7YXmPeZUyU9PgCKsa7yFYxkPIeR/A7gFwJfpm+6hlL7u1PEE1jCSKd370VFM37cpa7H85IEe1Ly6BL7RqcnbyEK9kTYKrY1h28pgWtjVlVXzSj+ZexDQy0YLcU3F7gV3vayu0LI+gT5Ol9oeppQ+6PAxBC7RFUkwLforEzEk01mI3noEwA8c8jYKgH1lMD2UQVfK2swEIq2Nljw7GaNjA/iBsRTWVPSUcm5Z16zuiODffleO7k17C2bzuFcRazwCw9SF/ZoW/TxIbzfCsy/JTJpGA4qZtSe74PnWAfoCBqektfHmaXihtnkoqHX5sagjklHTFfuail5W50bWN5RVpVYnRFZlDacDzz8TQm4CsA3Av1JKDzl8PEEOsPaJAOqr7EUTmtH1Vg1OZljpSEGDFVBSImOamcQHp1yuabwpYVcZDNBoI614TSMqMrfRKiUBQIUf6E/PsaNCBA+cXwUAOKdtf0FZO+WKXlbnRtYnDEHthVCq3WtD88mEvAmApbddCGArgB6kNpf/DMBJlNLvyh/U19eXOXi7lXqEIGdG7tiK+nXPwT84VD5L+vwASXcITZMoC6Lzyhuxvc+Hb/zfZ1ApK7dJ9wFQvRaLWNUo7LvkWox952WU9R3EYPrvQ2dPtvnsUqw/4Md9u4MYSJJUOwTFhk5p/IfOnozz/v0WdldTANt//KQj49Pj6t+VY39UrQOqClBEk8BAcmjE5T6KK0fHse5AQHX7PafFcMVo7+2nkX9+EvLz0bvfDiZtqmC2PCSg+HAKux1HqSN32K+qqsp68ywFHqMQQk4BsJZSerb8dnngMYqyZUAx4uY5Vt55HXw6LQckkum2B1o76ZUtsdmTOEHkmXdsOU8j8u5z2vZnrog73l/ANPSUzo33fiQ1Wj7IYb03O088xdJ5jlyxl2sNxMJPst3EJerDfsf2ZTn9ndX7nO2S+fOQf4fkOPme5hO7P09l4HFS1XYSpfSv6T+vBbDDqWMJcseMQkx6LEsBp5xwo/Pv4fbRsUscYETNtLojkjVh6K1RWVlX4q0PjfzmdwALP2KtnjgseA1DveweoLd+5rR1TSmspbmJk/t4/g8h5GNCyB8BXALgDr0nCNzHTBDgPZZn6xI/d7Ite2R46O1hkQKTHD3TTqN7mljw1ofGvvOy4XNiweviOSrI9utmOPcAKC6lm9tI+5PGhJLCENQGHMt4KKX6jVYEeYd5he/3A4SAxONDt2kEDN6EG/jDVkTn/tAWcQALPTUTKzAZMe3Mxf0B4GePZX0Hob3qpQ1Pqg2AeRX+7dPK8evdA+Lq3GZaG8M4L7mv6Ev9biDk1CUOTznGuo03GWtt9Mx1EjeCnpqJFZgkW57/2reaqWqzAk8mPlg1yvJra5WSWGsbk2udXfMQCKwgAo+AGxyMTsZWN3rmujCsV3fnBabNDVOBu2fmZNqpBW99aN8l16IGxtobmIUXkIRdv6CQEV5tJYIRD7VcseJ3ZrZxmBw9XzDe2ohTJSfe+tChsycbbm9gB8K3TFDoiIynBLCyG98IVjZ6Wt2Yp3Vlr2Vj4xTM7LG93bWNqcK3TOAFROApAdyY9HJdy3Ha7qRQSk5uGJ4CYoe9wBuIUlsJ4NaklwtOtRMoNPTaG9hFIbpVCwRKROApAdya9HLB7XUYK1hZJ7Or748evIA9MsTZ3CMQ5AEReEoAvUnPSeGBHk41DrObba+sRXL5L3IWB1jZmGqGRROGI8j4VR+NUSEyEBQMYo2nBNBa/GcKD556AMHnHwWJHM1Z9mtGIm11HcYJmbKc1R0RXLZ+JSoT2dtAza6TObmnSaK1MYy7tvbhYCx7nWeQAt977zBu3XhY7OsR5B1PBx6nJ5xigjfpMYUH8ThI/Ejq3zko4NxUVjmt2ANSC/azGa0gUsfL/zqZkkMxtlmbpDkQSrcUThuLCvh4ttTm5r6IYsbIxCld2RtFz0PNTrQUe3bRFUnoerwVEkaEGU59Hl7Byv4xgXU8G3jcmHBKAaMTp5krezeVVW4o9urCfixsmImIL5h1+3F/0NFuqLnCEmywKGWlm5sXRwI1ng08hSwR9hIs4QETHzGcTbopkXZDsbdownCsGTsF88+Yh89DNUgC+KK8Bh9+6/aCLO0qBRtW3aqL0QlByM7zi2fXeKz6gwlSqIQHw4aD9EdAEtk/QJJMGl47cbN3iZX+OUYZckCYitNqWzyxHiAXbCjX3ADjn0exOiG40S5bwMezgceNCcdOClkIoRQeBLZsQOjJ+0GSyazHGVVxuWlVY8WuxwyF4oCQC1Y+D685IRgVDIjGbvnFs4HHrQnHDtxQXhkZg9H3Kt48DaGl9zHvM1rKZE3UyknhlnF+K405M7ghU/Y6uQZOL5WkzGRn+fDxEwzh2cADeGfCccsgkkcugc/uUiZrUrhvdxAnjYmofuy5ZIdCGusMXipJmc3OvJzFeh3Pigu8RL6FELkoAO22eGFNCgNJolIR5SKTF9JY5/CSpZGXsrNSx9MZj1dwWgihlyHkEvjsLmUanRRyyQ7lQW1W92Ys3tOG+mgP9m2pQeCm+a5nxU6VFJ1AL1P0UknKS9lZqSMCjws4KYQwUkbLNfDZWco0Oilwg2FvN3B7K7NdtRS8ZnVvxtJdyxFOpqxt6vp7QBklRSeFHmZKivmm2BRrQjDgHUSpzQWcNIg0UkZzyxlZC1bJptxHVZOCVjAc1vcls/wmBa/Fe9oyQUdC+V447XhhtKRYCBjZROmlMqZXDGcFIuNxDaeEEEbKaIWgAGSVbG4Zd1w1KbCywyTUV0jy8pt0pVsf1fdTc1ro4aV1BiNj9ZqcWggGvIEIPB7HaBmtEBSAykmhvb1d9RhpjN0rl6I+2oPOUA1O1gko0mvu21KDun71Y+XvhdNCDy+tMxgZq5cCqcA7iFKbxymEMprdxJun4eKv/xLBi3+FxgsewRcGDDpbG8OouWm+7nvhtMWO0ZJiIWBEsVYqHWIF7iICj4J8NkXLBbcajLmNfFJkGXSygquR98JIoLbiTcZaZ7jntFhBln+MrIl4SU4t8A6i1CajEBwGcqEQymh2I18TerG2BdUhHxZ/+iJT1SZH773QW++yQ+llpKRYKOitiXhJTi3wDiLwyMi3w4Agm+xJcQaAGbBDS6UVnLy2mO4GYsFeYDei1CYj3w4DgvwjFtMFAucRgUeGG71dBIWNWEwXCJxHBB4ZxagQE5hDLKYLBM4j1nhkFMJGSy9RDI7QSvuc62fMA1qaPX9eAkEhIwKPAqsKsUJu+GYnxeDzxVMxXj/3h2idWfifWTEEfkFpIkptNuK0D1ghYcTnq9DJpV1EoeAlDzWBQIkIPDbi5YnMLMWg/vKyirEYAr+gdBGBx0a8PJGZpRjUX15WMRZD4BeULiLw2IiXJzKzFIP6y8sqxmII/ILSRQQeG/HyRGaWYuh94mWfu2II/ILSRajabKTU5NjFYKXiVZ874aEm8DKWAg8hpBXA/wZwJoBJlNJtsvt+BOBmAAkAt1NK/8fKsbyCVycygfcohsAvKE2sZjw7APwDgKXyGwkhZwGYBeArAMYCeJMQcjqlVKx8CgQCQYljaY2HUrqTUrqLcdd0AKsopVFK6acAdgOYZOVYAoFAICgOCKVU/1F6L0LIuwB+KJXaCCG/BLCVUvp8+u+nAKynlL4kf15fX1/m4IXcs8TrjNyxFWPfeRllfQcxWDUK+y65FofOnpzvYQkcZv0BPx7/vAzdUYLaEMVt4wdxxWhRdBC4Q1NTU+bfVVVVWUoY3VIbIeRNAGMYdy2klK7hPY1xm2aEkw9Si/b2dsOP9Sp2nmNgywaEXn8+s7E12HcQ419/HmNqx+R9LaoUPksgP+e5uiOC+zuGLI32Rwnu7yjHSWOcUR6Kz7K4cPo8dQMPpfTrObxuF4B62d91APbl8DoCi4jmdqWJaGgnKGSc2sfzKoBZhJAQIeRUAE0APnToWAINSslNQTCEcDYQFDKWAg8h5FpCSBeACwCsI4T8DwBQSv8EoA3AnwG8AeD7QtGWH0rJTUEwhHA2EBQyVlVtL1NK6yilIUppLaX072X3LaaUNlJKz6CUrrc+VEEulJKbgmAI4WwgKGSEc0GRU2puCoIUwtlAUMiIwFMCCDeF0kQ4GwgKFWESKhAIBAJXEYFHIBAIBK4iAo9AIBAIXEUEHoFAIBC4ii1ebbki92oTCAQCQXGi9GoTGY9AIBAIXEUEHoFAIBC4Sl5LbQKBQCAoPUTGIxAIBAJX8WTgIYT8jBDyR0LIdkLIbwkhY/M9JrshhPyCEPKX9Hm+TAgZke8xOQEhpJUQ8idCSJIQMjHf47ETQsjlhJBdhJDdhJC78z0eJyCEPE0IOUAI2ZHvsTgJIaSeEPIOIWRn+vu6IN9jshtCSDkh5ENCyB/S5/hTx47lxVIbIeQESumR9L9vB3AWpfR7eR6WrRBCvgHgbUppnBDyAABQSu/K87BshxByJoAkgKWQdbH1OoQQP4BPAExDqj/V7wBcTyn9c14HZjOEkKkAjgF4llJ6dr7H4xSEkJMAnEQp/T0hZDiAjwBcU0yfJyGEAAhTSo8RQsoAbAKwgFK61e5jeTLjkYJOmjB0upt6EUrpbyml8fSfW5Fqpld0UEp3Ukp35XscDjAJwG5K6R5KaQzAKgDT8zwm26GUbgRwMN/jcBpK6V8ppb9P//sogJ0AxuV3VPZCUxxL/1mW/s+RudWTgQcACCGLCSGdAG4AsCjf43GY7wIQrSW8xTgAnbK/u1BkE1WpQgg5BcDXAHyQ35HYDyHETwjZDuAAgA2UUkfOsWADDyHkTULIDsZ/0wGAUrqQUloP4FcA/jm/o80NvXNMP2YhgDhS5+lJjJxnEUIYtxVdZl5qEEKGAfgNgH9RVF6KAkppglJ6HlIVlkmEEEfKpwXbFoFS+nWDD/01gHUAfuLgcBxB7xwJIbMBXAXgMurFxbg0Jj7LYqILQL3s7zoA+/I0FoENpNc9fgPgV5TS/873eJyEUnqYEPIugMsB2C4cKdiMRwtCSJPsz28B+Eu+xuIUhJDLAdwF4FuU0uP5Ho/ANL8D0EQIOZUQEgQwC8CreR6TIEfSC+9PAdhJKX0o3+NxAkLIiZJ6lhBSAeDrcGhu9aqq7TcAzkBKDfU5gO9RSvfmd1T2QgjZDSAEoDd909ZiU+4BACHkWgCPAjgRwGEA2+Ut1L0MIeSbAP4TgB/A05TSxXkeku0QQl4AcDGAGgDdAH5CKX0qr4NyAELIFADvAfgYqXkHAO6hlL6ev1HZCyHkqwCeQer76gPQRim915FjeTHwCAQCgcC7eLLUJhAIBALvIgKPQCAQCFxFBB6BQCAQuIoIPAKBQCBwFRF4BAKBQOAqIvAIBAKBwFVE4BEIBAKBq4jAIxAIBAJX+f9c8zLYtQ9VJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 250\n",
    "m, ep = np.array([1, 3]), 3\n",
    "\n",
    "x = stats.norm(0, 1).rvs(size=(n, 2))\n",
    "y = x @ m + stats.norm(0, ep).rvs(size=n)\n",
    "\n",
    "print(f'Rsq = {r2_score(y, x @ m):0.3f}')\n",
    "print(f'MSE = {mean_squared_error(y, x @ m):0.3f}')\n",
    "print(f'Mean loglikelihood = {stats.norm(x @ m, ep).logpdf(y).mean()}')\n",
    "\n",
    "graph.plot(x, y, 'o')\n",
    "graph.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.408343593549797 4.196660909092196\n",
      "-2.853227720651972\n"
     ]
    }
   ],
   "source": [
    "print(y.mean(), y.std())\n",
    "print(stats.norm(y.mean(), y.std()).logpdf(y).mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statsmodels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>          <td>y</td>        <th>  No. Observations:  </th>      <td>   250</td>      \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>GLM</td>       <th>  Df Residuals:      </th>      <td>   248</td>      \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>       <td>Gaussian</td>     <th>  Df Model:          </th>      <td>     1</td>      \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>      <td>identity</td>     <th>  Scale:             </th> <td>7.810402211331509</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>IRLS</td>       <th>  Log-Likelihood:    </th>     <td> -610.66</td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Tue, 24 Jul 2018</td> <th>  Deviance:          </th>     <td>  1937.0</td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>15:02:50</td>     <th>  Pearson chi2:      </th>     <td>1.94e+03</td>     \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>         <td>2</td>        <th>                     </th>         <td> </td>        \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "   <td></td>     <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th> <td>    0.8828</td> <td>    0.175</td> <td>    5.057</td> <td> 0.000</td> <td>    0.541</td> <td>    1.225</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x2</th> <td>    3.1216</td> <td>    0.177</td> <td>   17.650</td> <td> 0.000</td> <td>    2.775</td> <td>    3.468</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   No. Observations:                  250\n",
       "Model:                            GLM   Df Residuals:                      248\n",
       "Model Family:                Gaussian   Df Model:                            1\n",
       "Link Function:               identity   Scale:               7.810402211331509\n",
       "Method:                          IRLS   Log-Likelihood:                -610.66\n",
       "Date:                Tue, 24 Jul 2018   Deviance:                       1937.0\n",
       "Time:                        15:02:50   Pearson chi2:                 1.94e+03\n",
       "No. Iterations:                     2                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "x1             0.8828      0.175      5.057      0.000       0.541       1.225\n",
       "x2             3.1216      0.177     17.650      0.000       2.775       3.468\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "glm = sm.GLM(y, x).fit()\n",
    "display(glm.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Keras OLS Regression\n",
    "\n",
    "Learning the conditional mean given data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_layer = Input(shape=(x.shape[1],))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_1 (InputLayer)         (None, 2)                 0         \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 1)                 3         \n",
      "=================================================================\n",
      "Total params: 3\n",
      "Trainable params: 3\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n",
      "None\n",
      "Epoch 1/100\n",
      "250/250 [==============================] - 0s 283us/step - loss: 16.5699 - keras_r2: 0.0285\n",
      "Epoch 2/100\n",
      "250/250 [==============================] - 0s 63us/step - loss: 14.1769 - keras_r2: 0.1661\n",
      "Epoch 3/100\n",
      "250/250 [==============================] - 0s 46us/step - loss: 12.4672 - keras_r2: 0.2681\n",
      "Epoch 4/100\n",
      "250/250 [==============================] - 0s 51us/step - loss: 11.2100 - keras_r2: 0.3379\n",
      "Epoch 5/100\n",
      "250/250 [==============================] - 0s 38us/step - loss: 10.2817 - keras_r2: 0.3961\n",
      "Epoch 6/100\n",
      "250/250 [==============================] - 0s 50us/step - loss: 9.5937 - keras_r2: 0.4289\n",
      "Epoch 7/100\n",
      "250/250 [==============================] - 0s 37us/step - loss: 9.0927 - keras_r2: 0.4746\n",
      "Epoch 8/100\n",
      "250/250 [==============================] - 0s 39us/step - loss: 8.7414 - keras_r2: 0.4702\n",
      "Epoch 9/100\n",
      "250/250 [==============================] - 0s 45us/step - loss: 8.4720 - keras_r2: 0.5026\n",
      "Epoch 10/100\n",
      "250/250 [==============================] - 0s 44us/step - loss: 8.2711 - keras_r2: 0.5052\n",
      "Epoch 11/100\n",
      "250/250 [==============================] - 0s 44us/step - loss: 8.1312 - keras_r2: 0.5217\n",
      "Epoch 12/100\n",
      "250/250 [==============================] - 0s 41us/step - loss: 8.0246 - keras_r2: 0.5169\n",
      "Epoch 13/100\n",
      "250/250 [==============================] - 0s 32us/step - loss: 7.9453 - keras_r2: 0.5448\n",
      "Epoch 14/100\n",
      "250/250 [==============================] - 0s 37us/step - loss: 7.8857 - keras_r2: 0.5260\n",
      "Epoch 15/100\n",
      "250/250 [==============================] - 0s 44us/step - loss: 7.8456 - keras_r2: 0.5255\n",
      "Epoch 16/100\n",
      "250/250 [==============================] - 0s 39us/step - loss: 7.8226 - keras_r2: 0.5261\n",
      "Epoch 17/100\n",
      "250/250 [==============================] - 0s 49us/step - loss: 7.7850 - keras_r2: 0.5437\n",
      "Epoch 18/100\n",
      "250/250 [==============================] - 0s 45us/step - loss: 7.7796 - keras_r2: 0.5272\n",
      "Epoch 19/100\n",
      "250/250 [==============================] - 0s 33us/step - loss: 7.7547 - keras_r2: 0.5445\n",
      "Epoch 20/100\n",
      "250/250 [==============================] - 0s 40us/step - loss: 7.7488 - keras_r2: 0.5245\n",
      "Epoch 21/100\n",
      "250/250 [==============================] - 0s 56us/step - loss: 7.7441 - keras_r2: 0.5287\n",
      "Epoch 22/100\n",
      "250/250 [==============================] - 0s 37us/step - loss: 7.7416 - keras_r2: 0.5287\n",
      "Epoch 23/100\n",
      "250/250 [==============================] - 0s 40us/step - loss: 7.7313 - keras_r2: 0.5209\n",
      "Epoch 24/100\n",
      "250/250 [==============================] - 0s 45us/step - loss: 7.7355 - keras_r2: 0.5392\n",
      "Epoch 25/100\n",
      "250/250 [==============================] - 0s 47us/step - loss: 7.7346 - keras_r2: 0.5247\n",
      "Epoch 26/100\n",
      "250/250 [==============================] - 0s 37us/step - loss: 7.7283 - keras_r2: 0.5297\n",
      "Epoch 27/100\n",
      "250/250 [==============================] - 0s 42us/step - loss: 7.7262 - keras_r2: 0.5418\n",
      "Epoch 28/100\n",
      "250/250 [==============================] - 0s 29us/step - loss: 7.7281 - keras_r2: 0.5371\n",
      "Epoch 29/100\n",
      "250/250 [==============================] - 0s 49us/step - loss: 7.7195 - keras_r2: 0.5454\n",
      "Epoch 30/100\n",
      "250/250 [==============================] - 0s 55us/step - loss: 7.7217 - keras_r2: 0.5317\n",
      "Epoch 31/100\n",
      "250/250 [==============================] - 0s 43us/step - loss: 7.7227 - keras_r2: 0.5342\n",
      "Epoch 32/100\n",
      "250/250 [==============================] - 0s 24us/step - loss: 7.7180 - keras_r2: 0.5438\n",
      "Epoch 33/100\n",
      "250/250 [==============================] - 0s 37us/step - loss: 7.7210 - keras_r2: 0.5403\n",
      "Epoch 34/100\n",
      "250/250 [==============================] - 0s 41us/step - loss: 7.7212 - keras_r2: 0.5432\n",
      "Epoch 35/100\n",
      "250/250 [==============================] - 0s 39us/step - loss: 7.7232 - keras_r2: 0.5368\n",
      "Epoch 36/100\n",
      "250/250 [==============================] - 0s 41us/step - loss: 7.7266 - keras_r2: 0.5275\n",
      "Epoch 37/100\n",
      "250/250 [==============================] - 0s 57us/step - loss: 7.7257 - keras_r2: 0.5546\n",
      "[array([[0.8484873],\n",
      "       [3.0998402]], dtype=float32), array([0.19270161], dtype=float32)]\n"
     ]
    }
   ],
   "source": [
    "# Keras OLS\n",
    "ols_layer = Dense(1, kernel_initializer='zeros')(input_layer)\n",
    "ols_model = Model(input_layer, ols_layer)\n",
    "\n",
    "print(ols_model.summary())\n",
    "\n",
    "# Fit OLS model\n",
    "ols_model.compile(loss='mse', metrics=[keras_r2], optimizer=kopt.SGD(0.01))\n",
    "ols_model.fit(\n",
    "    x, y, \n",
    "    epochs=100, verbose=1, \n",
    "    callbacks=[keras.callbacks.EarlyStopping(monitor='loss', patience=5)]\n",
    ")\n",
    "\n",
    "print(ols_model.layers[1].get_weights())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Keras Distribution as Loss Function\n",
    "\n",
    "Here I fit a linear model where the loss function is a the mean negative log likelihood of observing `y_true` given proposed model parameters.\n",
    "\n",
    "I am using the **mean** rather than the **sum** of the negative log likelihood because it is has identical optima while being invariant to batch size.\n",
    "\n",
    "Finding the optimal value for likelihood functions are famously hard so I use stochastic gradient descent with a learning rate drop off. Needless to say, it takes a lot longer to converge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normal_loglikelihood(y_true, parameters):\n",
    "    m, s = parameters[:, 0], parameters[:, 1]\n",
    "    m, s = K.expand_dims(m), K.expand_dims(s)\n",
    "    \n",
    "    norm = tfdist.Normal(m, s)\n",
    "    return K.mean(-norm.log_prob(y_true))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "__________________________________________________________________________________________________\n",
      "Layer (type)                    Output Shape         Param #     Connected to                     \n",
      "==================================================================================================\n",
      "input_1 (InputLayer)            (None, 2)            0                                            \n",
      "__________________________________________________________________________________________________\n",
      "dense_2 (Dense)                 (None, 1)            3           input_1[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "dense_3 (Dense)                 (None, 1)            3           input_1[0][0]                    \n",
      "__________________________________________________________________________________________________\n",
      "concatenate_1 (Concatenate)     (None, 2)            0           dense_2[0][0]                    \n",
      "                                                                 dense_3[0][0]                    \n",
      "==================================================================================================\n",
      "Total params: 6\n",
      "Trainable params: 6\n",
      "Non-trainable params: 0\n",
      "__________________________________________________________________________________________________\n",
      "None\n",
      "Tensor(\"loss_1/concatenate_1_loss/ExpandDims:0\", shape=(?, 1), dtype=float32) Tensor(\"loss_1/concatenate_1_loss/ExpandDims_1:0\", shape=(?, 1), dtype=float32)\n",
      "Epoch 1/5000\n",
      "250/250 [==============================] - 0s 379us/step - loss: 10.2526\n",
      "Epoch 2/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 5.4230\n",
      "Epoch 3/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 4.4141\n",
      "Epoch 4/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 3.9416\n",
      "Epoch 5/5000\n",
      "250/250 [==============================] - 0s 54us/step - loss: 3.6498\n",
      "Epoch 6/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 3.4604\n",
      "Epoch 7/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 3.3254\n",
      "Epoch 8/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 3.2217\n",
      "Epoch 9/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 3.1414\n",
      "Epoch 10/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 3.0759\n",
      "Epoch 11/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 3.0221\n",
      "Epoch 12/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.9774\n",
      "Epoch 13/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.9382\n",
      "Epoch 14/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.9050\n",
      "Epoch 15/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.8759\n",
      "Epoch 16/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 2.8506\n",
      "Epoch 17/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.8277\n",
      "Epoch 18/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.8068\n",
      "Epoch 19/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.7884\n",
      "Epoch 20/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.7718\n",
      "Epoch 21/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.7569\n",
      "Epoch 22/5000\n",
      "250/250 [==============================] - 0s 35us/step - loss: 2.7436\n",
      "Epoch 23/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.7307\n",
      "Epoch 24/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.7193\n",
      "Epoch 25/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.7086\n",
      "Epoch 26/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.6986\n",
      "Epoch 27/5000\n",
      "250/250 [==============================] - 0s 42us/step - loss: 2.6895\n",
      "Epoch 28/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.6808\n",
      "Epoch 29/5000\n",
      "250/250 [==============================] - 0s 33us/step - loss: 2.6728\n",
      "Epoch 30/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.6653\n",
      "Epoch 31/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.6581\n",
      "Epoch 32/5000\n",
      "250/250 [==============================] - 0s 38us/step - loss: 2.6513\n",
      "Epoch 33/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.6451\n",
      "Epoch 34/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.6391\n",
      "Epoch 35/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.6333\n",
      "Epoch 36/5000\n",
      "250/250 [==============================] - 0s 36us/step - loss: 2.6279\n",
      "Epoch 37/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.6228\n",
      "Epoch 38/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.6180\n",
      "Epoch 39/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.6133\n",
      "Epoch 40/5000\n",
      "250/250 [==============================] - 0s 42us/step - loss: 2.6090\n",
      "Epoch 41/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.6047\n",
      "Epoch 42/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.6006\n",
      "Epoch 43/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.5968\n",
      "Epoch 44/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.5930\n",
      "Epoch 45/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.5893\n",
      "Epoch 46/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.5858\n",
      "Epoch 47/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.5824\n",
      "Epoch 48/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.5793\n",
      "Epoch 49/5000\n",
      "250/250 [==============================] - 0s 38us/step - loss: 2.5760\n",
      "Epoch 50/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.5730\n",
      "Epoch 51/5000\n",
      "250/250 [==============================] - 0s 33us/step - loss: 2.5701\n",
      "Epoch 52/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.5671\n",
      "Epoch 53/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.5644\n",
      "Epoch 54/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.5616\n",
      "Epoch 55/5000\n",
      "250/250 [==============================] - 0s 34us/step - loss: 2.5590\n",
      "Epoch 56/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.5566\n",
      "Epoch 57/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.5542\n",
      "Epoch 58/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.5517\n",
      "Epoch 59/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.5495\n",
      "Epoch 60/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.5471\n",
      "Epoch 61/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.5449\n",
      "Epoch 62/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.5428\n",
      "Epoch 63/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.5407\n",
      "Epoch 64/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.5387\n",
      "Epoch 65/5000\n",
      "250/250 [==============================] - 0s 28us/step - loss: 2.5367\n",
      "Epoch 66/5000\n",
      "250/250 [==============================] - 0s 67us/step - loss: 2.5347\n",
      "Epoch 67/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.5329\n",
      "Epoch 68/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.5310\n",
      "Epoch 69/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.5292\n",
      "Epoch 70/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.5274\n",
      "Epoch 71/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.5256\n",
      "Epoch 72/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.5239\n",
      "Epoch 73/5000\n",
      "250/250 [==============================] - 0s 33us/step - loss: 2.5222\n",
      "Epoch 74/5000\n",
      "250/250 [==============================] - 0s 65us/step - loss: 2.5206\n",
      "Epoch 75/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.5191\n",
      "Epoch 76/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.5175\n",
      "Epoch 77/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.5160\n",
      "Epoch 78/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.5144\n",
      "Epoch 79/5000\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "250/250 [==============================] - 0s 48us/step - loss: 2.5130\n",
      "Epoch 80/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.5115\n",
      "Epoch 81/5000\n",
      "250/250 [==============================] - 0s 39us/step - loss: 2.5101\n",
      "Epoch 82/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.5088\n",
      "Epoch 83/5000\n",
      "250/250 [==============================] - 0s 54us/step - loss: 2.5073\n",
      "Epoch 84/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 2.5060\n",
      "Epoch 85/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.5047\n",
      "Epoch 86/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.5034\n",
      "Epoch 87/5000\n",
      "250/250 [==============================] - 0s 65us/step - loss: 2.5021\n",
      "Epoch 88/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 2.5008\n",
      "Epoch 89/5000\n",
      "250/250 [==============================] - 0s 69us/step - loss: 2.4997\n",
      "Epoch 90/5000\n",
      "250/250 [==============================] - 0s 69us/step - loss: 2.4984\n",
      "Epoch 91/5000\n",
      "250/250 [==============================] - 0s 73us/step - loss: 2.4972\n",
      "Epoch 92/5000\n",
      "250/250 [==============================] - 0s 70us/step - loss: 2.4961\n",
      "Epoch 93/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4949\n",
      "Epoch 94/5000\n",
      "250/250 [==============================] - 0s 54us/step - loss: 2.4938\n",
      "Epoch 95/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4927\n",
      "Epoch 96/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4916\n",
      "Epoch 97/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.4905\n",
      "Epoch 98/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4894\n",
      "Epoch 99/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4884\n",
      "Epoch 100/5000\n",
      "250/250 [==============================] - 0s 65us/step - loss: 2.4874\n",
      "Epoch 101/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4863\n",
      "Epoch 102/5000\n",
      "250/250 [==============================] - 0s 94us/step - loss: 2.4854\n",
      "Epoch 103/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4843\n",
      "Epoch 104/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4834\n",
      "Epoch 105/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4824\n",
      "Epoch 106/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4815\n",
      "Epoch 107/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4805\n",
      "Epoch 108/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4797\n",
      "Epoch 109/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.4788\n",
      "Epoch 110/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4779\n",
      "Epoch 111/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4770\n",
      "Epoch 112/5000\n",
      "250/250 [==============================] - 0s 69us/step - loss: 2.4762\n",
      "Epoch 113/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.4753\n",
      "Epoch 114/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4744\n",
      "Epoch 115/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.4736\n",
      "Epoch 116/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.4729\n",
      "Epoch 117/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4720\n",
      "Epoch 118/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4713\n",
      "Epoch 119/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4705\n",
      "Epoch 120/5000\n",
      "250/250 [==============================] - 0s 36us/step - loss: 2.4697\n",
      "Epoch 121/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4690\n",
      "Epoch 122/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4682\n",
      "Epoch 123/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.4674\n",
      "Epoch 124/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.4668\n",
      "Epoch 125/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4660\n",
      "Epoch 126/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4653\n",
      "Epoch 127/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.4646\n",
      "Epoch 128/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.4639\n",
      "Epoch 129/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4632\n",
      "Epoch 130/5000\n",
      "250/250 [==============================] - 0s 30us/step - loss: 2.4626\n",
      "Epoch 131/5000\n",
      "250/250 [==============================] - 0s 41us/step - loss: 2.4620\n",
      "Epoch 132/5000\n",
      "250/250 [==============================] - 0s 39us/step - loss: 2.4613\n",
      "Epoch 133/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4607\n",
      "Epoch 134/5000\n",
      "250/250 [==============================] - 0s 31us/step - loss: 2.4601\n",
      "Epoch 135/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4594\n",
      "Epoch 136/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4589\n",
      "Epoch 137/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4583\n",
      "Epoch 138/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4577\n",
      "Epoch 139/5000\n",
      "250/250 [==============================] - 0s 84us/step - loss: 2.4571\n",
      "Epoch 140/5000\n",
      "250/250 [==============================] - 0s 28us/step - loss: 2.4566\n",
      "Epoch 141/5000\n",
      "250/250 [==============================] - 0s 32us/step - loss: 2.4559\n",
      "Epoch 142/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4554\n",
      "Epoch 143/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 2.4548\n",
      "Epoch 144/5000\n",
      "250/250 [==============================] - 0s 38us/step - loss: 2.4543\n",
      "Epoch 145/5000\n",
      "250/250 [==============================] - 0s 39us/step - loss: 2.4538\n",
      "Epoch 146/5000\n",
      "250/250 [==============================] - 0s 36us/step - loss: 2.4532\n",
      "Epoch 147/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.4527\n",
      "Epoch 148/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.4522\n",
      "Epoch 149/5000\n",
      "250/250 [==============================] - 0s 54us/step - loss: 2.4518\n",
      "Epoch 150/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4512\n",
      "Epoch 151/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4508\n",
      "Epoch 152/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4503\n",
      "Epoch 153/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4498\n",
      "Epoch 154/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4494\n",
      "Epoch 155/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4489\n",
      "Epoch 156/5000\n",
      "250/250 [==============================] - 0s 35us/step - loss: 2.4485\n",
      "Epoch 157/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.4480\n",
      "Epoch 158/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4476\n",
      "Epoch 159/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4472\n",
      "Epoch 160/5000\n",
      "250/250 [==============================] - 0s 42us/step - loss: 2.4468\n",
      "Epoch 161/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.4463\n",
      "Epoch 162/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4460\n",
      "Epoch 163/5000\n",
      "250/250 [==============================] - 0s 42us/step - loss: 2.4456\n",
      "Epoch 164/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.4451\n",
      "Epoch 165/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4448\n",
      "Epoch 166/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4444\n",
      "Epoch 167/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4440\n",
      "Epoch 168/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4437\n",
      "Epoch 169/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4432\n",
      "Epoch 170/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4429\n",
      "Epoch 171/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4425\n",
      "Epoch 172/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.4422\n",
      "Epoch 173/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.4419\n",
      "Epoch 174/5000\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "250/250 [==============================] - 0s 50us/step - loss: 2.4415\n",
      "Epoch 175/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4412\n",
      "Epoch 176/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4409\n",
      "Epoch 177/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4406\n",
      "Epoch 178/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.4403\n",
      "Epoch 179/5000\n",
      "250/250 [==============================] - 0s 66us/step - loss: 2.4400\n",
      "Epoch 180/5000\n",
      "250/250 [==============================] - 0s 67us/step - loss: 2.4396\n",
      "Epoch 181/5000\n",
      "250/250 [==============================] - 0s 67us/step - loss: 2.4393\n",
      "Epoch 182/5000\n",
      "250/250 [==============================] - 0s 66us/step - loss: 2.4390\n",
      "Epoch 183/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.4387\n",
      "Epoch 184/5000\n",
      "250/250 [==============================] - 0s 74us/step - loss: 2.4385\n",
      "Epoch 185/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4381\n",
      "Epoch 186/5000\n",
      "250/250 [==============================] - 0s 67us/step - loss: 2.4379\n",
      "Epoch 187/5000\n",
      "250/250 [==============================] - 0s 72us/step - loss: 2.4376\n",
      "Epoch 188/5000\n",
      "250/250 [==============================] - 0s 68us/step - loss: 2.4374\n",
      "Epoch 189/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4371\n",
      "Epoch 190/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4369\n",
      "Epoch 191/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4366\n",
      "Epoch 192/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.4364\n",
      "Epoch 193/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4362\n",
      "Epoch 194/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4359\n",
      "Epoch 195/5000\n",
      "250/250 [==============================] - 0s 69us/step - loss: 2.4356\n",
      "Epoch 196/5000\n",
      "250/250 [==============================] - 0s 39us/step - loss: 2.4354\n",
      "Epoch 197/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4352\n",
      "Epoch 198/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4349\n",
      "Epoch 199/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4348\n",
      "Epoch 200/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4345\n",
      "Epoch 201/5000\n",
      "250/250 [==============================] - 0s 69us/step - loss: 2.4343\n",
      "Epoch 202/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4341\n",
      "Epoch 203/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4339\n",
      "Epoch 204/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.4337\n",
      "Epoch 205/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4335\n",
      "Epoch 206/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.4334\n",
      "Epoch 207/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4331\n",
      "Epoch 208/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4329\n",
      "Epoch 209/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4327\n",
      "Epoch 210/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4326\n",
      "Epoch 211/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4324\n",
      "Epoch 212/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4323\n",
      "Epoch 213/5000\n",
      "250/250 [==============================] - 0s 42us/step - loss: 2.4320\n",
      "Epoch 214/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4318\n",
      "Epoch 215/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4317\n",
      "Epoch 216/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4316\n",
      "Epoch 217/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4314\n",
      "Epoch 218/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4312\n",
      "Epoch 219/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4311\n",
      "Epoch 220/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4309\n",
      "Epoch 221/5000\n",
      "250/250 [==============================] - 0s 41us/step - loss: 2.4308\n",
      "Epoch 222/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4306\n",
      "Epoch 223/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4304\n",
      "Epoch 224/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4303\n",
      "Epoch 225/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4302\n",
      "Epoch 226/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4301\n",
      "Epoch 227/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4299\n",
      "Epoch 228/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4297\n",
      "Epoch 229/5000\n",
      "250/250 [==============================] - 0s 33us/step - loss: 2.4297\n",
      "Epoch 230/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4296\n",
      "Epoch 231/5000\n",
      "250/250 [==============================] - 0s 41us/step - loss: 2.4295\n",
      "Epoch 232/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4293\n",
      "Epoch 233/5000\n",
      "250/250 [==============================] - 0s 33us/step - loss: 2.4292\n",
      "Epoch 234/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4290\n",
      "Epoch 235/5000\n",
      "250/250 [==============================] - 0s 30us/step - loss: 2.4289\n",
      "Epoch 236/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.4289\n",
      "Epoch 237/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4288\n",
      "Epoch 238/5000\n",
      "250/250 [==============================] - 0s 39us/step - loss: 2.4286\n",
      "Epoch 239/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4286\n",
      "Epoch 240/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4284\n",
      "Epoch 241/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4283\n",
      "Epoch 242/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4282\n",
      "Epoch 243/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4281\n",
      "Epoch 244/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4280\n",
      "Epoch 245/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4279\n",
      "Epoch 246/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4278\n",
      "Epoch 247/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4277\n",
      "Epoch 248/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4277\n",
      "Epoch 249/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4275\n",
      "Epoch 250/5000\n",
      "250/250 [==============================] - 0s 42us/step - loss: 2.4274\n",
      "Epoch 251/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.4273\n",
      "Epoch 252/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4272\n",
      "Epoch 253/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4272\n",
      "Epoch 254/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4271\n",
      "Epoch 255/5000\n",
      "250/250 [==============================] - 0s 38us/step - loss: 2.4270\n",
      "Epoch 256/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4270\n",
      "Epoch 257/5000\n",
      "250/250 [==============================] - 0s 65us/step - loss: 2.4269\n",
      "Epoch 258/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.4268\n",
      "Epoch 259/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4267\n",
      "Epoch 260/5000\n",
      "250/250 [==============================] - 0s 65us/step - loss: 2.4267\n",
      "Epoch 261/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4266\n",
      "Epoch 262/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4265\n",
      "Epoch 263/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.4264\n",
      "Epoch 264/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4264\n",
      "Epoch 265/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4263\n",
      "Epoch 266/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4263\n",
      "Epoch 267/5000\n",
      "250/250 [==============================] - 0s 36us/step - loss: 2.4262\n",
      "Epoch 268/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4261\n",
      "Epoch 269/5000\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "250/250 [==============================] - 0s 37us/step - loss: 2.4260\n",
      "Epoch 270/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4260\n",
      "Epoch 271/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4260\n",
      "Epoch 272/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4258\n",
      "Epoch 273/5000\n",
      "250/250 [==============================] - 0s 68us/step - loss: 2.4258\n",
      "Epoch 274/5000\n",
      "250/250 [==============================] - 0s 67us/step - loss: 2.4258\n",
      "Epoch 275/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4257\n",
      "Epoch 276/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.4256\n",
      "Epoch 277/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4256\n",
      "Epoch 278/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4256\n",
      "Epoch 279/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4255\n",
      "Epoch 280/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4255\n",
      "Epoch 281/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4254\n",
      "Epoch 282/5000\n",
      "250/250 [==============================] - 0s 54us/step - loss: 2.4253\n",
      "Epoch 283/5000\n",
      "250/250 [==============================] - 0s 55us/step - loss: 2.4253\n",
      "Epoch 284/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4252\n",
      "Epoch 285/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4252\n",
      "Epoch 286/5000\n",
      "250/250 [==============================] - 0s 41us/step - loss: 2.4252\n",
      "Epoch 287/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4252\n",
      "Epoch 288/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4251\n",
      "Epoch 289/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4250\n",
      "Epoch 290/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4250\n",
      "Epoch 291/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.4250\n",
      "Epoch 292/5000\n",
      "250/250 [==============================] - 0s 41us/step - loss: 2.4249\n",
      "Epoch 293/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4249\n",
      "Epoch 294/5000\n",
      "250/250 [==============================] - 0s 39us/step - loss: 2.4248\n",
      "Epoch 295/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4248\n",
      "Epoch 296/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4248\n",
      "Epoch 297/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.4247\n",
      "Epoch 298/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4247\n",
      "Epoch 299/5000\n",
      "250/250 [==============================] - 0s 47us/step - loss: 2.4246\n",
      "Epoch 300/5000\n",
      "250/250 [==============================] - 0s 33us/step - loss: 2.4246\n",
      "Epoch 301/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.4246\n",
      "\n",
      "Epoch 00301: ReduceLROnPlateau reducing learning rate to 0.004999999888241291.\n",
      "Epoch 302/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4245\n",
      "Epoch 303/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4245\n",
      "Epoch 304/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4245\n",
      "Epoch 305/5000\n",
      "250/250 [==============================] - 0s 42us/step - loss: 2.4244\n",
      "Epoch 306/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4244\n",
      "\n",
      "Epoch 00306: ReduceLROnPlateau reducing learning rate to 0.0024999999441206455.\n",
      "Epoch 307/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4244\n",
      "Epoch 308/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4244\n",
      "Epoch 309/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4244\n",
      "Epoch 310/5000\n",
      "250/250 [==============================] - 0s 41us/step - loss: 2.4243\n",
      "Epoch 311/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00311: ReduceLROnPlateau reducing learning rate to 0.0012499999720603228.\n",
      "Epoch 312/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.4243\n",
      "Epoch 313/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4243\n",
      "Epoch 314/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00314: ReduceLROnPlateau reducing learning rate to 0.0006249999860301614.\n",
      "Epoch 315/5000\n",
      "250/250 [==============================] - 0s 58us/step - loss: 2.4243\n",
      "Epoch 316/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.4243\n",
      "Epoch 317/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00317: ReduceLROnPlateau reducing learning rate to 0.0003124999930150807.\n",
      "Epoch 318/5000\n",
      "250/250 [==============================] - 0s 54us/step - loss: 2.4243\n",
      "Epoch 319/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4243\n",
      "Epoch 320/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00320: ReduceLROnPlateau reducing learning rate to 0.00015624999650754035.\n",
      "Epoch 321/5000\n",
      "250/250 [==============================] - 0s 71us/step - loss: 2.4243\n",
      "Epoch 322/5000\n",
      "250/250 [==============================] - 0s 50us/step - loss: 2.4243\n",
      "Epoch 323/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00323: ReduceLROnPlateau reducing learning rate to 7.812499825377017e-05.\n",
      "Epoch 324/5000\n",
      "250/250 [==============================] - 0s 49us/step - loss: 2.4243\n",
      "Epoch 325/5000\n",
      "250/250 [==============================] - 0s 38us/step - loss: 2.4243\n",
      "Epoch 326/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00326: ReduceLROnPlateau reducing learning rate to 3.9062499126885086e-05.\n",
      "Epoch 327/5000\n",
      "250/250 [==============================] - 0s 40us/step - loss: 2.4243\n",
      "Epoch 328/5000\n",
      "250/250 [==============================] - 0s 37us/step - loss: 2.4243\n",
      "Epoch 329/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00329: ReduceLROnPlateau reducing learning rate to 1.9531249563442543e-05.\n",
      "Epoch 330/5000\n",
      "250/250 [==============================] - 0s 66us/step - loss: 2.4243\n",
      "Epoch 331/5000\n",
      "250/250 [==============================] - 0s 38us/step - loss: 2.4243\n",
      "Epoch 332/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00332: ReduceLROnPlateau reducing learning rate to 9.765624781721272e-06.\n",
      "Epoch 333/5000\n",
      "250/250 [==============================] - 0s 39us/step - loss: 2.4243\n",
      "Epoch 334/5000\n",
      "250/250 [==============================] - 0s 44us/step - loss: 2.4243\n",
      "Epoch 335/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00335: ReduceLROnPlateau reducing learning rate to 4.882812390860636e-06.\n",
      "Epoch 336/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.4243\n",
      "Epoch 337/5000\n",
      "250/250 [==============================] - 0s 62us/step - loss: 2.4243\n",
      "Epoch 338/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00338: ReduceLROnPlateau reducing learning rate to 2.441406195430318e-06.\n",
      "Epoch 339/5000\n",
      "250/250 [==============================] - 0s 59us/step - loss: 2.4243\n",
      "Epoch 340/5000\n",
      "250/250 [==============================] - 0s 46us/step - loss: 2.4243\n",
      "Epoch 341/5000\n",
      "250/250 [==============================] - 0s 45us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00341: ReduceLROnPlateau reducing learning rate to 1.220703097715159e-06.\n",
      "Epoch 342/5000\n",
      "250/250 [==============================] - 0s 54us/step - loss: 2.4243\n",
      "Epoch 343/5000\n",
      "250/250 [==============================] - 0s 63us/step - loss: 2.4243\n",
      "Epoch 344/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00344: ReduceLROnPlateau reducing learning rate to 6.103515488575795e-07.\n",
      "Epoch 345/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4243\n",
      "Epoch 346/5000\n",
      "250/250 [==============================] - 0s 48us/step - loss: 2.4243\n",
      "Epoch 347/5000\n",
      "250/250 [==============================] - 0s 52us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00347: ReduceLROnPlateau reducing learning rate to 3.0517577442878974e-07.\n",
      "Epoch 348/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4243\n",
      "Epoch 349/5000\n",
      "250/250 [==============================] - 0s 66us/step - loss: 2.4243\n",
      "Epoch 350/5000\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "250/250 [==============================] - 0s 53us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00350: ReduceLROnPlateau reducing learning rate to 1.5258788721439487e-07.\n",
      "Epoch 351/5000\n",
      "250/250 [==============================] - 0s 57us/step - loss: 2.4243\n",
      "Epoch 352/5000\n",
      "250/250 [==============================] - 0s 53us/step - loss: 2.4243\n",
      "Epoch 353/5000\n",
      "250/250 [==============================] - 0s 43us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00353: ReduceLROnPlateau reducing learning rate to 7.629394360719743e-08.\n",
      "Epoch 354/5000\n",
      "250/250 [==============================] - 0s 56us/step - loss: 2.4243\n",
      "Epoch 355/5000\n",
      "250/250 [==============================] - 0s 51us/step - loss: 2.4243\n",
      "Epoch 356/5000\n",
      "250/250 [==============================] - 0s 72us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00356: ReduceLROnPlateau reducing learning rate to 3.814697180359872e-08.\n",
      "Epoch 357/5000\n",
      "250/250 [==============================] - 0s 64us/step - loss: 2.4243\n",
      "Epoch 358/5000\n",
      "250/250 [==============================] - 0s 60us/step - loss: 2.4243\n",
      "Epoch 359/5000\n",
      "250/250 [==============================] - 0s 71us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00359: ReduceLROnPlateau reducing learning rate to 1.907348590179936e-08.\n",
      "Epoch 360/5000\n",
      "250/250 [==============================] - 0s 72us/step - loss: 2.4243\n",
      "Epoch 361/5000\n",
      "250/250 [==============================] - 0s 68us/step - loss: 2.4243\n",
      "Epoch 362/5000\n",
      "250/250 [==============================] - 0s 61us/step - loss: 2.4243\n",
      "\n",
      "Epoch 00362: ReduceLROnPlateau reducing learning rate to 9.53674295089968e-09.\n",
      "Epoch 363/5000\n",
      "250/250 [==============================] - 0s 70us/step - loss: 2.4243\n",
      "[array([[0.8283295],\n",
      "       [2.9777029]], dtype=float32), array([0.19800474], dtype=float32)]\n",
      "[array([[ 0.38376114],\n",
      "       [-0.12282857]], dtype=float32), array([2.6940897], dtype=float32)]\n"
     ]
    }
   ],
   "source": [
    "mu = Dense(1, kernel_initializer='zeros')(input_layer)\n",
    "sigma = Dense(1, kernel_initializer='zeros', activation='softplus')(input_layer)\n",
    "\n",
    "params = concatenate([mu, sigma])\n",
    "\n",
    "linear_model = Model(input_layer, params)\n",
    "\n",
    "print(linear_model.summary())\n",
    "\n",
    "# Fit sigma conditional model\n",
    "linear_model.compile(loss=normal_loglikelihood, optimizer=kopt.SGD(0.01))\n",
    "hist = linear_model.fit(\n",
    "    x, y, \n",
    "    epochs=5000, verbose=1,\n",
    "    callbacks=[\n",
    "        keras.callbacks.EarlyStopping(monitor='loss', patience=10),\n",
    "        keras.callbacks.ReduceLROnPlateau(monitor='loss', factor=0.5, patience=3, verbose=1),\n",
    "        keras.callbacks.ModelCheckpoint('temp', monitor='loss', save_best_only=True)\n",
    "    ]\n",
    ")\n",
    "\n",
    "# Check coefficents\n",
    "print(linear_model.layers[1].get_weights())\n",
    "print(linear_model.layers[2].get_weights())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparsion\n",
    "\n",
    "Check the results from a single point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 4.]]\n",
      "(1, 2)\n",
      "\n",
      "True value y     = [[13.]]\n",
      "OLS value E(y)   = [[13.440549]]\n",
      "y|x ~ N(mu, sig) = [[12.937145  2.659118]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAEGCAYAAADIRPqpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAGVVJREFUeJzt3X2QXGW94PHvb8ZLYhIIb1dvSKAIW5OVwLq5YJCVxJuVFVBURIq7BAII1qJXUCndNaIIlIBSgvheWGJQ2ABKyVVSFHcBXyJJKULgpoAQ3E6AC5NEstdIYkxCTPLbP+ZM7EzmLdM9c3q6v5+qqT79Oy/9nKdO55fnOU8/JzITSZLK1FZ2ASRJMhlJkkpnMpIklc5kJEkqnclIklS615VdgGobN250aJ8kNbmJEydGz5gtI0lS6UxGkqTSNVQ3nbpUKhU6OjrKLkbDGi318573vGeP9/fffz8A43vE/1zE62m01FGZrKP+jXT92DKSJJXOZCRJKp3ddJJUB5nJ5s2b2bVrV9lFqYuxY8eycePGIe3b1tbGhAkTiNhr0FyfTEaSVAebN29mzJgx7LfffmUXpS7GjBnD2LFjh7Tv9u3b2bx5M/vvv/+g97GbTpLqYNeuXU2TiGq133777XML0WQkSSqd3XSS1AQ2bNjA+973PgDWr19Pe3s7hxxyCAC/+MUvGr7VZjKShtnS378GwIHfXwPAX8osjJrWwQcfzNKlSwH40pe+xIQJE/jYxz62xzaZSWbS1tZ4nWImI0kaBj1/9Fyr+4f44+jnn3+e8847jxNPPJFly5Zx5513MmvWLF566SUA7r33XhYvXsw3v/lN1q9fzyc/+Uk6OzsBuPHGG5k5c+YexzvllFP42te+xvTp0wE4+eST+da3vsXRRx9dw9kN4p5RRNwWEesj4pmq2I8iYnnx92JELC/iR0bE1qp136na5/iIeDoiVkXEN2JfxvxJkobsueee4/zzz2fJkiUcdthhfW43f/58PvGJT7B48WK++93v8vGPf3yvbebNm8ddd921+7hAzYkIBtcy+gHwLeCO7kBm/vfu5Yj4ClA9GH11Zs7o5Ti3AJcAjwIPAKcB/7LvRZYk7YupU6dy3HHHDbjd4sWLqVQqQFeX3quvvsrWrVt5/etfv3ubs846i9mzZ3PNNdewcOFCzjvvvLqUccBklJmPRMSRva0rWjf/CLyjv2NExCTggMz8TfH+DuD9mIwkadiNGzdu93LP+0Xbtm3bvZyZuwc7bNu2rdffGY0fP55Zs2bxwAMPsGjRIpYsWVKXMtZ6z2g28EpmVqpiUyPiX4FNwJWZuQSYDHRWbdNZxKRRq3tAAsCrF3k5a09Dvccz3Nra2pg4cSKrV69m6tSp3H///Rx66KEAzJkzh1tvvZVLL70UgKeeeoo3v/nNex3jggsuYN68ecyePZuJEyfWpVy1JqO5wN1V79cBR2TmHyLieOCnEXEM0Nv9oX4fpNfdVGxVrX7+A2mM+vnr/zZ7K8/WrVuLpT3/J7pld7zvfeuhMeqosdWzjsaOHcuYMWPqdrxa7Nixg7/85S9s27aN1157jczcowX0uc99jg984ANMnjyZadOmsX37drZt28Z1113H/PnzWbhwITt27OCkk07ihhtu2Ov4xx57LGPHjuXss8/e47jVNm3axPr163e/H2gG8Mgc+OGqRTfd/Zl5bFXsdcAa4PjM7Oxjv8XA/yy2+2VmvqmIzwXmZOaHq7f3Sa9dnNq+f41SP9Uto55evWjy7tFU3UO7ufQHAPzi2x/cvd3svxvjIyRKUu862rhxY91aCY2gr246gM7OTs4880wee+yxPuef668+6v2k1/8GPFediCLibyOivVg+CugAns/MdcCfIuLE4j7TBcB9NXy2JKkECxcu5NRTT+Xzn//8Pk2EOpABu+ki4m5gDnBoRHQCV2fmAuAc9uyiA3g78IWI2AHsBD6SmRuKdf9E18i819M1cMHBC2paB35/DXS3iKQmMm/ePObNm1f34w5mNN3cPuIf7CV2L3BvH9svA47tbZ0kqbU13pwQkqSWYzKSSrTk969x4PfX9DsYQmoFJiNJUulMRpKk0pmMJEmDtmLFCqZNm8azzz5b1+OajCRJg3bzzTfz0EMPcfPNN9f1uCYjSWoia9asYe7cuRx33HHMmDGD+fPns337dgAmT+59DsWbbrqJE088kbe97W3MmjWLZcuW9Xn8BQsWcOSRR/K9732vruU2GUn96B7p5mg3jQaZyfnnn8/pp5/Ok08+yRNPPMGf//xnrr322j73eeyxx3jwwQf51a9+xa9//Wvuu+++PpPWcDIZSVKTeOSRRxgzZszuGRLa29v54he/yMKFC9myZUuv+/z+97/n4IMP3j3J6yGHHMKkSZP22GbFihWceuqpu98vX76c9773vXUtu48dl3ooqxXU83N9LIX21cqVK5kxY89nmx5wwAFMmTKF559/vtd93vGOd/DlL3+Z448/njlz5nDmmWcya9asPbY5+uijeeGFF9i5cyft7e1ceeWVXHfddXUtu8lIGqSRTlI+L2l0G1/M2l4vg5ndPTN7nby0rzjAhAkTdnfRLVmyhIsvvpirr76as846a/c2bW1tvOlNb2LlypWsXr2aKVOm7JX0amUykqQmcfTRR7No0aI9Yps2bWLNmjVMnTq1z/3a29uZPXs2s2fP5phjjuGuu+7aIxkBzJw5k9/+9rcsWLCAH//4x3Uvu/eMJKlJ/MM//ANbt27l7ru7Hqiwc+dOrrzySs4999w9Hj1erVKpsHr16t3vn376aY444oi9tnvLW97Cddddx+mnn85hhx1W97LbMpKkJhERLFy4kE996lPceOON7Nq1i3e+851cddVVAGzZsoXp06fv3v6jH/0os2bN4tOf/jQbN26kvb2do446iq9//et7HXvatGmMGTOGyy+/fFjKbjKSpGEwHE/wHYwpU6bwox/9qNd1f/zjH3uNP/TQQ3vFej5O/Dvf+Q5XXXUV48ePr72QvbCbTpLUpxdeeIGZM2eydetWzj333GH7HFtGkqQ+TZ06lccff3zYP8eWkSSpdCYjSVLpTEaSpNINeM8oIm4D3gOsz8xji9g1wP8A/l+x2Wcz84Fi3RXAh4CdwMcz88EifhrwdaAd+F5m3lDfU5GaV3+zPzg7g5rBYAYw/AD4FnBHj/hXM/Om6kBETAfOAY4BDgN+FhHTitXfBt4JdAKPR8SizKzv05mkIXBGbql8A3bTZeYjwIZBHu8M4IeZ+VpmvgCsAk4o/lZl5vOZuR34YbGtJDWFtra23c8NanXbt2+nrW3f7gLVMrT7soi4AFgGfCoz/whMBh6t2qaziAG83CP+1ho+W5IayoQJE9i8eTNbt24tuyh1sWnTJg444IAh7dvW1saECRP2aZ+hJqNbgGuBLF6/AlwM9DYtbNJ7Cyz7+4BKpTLEojWHVj//gdS3fnqfs2u06KsuvIYGZh31r+csDLXo6Ojod/2QklFmvtK9HBG3At3zXnQCh1dtOgVYWyz3Fe/VQAVvZpVKpaXPfyB1r5+lo/ueUW914TU0MOuofyNdP0Ma2h0R1Y8BPBN4plheBJwTEWMiYirQATwGPA50RMTUiNiPrkEOe85zLklqWYMZ2n03MAc4NCI6gauBORExg66utheBDwNk5oqIuAd4FtgBXJqZO4vjXAY8SNfQ7tsyc0Xdz0aSNCoNmIwyc24v4QX9bH89cH0v8QeAB/apdJKkluAMDJKk0pmMJEmlMxlJkkrn84zUkpwCSGostowkSaUzGUmSSmcykiSVzmQkSSqdAxikUa56MIYP2tNoZctIklQ6k5EkqXQmI0lS6UxGkqTSmYwkSaUzGUmSSufQbqmJ/HWY9zhYusah3ho1TEZqCU6MKjU2u+kkSaUzGUmSSmcykiSVbsBkFBG3RcT6iHimKnZjRDwXEU9FxE8i4sAifmREbI2I5cXfd6r2OT4ino6IVRHxjYiI4TklSdJoM5iW0Q+A03rEHgaOzcw3A/8XuKJq3erMnFH8faQqfgtwCdBR/PU8piSpRQ2YjDLzEWBDj9hDmbmjePsoMKW/Y0TEJOCAzPxNZiZwB/D+oRVZktRs6jG0+2LgR1Xvp0bEvwKbgCszcwkwGeis2qaziPWpUqnUoWijV6uf/0AGUz8zl44bgZI0Nq+j/lk//atn/XR0dPS7vqZkFBGfA3YAdxahdcARmfmHiDge+GlEHAP0dn8o+zv2QAVvZpVKpaXPfyCDrp+l/rbI66hvfs/6N9L1M+RkFBEXAu8BTi663sjM14DXiuUnImI1MI2ullB1V94UYO1QP1uS1FyGNLQ7Ik4D5gPvy8wtVfG/jYj2YvkougYqPJ+Z64A/RcSJxSi6C4D7ai69JKkpDNgyioi7gTnAoRHRCVxN1+i5McDDxQjtR4uRc28HvhARO4CdwEcys3vwwz/RNTLv9cC/FH+SJA2cjDJzbi/hBX1sey9wbx/rlgHH7lPpJNWkek4+J01VI3MGBklS6UxGkqTS+QgJNRUfFSGNTraMJEmlMxlJkkpnMpIklc5kJEkqnQMYpBbRc3CHvztSI7FlJEkqnclIklQ6k5EkqXQmI0lS6RzAoFHNGRek5mDLSJJUOpORJKl0JiNJUulMRpKk0jmAQWpRPgVWjcSWkSSpdCYjSVLpBpWMIuK2iFgfEc9UxQ6OiIcjolK8HlTEIyK+ERGrIuKpiDiuap8Li+0rEXFh/U9HkjQaDbZl9APgtB6xzwA/z8wO4OfFe4B3AR3F3yXALdCVvICrgbcCJwBXdycwSVJrG1QyysxHgA09wmcAtxfLtwPvr4rfkV0eBQ6MiEnAqcDDmbkhM/8IPMzeCU6S1IJqGU33xsxcB5CZ6yLiDUV8MvBy1XadRayveK8qlUoNRRv9Wv38B/LX+hlXajmaRateb6163oNVz/rp6Ojod/1wDO2OXmLZT7xXAxW8mVUqlZY+/4E4H139teL15vesfyNdP7WMpnul6H6jeF1fxDuBw6u2mwKs7ScuSWpxtSSjRUD3iLgLgfuq4hcUo+pOBDYW3XkPAqdExEHFwIVTipgkqcUNqpsuIu4G5gCHRkQnXaPibgDuiYgPAS8BZxebPwC8G1gFbAEuAsjMDRFxLfB4sd0XMrPnoAhJUgsaVDLKzLl9rDq5l20TuLSP49wG3Dbo0kmSWoIzMEiSSmcykiSVzmQkSSqdyUiSVDqTkSSpdCYjSVLpfNKrGp7T/0jNz5aRJKl0JiNJUulMRpKk0nnPSNJe9+VevajPR41Jw8KWkSSpdCYjSVLpTEaSpNJ5z0jSXqrvIXn/SCPBlpEkqXS2jNSQnHVBai22jCRJpTMZSZJKN+RkFBH/MSKWV/1tiojLI+KaiFhTFX931T5XRMSqiPhdRJxan1OQJI12Q75nlJm/A2YAREQ7sAb4CXAR8NXMvKl6+4iYDpwDHAMcBvwsIqZl5s6hlkGS1Bzq1U13MrA6M/+tn23OAH6Yma9l5gvAKuCEOn2+JGkUq1cyOge4u+r9ZRHxVETcFhEHFbHJwMtV23QWMUlSi6t5aHdE7Ae8D7iiCN0CXAtk8foV4GIgetk9+zpupVKptWijWqufP4wruwAqNPO12MznVg/1rJ+Ojo5+19fjd0bvAp7MzFcAul8BIuJW4P7ibSdweNV+U4C1fR10oII3s0ql0tLnD8BSf2fUKJr1WvR71r+Rrp96JKO5VHXRRcSkzFxXvD0TeKZYXgTcFRE30zWAoQN4rA6fL2kYOTWQRkJNySgixgHvBD5cFf5yRMygqwvuxe51mbkiIu4BngV2AJc6kk6SBDUmo8zcAhzSI3Z+P9tfD1xfy2dKkpqPMzBIkkrnRKlqCE6MKrU2W0aSpNKZjCRJpTMZSZJKZzKSJJXOAQySBq3nQBN/BKt6sWUkSSqdyUiSVDqTkSSpdN4zUmn8oaukbraMJEmlMxlJkkpnMpIklc5kJEkqnclIklQ6k5EkqXQO7ZY0ZNXD850aSLUwGUmqC+etUy1MRhpR/tBVUm9qvmcUES9GxNMRsTwilhWxgyPi4YioFK8HFfGIiG9ExKqIeCoijqv18yVJo1+9BjD818yckZlvKd5/Bvh5ZnYAPy/eA7wL6Cj+LgFuqdPnS5JGseHqpjsDmFMs3w4sBuYX8TsyM4FHI+LAiJiUmeuGqRwqmd1ykgajHi2jBB6KiCci4pIi9sbuBFO8vqGITwZertq3s4hJklpYPVpGJ2Xm2oh4A/BwRDzXz7bRSyx727BSqdShaKNX85z/uLILoJKMhmt4NJSxTPWsn46Ojn7X15yMMnNt8bo+In4CnAC80t39FhGTgPXF5p3A4VW7TwHWDqXgzaxSqTTP+S+1m65VNfo13FTfs2Ew0vVTUzddRIyPiP27l4FTgGeARcCFxWYXAvcVy4uAC4pRdScCG71fJEmqtWX0RuAnEdF9rLsy8/9ExOPAPRHxIeAl4Oxi+weAdwOrgC3ARTV+viSpCdSUjDLzeeA/9xL/A3ByL/EELq3lMyVJzceJUiVJpXM6IEnDwklUtS9sGUmSSmcykiSVzmQkSSqdyUiSVDoHMKjunBxV0r6yZSRJKp0tI0nDzmHeGojJSDWzW05SreymkySVzmQkSSqdyUiSVDqTkSSpdA5gkDSieg54cXSdwGSkIXIEnaR6sptOklQ6W0aSSuUPYgW2jCRJDcBkJEkq3ZCTUUQcHhG/jIiVEbEiIj5RxK+JiDURsbz4e3fVPldExKqI+F1EnFqPE5AkjX613DPaAXwqM5+MiP2BJyLi4WLdVzPzpuqNI2I6cA5wDHAY8LOImJaZO2sogySpCQw5GWXmOmBdsfyniFgJ9Hf38Qzgh5n5GvBCRKwCTgB+M9QyaGQ5nFvScKnLPaOIOBL4e+C3ReiyiHgqIm6LiIOK2GTg5ardOuk/eUmSWkTNQ7sjYgJwL3B5Zm6KiFuAa4EsXr8CXAxEL7tnX8etVCq1Fm1Ua8zzH1d2AdTkera+H5+1ZVg/rzG/Z42jnvXT0dHR7/qaklFE/A1diejOzPxngMx8pWr9rcD9xdtO4PCq3acAa/s69kAFb2aVSqUxz3+p3XQaWcP5PWjY71mDGOn6qWU0XQALgJWZeXNVfFLVZmcCzxTLi4BzImJMREwFOoDHhvr5kqTmUUvL6CTgfODpiFhexD4LzI2IGXR1wb0IfBggM1dExD3As3SNxLvUkXSSJKhtNN1Ser8P9EA/+1wPXD/Uz5QkNSfnplOfHMotaaSYjCQ1LCdRbR3OTSdJKp3JSJJUOrvptAfvE0kqgy0jSVLpbBlJGhV6ttod0NBcTEYtzm45SY3AbjpJUulMRpKk0tlNJ2lU8gexzcVk1IK8TySp0ZiMmpQJR63KUXejk8lI0qjnf75GP5ORpKbmvaXRwdF0kqTS2TJqInZVSP3b8zsyjlc7SiuKejAZjTImHEnNyGQkqWV5P6lxmIwkCRNT2UY8GUXEacDXgXbge5l5w0iXodHNXDoOltodJ5Wlv+5wE9XwGNHRdBHRDnwbeBcwHZgbEdNHsgySpMYTmTlyHxbxX4BrMvPU4v0VAJn5JYCNGzeOXGEkSaWYOHFi9IyN9O+MJgMvV73vLGKSpBY20slor2wI2BqSpBY30gMYOoHDq95PAdZ2v+mt6SZJan4j3TJ6HOiIiKkRsR9wDrBohMsgSWowI5qMMnMHcBnwILASuCczV4xkGRpZRLwYEU9HxPKIWFZ2eRpBRNwWEesj4pmq2MER8XBEVIrXg8osY9n6qKNrImJNcS0tj4h3l1nGMkXE4RHxy4hYGRErIuITRdzrqNBPHY3YdTSio+nUv4h4EXhLZv572WVpFBHxdmAzcEdmHlvEvgxsyMwbIuIzwEGZOb/Mcpapjzq6BticmTeVWbZGEBGTgEmZ+WRE7A88Abwf+CBeR0C/dfSPjNB15KzdamiZ+QiwoUf4DOD2Yvl2ur40LauPOlIhM9dl5pPF8p/o6pWZjNfRbv3U0YgxGTWWBB6KiCci4pKyC9PA3piZ66DrSwS8oeTyNKrLIuKpohuvZbugqkXEkcDfA7/F66hXPeoIRug6Mhk1lpMy8zi6Zqi4tOh+kYbiFuA/ADOAdcBXyi1O+SJiAnAvcHlmbiq7PI2olzoasevIZNRAMnNt8boe+AlwQrklalivFH3c3X3d60suT8PJzFcyc2dm7gJupcWvpYj4G7r+kb0zM/+5CHsdVemtjkbyOjIZNYiIGF/cOCQixgOnAM/0v1fLWgRcWCxfCNxXYlkaUvc/soUzaeFrKSICWACszMybq1Z5HRX6qqORvI4cTdcgIuIoulpD0PVj5Lsy8/oSi9QQIuJuYA5wKPAKcDXwU+Ae4AjgJeDszGzZG/h91NEcurpWEngR+HD3/ZFWExGzgCXA08CuIvxZuu6JeB3Rbx3NZYSuI5ORJKl0dtNJkkpnMpIklc5kJEkqnclIklQ6k5EkqXQmI0lS6UxGkqTSmYwkSaUzGUkjLCL+V0Tc2yP2zYj4WlllksrmDAzSCCvm+1oFTM7MVyPidcBa4F2Z+US5pZPKYctIGmHF3F6PAGcXodOAfzcRqZWZjKRy3A7MK5bnAf+7xLJIpbObTipBRIyl62Fls4FHgemZ+VK5pZLKYzKSShIRtwJvpauL7h1ll0cqk910UnluB/4TdtFJtoykskTEEcBzwN9l5qayyyOVyZaRVIKIaAM+CfzQRCR1Pd5a0giKiPF0PR783+ga1i21PLvpJEmls5tOklQ6k5EkqXQmI0lS6UxGkqTSmYwkSaUzGUmSSvf/AbKz7mK+0+rhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "query = np.asmatrix([1, 4], dtype=np.float)\n",
    "print(query)\n",
    "print(query.shape)\n",
    "print()\n",
    "\n",
    "print('True value y     =', query @ m)\n",
    "print('OLS value E(y)   =', ols_model.predict(query))\n",
    "print('y|x ~ N(mu, sig) =', linear_model.predict(query))\n",
    "\n",
    "lm_y_params = linear_model.predict(query)[0]\n",
    "\n",
    "graph.hist(stats.norm(*lm_y_params).rvs(50000), bins=100, label='')\n",
    "graph.axvline(query @ m, color='black', alpha=0.75, label='True y')\n",
    "graph.axvline(ols_model.predict(query)[0], color='red', alpha=0.75, label=r'OLS $\\hat{y}$')\n",
    "graph.xlabel('y')\n",
    "graph.legend()\n",
    "graph.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Stephen Anthony Rose_"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}