Linear_Mixed_Effects_Model_Variational_Inference_SSTest.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Copy of Linear_Mixed_Effects_Model_Variational_Inference.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/shihaoshen/TensorflowTest/blob/master/Linear_Mixed_Effects_Model_Variational_Inference_SSTest.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "u3Zq5VrfiDqB"
      },
      "cell_type": "markdown",
      "source": [
        "##### Copyright 2018 The TensorFlow Authors.\n",
        "\n",
        "Licensed under the Apache License, Version 2.0 (the \"License\");"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "3jTEqPzFiHQ0",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "#@title Licensed under the Apache License, Version 2.0 (the \"License\"); { display-mode: \"form\" }\n",
        "# you may not use this file except in compliance with the License.\n",
        "# You may obtain a copy of the License at\n",
        "#\n",
        "# https://www.apache.org/licenses/LICENSE-2.0\n",
        "#\n",
        "# Unless required by applicable law or agreed to in writing, software\n",
        "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
        "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
        "# See the License for the specific language governing permissions and\n",
        "# limitations under the License."
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "x97n3SaNmNpB"
      },
      "cell_type": "markdown",
      "source": [
        "# Fitting Generalized Linear Mixed-effects Models Using Variational Inference\n",
        "\n",
        "<table class=\"tfo-notebook-buttons\" align=\"left\">\n",
        "  <td>\n",
        "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/tensorflow/probability/blob/master/tensorflow_probability/examples/jupyter_notebooks/Linear_Mixed_Effects_Model_Variational_Inference.ipynb\"><img src=\"https://www.tensorflow.org/images/colab_logo_32px.png\" />Run in Google Colab</a>\n",
        "  </td>\n",
        "  <td>\n",
        "    <a target=\"_blank\" href=\"https://github.com/tensorflow/probability/blob/master/tensorflow_probability/examples/jupyter_notebooks/Linear_Mixed_Effects_Model_Variational_Inference.ipynb\"><img src=\"https://www.tensorflow.org/images/GitHub-Mark-32px.png\" />View source on GitHub</a>\n",
        "  </td>\n",
        "</table>"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "0GVst7yy6Aww"
      },
      "cell_type": "markdown",
      "source": [
        "## Abstract\n"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "Lt4RS9whJhQh"
      },
      "cell_type": "markdown",
      "source": [
        "\n",
        "In this colab we demonstrate how to fit a Generalized Linear Mixed-effects Model using Variational Inference and TensorFlow.\n"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "0-lfIBVAzi7D"
      },
      "cell_type": "markdown",
      "source": [
        "## Model Family"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "ljSRsKrXwqb6"
      },
      "cell_type": "markdown",
      "source": [
        "[Generalized linear mixed-effect models](https://en.wikipedia.org/wiki/Generalized_linear_mixed_model) (GLMM) are similar to [generalized linear models](https://en.wikipedia.org/wiki/Generalized_linear_model) (GLM) except that they incorporate a sample specific noise into the predicted linear response.  This is useful in part because it allows rarely seen features to share information with more commonly seen features.\n",
        "\n"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "H-B38entvltq"
      },
      "cell_type": "markdown",
      "source": [
        "As a generative process, a Generalized Linear Mixed-effects Model (GLMM) is characterized by:\n",
        "\n",
        "\\begin{align}\n",
        "\\text{for } & r = 1\\ldots R:  \\hspace{2.45cm}\\text{# for each random-effect group}\\\\\n",
        " &\\begin{aligned}\n",
        "  \\text{for } &c = 1\\ldots |C_r|:  \\hspace{1.3cm}\\text{# for each category (\"level\") of group $r$}\\\\\n",
        "  &\\begin{aligned}\n",
        "    \\beta_{rc}\n",
        "    &\\sim \\text{MultivariateNormal}(\\text{loc}=0_{D_r}, \\text{scale}=\\Sigma_r^{1/2})\n",
        "  \\end{aligned}\n",
        "\\end{aligned}\\\\\\\\\n",
        "\\text{for } & i = 1 \\ldots N:  \\hspace{2.45cm}\\text{# for each sample}\\\\\n",
        "&\\begin{aligned}\n",
        "  &\\eta_i = \\underbrace{\\vphantom{\\sum_{r=1}^R}x_i^\\top\\omega}_\\text{fixed-effects} + \\underbrace{\\sum_{r=1}^R z_{r,i}^\\top \\beta_{r,C_r(i) }}_\\text{random-effects} \\\\\n",
        "  &Y_i|x_i,\\omega,\\{z_{r,i} , \\beta_r\\}_{r=1}^R \\sim \\text{Distribution}(\\text{mean}= g^{-1}(\\eta_i))\n",
        "\\end{aligned}\n",
        "\\end{align}\n"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "3gZmFJXAHwfy"
      },
      "cell_type": "markdown",
      "source": [
        "where:\n",
        "\n",
        "\\begin{align}\n",
        "R &= \\text{number of random-effect groups}\\\\\n",
        "|C_r| &= \\text{number of categories for group $r$}\\\\\n",
        "N &= \\text{number of training samples}\\\\\n",
        "x_i,\\omega &\\in \\mathbb{R}^{D_0}\\\\\n",
        "D_0 &= \\text{number of fixed-effects}\\\\\n",
        "C_r(i) &= \\text{category (under group $r$) of the $i$th sample}\\\\\n",
        "z_{r,i} &\\in \\mathbb{R}^{D_r}\\\\\n",
        "D_r &= \\text{number of random-effects associated with group $r$}\\\\\n",
        "\\Sigma_{r} &\\in \\{S\\in\\mathbb{R}^{D_r \\times D_r} : S \\succ 0 \\}\\\\\n",
        "\\eta_i\\mapsto g^{-1}(\\eta_i) &= \\mu_i, \\text{inverse link function}\\\\\n",
        "\\text{Distribution} &=\\text{some distribution parameterizable solely by its mean}\n",
        "\\end{align}"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "5AYonR45P1Hr"
      },
      "cell_type": "markdown",
      "source": [
        "In words, this says that every category of each group is associated with an iid MVN, $\\beta_{rc}$. Although the $\\beta_{rc}$ draws are always independent, they are only indentically distributed for a group $r$; notice there is exactly one $\\Sigma_r$ for each $r\\in\\{1,\\ldots,R\\}$.\n",
        "\n",
        "When affinely combined with a sample's group's features ($z_{r,i}$), the result is sample-specific noise on the $i$-th predicted linear response (which is otherwise $x_i^\\top\\omega$)."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "__dP1MdYKda0"
      },
      "cell_type": "markdown",
      "source": [
        "When we estimate $\\{\\Sigma_r:r\\in\\{1,\\ldots,R\\}\\}$ we're essentially estimating the amount of noise a random-effect group carries which would otherwise drown out the signal present in $x_i^\\top\\omega$."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "0EZXZzlYSbM7"
      },
      "cell_type": "markdown",
      "source": [
        "There are a variety of options for the $\\text{Distribution}$ and [inverse link function](https://en.wikipedia.org/wiki/Generalized_linear_model#Link_function), $g^{-1}$. Common choices are:\n",
        "- $Y_i\\sim\\text{Normal}(\\text{mean}=\\eta_i, \\text{scale}=\\sigma)$,\n",
        "- $Y_i\\sim\\text{Binomial}(\\text{mean}=n_i \\cdot \\text{sigmoid}(\\eta_i), \\text{total_count}=n_i)$, and, \n",
        "- $Y_i\\sim\\text{Poisson}(\\text{mean}=\\exp(\\eta_i))$.\n",
        "\n",
        "For more possibilities, see the [`tfp.glm`](https://github.com/tensorflow/probability/tree/master/tensorflow_probability/python/glm) module."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "XajrojApx5cR"
      },
      "cell_type": "markdown",
      "source": [
        "## Variational Inference"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "fIQn1mlYAUzx"
      },
      "cell_type": "markdown",
      "source": [
        "Unfortunately, finding the maximum likelihood estimates of the parameters $\\beta,\\{\\Sigma_r\\}_r^R$ entails a non-analytical integral. To circumvent this problem, we instead find the parameters which minimize an upper bound.  Writing $\\phi(x)=\\exp(-x^2/2)/\\sqrt{2\\pi}$ we find:"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "9epWkJw0x8hw"
      },
      "cell_type": "markdown",
      "source": [
        "$\\begin{align}\n",
        "-\\log \\overbrace{p(\\{y\\}_i^N|\\{x_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R)}^{\\text{evidence}}\n",
        "& = -\\log \\int_{\\mathbb{R}^{\\sum_r |C_r|D_r}} \\overbrace{\\left(\\prod_i^N p(y_i | \\eta_i(u)) \\right)}^{\\text{likelihood}} \\overbrace{\\left(\\prod_r^R |\\Sigma_r^{-1/2}| \\prod_c^{|C_r|} \\phi(\\Sigma_r^{-1/2} u_{rc}) \\right)}^{\\text{prior}} \\, du\\\\\n",
        "& = -\\log \\int_{\\mathbb{R}^{\\sum_r |C_r|D_r}} \\frac{\n",
        "  q_\\lambda(u|\\{x_i,y_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R)\n",
        " }{\n",
        "  q_\\lambda(u|\\{x_i,y_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R)} \\left( \\prod_i^N p(y_i | \\eta_i(u)) \\right) \\left(\\prod_r^R |\\Sigma_r^{-1/2}| \\prod_c^{|C_r|} \\phi(\\Sigma_r^{-1/2} u_{rc}) \\right) \\, du\\\\\n",
        "& \\le \\text{E}_{q_\\lambda(U|\\{x_i,y_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R)} \\left[ -\\log\n",
        "\\frac{\n",
        "  \\left( \\prod_i^N p(y_i | \\eta_i(U)) \\right) \\left(\\prod_r^R |\\Sigma_r^{-1/2}| \\prod_c^{|C_r|} \\phi(\\Sigma_r^{-1/2} U_{rc}) \\right)\n",
        "  }{\n",
        "    q_\\lambda(U|\\{x_i,y_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R)\n",
        "  }  \\right]\\\\\n",
        "&= \\text{KL}\\left[q_\\lambda(U|\\{x_i,y_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R) \\Bigg| \\left( \\prod_i^N p(y_i | \\eta_i(U)) \\right) \\left(\\prod_r^R |\\Sigma_r^{-1/2}| \\prod_c^{|C_r|} \\phi(\\Sigma_r^{-1/2} U_{rc}) \\right) \\right]\n",
        "\\end{align}$"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "ecNUsIBQN-zz"
      },
      "cell_type": "markdown",
      "source": [
        "The inequality follows from [Jensen's Inequality](https://en.wikipedia.org/wiki/Jensen%27s_inequality)."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "GCDXHiEx-dY5"
      },
      "cell_type": "markdown",
      "source": [
        "So, instead of solving:\n",
        "\n",
        "$$\\begin{align}\n",
        "\\{\\beta^*, \\{\\Sigma_r^*\\}_r^R\\} = \\operatorname{\\arg\\min}_{\\beta,\\{\\Sigma_r\\}_r^R} \\left\\{\n",
        " -\\log  p(\\{y\\}_i^N|\\{x_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R)\n",
        "\\right\\}\n",
        "\\end{align}$$\n",
        "\n",
        "we solve:\n",
        "\n",
        "$$\\begin{align}\n",
        "\\{\\lambda^*, \\beta^*, \\{\\Sigma_r^*\\}_r^R\\} = \\operatorname{\\arg\\min}_{\\lambda,\\beta,\\{\\Sigma_r\\}_r^R} \\left\\{\n",
        "\\text{KL}\\left[q_\\lambda(U|\\{x_i,y_i,z_i\\}_i^N,\\beta,\\{\\Sigma_r\\}_r^R) \\Bigg| \\left( \\prod_i^N p(y_i | \\eta_i(U)) \\right) \\left(\\prod_r^R |\\Sigma_r^{-1/2}| \\prod_c^{|C_r|} \\phi(\\Sigma_r^{-1/2} U_{rc}) \\right) \\right]\n",
        "\\right\\}\n",
        "\\end{align}$$"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "Nu8B7ylx3UdL"
      },
      "cell_type": "markdown",
      "source": [
        "## Toy Problem"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "vDmAfghTJcLo"
      },
      "cell_type": "markdown",
      "source": [
        "[Gelman et al.'s (2007) \"radon dataset\"](http://www.stat.columbia.edu/~gelman/arm/) is a dataset sometimes used to demonstrate approaches for regression. (E.g., this closely related [PyMC3 blog post](http://twiecki.github.io/blog/2014/03/17/bayesian-glms-3/).) The radon dataset contains indoor measurements of Radon taken throughout the United States. [Radon](https://en.wikipedia.org/wiki/Radon) is naturally ocurring radioactive gas which is [toxic](http://www.radon.com/radon_facts/) in high concentrations.\n",
        "\n",
        "For our demonstration, let's suppose we're interested in validating the hypothesis that Radon levels are higher in households containing a basement. We also suspect Radon concentration is related to soil-type, i.e., geography matters.\n",
        "\n",
        "To frame this as an ML problem, we'll try to predict log-radon levels based on a linear function of the floor on which the reading was taken.  We'll also use the county as a random-effect and in so doing account for variances due to geography. In other words, we'll use a [generalized linear mixed-effect model](https://en.wikipedia.org/wiki/Generalized_linear_mixed_model)."
      ]
    },
    {
      "metadata": {
        "id": "ht5yNGWOMEPC",
        "colab_type": "code",
        "outputId": "ffb6e8f2-4828-468d-8d1e-e58f4876a123",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 68
        }
      },
      "cell_type": "code",
      "source": [
        "!pip install --upgrade tensorflow-probability\n"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already up-to-date: tensorflow-probability in /usr/local/lib/python2.7/dist-packages (0.5.0)\n",
            "Requirement already satisfied, skipping upgrade: six>=1.10.0 in /usr/local/lib/python2.7/dist-packages (from tensorflow-probability) (1.11.0)\n",
            "Requirement already satisfied, skipping upgrade: numpy>=1.13.3 in /usr/local/lib/python2.7/dist-packages (from tensorflow-probability) (1.14.6)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "_zr34b0IBqgY",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "%matplotlib inline\n",
        "\n",
        "from __future__ import absolute_import\n",
        "from __future__ import division\n",
        "from __future__ import print_function\n",
        "\n",
        "from pprint import pprint\n",
        "import collections\n",
        "import matplotlib.pyplot as plt; plt.style.use('ggplot')\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import seaborn as sns; sns.set_context('notebook')\n",
        "import tensorflow as tf\n",
        "import warnings\n",
        "import os\n",
        "from scipy.stats import norm"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "UzepxbtrpHpD"
      },
      "cell_type": "markdown",
      "source": [
        "### Obtain Dataset:"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "Z5aSvBj8v_vp",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "from six.moves import urllib\n",
        "CACHE_DIR = os.path.join(os.sep, 'tmp', 'radon')\n",
        "\n",
        "def cache_or_download_file(cache_dir, url_base, filename):\n",
        "  \"\"\"Read a cached file or download it.\"\"\"\n",
        "  filepath = os.path.join(cache_dir, filename)\n",
        "  if tf.gfile.Exists(filepath):\n",
        "    return filepath\n",
        "  if not tf.gfile.Exists(cache_dir):\n",
        "    tf.gfile.MakeDirs(cache_dir)\n",
        "  url = os.path.join(url_base, filename)\n",
        "  print(\"Downloading {url} to {filepath}.\".format(url=url, filepath=filepath))\n",
        "  urllib.request.urlretrieve(url, filepath)\n",
        "  return filepath\n",
        "\n",
        "\n",
        "def download_radon_dataset(cache_dir=CACHE_DIR):\n",
        "  \"\"\"Download the radon dataset and read as Pandas dataframe.\"\"\"\n",
        "  url_base = 'http://www.stat.columbia.edu/~gelman/arm/examples/radon/'\n",
        "  # Alternative source:\n",
        "  # url_base = ('https://raw.githubusercontent.com/pymc-devs/uq_chapter/'\n",
        "  #             'master/reference/data/')\n",
        "  srrs2 = pd.read_csv(cache_or_download_file(cache_dir, url_base, 'srrs2.dat'))\n",
        "  srrs2.rename(columns=str.strip, inplace=True)\n",
        "  cty = pd.read_csv(cache_or_download_file(cache_dir, url_base, 'cty.dat'))\n",
        "  cty.rename(columns=str.strip, inplace=True)\n",
        "  return srrs2, cty\n",
        "\n",
        "\n",
        "def preprocess_radon_dataset(srrs2, cty, state='MN'):\n",
        "  \"\"\"Preprocess radon dataset as done in \"Bayesian Data Analysis\" book.\"\"\"\n",
        "  srrs2 = srrs2[srrs2.state==state].copy()\n",
        "  cty = cty[cty.st==state].copy()\n",
        "  \n",
        "  # We will now join datasets on Federal Information Processing Standards\n",
        "  # (FIPS) id, ie, codes that link geographic units, counties and county\n",
        "  # equivalents. http://jeffgill.org/Teaching/rpqm_9.pdf\n",
        "  srrs2['fips'] = 1000 * srrs2.stfips + srrs2.cntyfips\n",
        "  cty['fips'] = 1000 * cty.stfips + cty.ctfips\n",
        "\n",
        "  df = srrs2.merge(cty[['fips', 'Uppm']], on='fips')\n",
        "  df = df.drop_duplicates(subset='idnum')\n",
        "  df = df.rename(index=str, columns={'Uppm': 'uranium_ppm'})\n",
        "  \n",
        "  df['radon'] = df.activity.apply(lambda x: x if x > 0. else 0.1)\n",
        "  \n",
        "  # Remap categories to start from 0 and end at max(category).\n",
        "  county_name = sorted(df.county.unique())\n",
        "  df['county'] = df.county.astype(\n",
        "      pd.api.types.CategoricalDtype(categories=county_name)).cat.codes\n",
        "  county_name = map(str.strip, county_name)\n",
        "  \n",
        "  df['log_radon'] = df['radon'].apply(np.log)\n",
        "  df['log_uranium_ppm'] = df['uranium_ppm'].apply(np.log) \n",
        "  df = df[['log_radon', 'floor', 'county', 'log_uranium_ppm']]\n",
        " \n",
        "  return df, county_name"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "DLn6q9fSNizw",
        "outputId": "7c818c68-db85-4946-f4a1-846189f1745d",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 51
        }
      },
      "cell_type": "code",
      "source": [
        "df, counties = preprocess_radon_dataset(*download_radon_dataset())"
      ],
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Downloading http://www.stat.columbia.edu/~gelman/arm/examples/radon/srrs2.dat to /tmp/radon/srrs2.dat.\n",
            "Downloading http://www.stat.columbia.edu/~gelman/arm/examples/radon/cty.dat to /tmp/radon/cty.dat.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "id": "2rQFiLnvAhm3",
        "colab_type": "code",
        "outputId": "0f9160dc-1872-455a-b69f-7dd2c66d9b08",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "cell_type": "code",
      "source": [
        "type(df)\n"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "pandas.core.frame.DataFrame"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        }
      ]
    },
    {
      "metadata": {
        "id": "SrJ_QEE1a8zt",
        "colab_type": "code",
        "outputId": "8635977f-865c-4580-e043-aa435a465e81",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 102
        }
      },
      "cell_type": "code",
      "source": [
        "np.random.seed(1234)\n",
        "#add another variable for test\n",
        "df['log_u2'] = df['log_uranium_ppm']+norm.rvs(size=len(df['log_uranium_ppm']))\n",
        "df['floor2'] = np.array(df['floor']+np.rint(0.5*norm.rvs(size=len(df['log_uranium_ppm']))))\n",
        "df['floor3']=df['floor']\n",
        "df['log_uranium_ppm'][0]=-2"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py:6: SettingWithCopyWarning: \n",
            "A value is trying to be set on a copy of a slice from a DataFrame\n",
            "\n",
            "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
            "  \n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "5pXUQLVxi8E0",
        "outputId": "12a9e645-6b59-4f43-dad2-ba2f2c6ce6f3",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 385
        }
      },
      "cell_type": "code",
      "source": [
        "df.head(10)\n",
        "\n",
        "df['floor2'].value_counts().plot(kind='bar')\n",
        "plt.xlabel('Floor2')\n",
        "plt.ylabel('Count')\n"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0,0.5,u'Count')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 8
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFfCAYAAAC1E07gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGVtJREFUeJzt3X9s1IX9x/HXtb1Sbnj2oDo8Bwj1\nqHQyaFYIlQQca8REv8wNUDeyhIURGAYikDk7AUHoADEyt4IkZoEwRhyrm0PpNiAw10W7oWy2FFuh\nneNHBwWOEyhcC+19/3B26yqlQD/36fv6fCQm8Kn3+bz7+Vif/Xw+1089sVgsJgAAYFKS2wMAAIAb\nR8gBADCMkAMAYBghBwDAMEIOAIBhhBwAAMNS3B7gRtXV1bk9gmOCwWBCf36JjGNnG8fPtkQ+fsFg\n8Kof44wcAADDCDkAAIYRcgAADCPkAAAYRsgBADCMkAMAYBghBwDAMEIOAIBhhBwAAMMIOQAAhhFy\nAAAMI+QAABhGyAEAMMzsbz+Lp+aZk+K6vaNx3ZqU/Mr2OG8RANBVOCMHAMAwQg4AgGGEHAAAwwg5\nAACGEXIAAAwj5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCAYYQcAADDCDkAAIYRcgAADCPkAAAYRsgB\nADCMkAMAYBghBwDAMEIOAIBhhBwAAMMIOQAAhhFyAAAMS3FqxZWVlXrxxRc1YMAASdLAgQM1adIk\nFRUVqaWlRenp6Zo7d668Xq9KS0tVUlIij8ej/Px8TZgwwamxAABIKI6FXJKys7O1cOHC1r+vX79e\nEydOVF5enrZu3aq9e/dq3LhxKi4u1sqVK5WSkqKCggKNHj1affr0cXI0AAASQlwvrVdWVio3N1eS\nlJubq/Lych0+fFiZmZny+XxKTU1VVlaWqqqq4jkWAABmOXpGfuzYMa1evVoXLlzQ1KlT1djYKK/X\nK0ny+/2KRCKKRCLy+/2tr/l0OQAAuDbHQn7HHXdo6tSpysvL08mTJ7Vs2TI1Nzd32fqDwWCXreta\njsZtS+6I577sCdiftnH8bOuJx8+xkPft21f33XefJKl///5KT09XTU2NmpqalJqaqnA4rEAgoEAg\n0OYMPBwOKxQKXXP9dXV1To3e47Avu04wGGR/Gsbxsy2Rj19H36A4do+8tLRU27dvlyRFIhF9/PHH\nuv/++1VWViZJKisr08iRIxUKhVRTU6OGhgZFo1FVV1dr2LBhTo0FAEBCceyMPDc3Vy+99JLeffdd\nXblyRd/97nc1ePBgFRUVaffu3crIyND48eOVkpKiadOmqbCwUB6PR1OmTJHP53NqLAAAEoonFovF\n3B7iRsTz8knzzElx25Ybkl/Z7vYICSORL+31BBw/2xL5+LlyaR0AADiPkAMAYBghBwDAMEIOAIBh\nhBwAAMMIOQAAhhFyAAAMI+QAABhGyAEAMIyQAwBgGCEHAMAwQg4AgGGEHAAAwwg5AACGEXIAAAwj\n5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCAYYQcAADDCDkAAIYRcgAADCPkAAAYRsgBADCMkAMAYBgh\nBwDAMEIOAIBhhBwAAMMIOQAAhhFyAAAMI+QAABhGyAEAMIyQAwBgGCEHAMAwQg4AgGGEHAAAwwg5\nAACGEXIAAAwj5AAAGEbIAQAwjJADAGAYIQcAwLAUJ1fe1NSkhQsXavLkybr33ntVVFSklpYWpaen\na+7cufJ6vSotLVVJSYk8Ho/y8/M1YcIEJ0cCACChOHpG/tprr6lPnz6SpG3btmnixIl67rnn1L9/\nf+3du1fRaFTFxcVavHixli5dqh07dujChQtOjgQAQEJxLOTHjx/XsWPHlJOTI0mqrKxUbm6uJCk3\nN1fl5eU6fPiwMjMz5fP5lJqaqqysLFVVVTk1EgAACcexS+ubN2/WjBkz9Mc//lGS1NjYKK/XK0ny\n+/2KRCKKRCLy+/2tr/l0eWcEg8Eun/lqjsZtS+6I577sCdiftnH8bOuJx8+RkL/11lsaOnSobr/9\ndidWL0mqq6tzbN09Dfuy6wSDQfanYRw/2xL5+HX0DYojId+/f7/q6+u1f/9+nTlzRl6vV2lpaWpq\nalJqaqrC4bACgYACgUCbM/BwOKxQKOTESAAAJCRHQj5//vzWP2/btk233367qqurVVZWpnHjxqms\nrEwjR45UKBTShg0b1NDQoOTkZFVXV2v69OlOjAQAQEJy9MfP/tujjz6qoqIi7d69WxkZGRo/frxS\nUlI0bdo0FRYWyuPxaMqUKfL5fPEaCQAA8xwP+aOPPtr658WLF7f7+JgxYzRmzBinxwAAICHxZDcA\nAAwj5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCAYYQcAADDCDkAAIYRcgAADCPkAAAYRsgBADCMkAMA\nYBghBwDAMEIOAIBhhBwAAMMIOQAAhhFyAAAMI+QAABhGyAEAMIyQAwBgGCEHAMAwQg4AgGGEHAAA\nwwg5AACGEXIAAAwj5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCAYYQcAADDCDkAAIYRcgAADCPkAAAY\nRsgBADCMkAMAYBghBwDAMEIOAIBhhBwAAMM6FfLjx4+3W/bhhx92+TAAAOD6dBjyhoYGnThxQuvX\nr1d9fb1OnjypkydPqq6uTuvWrYvXjAAA4CpSOvrghx9+qB07duijjz7SsmXLWpcnJSVpxIgRjg8H\nAAA61mHIc3JylJOTo507d+qBBx6I10wAAKCTOgz5p0aPHq2SkhJduHBBsVisdfljjz121dc0NjZq\n3bp1+vjjj3X58mVNnjxZgwYNUlFRkVpaWpSenq65c+fK6/WqtLRUJSUl8ng8ys/P14QJE27+MwMA\noAfoVMhXrlypQYMG6bbbbpPH4+nUit977z1lZmbqa1/7mk6dOqUVK1YoKytLEydOVF5enrZu3aq9\ne/dq3LhxKi4u1sqVK5WSkqKCggKNHj1affr0ualPDACAnqBTIU9LS9OcOXOua8X33Xdf65/PnDmj\nvn37qrKyUjNnzpQk5ebmavv27QoGg8rMzJTP55MkZWVlqaqqSrm5ude1PQAAeqJOhTwUCun48eO6\n8847r3sDixYt0pkzZ/T0009r+fLl8nq9kiS/369IJKJIJCK/39/673+6/FqCweB1z3KjjsZtS+6I\n577sCdiftnH8bOuJx69TIX///fe1Y8cO+f1+JSX95yfWXn755Wu+dsWKFfroo4/005/+tM399ZtV\nV1fXZevq6diXXScYDLI/DeP42ZbIx6+jb1A6FfKnnnrqujdaW1srv9+vjIwM3XXXXWpublbv3r3V\n1NSk1NRUhcNhBQIBBQKBNmfg4XBYoVDourcHAEBP1Kknu1VUVHzmPx05ePCg3nzzTUlSJBJRNBrV\n8OHDVVZWJkkqKyvTyJEjFQqFVFNTo4aGBkWjUVVXV2vYsGE3+WkBANAzdOqMvKqqqvXPV65c0aFD\nh5SVldXhj4k98MADevnll7VkyRI1NTVpxowZyszMVFFRkXbv3q2MjAyNHz9eKSkpmjZtmgoLC+Xx\neDRlypTWN74BAICOeWI3cOO6sbFR69ev1/z5852YqVPieR+keeakuG3LDcmvbHd7hISRyPfoegKO\nn22JfPw6ukd+Q7/9rFevXjpx4sQNDwQAALpGpy6tL1mypM2DYMLhsAYOHOjYUAAAoHM6FfLHH3+8\nzd99Pp8GDRrkyEAAAKDzOnVpPTs7Wx6PR7W1taqtrVVTU1OnH9UKAACc06mQ//KXv9SWLVt09uxZ\nhcNhbdy4Ub/5zW+cng0AAFxDpy6tV1ZWavny5a1PdWtubtazzz6rr3/9644OBwAAOtapM/JYLNbm\n0azJyclcWgcAoBvo1Bn54MGDtXr1ag0fPlySVF5erszMTEcHAwAA13bNkNfX12v69Ol65513dOjQ\nIXk8HmVnZ2vSpMR+SAoAABZ0eGm9oqJCixcvVjQa1dixYzV9+nR95Stf0c6dO1VbWxuvGQEAwFV0\nGPLi4mI988wzbZ59PnDgQP3gBz/Qq6++6vhwAACgYx2GPBaLfeYT3AYMGKDLly87NhQAAOicDkMe\njUav+rHz5893+TAAAOD6dBjyAQMGaOfOne2W//a3v9Xdd9/t2FAAAKBzOnzX+re//W2tWbNGf/rT\nn5SZmamWlhZVV1erd+/eKigoiNeMAADgKjoMeXp6ugoLC1VRUaGjR48qKSlJeXl5ys7Ojtd8AACg\nA516IMzw4cNbHwYDAAC6j049ohUAAHRPhBwAAMMIOQAAhhFyAAAMI+QAABhGyAEAMIyQAwBgGCEH\nAMAwQg4AgGGderIbYNnXflEVx63Fc1vSb6fdE9ftAeh+OCMHAMAwQg4AgGGEHAAAwwg5AACGEXIA\nAAwj5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCAYYQcAADDCDkAAIYRcgAADCPkAAAYRsgBADCMkAMA\nYFiKkyvfsmWLPvjgA7W0tOiRRx5RZmamioqK1NLSovT0dM2dO1der1elpaUqKSmRx+NRfn6+JkyY\n4ORYAAAkDMdCfuDAAR09elSFhYU6f/68nnrqKQ0fPlwTJ05UXl6etm7dqr1792rcuHEqLi7WypUr\nlZKSooKCAo0ePVp9+vRxajQAABKGY5fWs7OzNX/+fEnS5z73OTU2NqqyslK5ubmSpNzcXJWXl+vw\n4cPKzMyUz+dTamqqsrKyVFVV5dRYAAAkFMfOyJOSkpSWliZJ2rNnj3JycvT+++/L6/VKkvx+vyKR\niCKRiPx+f+vrPl1+LcFg0JnBP8PRuG3JHfHcl+5I3G8ME//YxR/71LaeePwcvUcuSfv27dOePXu0\naNEizZs3r8vWW1dX12Xr6unYl3Zx7LpWMBhknxqWyMevo29QHH3X+t///nf9+te/1g9/+EP5fD6l\npaWpqalJkhQOhxUIBBQIBNqcgX+6HAAAXJtjIb948aK2bNmip59+uvWNa8OHD1dZWZkkqaysTCNH\njlQoFFJNTY0aGhoUjUZVXV2tYcOGOTUWAAAJxbFL62+//bbOnz+vtWvXti574okntGHDBu3evVsZ\nGRkaP368UlJSNG3aNBUWFsrj8WjKlCny+XxOjQUAQEJxLOT5+fnKz89vt3zx4sXtlo0ZM0Zjxoxx\nahQAABIWT3YDAMAwQg4AgGGEHAAAwwg5AACGEXIAAAwj5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCA\nYYQcAADDCDkAAIYRcgAADCPkAAAYRsgBADCMkAMAYBghBwDAMEIOAIBhhBwAAMMIOQAAhhFyAAAM\nI+QAABhGyAEAMIyQAwBgGCEHAMAwQg4AgGGEHAAAwwg5AACGEXIAAAwj5AAAGEbIAQAwjJADAGAY\nIQcAwDBCDgCAYYQcAADDCDkAAIYRcgAADCPkAAAYRsgBADCMkAMAYBghBwDAMEIOAIBhKU6u/MiR\nI1qzZo0eeughPfjggzp9+rSKiorU0tKi9PR0zZ07V16vV6WlpSopKZHH41F+fr4mTJjg5FgAjHjj\nl5E4bzG+2/u/x9Ljuj0kJsfOyKPRqDZu3Kh77723ddm2bds0ceJEPffcc+rfv7/27t2raDSq4uJi\nLV68WEuXLtWOHTt04cIFp8YCACChOBZyr9ergoICBQKB1mWVlZXKzc2VJOXm5qq8vFyHDx9WZmam\nfD6fUlNTlZWVpaqqKqfGAgAgoTh2aT05OVnJycltljU2Nsrr9UqS/H6/IpGIIpGI/H5/67/z6XIA\nAHBtjt4jd1IwGIzbto7GbUvuiOe+dEfiXuFJ/GOX2N/UJ/7xi7+euE/jGvK0tDQ1NTUpNTVV4XBY\ngUBAgUCgzRl4OBxWKBS65rrq6uqcHLVHYV/axbGzjePXtYLBYMLu046+QYnrj58NHz5cZWVlkqSy\nsjKNHDlSoVBINTU1amhoUDQaVXV1tYYNGxbPsQAAMMuxM/La2lpt3rxZp06dUnJyssrKyjRv3jyt\nW7dOu3fvVkZGhsaPH6+UlBRNmzZNhYWF8ng8mjJlinw+n1NjAQCQUBwL+ZAhQ7R06dJ2yxcvXtxu\n2ZgxYzRmzBinRgEAIGHxZDcAAAwj5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCAYYQcAADDCDkAAIYR\ncgAADCPkAAAYRsgBADCMkAMAYBghBwDAMEIOAIBhhBwAAMMIOQAAhhFyAAAMI+QAABhGyAEAMIyQ\nAwBgGCEHAMAwQg4AgGGEHAAAwwg5AACGEXIAAAwj5AAAGEbIAQAwjJADAGAYIQcAwDBCDgCAYYQc\nAADDCDkAAIYRcgAADCPkAAAYluL2AACAxPSTn/zE7REcNW/ePLdHkMQZOQAAphFyAAAMI+QAABhG\nyAEAMIyQAwBgGCEHAMAwQg4AgGGEHAAAw7rNA2E2bdqkQ4cOyePxaPr06br77rvdHgkAgG6vW5yR\nHzx4UCdOnFBhYaFmz56tjRs3uj0SAAAmdIuQV1RUaNSoUZKkL3zhC2poaNDFixddngoAgO6vW1xa\nj0QiGjJkSOvf/X6/IpGIfD7fVV8TDAbjMdondrwbv22hy+37fhz/W0GXmjWfY2fZqlWr3B6hR+gW\nZ+T/KxaLuT0CAAAmdIuQBwIBRSKR1r+fPXtWgUDAxYkAALChW4R8xIgRKisrkyTV1tYqEAiod+/e\nLk8FAED354l1k+vYv/jFL/TBBx/I4/FoxowZuuuuu9weCQCAbq/bhBwAAFy/bnFpHQAA3BhCDgCA\nYYQcAADDCDkAAIZ1iye79WTHjx9XRUVF68/RBwIBjRgxQv3793d5MnQGxw9wB197/8G71l302muv\nqby8XDk5OfL7/ZKkcDis9957T2PHjtXDDz/s8oToCMfPtoaGBr3++uvtYjBy5EhNmjSJZ1l0Y3zt\n/Y8YXLNo0aJYS0tLu+VXrlyJLVq0yIWJcD04frb96Ec/iu3atSt29uzZWHNzc6y5uTl26tSp2Btv\nvBFbvXq12+OhA3zttcU9chc1Nzfr7Nmz7ZafPXuW580bwPGz7dKlS8rPz1d6erqSkpKUlJSkjIwM\nPfzww2poaHB7PHSAr722uEfuoscff1wrVqzQLbfcoltuuUXSJ78J7tKlS5o5c6bL0+FaOH62+Xw+\nvfHGGxo9enTr5dlIJKJ33nlHffr0cXk6dISvvba4R94N1NfXt7lHd9ttt7k8Ea4Hx8+mixcvqqSk\npPUeucfjaX3D1IMPPqi0tDS3R8Q18LX3CULeTe3bt0+jRo1yewzcII4f4I6e+LXHPfJuIBqN6sSJ\nEzpx4oSi0agkcY/OuPr6erdHwE145ZVX3B4BnfBZ56FHjhxxYRJ3cY/cRTU1Ndq4caMaGhrk9/sV\ni8V09uxZ9e3bVzNmzHB7PNyEd999Vw899JDbY+AGTZ061e0R0IG//vWv2rRpkxobG5WTk6MZM2a0\n/rjggQMHNHnyZJcnjC9C7qJNmzbpe9/7nu688842y2tra/Wzn/1My5Ytc2kydMYf/vCHz1wei8UU\nDofjPA26UmVlpcaOHev2GLiK119/Xc8//7x8Pp/27NmjFStW6JlnnpHP53N7NFdwad1FsVisXcQl\naciQIWppaXFhIlyPN998U//85z917ty5Nv+cP39ezc3Nbo+Hm7B79263R0AHkpKS1KdPHyUlJSk/\nP1+PPPKICgsLde7cObdHcwVn5C4KhUJavXq1Ro0a1ebHX/7yl79o2LBhLk+Ha/n+97+vjRs36jvf\n+Y68Xm+bjx08eNClqdBZBQUFV/1YXV1dHCfB9brnnnu0atUqLViwQKmpqRo1apS8Xq+WL1+u8+fP\nuz1e3PGudZcdPHhQBw4caPe84KFDh7o8GTqjsbFRXq9XSUltL27V1tZqyJAhLk2Fzli3bp0GDx7c\n7h3OsVhML730kgoLC12aDJ1RWVmp7OxseTye1mUXL17U22+/rfz8fBcniz9CDqBHunz5sn7+85/r\nW9/6VrufGX/xxRe1YMEClyYDrg8hB4B/W7t2rebPn+/2GMB14c1uAPBvPfXNUrCNkAPAv+Xl5bk9\nAnDduLQOAIBhnJEDAGAYIQcAwDAeCAP0QPX19XryySfbPa/A5/Ppy1/+sr761a/e9DZ+//vf6623\n3lJSUpI+//nPa86cOUpJ4X85QFfjqwroofx+v5YuXdpm2bp167pk3UeOHNHvfvc7vfDCC/J6vXrh\nhRf05z//Wffff3+XrB/AfxByAJ9pz5492rVrl3r16qVbb71Vs2bNks/n0/79+1VcXKxevXopNTVV\ns2bNUt++ffXEE08oLy+v9Wx/1apVrY+u9fv9PfLRmUA8cI8cQDunT5/Wr371Ky1ZskRLly5Vv379\ntGPHDjU2NmrDhg1auHChnn32WeXk5OjVV19tfd0dd9yhBQsWKCkpqfXXStbX1+tvf/sbP9oFOIQz\ncqCHOnfuXLtL67169ZL0ybPiBw8e3BrjL37xi9q1a5f+9a9/6dZbb1W/fv0kSdnZ2dq1a1fr67Oy\nstqs79ixY1qzZo1mzZqljIwMBz8boOci5EAP1dE98v/+RRTSJ79IpDP++81sx44d0/PPP685c+bo\nnnvuublhAVwVl9YBtDNkyBD94x//0KVLlyRJFRUVGjp0qILBoM6dO6fTp09LksrLyxUKhdq9/sqV\nK/rxj3+sJ598kogDDuOMHEA7/fr102OPPably5fL6/Wqb9+++uY3v6nU1FTNnj1ba9euldfrVVpa\nmmbPnt3u9fv27dPp06e1efPm1mVf+tKX9I1vfCOenwbQI/CIVgAADOPSOgAAhhFyAAAMI+QAABhG\nyAEAMIyQAwBgGCEHAMAwQg4AgGH/D8OlYx8BWYnuAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fc8478cb050>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "ZSPUVvLVnEZs",
        "colab_type": "code",
        "outputId": "f0b2bf0e-5209-4698-8d2a-a99b49df1a8e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 359
        }
      },
      "cell_type": "code",
      "source": [
        "df.head(10)"
      ],
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>log_radon</th>\n",
              "      <th>floor</th>\n",
              "      <th>county</th>\n",
              "      <th>log_uranium_ppm</th>\n",
              "      <th>log_u2</th>\n",
              "      <th>floor2</th>\n",
              "      <th>floor3</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>0.788457</td>\n",
              "      <td>1</td>\n",
              "      <td>0</td>\n",
              "      <td>-2.000000</td>\n",
              "      <td>-0.217612</td>\n",
              "      <td>1.0</td>\n",
              "      <td>1</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.788457</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>-0.689048</td>\n",
              "      <td>-1.880023</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>1.064711</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>-0.689048</td>\n",
              "      <td>0.743659</td>\n",
              "      <td>-1.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0</td>\n",
              "      <td>0</td>\n",
              "      <td>-0.689048</td>\n",
              "      <td>-1.001699</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>1.131402</td>\n",
              "      <td>0</td>\n",
              "      <td>1</td>\n",
              "      <td>-0.847313</td>\n",
              "      <td>-1.567902</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>0.916291</td>\n",
              "      <td>0</td>\n",
              "      <td>1</td>\n",
              "      <td>-0.847313</td>\n",
              "      <td>0.039850</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>0.405465</td>\n",
              "      <td>0</td>\n",
              "      <td>1</td>\n",
              "      <td>-0.847313</td>\n",
              "      <td>0.012276</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>0.000000</td>\n",
              "      <td>0</td>\n",
              "      <td>1</td>\n",
              "      <td>-0.847313</td>\n",
              "      <td>-1.483836</td>\n",
              "      <td>1.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>-0.356675</td>\n",
              "      <td>0</td>\n",
              "      <td>1</td>\n",
              "      <td>-0.847313</td>\n",
              "      <td>-0.831616</td>\n",
              "      <td>0.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>0.182322</td>\n",
              "      <td>0</td>\n",
              "      <td>1</td>\n",
              "      <td>-0.847313</td>\n",
              "      <td>-3.089998</td>\n",
              "      <td>-1.0</td>\n",
              "      <td>0</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   log_radon  floor  county  log_uranium_ppm    log_u2  floor2  floor3\n",
              "0   0.788457      1       0        -2.000000 -0.217612     1.0       1\n",
              "1   0.788457      0       0        -0.689048 -1.880023     0.0       0\n",
              "2   1.064711      0       0        -0.689048  0.743659    -1.0       0\n",
              "3   0.000000      0       0        -0.689048 -1.001699     0.0       0\n",
              "4   1.131402      0       1        -0.847313 -1.567902     1.0       0\n",
              "5   0.916291      0       1        -0.847313  0.039850     0.0       0\n",
              "6   0.405465      0       1        -0.847313  0.012276     0.0       0\n",
              "7   0.000000      0       1        -0.847313 -1.483836     1.0       0\n",
              "8  -0.356675      0       1        -0.847313 -0.831616     0.0       0\n",
              "9   0.182322      0       1        -0.847313 -3.089998    -1.0       0"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 9
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "_OxaVNnjYZyL"
      },
      "cell_type": "markdown",
      "source": [
        "### Specializing the GLMM Family"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "niha5M54Yjf-"
      },
      "cell_type": "markdown",
      "source": [
        "In this section, we specialize the GLMM family to the task of predicting radon levels. To do this, we first consider the fixed-effect special case of a GLMM.\n",
        "\n",
        "Using R notation, we might consider the following GLM:\n",
        "\n",
        "    log(radon) ~ 1 + floor\n",
        "\n",
        "This model posits that the `radon` response is governed by the floor of a building (e.g., \"ground floor\" of a \"two story\" home).  More concretely it states that the `log(radon)` reading is explainable (in expectation) by the formula `offset + weight_floor[floor]`, i.e., there's a weight learned for every floor and a universal `intercept` term. In R's \"tilde notation\", the `1` indicates there's a weight associated with the \"null feature\", i.e., a weight associated with every sample.\n",
        "\n",
        "Given our data, seems like it might be a good start. I.e.,"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "YwzykNvJgfJo",
        "outputId": "a5c7372c-2d64-47ed-c117-62eb12e7a551",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 720
        }
      },
      "cell_type": "code",
      "source": [
        "df['log_radon'].plot(kind='density');\n",
        "plt.xlabel('log(radon)')\n",
        "plt.figure()\n",
        "df['floor'].value_counts().plot(kind='bar')\n",
        "plt.xlabel('Floor')\n",
        "plt.ylabel('Count')\n"
      ],
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0,0.5,u'Count')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 10
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAFYCAYAAAB6RnQAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VeWB7vHfu3MPEEggF5IAIQQS\nghBAQTzF1nrpdDqnR+fUqh06U45oh+nYnurUU62oUA9eznQcO606rT2FOta22su0pzPTjhVRW4oi\n13APSQiQEAIhG0ICue33/LGSDZEEk7D3Xnvt/Xw/+mHv7GSvh/VJeLLWetf7GmutRURERDzD53YA\nERERGR6Vt4iIiMeovEVERDxG5S0iIuIxKm8RERGPUXmLiIh4jMpbRETEYxLdDjBcDQ0NbkdwVX5+\nftzvA9B+6KP9oH3QR/vBEWv7IT8/f8CP68hbRETEY1TeIiIiHqPyFhER8RiVt4iIiMeovEVERDxG\n5S0iIuIxKm8RERGPUXmLiIh4jMpbRETEY1TeIiIiHuO56VFFxJtsTw9sfwd7qgUzaz4mZ6LbkUQ8\nS+UtImFnW5oJ/NPX4Uit89znw3zui/j+yw0uJxPxJp02F5GwsufaCfzjI3CkFrPoo5i/WA6p6di1\n38LW7HM7nognqbxFJKzsj78HRw9jbvgk5s4v4/voJ/B94WtgAwRe/DY2EHA7oojnqLxFJGzsgT3Y\nP/wOCqdibv0fGGMAMKVXYBZ9FOrrYNs7LqcU8R6Vt4iEhbWWwI9fAMC35K8xif2H2JhPfBqAwOu/\ning2Ea9TeYtIeOzaCnUH4Mr/gikpv+hlM7EQSmfD/l3YpgYXAop4l8pbRMIi8B+vAuD7xG2Dfo5Z\nfCMAduObEckkEitU3iIScvbAbti/C2ZfhZlcPOjnmTkLISERu2NTBNOJeJ/KW0RCzr7+awB8H//U\nJT/PpI+CGbOg7gD25IlIRBOJCSpvEQkpe9qP3boRCqbA9Iuvdb+fqbja+TodfYsMmcpbRELKbngd\nerox1/5J8NawSzEVC5yv27k53NFEYobKW0RCxgYC2Ld+C0nJmEXXDelrzIRcyM6Dql3YQE94A4rE\nCJW3iIRO9V443oi58kOYUaOH/GVmxhXQ3gZHDoYvm0gMUXmLSMjYd98CwFz9keF94YwrnK/fvzPU\nkURikspbRELC9vRgN/8BxoyFmRXD+lpT2lve+3aFI5pIzFF5i0ho7NkOraecU+YJCcP6UjM+B8bn\nwIFdWGvDFFAkdqi8RSQkgqfMF354RF9vikvhTCucOBbKWCIxSeUtIpfNdnVit22ErAkwrWxkbzJ1\nhvNetftDmEwkNqm8ReTy7d0BZ9udU+a+kf2zYqZOdx6ovEU+kMpbRC6b7V2T28xdNPI3mTQNfD4d\neYsMgcpbRC6LDQSw29+F0RlQMsJT5oBJSXGmVD1Ug+3uDmFCkdij8haRy3OwCk61YOYswPiGN8r8\n/czUGdDVCQ11IQonEptU3iJyWez2dwEwc6++/Dcrcq572xqdOhe5FJW3iFwWu+0dSEqG8rmX/V6m\nt7w5VH3Z7yUSy1TeIjJi9ngjNByCmRWYlNTLf8OJhZCYiD1ce/nvJRLDVN4iMmJ21xYAzOwrQ/J+\nJjEJ8idDfR22RyuMiQxG5S0iI2Z3bQPAzJofsvc0k6Y6g9Ya60P2niKxRuUtIiNiu7th73bImYjJ\nzgvdG08qdt7/cE3o3lMkxqi8RWRkavbBubOY8nkhfVvTW97ourfIoFTeIjIidvdWAMys0JY3hUXO\n++vIW2RQKm8RGRG7ayskJEDp7JC+r0kfBdl5cLhWy4OKDELlLSLDZs+chroDMK0Mk5Ye+g1Mmgpn\nToP/ZOjfWyQGqLxFZNjsnu1gbcivd/cJXvc+pFPnIgNReYvI8PXd3x3q6929jEaci1ySyltEhsVa\n69zfPXoMTJ4Wno1MmupsSyPORQak8haR4Wk4DP5mzMy5GF+Y/gnJHO/8cqAjb5EBqbxFZFj6bhEj\nhLOqvZ8xxpms5Xgj9mx72LYj4lUqbxEZluB85rMufxWxSzG9p845cjCs2xHxIpW3iAyZ7eyA/bug\nYApm3Pjwbix43VunzkXeLzGcb7527VqqqqowxrB06VJKSkou+pyXX36Z/fv3s3LlynBGEZFQOLAb\nujrDNsr8QmZSMRY0TarIAMJ25L17924aGxtZvXo1y5cvZ82aNRd9zpEjR9izZ0+4IohIiNldvVOi\nhun+7n5yCyAxSSPORQYQtvKurKxkwYIFABQWFtLW1kZ7e/+BJy+++CJ33HFHuCKISIjZXVshKRmm\nl4d9WyYxEQqmOGt7d3eHfXsiXhK20+Z+v5/i4uLg84yMDPx+P+npzlSK69evp7y8nOzs7GG9b35+\nfkhzepH2gUP7wRGp/dB9oomj9XWkzr+G7KKpEdnmydJZtNUdIDvQSXL+5EE/T98LDu0HRzzsh7Be\n877QhQsMnDlzhjfeeIOHH36YkyeHN3dxQ0NDqKN5Sn5+ftzvA9B+6BPJ/RD4w+sAdJSUR26b43MB\naNryDr7kgedQ1/eCQ/vBEWv7YbBfRMJ22jwzMxO/3x983tLSQmZmJgA7d+7k9OnTPProo3zjG9+g\ntraWtWvXhiuKiIRC3y1ikbje3Utre4sMLGxH3hUVFbzyyivcdNNN1NTUkJmZSVpaGgCLFi1i0aJF\nADQ1NfHcc8+xdOnScEURkctkAz3YPdtg3HjInxS5DQfX9lZ5i1wobOVdWlpKcXExK1aswBjDsmXL\nWL9+Penp6SxcuDBcmxWRcKirgTOtmA/d6Mx+FiEmLb13be8arLUR3bZINAvrNe8lS5b0e15UVHTR\n5+Tk5Ogeb5Eo1zerWjinRB3UpGLYsgFamiFrQuS3LxKFNMOaiHwgu2srGIMpr4j4toPTpOrUuUiQ\nyltELsm2t0HNXiiajhk1JuLbP7+2d3XEty0SrVTeInJpe3dAIBCRKVEHpLW9RS6i8haRS+pbAtS1\n8g6u7a3yFumj8haRQVlrsTu3QFo6TC11JUO/tb3b21zJIBJtVN4iMriGQ9DchJk1H5OQ4FoMre0t\n0p/KW0QGZbe94zyocHluBl33FulH5S0ig7Lb3wWfDzP7SldznJ8mtcbVHCLRQuUtIgOyp1rgYBWU\nlLtyi1g/WttbpB+Vt4gMyO7YBNZi3D5lzgVrezfUYbu73I4j4jqVt4gMyO7YBBAV5Q1gpkyD7m6o\nP+R2FBHXqbxF5CK2swN2b4W8QkzuwOsJR1zRdADswSqXg4i4T+UtIhfbuwM6O6PmqBvATHXKG5W3\niMpbRC5mt78LRM8pcwAmTobkZB15i6DyFpH3sYEAdvsmZ0rSae7MqjYQk5AAk6dBwyFsR4fbcURc\npfIWkf4OVcOpk5jZV2F87s2qNhBTNB0CAdAKYxLnVN4i0s/5U+ZXu5xkAH2D1mp16lzim8pbRPqx\n29+FxESYNdftKBcxRRq0JgIqbxG5gG0+7iy9WTobk5rudpyL5UyE9FEatCZxT+UtIkHRNjHL+xlj\nnFPnTUexbWfcjiPiGpW3iAQFy3vOApeTDC546rxOR98Sv1TeIgKA7TjnTM5SWIQZn+N2nEEZDVoT\nUXmLSK8926G7K6qPugGYOgMAW7vf5SAi7lF5iwjgjVPmAGZcFmRNgJp9WGvdjiPiCpW3iDizqu14\nD0ZnQN8c4tFs6gxoPQXNTW4nEXGFyltE4HBN76xqV0bdrGoDMVOdaVt16lzilcpbRJy5zIneW8Te\nzxT3zrleo/KW+KTyFhHszs2QkAAzo29WtQFNngY+H7Z2n9tJRFyh8haJc7btDBw8AMWlmPRRbscZ\nEpOSAoVFUFeN7epyO45IxKm8ReLdvkqwAUxZhdtJhsVMnQHdXXRpqlSJQypvkThn924HwJR7q7zp\nHbTWsW+ny0FEIk/lLRLn7J7tkJIGRTPcjjIsptjJ26nyljik8haJY/bkCWishxmzMImJbscZntwC\nSBul8pa4pPIWiWPBU+YzPXbKHDA+H0ydTnf9IWxbq9txRCJK5S0Sz/Z4t7yhd9AagCZrkTij8haJ\nY3bfThgzFgqmuB1lRIIzrWmyFokzKm+ROGWbj0PLCZg2E2OM23FGplgrjEl8UnmLxCl7YDcApmSm\ny0lGzowZS0JuARzcrxXGJK6ovEXiVfUewNvlDZBSdgWcaYXjR92OIhIxKm+ROGWr90JikjNPuIcl\nl14BgK3RPOcSP1TeInHInmuHwwehaDomKcntOJelr7y1wpjEE5W3SDyq2e/MZ+7xU+YAycUzICFR\ng9Ykrqi8ReJQ3ylmM63M5SSXzySnwKSpcLgW29XpdhyRiFB5i8QhW3fAeVBU4m6QEDHFpdDTDYdq\n3I4iEhEqb5F4VFcNY7Mw48a7nSQ0pup+b4kvKm+ROGNPtTiTs0zx9ijzC/WtMIZGnEucUHmLxJve\nU+ZmSmycMgcgeyKMHqMjb4kbKm+ROGPrqoHYKm9jjLMe+Ylj2NZTbscRCbuwLuC7du1aqqqqMMaw\ndOlSSkrO/2Pxu9/9jjfeeAOfz8eUKVNYtmyZd+dXFvGQ4GC1GDptDs4KY3bnZuc2uIoFbscRCauw\nHXnv3r2bxsZGVq9ezfLly1mzZk3wtY6ODjZs2MCqVat47LHHqK+vZ/9+ne4SiYi6AzAuCzMuy+0k\nIWWKe1cYq9V1b4l9YSvvyspKFixwfvstLCykra2N9vZ2AFJSUnjkkUdITEyko6OD9vZ2xo0bF64o\nItLL+k+C/yTE0CnzoL4R5xq0JnEgbOXt9/vJyMgIPs/IyMDv9/f7nH/913/li1/8Itdccw25ubnh\niiIifWLwencfM2o05BbAwSpsIOB2HJGwCus17wsNtFzfLbfcwic+8QmeeOIJysrKKCv74Nme8vPz\nwxHPU7QPHNoPjuHsh1PrmzgNjJ+3kLQY2n99+6B5VgXt6/6dnEAXSYVTXU4VefqZcMTDfghbeWdm\nZvY70m5paSEzMxOAM2fOcOjQIcrLy0lOTmbu3Lns27dvSOXd0NAQrsiekJ+fH/f7ALQf+gx3P/Ts\nqQTgZPoYTIzsvwv3QSBvEgDH3vk9vsQUN2NFnH4mHLG2Hwb7RSRsp80rKirYuHEjADU1NWRmZpKW\nlgZAd3c3zz33HOfOnQPgwIEDcfGbkojr6usgfRRkTnA7SVj0DVpDg9YkxoXtyLu0tJTi4mJWrFiB\nMYZly5axfv160tPTWbhwIbfeeiurVq0K3ip21VVXhSuKiAC2owOajkJJWezelllQBEnJmqxFYl5Y\nr3kvWbKk3/OioqLg4+uuu47rrrsunJsXkQsdPeQsA1pY5HaSsDGJiTC5GGr3Yzs6MCnxdepc4odm\nWBOJE7a+znlQUORqjnAzU2dAIACHtcKYxC6Vt0i8OOKUdywfeQPO2t6APVLrchCR8FF5i8QJW3/Q\neVAw2dUc4Wb6bhE7fNDVHCLhpPIWiRdHDsKEXExquttJwmviJEhI0JG3xDSVt0gcsKdboPUUFExx\nO0rYmaQkyCuE+jpsoMftOCJhofIWiQfxcr27l5lUDB3noKnR7SgiYaHyFokD8TLSPGhSkfOnTp1L\njFJ5i8SDo4cBMPmxPVitT9+gNatBaxKjVN4iccAePQw+H+ROdDtKZOh2MYlxKm+RGGethYbDkD0R\nk5jkdpyIMGPGwtgsOKzyltik8haJda1+aD/j3EIVTyZNhZYT2LZWt5OIhJzKWyTWHT0CgJlY6HKQ\nyAqOrO8daS8SS1TeIjHO9g5Wi7sj797BebbhkMtBREJP5S0S6+L1yLtvZL3KW2KQylskxgWPvPPi\nq7zJKwRjdOQtMWlI5f3DH/6Qo0ePhjuLiITD0cOQNQGTmuZ2kogyKSkwITd4j7tILEkcyieNHj2a\nf/qnfyIlJYXrr7+eRYsWkZycHO5sInKZbHsb+E9C+Ty3o7gjfzJsfxfbesq5fUwkRgypvG+++WZu\nvvlmjh07xh//+Ecee+wxpkyZwp/+6Z9SUFAQ7owiMlKNvde78+NssFovkz8Ju/1d57p36Wy344iE\nzLCueTc3N9PY2MjZs2dJTU3l2Wef5T//8z/DlU1ELpPtHaxGnA1WC9KIc4lRQzryfvXVV3n77beZ\nOHEiN910E5///Ofx+Xx0d3fz4IMP8rGPfSzcOUVkJI46pWXy4vTIe+JkLGjEucScIZX3qVOnePjh\nh8nOzg5+rKmpiZycHJYsWRK2cCJyec4fecdneZ8fca5BaxJbPvC0eSAQoL6+ngkTJhAIBAgEAnR1\ndfHUU08BMHfu3LCHFJEROnoYRmdgxmS4ncQVwRHnOvKWGHPJI+/f//73vPrqqzQ2NnLHHXcEP26M\noaKiIuzhRGTkbFcnnGiCkjK3o7irb8T5aT8mY5zbaURC4pLlvXjxYhYvXswrr7zCbbfdFqlMIhIK\nx+rBBjDxesq8V78R5ypviRGXLO+tW7cyb948JkyYwLp16y56/frrrw9bMBG5PMHrvHFe3sER50cP\nY8rmuBxGJDQuWd51dXXMmzePvXv3Dvi6ylskigXnNI/v8jb5GnEuseeS5X3LLbcA8IUvfAFrLcYY\nurq6OHXqFBMmTIhIQBEZoeBqYnF6j3efXGciqeDIe5EYMKRbxX7xi1+QmprKDTfcwFe/+lXS0tKo\nqKjg9ttvD3c+ERkh23gEUtIgM75/0TYpqZA1wRkDIBIjhjTD2ubNm/n4xz/Ohg0buPLKK3n88ccH\nPZUuIu6zgR6nrPIKMMa4Hcd9eYXgP4k91+52EpGQGFJ5JyQkYIxh27ZtLFy4EHDu/xaRKHXiGHR3\nx90a3oMxvafOOdbgbhCREBlSeY8aNYonnniC+vp6ZsyYwebNm/H5tBS4SNQ62nuKON7W8B5MXu91\n70adOpfYMKRr3l/60pfYsWMHZWXOZA9JSUn87d/+bViDicjI2UZnsJqOvB0mr8AZca7ylhgxpPL2\n+XwYY9i8eTPWWgBOnDihW8VEotVR3ePdT27vLzEatCYxYkjlvXr1anw+X7+FSUD3eYtEK9tYDz4f\nZOe5HSU6ZI6H5GTdLiYxY0jl3dPTw6pVq8KdRURCwFrrHHnnTMQkJrkdJyoYnw9yCqCpHhsIOM9F\nPGxI38GFhYW0traGO4uIhEKrH9rbIE7X8B6MySuAzk5oaXY7ishlG9KR98mTJ/nSl75EQUEBCQkJ\nwY/raFwkCgWnRS1wOUiUyeu7XewIjM++9OeKRLkhlXffNKkiEv2C13V15N1f7vnbxUz5PJfDiFye\nIZ02Ly8v59y5cxw6dIjy8nKysrKYOXNmuLOJyEg09h156zaxCwX3h24XkxgwpPJ+6aWXWLduHW+8\n8QYAv//971mzZk1Yg4nIyJw/8lZ595ObD4DV7WISA4ZU3nv27OErX/kK6enpANx6663U1taGNZiI\njFDjYRg3HpOW7naSqGJS02Fclo68JSYMqbyTk5P7PQ8EAvT09IQlkIiMnD13Fk6e0DKgg8ktgJPH\nsR0dbicRuSxDGrA2Y8YMnn32WVpaWvj1r3/NO++8Q3l5ebizichw9Z4SNjplPiCTV4DdVwlNDTBp\nqttxREZsSOX94Q9/mLq6Oqqrq9m7dy+f/OQng6uLiUj0sMFpUVXeA7pggRKj8hYPu2R5d3Z28s1v\nfpO6ujqKi4vJyspi7969JCcnM3/+fBITh9T9IhIpR3XkfSkmt9BZoOSYpkkVb7tk+/70pz8lKyuL\n++67Lzg5S0dHB2vWrOHll1/mr/7qryISUkSGpm81MR15DyJP63pLbLjkgLW9e/fyuc99rt+saikp\nKdx1113s2LEj7OFEZJiOHoG0dBib5XaS6DQ+GxISsSpv8bhLlndCQsKAp8YTExMZNWpU2EKJyPDZ\nnh5oOgp5hRhj3I4TlYwvAXImwrH64PLGIl404qV1fFqVRyS6HG+Enm5d7/4gufnOwi2tp9xOIjJi\nl7zmvX//fv7mb/5mwNdOnz79gW++du1aqqqqMMawdOlSSkpKgq/t3LmTH/3oR/h8PiZOnMjy5cv1\nC4HI5Qhe79ac5pdicgt6B601QMY4t+OIjMgly/uZZ54Z8Rvv3r2bxsZGVq9ezZEjR3j++edZvXp1\n8PXvfve7PProo4wfP56nn36abdu2MX/+/BFvTyTe2QanvDWn+Qe4YJpUM13zVYg3XbK8s7NHvmxe\nZWUlCxYsAJz1wNva2mhvbw9Osfrkk08GH2dkZHDmzJkRb0tEgPpDzp8FU9zNEeX6HXmLeFTYbtT2\n+/0UFxcHn2dkZOD3+4OF3fdnS0sL27dv5/bbbx/S++bn54c+rMdoHzi0Hxx9+6Hx+FG6U1LJnz0X\nE2eXoIbzvdCTlkIDkHqqmQkx9j2knwlHPOyHiM2yMtDIzlOnTvHUU09x1113MWbMmCG9T0NDfP+2\nnJ+fH/f7ALQf+vTtB9vTQ+BwDRQUcbSx0e1YETXc7wVrLaSlc/ZgdUx9D+lnwhFr+2GwX0TC9ut5\nZmYmfr8/+LylpYXMzMzg8/b2dh5//HHuuOMOKioqwhVDJD40HYXuboxOmX8gY4yzQMnxo9iAFlgS\nbwpbeVdUVLBx40YAampqyMzMJC0tLfj6iy++yJ/92Z8xd+7ccEUQiR8Nvde78ye7m8MjTG4+dHdD\n83G3o4iMSNhOm5eWllJcXMyKFSswxrBs2TLWr19Peno6FRUVvPXWWzQ2NrJu3ToAFi9ezI033hiu\nOCIxzdbXAWAKVN5DknvBNKnZee5mERmBsF7zXrJkSb/nRUVFwccvv/xyODctEldsg1Pe5Ou0+ZBc\neLvYFbpFVbwnvoakisSqhsPOnOaZ491O4gkmeORd724QkRFSeYt4nO3qckoof7LmNB+q4JF37IxK\nlvii8hbxumP1EAhopPkwmNQ0GJeliVrEs1TeIh7XN1hNI82HKbcATh7Hdna4nURk2FTeIl7Xe5uY\nUXkPi8nNB2ude+RFPEblLeJxtkFzmo/IhbeLiXiMylvE6+rrYHQGRstbDkvfiHOrEefiQSpvEQ8L\ntJ2B440waarbUbynd8S5jrzFi1TeIh7WVVsFgFF5D9+EXPD5dOQtnqTyFvGwztr9zgOV97CZxESY\nkKeJWsSTVN4iHtZV45S3mVTschKPyiuAM63YM6fdTiIyLCpvEQ/rrNkPiUnnR07LsBhd9xaPUnmL\neJTt7qarrhoKpjingGX4NOJcPErlLeJVx+qhqxNTWOR2Es/Skbd4lcpbxKPs4Vrnga53j5yOvMWj\nVN4iXtVb3rpN7DKMy4KUVB15i+eovEU8yh6ucR7otPmIGWOcyVqaGrCBgNtxRIZM5S3iQdZaOFxL\nQl4BJn2U23E8zeQWQGcntDS7HUVkyFTeIl7kPwlnTpM8dYbbSbwvOGhN173FO1TeIl50yDllnlSs\n8r5sveVtdd1bPETlLeJB9qAzp3ny9JkuJ/E+k1voPNCRt3iIylvEg2zdAQCSp5e7nCQG5E4EdLuY\neIvKW8RjrLVwsAqyskkYl+V2HM8z6aNhzFjdLiaeovIW8ZqTx6H1FBRNdztJ7MgtgBNN2K4ut5OI\nDInKW8Rreq93G5V3yJi8ArABONHodhSRIVF5i3iMPehc7zZFJe4GiSW6XUw8RuUt4jF9I82ZovIO\nFdM3x3mjylu8QeUt4iE2EIC6A5CrmdVCSquLiceovEW8pKkBzrajU+Yhlj0RjE+3i4lnqLxFPKTv\nerdGmoeWSUqCCTk68hbPUHmLeIlGmodPbj6c9mPb29xOIvKBVN4iHmJr94PPB5OK3Y4Sc/oGreno\nW7xA5S3iEbarEw5Vw6RiTEqK23FiT3CBEl33luin8hbxirpq6O7GTCtzO0lMMhpxLh6i8hbxCFuz\n13mg8g4PrS4mHqLyFvEIW+2Ut5mmZUDDInM8JCdjG4+4nUTkA6m8RTzAWgvVe2FcFmRNcDtOTDI+\nH+RNgqNHsIEet+OIXJLKW8QLmpvgVAtMK8MY43aamGUKJkN3FzRpgRKJbipvEQ/QKfMIyZ/s/NlQ\n524OkQ+g8hbxguo9AJjiUpeDxDZTMAUA23DI5SQil6byFvEAW70PEpNg8jS3o8S2viPvepW3RDeV\nt0iUsx3n4EgtFJU4c3BL+GRlQ0qajrwl6qm8RaJd9V4IBDQ5SwQYY6BgMhyrx3Z3uR1HZFAqb5Eo\nZ6t2AWBmXOFykvhg8idDT49mWpOopvIWiXJ2/04wBko00jwiCpzr3jp1LtFM5S0SxWxXJ9Tsh0lT\nMemj3Y4TF0y+M+Kcet0uJtFL5S0SzWr3Q3eXTplHUr6OvCX6qbxFopjdvxMAM32Wy0niyNhMGD0G\njhx0O4nIoBLD+eZr166lqqoKYwxLly6lpKQk+FpnZyff/e53OXLkCE8++WQ4Y4h4lt3vDFZD5R0x\nxhiYVAx7tmPb2zDpo9yOJHKRsB157969m8bGRlavXs3y5ctZs2ZNv9dfeuklioqKwrV5Ec+z3V3O\nzGoFUzBjMtyOE1fMpGLnweFad4OIDCJs5V1ZWcmCBQsAKCwspK2tjfb29uDrn/nMZ1i4cGG4Ni/i\nfXXV0NmpU+ZumOyUtz1c7XIQkYGFrbz9fj8ZGeePFjIyMvD7/cHnaWlp4dq0SEzou96NBqtFnOmb\nhvZQjbtBRAYR1mveF7LWhuR98vPzQ/I+XqZ94Ij1/dBUu48OYOK115NwiTW8Y30/DEWo94HNzaU+\nNY3EhkPkeWj/6nvBEQ/7IWzlnZmZ2e9Iu6WlhczMzMt+34aG+J71KD8/P+73AcT+frBdnQR2boWC\nKRw71wmD/F1jfT8MRbj2gS2YQlftfuoP1mKSU0L+/qGm7wVHrO2HwX4RCdtp84qKCjZu3AhATU0N\nmZmZOlUuMlQH9kBXJ2bmXLeTxC0zuRgCAU3WIlEpbEfepaWlFBcXs2LFCowxLFu2jPXr15Oens7C\nhQt5+umnaW5upqGhgZUrV3LjjTeyePHicMUR8RS7ZxsAplzl7Zre6972UA1m6gyXw4j0F9Zr3kuW\nLOn3/MJbw+67775wblrE0+zu7ZCQCNPL3Y4St8zkYixo0JpEJc2wJhJl7JnTcKgappViUnWpyTX5\nkyEhEVt3wO0kIhdReYtEm33TqcDqAAAUkElEQVSVYC1mZoXbSeKaSUyCwiI4ctBZIEYkiqi8RaKM\n3bUVQIPVooCZVgY93c6EOSJRROUtEkWstdjK92DUGJg63e04Mq0MAFuz1+UgIv2pvEWiyeEa8J/E\nzL4S40twO03cM33lXb3P5SQi/am8RaKI3bHJeTBngbtBxJGVDWOzoGZvyGaJFAkFlbdIFLE73oOE\nBMyseW5HEXqXB51WCv6TcPK423FEglTeIlHCnm6B2v1QUo5JH+12HOllivtOneu6t0QPlbdIlLCV\nmwEwOmUeVUzJTOdB1S53g4hcQOUtEiXsVmctAJV3lJlSAilp2L073E4iEqTyFokCtv0M7NoChUWY\nvAK348gFTGIizJgFjfXYlma344gAKm+RqGC3vgPd3ZirtDhPNDJlcwCwe7a7nETEofIWiQL2vbcB\nMAuudTmJDCQ4Ve1elbdEB5W3iMvsmdOwZztMKcHkTHQ7jgykYAqMzsDu2aH7vSUqqLxFXGbfeRN6\nejALddQdrYzP55w69zdDw2G344iovEXcZK3FvvVbSEjELPqo23HkUiqcuwDsto0uBxFReYu4q3oP\nNBzCzL8GkzHO7TRyCWb2AvD5sNvecTuKiMpbxE32rd8CYK79mMtJ5IOYUaOhdDYcrNItY+I6lbeI\nS6z/JHbT25BbAL23Ikl0M3OvBsBu19G3uEvlLeIS+9ovnXu7P3aLswCGRL1geW/e4HISiXcqbxEX\n2NZT2Df/A8ZlYa653u04MkQmKxtKZsK+SqxWGRMXqbxFXGB/9SPoOIf5+K2YpCS348gwmGuuB2ux\nf3zD7SgSx1TeIhFm6w9h3/oN5BZgPvJxt+PIMJmrFkNSMnbDOmwg4HYciVMqb5EIsj09BNZ+EwIB\nfJ++01n0QjzFpI/CXPkhaGqA3dvcjiNxSuUtEkH2334CB6swi67DVGjpT68yN/xXAAKv/z+Xk0i8\nUnmLRIjdthH7/34M43Mwd9ztdhy5DKZoujNwbedmbH2d23EkDqm8RSLAbt1I4Dv/B5KS8X3hQcyo\nMW5Hksvk+/itAAR++UOXk0g8UnmLhFlgwzoC//wkJCTiu2cFZvI0tyNJKMy5CqaVwdaN2Noqt9NI\nnFF5i4RR4PVfY9c8A6np+O79OqZ8rtuRJESMMfhu+SwAgZ+t1VKhElEqb5EwsNYS+PVPsD/+LozN\nxHf/45hpZW7HkhAzZXNgzgJn0pY/rnM7jsQRlbdIiFlrsa9+H/vLH8L4HHz/6wlMYZHbsSRMfH+x\nHFLSsK98H3va73YciRMqb5EQsoEA9ofPO/OWT5yE76tPYXLy3Y4lYWTGZ2P+/LPQ1kpgzTOauEUi\nQuUtEiLB4n7zNzBpKr77n8Bkjnc7lkSA+eifwRXzYecW7H/81O04EgdU3iIhYK3FvvzPzvrck4vx\n3fcYZkyG27EkQozPh+/O+yBzAvaXP8Ru2+h2JIlxKm+RELD/+QvniLtwqlPco1Xc8caMycD3hQch\nKZnAC9/A1u53O5LEMJW3yGWyWzdif/YDGDce3/98RBOwxDFTNB3f5++Hrm4C33oMe/SI25EkRqm8\nRS6Drasm8L1/cGZO++IKzDhd4453pmIhZslyaD1F4B9WYI81uB1JYpDKW2SEbEszgW8/Bl2d+O76\nO82cJkG+j3wcc/tdcOqkU+DHG92OJDFG5S0yArbjHIFv/2/wn8R86nOYeYvcjiRRxnfjf8Pc+j+g\n5QSBb3wN26hT6BI6Km+RYbKBAIH/+zQcqsYsvgnzsT93O5JEKd+f/Dnm1qVw8gSB//Mgtq7a7UgS\nI1TeIsNkf/4D2LoRSmdjlizHGON2JIlivj/575jPfgHOnCbwDw9h9+90O5LEAJW3yDAE3vwN9re/\ngLwCfH/zACYxye1I4gG+j3wcc/dXoLODwDMrCbz7ltuRxONU3iJDZCvfw778zzA6A98XdUuYDI9v\nwbX47lkBCQnYF75B4Bf/oqlUZcRU3iJDYHduIfDcE+fX5M6Z6HYk8SBzxZX4Hvx7yM7D/vurBL65\nEutvdjuWeJDKW+QD2Mr3CDy7GozB97cPaWlPuSwmfzK+h/4BrrgSdm8jsPJLBDas01G4DIvKW2QQ\n1loCr/2SwLf+9/ninjXP7VgSA8yoMfi+9AjmL5ZDVwd2zTMEnvxf2D3bsda6HU88INHtACK2pweO\nH4UzrYCFURmQnevqYDDbfJzAS8/Czi0wNtMp7qkzXMsjsccYg/noJ7BzrsL+7AfYTW8TePphKJyK\n+fCfYOYtwozLcjumRCmVt7jC+pux77yF3f4O1B2Azs7+n5CYCJOKMTMrMHMXwZRpGF94TxRZa+FQ\nDfat32I3/A66u6F8Lr6l/1NLe0rYmPE5mM/fj73pFuxr/4rd/AdnhboffQeKSzEzZmGmlUNJmQZJ\nSpDKWyLGnmvHbvkjduN62LsDrAXjg4LJztSiY8c5n3j6FPbIQadIa/dj//1VGJeFmbMQM2cBzJwz\n9G12d8EpP5xugXNnobsLujqxXV3Q2QHtZ+BMK/boYThUAy0nnC/MzsN88jOYRdfpPm6JCDN1ulPi\nLXdit2zAbtkAVXuw1Xux/Mz5pJyJmCklzi+zk6c5f6aPdje4uMLYMF5gWbt2LVVVVRhjWLp0KSUl\nJcHXduzYwY9+9CN8Ph/z5s3j1ltvHdJ7NjTE9yT/+fn5ntoHtrMDKjdjN72Nrdx0/gh7WplTjFct\nHnT5TNtxDnZtxW57B7tjE7S1Oi8kJpI0dQbdOfkwNgvSR0Ggxylj/0lsywnwn4RTJ3tPxQ9RxjhM\n2RzMlR+CioWYhITL/NuHn9e+H8IhlveBPdcONfuwVXuw1Xugrtr5hfNC2XmYoulkzJ5H67jsuC/0\nWPt+yM/PH/DjYTvy3r17N42NjaxevZojR47w/PPPs3r16uDra9as4aGHHiIrK4uVK1eyaNEiCgsL\nwxVHIsB2d0NzEzQewR6qwVbtguq9TqkC5OQ7hX31R4Z0q5VJSYX512DmX+NcF6/Zh92xCbtnO121\nVVC1e/AvThsF47Kc64djM2FsJqSmQ1ISJCVDovOnGTUGRo1ysmWMC9GeEAkNk5oO5fMw5c5ASWst\nnDgGdQewddXYugNQV43d9DanNr19/guz884foU8pgQm5zi+nKaku/U0k1MJW3pWVlSxYsACAwsJC\n2traaG9vJz09nWPHjjF69GgmTJgAwLx586isrIxYedtTLdhfvIg9d7b3A/1eveDh+05KDPXzBvr4\npU5w9HttkK/pfXg8NYWec+eGt+0hf96ltz3AE+jqco6I287A2baL3z9/MqZiAeaqa2HS1BGfgjYJ\nCTC9HDO9HICJOTk0bN8Mraed7SYkQFKKU9KZ4/WPlMQkYwxk5znlfNVi4HyhZ54+ycltm84X+nu/\nh/d+3/8nNiUVxoyF5BTnl9jkZOfnpu8skzHO/xc+NgYwvf9F/yWkE2lpBM6edWfj47Iwt92J8YX/\nrF3Yytvv91NcXBx8npGRgd/vJz09Hb/fT0bG+VOlY8eOpbFxaEvmDXYKYTg6Wppo+uN651Srx5z7\n4E8JnQuL9qLSdZ6bpER8Y8bhy8vHNzqDhJyJJBVOIWnSVJLL55IwNnxHswVXXh229/aSUPxMeF3c\n74OCAgDSr70RcAq951gDnQf20Fm9j57mJgItJ+lpaSZw2o9tPYXtPIft6BjWZrxwE5tLtQ2ASUsn\n7+4vkzBmbNi3FbEBa5e6tD6cy+4huZaRmYPvmR86g5eCzIAPByutiz6PQYpuKJ8z5O0b8vMn0tBw\ndOjvO4TPC8WArEDv/91A8J+Dtnbn/zCItetaI6X9oH3Q5+L9YKC43Pn/fXqPpZ1/e7u7oKcHsE47\nW3v+MRYCfc+9UN2Ql5c35IPBkEtJ41hrG7QOcBZyhCJ+zTszMxO/3x983tLSQmZm5oCvnTx5kqys\nyN7PaNLSI7q9UDFJyZgkLYYhIpfPGOOcPo+hf1ISxmVh2iN6jtIVYbtxtqKigo0bNwJQU1NDZmYm\naWlpAOTk5HD27Fmampro6elhy5YtzJkz9Nt/RERE4lnYjrxLS0spLi5mxYoVGGNYtmwZ69evJz09\nnYULF3LXXXfxzW9+E4BrrrlG16xERESGKKzXvJcsWdLveVFRUfBxeXl5v1vHREREZGi0MImIiIjH\nqLxFREQ8RuUtIiLiMSpvERERj1F5i4iIeIzKW0RExGNU3iIiIh6j8hYREfEYY4ezKoiIiIi4Tkfe\nIiIiHqPyFhER8RiVt4iIiMeovEVERDxG5S0iIuIxKm8RERGPCet63qH2q1/9irfffpvExESWLVtG\nSUmJ25Fc4/f7uffee/nKV77CrFmz3I4TcT09PTz//PMcO3aMQCDAX/7lX1JWVuZ2rIhZu3YtVVVV\nGGNYunRp3P4svPTSS+zZs4dAIMAtt9zC1Vdf7XYkV3R2dvJ3f/d3fOpTn+K6665zO44r3n77bX71\nq1/h8/m4/fbbmT9/vtuRwsozR96HDx9mw4YNPPnkk9x9991s2bLF7Uiueumll8jJyXE7hmveeust\nUlNTeeyxx1i+fDk/+MEP3I4UMbt376axsZHVq1ezfPly1qxZ43YkV+zcuZPDhw+zevVqvva1r7F2\n7Vq3I7nmZz/7GaNHj3Y7hmtaW1v56U9/yte//nUeeOABNm3a5HaksPPMkffmzZu55pprSEhIoLi4\nmOLiYrcjuWbnzp2kpqYyefJkt6O45tprr+VDH/oQABkZGZw5c8blRJFTWVnJggULACgsLKStrY32\n9nbS09NdThZZ5eXlwTMOo0aNoqOjg0AggM/nmWOSkKivr+fIkSPMmzfP7SiuqaysZPbs2aSlpZGW\nlsZf//Vfux0p7DzzXX78+HFOnDjB6tWr+frXv87BgwfdjuSK7u5uXn31VT7zmc+4HcVViYmJJCcn\nA/Bv//ZvwSKPB36/n4yMjODzjIwM/H6/i4nc4fP5SE1NBWDdunXMmzcv7oob4MUXX+Rzn/uc2zFc\n1dTUREdHB0899RSPPPIIlZWVbkcKu6g88n799ddZt25dv4/5/X7mzp3L1772Nfbt28d3vvMdnnji\nCZcSRsZA+2Hu3LnccMMNjBo1yqVUkTfQfvj0pz/N3Llz+c1vfkNtbS1f/epXXUrnvnif4XjTpk2s\nW7eOFStWuB0l4t58801mzJgR15fQ+rS2tnL//fdz/PhxVq1axXPPPYcxxu1YYROV5X3DDTdwww03\n9PvYK6+8Qn5+PsYYysrKaGpqcild5Ay0Hx5++GECgQC//e1vaWxs5MCBA9x3331MmjTJpZThN9B+\nAOdoa/Pmzdx///0kJkblt3JYZGZm9jvSbmlpITMz08VE7tm2bRs///nPeeihh+LusgHAli1baGpq\nYsuWLTQ3N5OUlERWVhZz5sxxO1pEjR07ltLSUhISEsjLyyMtLY3Tp08zduxYt6OFjWf+xZs7dy6v\nvfYaixcvpr6+ngkTJrgdyRWPPfZY8PGzzz7LddddF9PFPZhjx47x2muvsXLlyuDp83hRUVHBK6+8\nwk033URNTQ2ZmZmkpaW5HSvi2tvbeemll3j44YfjdrDWvffeG3z8yiuvkJOTE3fFDc7PxLPPPsvN\nN99MW1sb586dY8yYMW7HCivPlPeMGTPYtm0bDz30EADLli1zOZG46fXXX6e1tbXfpZMVK1bExRF4\naWkpxcXFrFixAmNM3P4sbNiwgdbWVv7xH/8x+LF77rknbn+xj2dZWVksWrQo2A933nlnzI9/0JKg\nIiIiHhPbv5qIiIjEIJW3iIiIx6i8RUREPEblLSIi4jEqbxEREY9ReYt43K5du3j44YdH9LXnzp3j\nkUceobGxcURfv2PHDlauXDnsr+vu7mbVqlXU1dWNaLsi8U7lLRLHXnrpJa699lry8vIiut2+ZX2/\n/e1vEwgEIrptkVgQ+zNaiMSJhoYGXnjhBay19PT0sGTJEsrKyjh27Bjf/va3ASgpKWHr1q088MAD\npKWl8c4777B06VIAVq5cSVFREbW1tTz66KP87ne/48033wwuAvPlL3+ZUaNG8e677/LjH/+Y8ePH\n9yv9wbb/7LPPkpmZyaFDhzh69CjXX389N998M4WFhWRnZ/Pee++xcOFCN3aZiGfpyFskRqxZs4aP\nfexjrFy5krvvvjtY2D/5yU+45ppreOyxx6ioqODo0aOAs4zizJkz+81Kl5qayqpVq/D5fHR2drJi\nxQpWrVpFdnY2b7/9NgDf//73ue+++3jooYf6Lfww2PbBWfXpgQceYMWKFfz85z8PfnzOnDls27Yt\nrPtFJBapvEViRFVVVXBe68mTJ3P27FlOnz5NXV0ds2bNApw1AvqW0Wxubmb8+PH93mPGjBnBx6NH\nj+aJJ57g0UcfZdu2bZw+fZrW1lY6OzspLCwE4IorrvjA7YOz9jZAdnY2Z8+eDZ4qz87O5vjx4yHf\nFyKxTqfNRWLE+5c/tNZijCEQCPR77VLLJPYdhTc3N/Mv//IvPP3004wdO5YXX3yx33v2ufB69WDb\nB0hISLjoNREZOR15i8SI6dOns337dgBqa2sZM2YMY8aMoaCggP379wPO6PCzZ88CMH78eJqbmwd8\nr1OnTpGRkcHYsWM5c+YMO3bsoLu7mzFjxuDz+fqdev+g7V/K8ePHyc7Ovry/uEgc0pG3SIy48847\neeGFF3jttdfo7u7mnnvuAeC2227jW9/6Fn/4wx+YPn0648ePx+fzMXv2bF588UW6u7svWo2tqKiI\nvLw8HnzwQXJzc7ntttt44YUXmD9/PkuXLuXv//7vycnJYeLEiR+4/UuprKzkIx/5SGh3hEgc0Kpi\nIjGuurqarq4uysrK8Pv93HvvvbzwwgskJibyve99j6KiIm688caI56qvr+eZZ57hqaeeivnlG0VC\nTUfeIjEuNTWVtWvXAs7kKHfffXfwSPuzn/0sjz/+OFdccUVE7/Xu7u7me9/7Hvfcc4+KW2QEdOQt\nIiLiMfqVV0RExGNU3iIiIh6j8hYREfEYlbeIiIjHqLxFREQ8RuUtIiLiMf8fCbQxlI/F1UgAAAAA\nSUVORK5CYII=\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fc814cf27d0>"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAFWCAYAAACSHB8oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGvRJREFUeJzt3X9MVff9x/HXgXsRruzKRdrYq6Uq\nXKlOJtRbIm2mHSWzyRLbpVi7mWWs1uja0KwzbUesv2qJM23WmKE1/RG7pjWdY0tKKsvUqJO13rWb\n3UAUKjgzLVPUy7V43QWF+/1j8X5HW4EyD5fP9flImnjP5Z6+iQef93zOvRcrGo1GBQAAjJQU7wEA\nAMDwEXIAAAxGyAEAMBghBwDAYIQcAACDEXIAAAzmsGvHkUhE1dXVCofDunz5ssrKypSRkaHXXntN\nlmUpOztbS5culSTV1tbq4MGDsixLZWVluuOOO+waCwCAhGJbyPfv3y+v16vvf//7CgaDeu655+Tx\neFReXq7c3Fxt2rRJH3/8sSZOnKj3339fVVVVunTpklavXq2CggIlJbFYAADAYGyr5de+9jV1dXVJ\nksLhsNLT09XR0aHc3FxJ0uzZs9XY2KjDhw+rsLBQDodDbrdbN910k06dOmXXWAAAJBTbzsjvvvtu\n7d+/XxUVFQqHw3rmmWf0+uuvx+4fN26cOjs7lZ6eLrfbHdvudrvV2dmp7OzsAfff3t5u1+iwmdfr\n5e8PiAN+9szl9XqveZ9tIT9w4ICysrK0cuVKnThxQi+++KJcLlfs/mt9MuxQPzF2oG8Kox9/f0B8\n8LOXeGwLeUtLi2bNmiVJmjx5snp6etTb2xu7PxgMyuPxKDMzs98zxM7OTnk8nkH3z7NKc3FWAMQH\nP3vmGugJmG3XyCdMmKDW1lZJ0tmzZ5WWlqaJEyequblZkvThhx+qoKBAM2fO1KFDh3TlyhUFg0EF\ng0FNmjTJrrEAAEgoll2//SwSiWjLli26cOGC+vr6tGjRImVkZOiVV15RNBpVbm6ufvjDH0qSfv/7\n3+tPf/qTJOnhhx9Wfn7+oPvnWaW5OCsA4oOfPXMNdEZuW8jtxsFoLv4xAeKDnz1zxWVpHQAA2I+Q\nAwBgMEIOAIDBCDkAAAYj5AAAGIyQAwBgMEIOAIDBbPuIVgxf79IF8R7BVifjPYCNkl+tjfcIAG4w\nnJEDAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAA\nBiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYzGHXjvfu3asD\nBw7Ebre1tWn9+vV67bXXZFmWsrOztXTpUklSbW2tDh48KMuyVFZWpjvuuMOusQAASCi2hbykpEQl\nJSWSpCNHjuiDDz7Qr371K5WXlys3N1ebNm3Sxx9/rIkTJ+r9999XVVWVLl26pNWrV6ugoEBJSSwW\nAAAwmBGpZU1NjR544AF1dHQoNzdXkjR79mw1Njbq8OHDKiwslMPhkNvt1k033aRTp06NxFgAABjP\n9pC3trZq/PjxSkpK0tixY2Pbx40bp87OToVCIbnd7th2t9utzs5Ou8cCACAh2La0ftXevXt1zz33\nfGF7NBr90q+/1vbP83q9/8tYo9rJeA+AYUvk4xKJgWM08dge8qamJj3yyCOSpK6urtj2YDAoj8ej\nzMxMtbe3x7Z3dnbK4/EMut//fgwwWnBcYjTzer0co4Ya6AmYrUvrwWBQqampcjgccjgcmjhxopqb\nmyVJH374oQoKCjRz5kwdOnRIV65cUTAYVDAY1KRJk+wcCwCAhGHrGXkoFNK4ceNit8vLy/XKK68o\nGo0qNzdX3/jGNyRJ9957r9asWSNJevTRR3nFOgAAQ2RFh3pRepRJ5OWh3qUL4j0Chin51dp4jwBc\nE0vr5orb0joAALAXIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcA\nwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAAgxFyAAAMRsgB\nADAYIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAAgxFyAAAMRsgBADCYw86d19fXq7a2VklJ\nSVq0aJGys7NVXV2tvr4+ZWRkqKKiQk6nU/X19aqrq5NlWSotLVVJSYmdYwEAkDBsC3lXV5dqamr0\n85//XJFIRDt27FAgEND8+fNVXFys7du3a9++fZo7d65qamq0YcMGORwOVVZWqqioSOnp6XaNBgBA\nwrBtab2xsVH5+flKS0uTx+PRsmXL1NTUJL/fL0ny+/1qaGhQa2urcnJy5HK5lJKSory8PDU3N9s1\nFgAACcW2M/KOjg51d3dr48aNCofDWrhwobq7u+V0OiVJbrdboVBIoVBIbrc79rir2wfj9XrtGj3u\nTsZ7AAxbIh+XSAwco4nH1mvkXV1deuqpp3T27FmtW7dO0Wj0uu27vb39uu0LuF44LjGaeb1ejlFD\nDfQEzLal9XHjxikvL0/JycmaMGGC0tLSlJaWpp6eHklSMBiUx+ORx+PpdwZ+dTsAABicbSGfNWuW\nDh8+rL6+PnV1dSkSiSg/P1+BQECSFAgEVFBQIJ/Pp7a2NoXDYUUiEbW0tGj69Ol2jQUAQEKxbWk9\nMzNTc+bM0cqVKyVJjzzyiHJyclRdXa09e/YoKytL8+bNk8Ph0OLFi1VVVSXLslRWViaXy2XXWAAA\nJBQrej0vXI+gRL7O07t0QbxHwDAlv1ob7xGAa+Iaubnico0cAADYj5ADAGAwQg4AgMEIOQAABiPk\nAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEI\nOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAw\nQg4AgMEIOQAABnPYteOmpib94he/0K233ipJys7O1oIFC1RdXa2+vj5lZGSooqJCTqdT9fX1qqur\nk2VZKi0tVUlJiV1jAQCQUGwLuSTNmDFDK1asiN3esmWL5s+fr+LiYm3fvl379u3T3LlzVVNTow0b\nNsjhcKiyslJFRUVKT0+3czQAABLCiC6tNzU1ye/3S5L8fr8aGhrU2tqqnJwcuVwupaSkKC8vT83N\nzSM5FgAAxrL1jPzUqVPauHGjLl68qIULF6q7u1tOp1OS5Ha7FQqFFAqF5Ha7Y4+5un0wXq/Xtrnj\n7WS8B8CwJfJxicTAMZp4bAv5LbfcooULF6q4uFhnzpzRunXr1Nvbe932397eft32BVwvHJcYzbxe\nL8eooQZ6Ambb0npmZqbuuusuWZalCRMmKCMjQ+FwWD09PZKkYDAoj8cjj8fT7wz86nYAADA420Je\nX1+v2tpaSVIoFNKFCxd0zz33KBAISJICgYAKCgrk8/nU1tamcDisSCSilpYWTZ8+3a6xAABIKLYt\nrfv9fm3atEl/+ctfdOXKFT366KOaMmWKqqurtWfPHmVlZWnevHlyOBxavHixqqqqZFmWysrK5HK5\n7BoLAICEYkWj0Wi8hxiORL7O07t0QbxHwDAlv1ob7xGAa+Iaubnico0cAADYj5ADAGAwQg4AgMEI\nOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAw\nQg4AgMEIOQAABiPkAAAYjJADAGCwIYX8008//cK2Tz755LoPAwAAvpoBQx4Oh3X69Glt2bJFHR0d\nOnPmjM6cOaP29nZt3rx5pGYEAADX4Bjozk8++UQ7d+7UiRMntG7dutj2pKQkzZo1y/bhAADAwAYM\neWFhoQoLC7Vr1y59+9vfHqmZAADAEA0Y8quKiopUV1enixcvKhqNxrYvWrTItsEAAMDghvRitw0b\nNujEiROyLEtJSUmx/wAAQHwN6Yw8NTVVjz32mN2zAACAr2hIp9U+n+9L34IGAADia0hn5H//+9+1\nc+dOud3ufkvqL7/8sm2DAQCAwQ0p5E8//bTdcwAAgGEYUsgbGxu/dHtJScmAj+vp6dGKFSv04IMP\naubMmaqurlZfX58yMjJUUVEhp9Op+vp61dXVybIslZaWDrpPAADw/4YU8ubm5tifr1y5omPHjikv\nL2/Q6P72t79Venq6JGnHjh2aP3++iouLtX37du3bt09z585VTU2NNmzYIIfDocrKShUVFcUeAwAA\nBjakkH/+Fevd3d3asmXLgI/59NNPderUKRUWFkqSmpqatHTpUkmS3+9XbW2tvF6vcnJy5HK5JEl5\neXlqbm6W3+//yt8IAAA3oiGF/PPGjBmj06dPD/g1b775ppYsWaL9+/dL+k/8nU6nJMntdisUCikU\nCsntdscec3X7UHi93uGMboST8R4Aw5bIxyUSA8do4hlSyFevXi3LsmK3g8GgsrOzr/n1f/zjHzVt\n2jTdfPPN//uE19De3m7bvoHh4rjEaOb1ejlGDTXQE7Ahhfzhhx/ud9vlcum222675tcfOnRIHR0d\nOnTokM6fPy+n06nU1FT19PQoJSVFwWBQHo9HHo+n3xl4MBiUz+cbykgAAEBDDPmMGTN09OhRtbW1\nSZKmTZvW7wz985588snYn3fs2KGbb75ZLS0tCgQCmjt3rgKBgAoKCuTz+bR161aFw2ElJyerpaVF\n5eXl/9t3BADADWRIIf/1r3+thoYG3X777YpGo9q2bZuKior03e9+d8j/o4ceekjV1dXas2ePsrKy\nNG/ePDkcDi1evFhVVVWyLEtlZWWxF74BAIDBDSnkTU1NWr9+fexT3Xp7e7VmzZohhfyhhx6K/XnV\nqlVfuH/OnDmaM2fOUOcFAAD/ZUiftR6NRvt9NGtycvKAS+sAAGBkDOmMfMqUKdq4caPy8/MlSQ0N\nDcrJybF1MAAAMLhBQ97R0aHy8nIdPHhQx44dk2VZmjFjhhYsWDAS8wEAgAEMuLTe2NioVatWKRKJ\n6O6771Z5ebm+9a1vadeuXTp+/PhIzQgAAK5hwJDX1NRo5cqV/V5Jnp2drWeeeUbvvPOO7cMBAICB\nDRjyaDT6pZ/gduutt+ry5cu2DQUAAIZmwJBHIpFr3tfV1XXdhwEAAF/NgCG/9dZbtWvXri9sf/fd\nd5Wbm2vbUAAAYGgGfNX6D37wA73wwgs6cOCAcnJy1NfXp5aWFqWlpamysnKkZgQAANcwYMgzMjJU\nVVWlxsZGnTx5UklJSSouLtaMGTNGaj4AADCAIX0gTH5+fuzDYAAAwOgxpI9oBQAAoxMhBwDAYIQc\nAACDEXIAAAxGyAEAMBghBwDAYIQcAACDEXIAAAxGyAEAMBghBwDAYIQcAACDEXIAAAxGyAEAMBgh\nBwDAYIQcAACDEXIAAAxGyAEAMJjDrh13d3dr8+bNunDhgi5fvqwHH3xQt912m6qrq9XX16eMjAxV\nVFTI6XSqvr5edXV1sixLpaWlKikpsWssAAASim0h/+tf/6qcnBzdf//9Onv2rJ5//nnl5eVp/vz5\nKi4u1vbt27Vv3z7NnTtXNTU12rBhgxwOhyorK1VUVKT09HS7RgMAIGHYtrR+11136f7775cknT9/\nXpmZmWpqapLf75ck+f1+NTQ0qLW1VTk5OXK5XEpJSVFeXp6am5vtGgsAgIRi+zXyZ599Vps2bVJ5\nebm6u7vldDolSW63W6FQSKFQSG63O/b1V7cDAIDB2ba0ftXzzz+vEydO6Je//KWi0eh126/X671u\n+xptTsZ7AAxbIh+XSAwco4nHtpAfP35cbrdbWVlZmjx5snp7e5WWlqaenh6lpKQoGAzK4/HI4/H0\nOwMPBoPy+XyD7r+9vd2u0YFh47jEaOb1ejlGDTXQEzDbltaPHDmi9957T5IUCoUUiUSUn5+vQCAg\nSQoEAiooKJDP51NbW5vC4bAikYhaWlo0ffp0u8YCACChWNHrud79X3p6evTyyy/r/Pnz6unpUVlZ\nmXJyclRdXa3Lly8rKytLjz32mBwOhwKBgGpra2VZlu677z5985vfHHT/ifyssnfpgniPgGFKfrU2\n3iMA18QZubkGOiO3LeR2S+SDkZCbi5BjNCPk5orL0joAALAfIQcAwGCEHAAAgxFyAAAMRsgBADAY\nIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAAgxFyAAAM\nRsgBADAYIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAAgxFyAAAMRsgBADAYIQcAwGCEHAAA\ngxFyAAAM5rBz52+99ZaOHj2qvr4+PfDAA8rJyVF1dbX6+vqUkZGhiooKOZ1O1dfXq66uTpZlqbS0\nVCUlJXaOBQBAwrAt5IcPH9bJkydVVVWlrq4uPf3008rPz9f8+fNVXFys7du3a9++fZo7d65qamq0\nYcMGORwOVVZWqqioSOnp6XaNBgBAwrBtaX3GjBl68sknJUljx45Vd3e3mpqa5Pf7JUl+v18NDQ1q\nbW1VTk6OXC6XUlJSlJeXp+bmZrvGAgAgodgW8qSkJKWmpkqS9u7dq8LCQnV3d8vpdEqS3G63QqGQ\nQqGQ3G537HFXtwMAgMHZeo1ckj766CPt3btXzz77rJ544onrtl+v13vd9jXanIz3ABi2RD4ukRg4\nRhOPrSH/29/+pt/97ndauXKlXC6XUlNT1dPTo5SUFAWDQXk8Hnk8nn5n4MFgUD6fb9B9t7e32zk6\nMCwclxjNvF4vx6ihBnoCZtvS+qVLl/TWW2/pZz/7WeyFa/n5+QoEApKkQCCggoIC+Xw+tbW1KRwO\nKxKJqKWlRdOnT7drLAAAEoptZ+QffPCBurq69NJLL8W2Pf7449q6dav27NmjrKwszZs3Tw6HQ4sX\nL1ZVVZUsy1JZWZlcLpddYwEAkFCsaDQajfcQw5HIy0O9SxfEewQMU/KrtfEeAbgmltbNFZeldQAA\nYD9CDgCAwQg5AAAGI+QAABiMkAMAYDBCDgCAwQg5AAAGI+QAABiMkAMAYDBCDgCAwQg5AAAGI+QA\nABjM1t9HDgAmuf/t5niPYLPE/v7eXXx7vEeIC87IAQAwGCEHAMBghBwAAIMRcgAADEbIAQAwGCEH\nAMBghBwAAIMRcgAADEbIAQAwGCEHAMBghBwAAIMRcgAADEbIAQAwGCEHAMBghBwAAIPZ+vvI//nP\nf+qFF17Qd77zHd133306d+6cqqur1dfXp4yMDFVUVMjpdKq+vl51dXWyLEulpaUqKSmxcywAABKG\nbWfkkUhE27Zt08yZM2PbduzYofnz5+u5557ThAkTtG/fPkUiEdXU1GjVqlVau3atdu7cqYsXL9o1\nFgAACcW2kDudTlVWVsrj8cS2NTU1ye/3S5L8fr8aGhrU2tqqnJwcuVwupaSkKC8vT83NzXaNBQBA\nQrFtaT05OVnJycn9tnV3d8vpdEqS3G63QqGQQqGQ3G537GuubgcAAIOz9Rq5nbxeb7xHsM3JeA+A\nYUvk4/LGwGqgyW7Un78RDXlqaqp6enqUkpKiYDAoj8cjj8fT7ww8GAzK5/MNuq/29nY7RwWGheMS\niJ9E/vkb6EnKiL79LD8/X4FAQJIUCARUUFAgn8+ntrY2hcNhRSIRtbS0aPr06SM5FgAAxrLtjPz4\n8eN68803dfbsWSUnJysQCOiJJ57Q5s2btWfPHmVlZWnevHlyOBxavHixqqqqZFmWysrK5HK57BoL\nAICEYkWj0Wi8hxiORF5C6V26IN4jYJiSX62N9wj4H9z/NtfITfbu4tvjPYJtRs3SOgAAuL4IOQAA\nBiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4A\ngMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJAD\nAGAwQg4AgMEIOQAABiPkAAAYjJADAGAwQg4AgMEc8R7gqjfeeEPHjh2TZVkqLy9Xbm5uvEcCAGDU\nGxVn5EeOHNHp06dVVVWl5cuXa9u2bfEeCQAAI4yKkDc2NurOO++UJE2aNEnhcFiXLl2K81QAAIx+\no2JpPRQKaerUqbHbbrdboVBILpfrmo/xer0jMVp87PxLvCcAbkgfPZXA/64gYY2KM/LPi0aj8R4B\nAAAjjIqQezwehUKh2O3Ozk55PJ44TgQAgBlGRchnzZqlQCAgSTp+/Lg8Ho/S0tLiPBUAAKOfFR0l\n69hvv/22jh49KsuytGTJEk2ePDneIwEAMOqNmpADAICvblQsrQMAgOEh5AAAGGxUvI8ciS0SicTe\nlZCRkaHU1NQ4TwTc2MLhsMaOHRvvMXCdEHLYpq2tTdu2bVM4HJbb7VY0GlVnZ6cyMzO1ZMkSZWdn\nx3tE4Ib04osvas2aNfEeA9cJIYdt3njjDf34xz/WxIkT+20/fvy4Xn/9da1bty5OkwGJ7w9/+MOX\nbo9GowoGgyM8DexEyGGbaDT6hYhL0tSpU9XX1xeHiYAbx3vvvaf8/Pwv/XCt3t7eOEwEuxBy2Mbn\n82njxo2688475Xa7Jf3nc/X//Oc/a/r06XGeDkhsTz31lLZt26Yf/ehHcjqd/e47cuRInKaCHXgf\nOWx15MgRHT58OPZiN4/Ho1mzZmnatGlxngxIfN3d3XI6nUpK6v8GpePHj/f7RVUwGyEHAMBgvI8c\nAACDEXIAAAzGi92AG1xHR4d+8pOffOF1Cy6XS7Nnz9a9994bp8kADAUhByC32621a9f227Z58+b4\nDAPgKyHkAAa1d+9e7d69W2PGjNG4ceO0bNkyuVwuHTp0SDU1NRozZoxSUlK0bNkyZWZm6vHHH1dx\ncbE6Ojr005/+NN7jAwmNa+QABnTu3Dn95je/0erVq7V27VqNHz9eO3fuVHd3t7Zu3aoVK1ZozZo1\nKiws1DvvvBN73C233ELEgRHAGTkAffbZZ19YWh8zZoyk/7zneMqUKUpLS5Mkff3rX9fu3bv1r3/9\nS+PGjdP48eMlSTNmzNDu3btjj8/LyxuZ4YEbHCEHMOA1csuy+m0f6kdPOBz88wKMBJbWAQxo6tSp\n+sc//qF///vfkqTGxkZNmzZNXq9Xn332mc6dOydJamhokM/ni+eowA2Jp8wABjR+/HgtWrRI69ev\nl9PpVGZmpr73ve8pJSVFy5cv10svvSSn06nU1FQtX7483uMCNxw+ohUAAIOxtA4AgMEIOQAABiPk\nAAAYjJADAGAwQg4AgMEIOQAABiPkAAAYjJADAGCw/wOJ0A/8dBfxZAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fc814cf5090>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "id": "FAzfFiGnCyls",
        "colab_type": "code",
        "outputId": "5ce1d951-11e9-4e28-dd4d-8317019fad13",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 364
        }
      },
      "cell_type": "code",
      "source": [
        "df['log_u2'].plot(kind='density')"
      ],
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.axes._subplots.AxesSubplot at 0x7fc811411110>"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 11
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFKCAYAAAAE1MaxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt8VOWB//HPMwkhCTAwAUIIIcQQ\nCEYloISCUBer1lbbalupUtrKGral19+v2m7dChXqpsq+dtmlrvKru12opbqrtrV0txergKKIy0Uk\nGC6BcAshCZBMIAkJJPP8/phkBAkhQM6cuXzfr5cvZzJzZr4PuXznnPOcc4y11iIiIiJRz+N2ABER\nEekdKnUREZEYoVIXERGJESp1ERGRGKFSFxERiREqdRERkRihUhcREYkRiU6++IoVKygvL8cYw5w5\nc8jLyzvvOc899xy7d+9m4cKFPV7mw6qqqno7+mXJzMyMmCxu0Pg1fo1f449H4R57ZmbmBR9zbE29\nrKyM6upqSkpKmDdvHsuXLz/vOZWVlezYseOSlhEREZGuOVbqpaWlFBUVAZCVlUVTUxPNzc3nPOfZ\nZ5/lvvvuu6RlREREpGuOlbrf78fr9Ybue71e/H5/6P7atWspKChg6NChPV5GRERELszRfepnO/sU\n842NjaxZs4YFCxZQV1fXo2W6093+hXCLpCxu0Pg1/nim8cfv+CNl7I6Vus/nO2ctu76+Hp/PB8D2\n7ds5ceIEjz76KGfOnKGmpoYVK1Z0u0x3ImVyRjxPFAGNX+PX+DX++Bx/JE2Uc6zUCwsLeeGFF7jt\nttuoqKjA5/ORkpICwJQpU5gyZQoAtbW1PP3008yZM4ddu3ZdcBkRERHpnmOlnp+fT25uLvPnz8cY\nQ3FxMWvXriU1NZXJkyf3eBkRERHpGUf3qc+ePfuc+zk5Oec9Jz09PXSMelfLiIiISM/ojHIiIiIx\nQqUuIiISI8J2SJuIRDfbdgbKtmJrqjBpQ+Ca6zHJmsgqEklU6iJyUXb3dgL/8S9wvDZ4H6C/FzPr\nq3gm3+RqNhH5gEpdRLplSzcReKoErMXcfCcm/1rswX3Y136P/bd/JHCiHs+td7kdU0RQqYtIN+yB\nvQT+3xPgScDz7QWYqwsBMDdMw06ZQeCf5mP/6+cEBgyCz85yOa2IaKKciHTJtrYS+Pd/gtOn8Xzt\nb0OF3skMz8Lz3R9Dcgr2l09zpuqQS0lFpJNKXUS6ZFc9B9WVmFs+jSns+oRRZkQ25kvfgNZT1P3T\no9hAIMwpReRsKnUROY+tqcK+9nsYnI753Fe6fa7nI3+FuWEap3duw65/LUwJRaQrKnUROU/gpRXQ\n3oZn5l9jkvpe9PnmC8WY5BTsr3+BbWl2PqCIdEmlLiLnsIf2wdYNkJsP19/Yo2VM2hAG3HM/NJ7A\nvrrK4YQiciEqdRE5h/3DiwB4PnUfxpgeLzfgrlnQ34t95WVs00mn4olIN1TqIhJij1RiN78F2aPh\n2usvaVlPaj/MJz8Pp5qxr/zOoYQi0h2VuoiE2NW/B2vx3DHzktbSO5kZd8CAgdi1f8C2nHIgoYh0\nR6UuIgDYU83Yt9dC2hCY8JHLeg2T1DdY7M2N2Lc0E14k3FTqIgKA3bAGWk9hbvoEJiHhsl/H3HwH\n9EnCvvo7bKC9FxOKyMWo1EUEay12zR8gIRHz0duu6LXMgIGYGz8Gx2rg3Q29lFBEekKlLiKwezsc\nOYS5YRrG67vilzMdF3gJvPr7K34tEek5lbqIYN96FQDzV7f3yuuZjBFQMAH2lGEPH+yV1xSRi1Op\ni8Q523IKu+VtGDIMxlzTa6/ruekTwddf9+dee00R6Z5KXSTO2Xc3QGsLZurNl3UY2wUVToaBPuzb\nq7GnW3vvdUXkglTqInHOvr0aADPl5l59XZOYiJl2GzQ3YTe92auvLSJdU6mLxDFbfxx2boPR4zDp\nw3v99c1NHwdjsK//qddfW0TOp1IXiWP2f98Aa3t9Lb2TGZwO11wPFbs0YU4kDFTqInHMbn4LPB7M\nDdMcew/PtFuC77VhjWPvISJBKnWROGWP18K+3ZB/HWaA17k3KpwMKf2w77yuM8yJOCzRyRdfsWIF\n5eXlGGOYM2cOeXl5ocdeffVV1qxZg8fjYdSoURQXF1NWVsaSJUsYOXIkANnZ2TzwwANORhSJW3bL\n2wCOrqUDmD5JmEnTsOtegV3b4epCR99PJJ45VuplZWVUV1dTUlJCZWUly5Yto6SkBIDW1lbWr1/P\nokWLSExMZNGiRezevRuAgoICHnroIadiiUgHu/ktMB7MxCmOv5eZ+jHsulewb6/GqNRFHOPY5vfS\n0lKKiooAyMrKoqmpiebmZgD69u3Lj370IxITE2ltbaW5uZlBgwY5FUVEPsTWH4e9O2HsNRhvGH73\n8q6GIcOwW97WJVlFHORYqfv9frzeD/bTeb1e/H7/Oc95+eWX+fa3v83UqVMZNmwYAJWVlSxevJgF\nCxawbds2p+KJxLVwbXrvZIwJzrBvbQme7EZEHOHoPvWzWWvP+9rdd9/NHXfcweOPP864ceMYPnw4\nM2fOZOrUqdTU1LBo0SKefPJJEhO7j5mZmelU7EsWSVncoPFHx/hrd71HK5Bx+2dIHJLea6/b3fjP\nfOYLVP/3f9K3dCNDPz+7194zkkTL998p8Tz+SBm7Y6Xu8/nOWTOvr6/H5wte/amxsZGDBw9SUFBA\nUlISEyZMYNeuXYwbN44bb7wRgIyMDAYNGkRdXR3p6d3/0amqqnJqGJckMzMzYrK4QeOPjvHblmYC\npZshO5fa023QS5kvOn6TCFlX0bJlA4f37Mak9u+V940U0fL9d0o8jz/cY+/uA4Rjm98LCwvZsCG4\nma2iogKfz0dKSgoAbW1tPP3007S0tACwZ88eMjMzWbduHatWrQKCm+8bGhpIS0tzKqJIfCp7D9ra\nMNdNCvtbm0nToL0Nu/WdsL+3SDxwbE09Pz+f3Nxc5s+fjzGG4uJi1q5dS2pqKpMnT+aee+5h0aJF\noUPaJk2aREtLC0uXLmXTpk20tbUxd+7ci256F5FLY0s3AWDGF4X9vc2k6diXV2I3vQU33hL29xeJ\ndY425uzZ5+43y8nJCd2eMWMGM2bMOOfxlJQUHn74YScjicQ1GwgES33AQMgZE/b3N8MyIesqKNuK\nbW6MuU3wIm7TGeVE4smhCmiox1x7A8bjzq//B5vg/9eV9xeJZSp1kThit3Vueg///vROnYfR2c1v\nuZZBJFap1EXiiN2xFYwHrp7gWgaTMQKycqDsXWxLs2s5RGKRSl0kTtiWU1CxG0aNxvRzd1+2KZwM\nbW3Bmfgi0mtU6iLxorwM2tswV493O0mw1AG7TfvVRXqTSl0kTtidwbViMy4CLqgyKg+8g7DbNmED\nAbfTiMQMlbpInLA73oPEPsGLq7jMeDzBk9+cbID95W7HEYkZKnWROGBPNsChfZB3NSapr9txgLM2\nwb+30eUkIrFDpS4SB+zOUgDMOPf3p4dcXQiJidqvLtKLVOoi8WBX8DLGkVTqJjkFxo2Hyv3Y47Vu\nxxGJCSp1kThg9+yApL7BCWoRxIzv2ATfcT56EbkyKnWRGGebGqHqIFw1FhNhF0gy10wEwJZtdTmJ\nSGxQqYvEuoqdYC0mAma9n2doBgxOh52l2EC722lEop5KXSTG2fIyAExegctJzmeMwRRMgFNNsH+P\n23FEop5KXSTG2b07gud7Hz3O7ShdMgXB89BrE7zIlVOpi8Qwe+YM7CuHEaMwKalux+nauPFgTPBi\nMyJyRVTqIrHs4F44cxozJgL3p3cw/b2QPRr27gpedEZELptKXSSG2T07gjdGR26pQ8cm+PY22L3d\n7SgiUU2lLhLDbMVOgMic+X6W0H71HboUq8iVUKmLxLL95TBgIKQNdTtJ90ZfDUlJKnWRK6RSF4lR\n9kQ91B2DnDEYY9yO0y3Tp0+w2A8fwJ484XYckailUheJVfuCx32bq8a6HKRnzJhrgjf2lLkbRCSK\nqdRFYpTdvxsAkzPG3SA9ZMZeC4DVZDmRy6ZSF4lRdn958EaUlDpXjQleinX3+24nEYlaKnWRGGSt\nDU6SG5yOGeB1O06PmKS+cNVYOLQP29zkdhyRqKRSF4lFx2qg8WTU7E/vZMZcCzYAe3e4HUUkKjl6\nHcYVK1ZQXl6OMYY5c+aQl/fBtZxfffVV1qxZg8fjYdSoURQXF2OM6XYZEemZqNv03sHkX4P9A9jd\n72Oum+R2HJGo49iaellZGdXV1ZSUlDBv3jyWL18eeqy1tZX169ezaNEiHnvsMQ4fPszu3bu7XUZE\nLkFHqUfLJLmQ3HHg8WiynMhlcqzUS0tLKSoqAiArK4umpiaam5sB6Nu3Lz/60Y9ITEyktbWV5uZm\nBg0a1O0yItJz9mBF8EZ2rrtBLpFJToFReXBgD7a1xe04IlHHsVL3+/14vR9M0PF6vfj9/nOe8/LL\nL/Ptb3+bqVOnMmzYsB4tIyLds9bCoX0wNCNyr8zWDTP2Gmhvh4pdbkcRiTqO7lM/m7X2vK/dfffd\n3HHHHTz++OOMG3f+tZ67WqYrmZmZV5yvt0RSFjdo/O6Pv+1oNUeaTpIyoYghYc7TG+M/VTSNY3/+\nLf2PHmbgLZ/shVThEwnffzfF8/gjZeyOlbrP5ztnLbu+vh6fzwdAY2MjBw8epKCggKSkJCZMmMCu\nXbu6XaY7VVVVvT+Ay5CZmRkxWdyg8UfG+O17GwFoGTI8rHl6a/zWlw7AiXc30nTTHVf8euESKd9/\nt8Tz+MM99u4+QDi2+b2wsJANGzYAUFFRgc/nIyUlBYC2tjaefvppWlqC+8z27NlDZmZmt8uISM/Y\nQ8H96WZkdO1P72S8g2BoBlTswgYCbscRiSqOrann5+eTm5vL/PnzMcZQXFzM2rVrSU1NZfLkydxz\nzz0sWrQodEjbpEmTMMact4yIXBp7aF/wxsir3A1yBczocdgNa6HmMAwf6XYckajh6D712bNnn3M/\nJycndHvGjBnMmDHjosuIyCU6VAH9B4BvsNtJLl/uONiwFrt3J0alLtJjOqOcSAyxp5rhaDWMzI34\ny612x4zumDirGfAil0SlLhJLKvcDYKJ40zsAI0ZB32TsHp0uVuRSqNRFYkjnJLlo3p8OYBISgqe4\nPXII29zodhyRqKFSF4klHZPkonXm+9nM6KuDN7QJXqTHVOoiMcQe2geJfWDYCLejXDEzOh8Au1el\nLtJTKnWRGGHb2+HwARgxCpMYtpNFOie3o9QrdrocRCR6qNRFYkXtEWg7gxkxyu0kvcL090LGiI6T\n0LS7HUckKqjURWJF1YHg/0dku5ujF5nccdByCo5Uuh1FJCqo1EVihD18EACTGTulzlXB68HbjuvD\ni0j3VOoisaIqWOrEUKmbnGCps2+3qzlEooVKXSRG2KqDkJIKviFuR+k9WTmQmIjdv8ftJCJRQaUu\nEgPsmTNQWwWZ2VF9etgPM4l9YGQuVO7DnjntdhyRiKdSF4kFNYehvT229qd3MDljoL0dDla4HUUk\n4qnURWKAjcH96SFXjQXQJniRHlCpi8SCqhic+d4hNFluvybLiVyMSl0kBnQezkaMnHjmHMMyISVV\nh7WJ9IBKXSQWVB2EfgPAO8jtJL3OeDzBK7ZVH9YV20QuQqUuEuXs6VY4Wg2ZI2Nq5vvZTE5e8Ib2\nq4t0S6UuEu2qK8EGYuac710xOZ2T5bQJXqQ7KnWRKBfTM987dc6A36dSF+mOSl0k2sXwzPdOxjcY\nBqaB1tRFuqVSF4lyoZnvMVzqQPDiLv7j2PrjbicRiVgqdZFoV3UQBgzEDBjodhJHfXC8utbWRS5E\npS4SxWzLKThWE/tr6YDRZVhFLkqlLhLNjlQCxPTM95BRHaWuy7CKXJBKXSSKxcXM9w6mX38YmgEH\n9mKtdTuOSERKdPLFV6xYQXl5OcYY5syZQ15eXuix7du38/zzz+PxeBg+fDjz5s1jx44dLFmyhJEj\nRwKQnZ3NAw884GREkehWdQCI7ZnvZzM5Y7Ab1wV3OQzNcDuOSMRxrNTLysqorq6mpKSEyspKli1b\nRklJSejxZ555hkcffZTBgwezZMkStm7dSt++fSkoKOChhx5yKpZITImnNXUARo2GjevgwB6VukgX\nHNv8XlpaSlFREQBZWVk0NTXR3NwcevyJJ55g8ODBAHi9XhobdU5nkUtWdRAGpQU3TccBkz0aAHtg\nr8tJRCKTY6Xu9/vxer2h+16vF7/fH7qfmpoKQH19Pe+99x4TJ04EoLKyksWLF7NgwQK2bdvmVDyR\nqGebm6DuWPyspUNwTR2wB3QOeJGuOLpP/WxdTWxpaGhg8eLFzJ07lwEDBjB8+HBmzpzJ1KlTqamp\nYdGiRTz55JMkJnYfMzMz06nYlyySsrhB4w/f+Ft3bKMW6D+2AF+E/LuHY/xHhmcROLSP4cOHR9wF\nbPTzH7/jj5SxO1bqPp/vnDXz+vp6fD5f6H5zczM/+clPmDVrFoWFhQCkpaVx4403ApCRkcGgQYOo\nq6sjPT292/eqqqpyYASXLjMzM2KyuEHjD+/4A9u2ANDkTeNUBPy7h2v8gRE52E1vUrXtXUwE7VfX\nz3/8jj/cY+/uA4Rjm98LCwvZsGEDABUVFfh8PlJSUkKPP/vss9x5551MmDAh9LV169axatUqILj5\nvqGhgbS0NKciikS3w/E18z2kYxM8B7VfXeTDHFtTz8/PJzc3l/nz52OMobi4mLVr15KamkphYSFv\nvPEG1dXVrF69GoDp06czbdo0li5dyqZNm2hra2Pu3LkX3fQuEq/ibuZ7BzMqD0twv7q5YZrbcUQi\niqONOXv27HPu5+TkhG4/99xzXS7z8MMPOxlJJHZUHYK0IZiUVLeThFfnZLn9miwn8mE6o5xIFLJN\njdBQF3dr6QAmVWeWE7kQlbpINIqDa6h3x4zKg+bG4JnlRCREpS4SheJ1f3qIJsuJdEmlLhKNtKYO\n6CQ0Ih+mUheJQqE19eEj3Q3iFp0uVqRLKnWRaFR1EAanY5JTLv7cGKTLsIp0TaUuEmVs4wk44Y/f\n/ekdzKg8aDqpyXIiZ1Gpi0Sb0P70ON303kmT5UTOo1IXiTJxP/O9gybLiZxPpS4SbeJ85nuIJsuJ\nnEelLhJlbNWh4I14nfneQZPlRM6nUheJNlUHYcgwTN9kt5O4zmSPDk6WO17rdhSRiKBSF4ki9mQD\nnGyI+/3pIR371dEmeBFApS4SXbQ//Rwmp3OyXLm7QUQihEpdJIpo5vuHhCbLVbgcRCQyqNRFoonW\n1M9h+vWHwelwqEKT5URQqYtEFVt1EIyBjCy3o0SO7NzgPAN/ndtJRFynUheJEtbas2a+93U7TsQw\n2bnBGwe1CV6kR6X+q1/9iiNHjjidRUS6c9IPjSe1P/1DzMiO/eqHNANeJLEnT+rfvz8//elP6du3\nLx/72MeYMmUKSUlJTmcTkbMd7tifPmKUy0EiTMeautWaukjPSv2uu+7irrvuoqamhrfffpvHHnuM\nUaNG8clPfpIRI0Y4nVFEAFu5P3gj6ypXc0ScQWkwYKA2v4twifvUjx8/TnV1NadOnSI5OZmnnnqK\nV155xalsInK2jlI3WTmuxog0xhgYmQvHa7FNJ92OI+KqHq2pv/jii6xbt47hw4dz22238dWvfhWP\nx0NbWxt/93d/x8c//nGnc4rEPVu5D/okQfpwt6NEHDMqF1v2bnBt/epCt+OIuKZHpd7Q0MCCBQsY\nOnRo6Gu1tbWkp6cze/Zsx8KJSJBtbw/OfB+Rg0lIcDtO5AlNlqvAqNQljl1083sgEODw4cMMGTKE\nQCBAIBDgzJkzLF68GIAJEyY4HlIk7tUchrY2TJYmyXVFh7WJBHW7pv7mm2/y4osvUl1dzX333Rf6\nujGGwkJ9GhYJlw8myeW4GSNyDc2A5BTNgJe4122pT58+nenTp/PCCy/whS984ZJffMWKFZSXl2OM\nYc6cOeTl5YUe2759O88//zwej4fhw4czb948PB5Pt8uIxK3OSXIjclyNEamMxxM8KmDvTmxrq07O\nI3Gr21J/9913mThxIkOGDGH16tXnPf6xj33sgsuWlZVRXV1NSUkJlZWVLFu2jJKSktDjzzzzDI8+\n+iiDBw9myZIlbN26leTk5G6XEYlXOpzt4kx2LnZPGVTug9Hj3I4j4opuS/3AgQNMnDiRnTt3dvl4\nd6VeWlpKUVERAFlZWTQ1NdHc3ExqaioATzzxROi21+ulsbGR8vLybpcRiVuH98OgNMwAr9tJIlf2\nWZPlVOoSp7ot9bvvvhuAb3zjG1hrMcZw5swZGhoaGDJkSLcv7Pf7yc3NDd33er34/f5QQXf+v76+\nnvfee497772XHTt2dLuMSDyyTY1Qdwyuvd7tKBHNZOdiQZPlJK716JC23/72tyQnJ3PLLbfwgx/8\ngJSUFAoLC7n33nt7/EZdXRaxoaGBxYsXM3fuXAYMGNCjZbqSmZnZ4xxOi6QsbtD4e3/8Ldu3cBQY\nkH8tgyL839fN779NT6cysQ+J1ZVkuJRDP//xO/5IGXuPSn3z5s089thjvP7669xwww186UtfYtGi\nRd0u4/P58Pv9ofv19fX4fL7Q/ebmZn7yk58wa9as0Ez6iy1zIVVVVT0ZhuMyMzMjJosbNH5nxh/Y\nugmAxkFDaI7gf9+I+P5nZnNmXzmHDx7EJPboz1vvvXUkjN9F8Tz+cI+9uw8QPTpNbEJCAsYYtm7d\nyuTJk4Hg8evdKSwsZMOGDQBUVFTg8/lISUkJPf7ss89y5513nnOc+8WWEYlLOj1sj5nsXGg7A9WV\nbkcRcUWPPsr269ePxx9/nLq6OsaOHcvmzZvxeLr/PJCfn09ubi7z58/HGENxcTFr164lNTWVwsJC\n3njjDaqrq0Oz6qdPn86tt9563jIi8c5W7oeERMjQxZMuKnTFtr36ECRxqUel/p3vfIdt27Yxblxw\nRmmfPn345je/edHlPnwK2ZycnNDt5557rkfLiMQzGwjA4QMwPAuT2MftOBHPZI/+YLLcjbe4HUck\n7HpU6h6PB2MMmzdvDk1eO3bsWLeHtIlIL6itgtOtGB2f3jNZOWAM9pBmwEt86lGpl5SU4PF4zrmg\nC3R/nLqIXDl7YG/wxqjR7gaJEqZvMgwbAYf2YQOB4JnmROJIj0q9vb39orPdRcQBHcdcm2yVek+Z\n7Fzs/1bCsRpdplbiTo8+xmZlZXHy5Emns4jIh9iDHWvq2dr83mOdV2zTJniJQz1aU6+rq+M73/kO\nI0aMIOGsazlr7V3EOdZaOLgXho3AJOusij1lRgbPLGcPVmBumOZ2HJGw6lGpd54uVkTC6FgNNDdh\nrtHpYS/JWYe1icSbHm1+LygooKWlhYMHD1JQUEBaWhpXX32109lE4lvnOcw1Se6SmP5eSBuqc8BL\nXOpRqa9cuZLVq1ezZs0aAN58802WL1/uaDCReNe5pqlJcpchOxdO+LH+OreTiIRVj0p9x44dfO97\n3wtdLe2ee+5h3759jgYTiXcfTJLL7f6Jch4zUpPlJD71qNSTkpLOuR8IBGhvb3ckkIh0TJI7sBcG\np2P6nX8FQ+me6dhlYbUJXuJMjybKjR07lqeeeor6+nr++7//m3feeYeCggKns4nEL38dnGyA66e6\nnSQ6jdRkOYlPPSr1m266iQMHDrB371527tzJpz/96dDV2kTEAdqffmV8g6G/V5PlJO50W+qnT59m\n6dKlHDhwgNzcXNLS0ti5cydJSUlcf/31JIb5esUi8cLu3wOo1C+XMQayR0PZu9imRky//m5HEgmL\nblv5pZdeIi0tjQcffDB00pnW1laWL1/Oc889x1e+8pWwhBSJN3b/7uCNnDHuBoliJjsXW/ZucLLc\nuPFuxxEJi24nyu3cuZP777//nLPI9e3bl7lz57Jt2zbHw4nEI2st7CuHoRmYAV6340SvbE2Wk/jT\nbaknJCR0uYk9MTGRfv36ORZKJK4dPQJNJzFaS78ipvNQQE2Wkzhy2dcl9OiShiKOsBUdm95zx7ob\nJNoNzYDkFK2pS1zpdp/67t27+frXv97lYydOnHAkkEjc218OgMlRqV8J4/HAyKtgz05sayumb1+3\nI4k4rttS/5d/+Zdw5RCRDnbfbkhI0JnkeoHJHo0tL4PKfTB6nNtxRBzXbakPHTo0XDlEBLBtZ4LH\nVo/IwSRpzfKKdV6x7VAFRqUucUA7xkUiSeV+aDuDuUqT5HrDB5PltF9d4oNKXSSC2H0dk+Suync3\nSKzIGAmJfTRZTuKGSl0kknSUutbUe4dJTIQRo+Dwfmxbm9txRBynUheJIHZfOSSnQMYIt6PEDJOd\nC21tcOSQ21FEHKdSF4kQtrkRqishZwzGk3DxBaRnQmeW00loJPap1EUiRefx6dr03qs0WU7iiaOX\nWVuxYgXl5eUYY5gzZw55eXmhx06fPs0zzzxDZWUlTzzxBADvv/8+S5YsYeTIkQBkZ2fzwAMPOBlR\nJGLYPTsBMLk69KpXZeWAx6PJchIXHCv1srIyqqurKSkpobKykmXLllFSUhJ6fOXKleTk5FBZWXnO\ncgUFBTz00ENOxRKJWHbvjuCN0Ve7GyTGmKS+kJEFh/ZhA4HgmeZEYpRjP92lpaUUFRUBkJWVRVNT\nE83NzaHHZ82axeTJk516e5GoYgPtULELho3QldkcYLJHQ+spqD3idhQRRzlW6n6/H6/3gz9OXq8X\nv98fup+SktLlcpWVlSxevJgFCxbo8q4SPw4fhJZTmDxtendE55nlNFlOYpyj+9TPZq296HOGDx/O\nzJkzmTp1KjU1NSxatIgnn3yyy8u/ni0zM7O3Yl6xSMriBo3/8sZ/cstb+IFBN0ylfxT/G0bq979l\nYhFHX/g5/etqGeRgxkgdf7jE8/gjZeyOlbrP5ztnzby+vh6fz9ftMmlpadx4440AZGRkMGjQIOrq\n6khPT+92uaqqqisP3AsyMzMjJosbNP7LH39g8wYAGoZkciJK/w0j+ftvU4NbDU/u2EazQxkjefzh\nEM/jD/fYu/sA4djm98LCQjZsCP6hqqiowOfzXXCTe6d169axatUqILj5vqGhgbS0NKciikQMu3cH\n9B+gk844xKT2D15f/eDeHm1P6t2UAAAabklEQVQ1FIlWjq2p5+fnk5uby/z58zHGUFxczNq1a0lN\nTWXy5MksWbKE48ePU1VVxcKFC7n11luZNGkSS5cuZdOmTbS1tTF37tyLbnoXiXa2/jgcr4XCyRhj\n3I4Tu7JzYfN6qD8GaboCpcQmRxtz9uzZ59zPyckJ3X7wwQe7XObhhx92MpJI5Ok4lE2XBnWWGZmL\n3bweDu5VqUvM0gGbIi6zezpLXcenO8mM6jxdrE5CI7FLpS7iMrtnByQkQk7exZ8sly90WJtKXWKX\nSl3ERba1BQ5VwKjRwTOfiWOM1wcD03QOeIlpKnURN+0vh0AAk6dN72ExajTUH8OeqHc7iYgjVOoi\nLtL+9PAKXQFv3x53g4g4RKUu4iK7+/3gDa2ph4XJGQuA3b/b3SAiDlGpi7jEtrUFD2cbPhLjHeR2\nnPjQMRnR7lOpS2xSqYu45cAeaG3BjL3G7SRxw/T3Bs8st69cZ5aTmKRSF3FJaNP72GvdDRJnzFVj\nobkRjuoyrBJ7VOoiLrG7twNgVOrh1TFZzu4rdzmISO9TqYu4wLa3w54ySM/EDNJFi8Kpc7Ic2q8u\nMUilLuKGQxXQcgqTr7X0sMvOBY8Hu19r6hJ7VOoiLujc9K796eFnkvpCVg4c2Bs8AkEkhqjURVzQ\nOUlO+9PdYXLGQtsZOHzA7SgivUqlLhJmNtAO5e/D0AxM2hC348Snzsly2gQvMUalLhJulQeguUnH\np7vIXKXJchKbVOoiYab96RFgeBb0TcZW7HI7iUivUqmLhJmOT3ef8STAVWPhyCFsU6PbcUR6jUpd\nJIxsoB12bYfB6Zghw9yOE9dMXkHwxt4d7gYR6UUqdZFwOlgBzY2YqwvdThL3zJjglfHsnjKXk4j0\nHpW6SBjZHe8Fb6jU3ZebD8YTuqa9SCxQqYuEUWepm3HjXU4iJjkVRuYEr9h25ozbcUR6hUpdJEzs\n6VYoL4Osq3T99Ahh8gqCJ6E5sMftKCK9QqUuEi57dkDbGUyBNr1HjDztV5fYolIXCZPQpnftT48Y\nZnRnqWu/usQGlbpImNgd70FCIozRmeQihUkbAoPTobwseLihSJRLdPLFV6xYQXl5OcYY5syZQ15e\nXuix06dP88wzz1BZWckTTzzRo2VEopVtPAEH98KYazB9k92OI2cx48Zj33oVDu2DUfp7I9HNsTX1\nsrIyqqurKSkpYd68eSxfvvycx1euXElOTs4lLSMStXaVgrXa9B6JOr4nocMNRaKYY6VeWlpKUVER\nAFlZWTQ1NdHc3Bx6fNasWUyePPmSlhGJVrZM+9Mjlbk6eHih3bHN5SQiV86xUvf7/Xi93tB9r9eL\n3+8P3U9JSbnkZUSikbUWu30TpPaHnDFux5EPMV4fjBgFe97X8eoS9Rzdp342a61jy2RmZl7yazsl\nkrK4QeM/f/yn9++hpu4YqTd9nMEjR7qQKnyi9ftfP+lGGn/3PIMbjpI8ftJlv060jr+3xPP4I2Xs\njpW6z+c7Zy27vr4en8/X68sAVFVVXX7QXpSZmRkxWdyg8Xc9/sDqPwJwKu+amP73iebvvx05GoBj\nb67GM+Ty/jhH8/h7QzyPP9xj7+4DhGOb3wsLC9mwYQMAFRUV+Hy+Lje5X+kyIpHOlm4GYzDXXu92\nFLmQsdeCx4Pdqf3qEt0cW1PPz88nNzeX+fPnY4yhuLiYtWvXkpqayuTJk1myZAnHjx+nqqqKhQsX\ncuuttzJ9+vTzlhGJZra5EfaUQc4YzICBbseRCzApqcELvOzdhW08genvvfhCIhHI0X3qs2fPPuf+\n2YewPfjggz1aRiSqlW2FQABz3eXvp5XwMOOLsHt2YLdvxky52e04IpdFZ5QTcZAt3QyAuYLJVxIe\nZnzwcFre2+huEJEroFIXcYgNBLDbN4N3EIzMdTuOXExmNgxOx76/BdumQ9skOqnURZyybzec8GOu\nuwHj0a9apDPGYAonw6nm4CVyRaKQ/tKIOMRueRsAM/FGl5NIT3Vugrdb33E5icjlUamLOMBai333\nbeibArp+evTIvw76DcBufktXbZOopFIXccKhfXC0GjN+EqZPkttppIdMYiLmhhuhoR52v+92HJFL\nplIXcYB9t2PT+/VTXU4il8oUfRQAu3Gdy0lELp1KXcQBdvN66JME197gdhS5VGOvgYE+7Ob1mgUv\nUUelLtLLbNVBOHIIrpmISdZpjqON8SRgJk2HppOwbZPbcUQuiUpdpJfZDWsAMJP/yuUkcrnM9NsA\nCLzxJ5eTiFwalbpIL7KBAPad1yElFVNY5HYcuUwmKwdGj4Oyrdij1W7HEekxlbpIb9q9HeqOYW6Y\nhknq63YauQLmptvBWuybf3E7ikiPqdRFelFo07suCBL1zKTpkNof+8afsa2tbscR6RGVukgvCbS0\nBGe9pw2FMQVux5ErZJL6Ym6+AxpPaG1dooZKXaSXNL/xZ2g5hZl6s871HiPMLZ+GpCTsK7/FtrW5\nHUfkovSXR6QXWGtp/O8XwXiC+2IlJpgBAzHTPw51R7Fvr3Y7jshFqdRFesO+3ZzZuxMKJ2PShrqd\nRnqRuf1z0CcJu+o5bGuL23FEuqVSF+kFds0fAPDcfIfLSaS3mbQhmNvuBn8d9i8vux1HpFsqdZEr\nZI8fxW58g8SRV8G48W7HEQeYT34OBgzE/uk32OO1bscRuSCVusgVsq/8Ftrb8c68XxPkYpRJTsXc\nMwdaWwj88imstW5HEumS/gKJXAF7wo9d9woMTif1rz7hdhxxkJn6MbhmIrz/Lna9Js1JZFKpi1wB\n+4cX4cxpzO2fxSQmuh1HHGSMwfPlb0JyCvb5n2GPVLodSeQ8KnWRy2Rrj2DX/hGGDAse9iQxzwxO\nx3zlW8HN8D9brDPNScRRqYtcJvubZ6G9DfO5+zF9+rgdR8LEU/RRzIxPwuED2JVPa/+6RBSVushl\nsO+/i938Flw1FjNpmttxJMzMF4ohZwx2w5rgLhiRCKFSF7lEtmMGNB4Pni99HWOM25EkzEyfJDzf\nmg9pQ7EvrySw8U23I4kAKnWRS2Z//Qs4Xov5+Gcx2aPdjiMuMQN9eL69IDhxbvm/0Lpzu9uRRHB0\nuu6KFSsoLy/HGMOcOXPIy8sLPbZt2zaef/55PB4PEydO5J577uH9999nyZIljBw5EoDs7GweeOAB\nJyOKXBK75W3smv+BjCzMp+9zO464zGTl4Pnq3xJ48jGOPfYg/GAxZsgwt2NJHHOs1MvKyqiurqak\npITKykqWLVtGSUlJ6PHly5fzyCOPkJaWxsKFC5kyZQoABQUFPPTQQ07FErls9mg1gV/8FJKS8Mz7\nASapr9uRJAKY627AzPobAs/9DJ58DM8PFmNS+7kdS+KUY5vfS0tLKSoqAiArK4umpiaam5sBqKmp\noX///gwZMiS0pl5aWupUFJErZluaCTz9ODQ3YWZ9DTNilNuRJIJ4br6T/p+5D6oOEvjZP2Db292O\nJHHKsTV1v99Pbm5u6L7X68Xv95Oamorf78fr9YYeGzhwINXV1WRnZ1NZWcnixYtpbGxk5syZjB9/\n8XNpZ2ZmOjKGyxFJWdwQi+O3bW0cW/RdWir30e+Tn8M38ysXnBwXi+O/FPE8fjv3u7QdqaRl45uk\nrFqJ7xsPx90kynj+/kfK2MN2CqzujuXsfGz48OHMnDmTqVOnUlNTw6JFi3jyySdJvMiZuqqqqno1\n6+XKzMyMmCxuiMXxW2uxv3wKu+VtuG4Sp+76Mi1HjnT53Fgc/6XQ+DM5/ZVvwZHDNP3h1zT3H4Tn\ntrvcjhU28fz9D/fYu/sA4djmd5/Ph9/vD92vr6/H5/N1+VhdXR1paWmkpaVx4403YowhIyODQYMG\nUVdX51REkYuyf3gxeG737NF4vvp9TEKC25EkgpnkVDzfng8D07AvLsfu3OZ2JIkzjpV6YWEhGzZs\nAKCiogKfz0dKSgoA6enpnDp1itraWtrb29myZQvjx49n3bp1rFq1Cghuvm9oaCAtLc2piCLdCrzz\nOvbllZA2FM+3F2CSU9yOJFHApA3F8/WHwWMI/Ns/Yhvq3Y4kccSxze/5+fnk5uYyf/58jDEUFxez\ndu1aUlNTmTx5MnPnzmXp0qUATJ06lczMTHw+H0uXLmXTpk20tbUxd+7ci256F3GC3bkNu3wppPTD\n851HMYP04VJ6zoweh/nc/dgX/4PAv/0jngd/jPFoK484z9HGnD179jn3c3JyQrcLCgrOOcQNICUl\nhYcfftjJSCIXZQ8fDM50Bzzf+DvMiGyXE0k0MrfdhS1/H7a+g131PObuL7kdSeKAzignchbrP07g\npwvhVBNmzncw4y5+9IVIV4wxeP76/8DgdOwfXsLu2eF2JIkDKnWRDralmcBPfwx1xzCf/TKeKTPc\njiRRzqT2x1P8IGAJ/Mc/Y1tOuR1JYpxKXYTgseiBn/0DHNqHuekTmE/e43YkiRFmTAHm43fD0Wrs\nr1e4HUdinEpd4p61Fvv8z2D7FrhuEuaLX4u7k4aIs8xds2HEKOzaP2K3b3E7jsQwlbrEPfvKb7Fv\n/Bmyc3UsujjC9EnC88B3ISGRwC+exDY3uh1JYpRKXeKa3fwW9qUV4BuiY9HFUSY7F/OpL4D/OPa/\nfu52HIlRKnWJW7ZiF4Gf/zP0TQkW+qDBbkeSGGc+cQ9kj8aufw373ka340gMUqlLXLK1VQT+9e+h\nrQ3P1/4WM/IqtyNJHDCJiXge+L/BzfC/fArbdNLtSBJjVOoSd+zxowT+aQGcbMB88auY625wO5LE\nETNiFOYzs6ChDvuf/+Z2HIkxKnWJK7ahnsCSBVB3FHP3l/DMuMPtSBKHzO2fg5wx2A1rsVs3uB1H\nYohKXeKGPV5L4B9/CLVVmE/eg+fOL7gdSeKUSUgInm0usQ+BXz6NbTzhdiSJESp1iQu2cj+BJ/4W\nqg9jbv8c5rNfdjuSxDmTmY2564twwo99/hm340iMUKlLzAu88zqBx78P/jrMzAfw3DNHJ5eRiGA+\nfjfk5mP/9w3s5vVux5EYoFKXmGVP1BNYvhT77/8EHg+eeT/A8/G73Y4lEmI8HZvh+yQR+NUy7MkG\ntyNJlFOpS8yxZ84Q+NOvCTwyD7v+NcjKwfPIEswN09yOJnIek5EVvCzryQbsr/6f23Ekyjl6PXWR\ncLLWwrsbCLy0HI5WQ/8BmC/Ow9x0u079KhHN3Ppp7LtvYze/ReDtNXim3ux2JIlSKnWJCfbQPgL/\n9e+wqxQSEjC3fgbzqfsw/fq7HU3kojo3wwf+/kHsL5/CjsjGZI92O5ZEIW1+l6hmT/gJPPuvBB77\nv8FCv24SnoVP4rl3rgpdoopJz8RT/BC0nSHw9OPavy6XRWvqEpVsoB275o/Y362EU80wfCSeLxRj\nrr3e7Wgil80UFmE+Mwv7u+cILF2E56G/x6Skuh1LoohKXaKO3V9OYOUyOLAHUvsHr39+0ye031xi\ngrnzXjhWi33rVQL/+vd4/s+jmKS+bseSKKFSl6hhm5uwv/sVds0fwAYwU27GzPxrjHeQ29FEeo0x\nBr7yTeypZtiynsDShXi++QgmVbuT5OJU6hLxrLWw+S0C//nv0FAHw0bg+dLXMePGux1NxBHGk4Bn\n7kMEfg5sXk9g8cN4vv4wJiPL7WgS4VTqEtFsTRWB/3wGtm+BxD6Yu76Iuf3zmD593I4m4ijTpw+e\nr34f+8J/YF/7PYG/fxBz71zMtFsxHs1xlq6p1CUi2dOt2D++hP3Tr6GtDQom4PniPMywTLejiYSN\n8SRg7vsbArn52Gefwj77r9g3/xKcFDp6nNvxJAKp1CWi2Pb24OUof/88HK+FQYPx3FsMN0zT+dol\nbnkm34TNuxr74nLspjeDFyfKK8B87FOY8UWYvppIJ0EqdYkItrUF+87r2L+8DNWHITERc/tnMZ+6\nF5OsQ3pETNpQzNf+FnvzHQT++GvYvhm7pwzbNxkzvgiumYgZNx4zON3tqOIiR0t9xYoVlJeXY4xh\nzpw55OXlhR7btm0bzz//PB6Ph4kTJ3LPPfdcdBmJLbbtDOwsxW7dgP3fN4LHmyckBE/reucXMGlD\n3Y4oEnHM2GtJGHsttupg8IPwxnXYjetg4zoswOB0GJWHyc7FjBoNI3PBO0hbuuKEY6VeVlZGdXU1\nJSUlVFZWsmzZMkpKSkKPL1++nEceeYS0tDQWLlzIlClTOHHiRLfLSHSy1sIJPxyrwR6rgcMHsPvL\noWI3tJ4KPmlgWvDUrh+9HeMb7G5gkShgMrMxn/0y9u4vBX+ndpVid7wHe3fClvXYLeuDJQ+Q0g+G\nZQbnpAwbAWlDMAN9MDANBvqgv1eT72KEY6VeWlpKUVERAFlZWTQ1NdHc3Exqaio1NTX079+fIUOG\nADBx4kRKS0s5ceLEBZcJB7ttI3b9amznr4I959Gzbtouv3wsuS/tLS3nP+ecl7nA61zoOT143ws+\nJ8zvW5OURPvp08H7rS3BNe9TzdDSDIHA+e+XkYW59nrMhI8E9w/q5DEil8wYA1k5mKwcuOXTwQ/R\n9cfg4F7sgQps5X6oOQyV+4IfpjvYc18EkpIhORn6pgT/3ycJPB7wJHT8/0O3P7Tmfyw5+YO/f+en\nvIQB9fypwae7vwXiWEoKgVOnLvi4KZoetqtEOlbqfr+f3Nzc0H2v14vf7yc1NRW/34/X6w09NnDg\nQKqrqzl58uQFl+lOZmbvzIiu+837NG1+67KXv/C3NAac/Qt8zi/zB7dPn/2U5BQ8qf3wDB2GSe1H\nQtoQEoeNIGHYcPqMyCZpzDV4+g9wPneY9dbPYrTS+CNk/CNGwLWF53zJtrfTfrSaM4cP0n68lkDd\nMdrrj9Fed4xAg59ASzP21CkCp5qwJ/3Y1lYItPf4Ld36+3eBVZSwutjYUwYMYPCnZ4YlS9gmytkL\nrR1281h3y5ytqqrqsjKd936fvR/PbZ/90FdNlze7KrbhGRkcqam+wLIXeJ0LPecS3rfHr9mD51zJ\nfrfMzMzzvhed6+ftwOmzHzhxMvhfDOlq/PFE44+G8RsYPir434Wfcc6fChsIBMs9EOj4r73LLW8Z\nGRlUV1ef9/VLa91LrOgedoTTLjj2Di0DBvbqz0Z3Hx4dK3Wfz4ff7w/dr6+vx+fzdflYXV0daWlp\nJCYmXnCZcDDGwICBl728Z4AXc7KxFxOJiLjLdG5uv4gE7yBMY3MYEkWehEFpmOYL7XoIL8dmRhQW\nFrJhwwYAKioq8Pl8pKSkAJCens6pU6eora2lvb2dLVu2MH78+G6XERERke45tqaen59Pbm4u8+fP\nxxhDcXExa9euJTU1lcmTJzN37lyWLl0KwNSpU0ObEz68jIiIiPSMo/vUZ8+efc79nJyc0O2CgoIu\nD1f78DIiIiLSMzowUUREJEao1EVERGKESl1ERCRGqNRFRERihEpdREQkRqjURUREYoRKXUREJEao\n1EVERGKEsT29aoqIiIhENK2pi4iIxAiVuoiISIxQqYuIiMQIlbqIiEiMUKmLiIjECJW6iIhIjHD0\neurxaNWqVaxbt47ExESKi4vJy8tzO1LY+f1+vvvd7/K9732Pa665xu04YdPe3s6yZcuoqakhEAjw\n5S9/mXHjxrkdy3ErVqygvLwcYwxz5syJu5/5lStXsmPHDgKBAHfffTcf+chH3I4UdqdPn+ahhx7i\n85//PDNmzHA7TlitW7eOVatW4fF4uPfee7n++utdzaM19V506NAh1q9fzxNPPMHf/M3fsGXLFrcj\nuWLlypWkp6e7HSPs3njjDZKTk3nssceYN28ev/jFL9yO5LiysjKqq6spKSlh3rx5LF++3O1IYbV9\n+3YOHTpESUkJP/zhD1mxYoXbkVzx61//mv79+7sdI+xOnjzJSy+9xI9//GMefvhhNm7c6HYkran3\nps2bNzN16lQSEhLIzc0lNzfX7Uhht337dpKTk8nOznY7Sth99KMfZdq0aQB4vV4aGxtdTuS80tJS\nioqKAMjKyqKpqYnm5mZSU1NdThYeBQUFoS0T/fr1o7W1lUAggMcTP+tLhw8fprKykokTJ7odJexK\nS0u57rrrSElJISUlha997WtuR9Kaem86evQox44do6SkhB//+Mfs37/f7Uhh1dbWxosvvsisWbPc\njuKKxMREkpKSAPif//mfUMHHMr/fj9frDd33er34/X4XE4WXx+MhOTkZgNWrVzNx4sS4KnSAZ599\nlvvvv9/tGK6ora2ltbWVxYsX86Mf/YjS0lK3I2lN/XK99tprrF69+pyv+f1+JkyYwA9/+EN27drF\nz372Mx5//HGXEjqrq/FPmDCBW265hX79+rmUKny6Gv/MmTOZMGECf/rTn9i3bx8/+MEPXErnnng9\n6/TGjRtZvXo18+fPdztKWL3++uuMHTs2Lne3dTp58iTf//73OXr0KIsWLeLpp5/GGONaHpX6Zbrl\nllu45ZZbzvnaCy+8QGZmJsYYxo0bR21trUvpnNfV+BcsWEAgEODPf/4z1dXV7NmzhwcffJCRI0e6\nlNI5XY0fgmtrmzdv5vvf/z6JibH/6+Xz+c5ZM6+vr8fn87mYKPy2bt3Kb37zGx555JG42e3QacuW\nLdTW1rJlyxaOHz9Onz59SEtLY/z48W5HC4uBAweSn59PQkICGRkZpKSkcOLECQYOHOhaptj/qxNG\nEyZM4C9/+QvTp0/n8OHDDBkyxO1IYfXYY4+Fbj/11FPMmDEjJgv9QmpqavjLX/7CwoULQ5vhY11h\nYSEvvPACt912GxUVFfh8PlJSUtyOFTbNzc2sXLmSBQsWxOVEse9+97uh2y+88ALp6elxU+gQ/Pl/\n6qmnuOuuu2hqaqKlpYUBAwa4mkml3ovGjh3L1q1beeSRRwAoLi52OZGE02uvvcbJkyfP2eUyf/78\nmF5jz8/PJzc3l/nz52OMibuf+fXr13Py5En++Z//OfS1b33rW3H3gT5epaWlMWXKlNDf/AceeMD1\nORW69KqIiEiMiK9pmiIiIjFMpS4iIhIjVOoiIiIxQqUuIiISI1TqIiIiMUKlLiIiEiNU6iIiIjFC\npS4iIhIj/j84hTtU01FsBAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fc811421350>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "8MqU1SefgRy5"
      },
      "cell_type": "markdown",
      "source": [
        "Although this seems like a good start, including something about geography is probably even better. I.e.,\n",
        "\n",
        "    radon ~ 1 + floor + county\n",
        "\n",
        "This model posits that the `radon` response is goverened by `offset + weight_floor[floor] + weight_county[county]`, i.e., the same as before except with a county-specific weight."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "rcddRr2Ug1cH"
      },
      "cell_type": "markdown",
      "source": [
        "In the absence of a training set, `radon ~ 1 + floor + county` feels right.  However, if we studied the training data, we'd discover that there's there's a large number of counties with a very small number of measurements. I.e.,"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "15f4k6gQg40_",
        "outputId": "f1f64c6f-001f-499f-dd5c-c48859190025",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 335
        }
      },
      "cell_type": "code",
      "source": [
        "fig, ax = plt.subplots(figsize=(22, 5));\n",
        "county_freq = df['county'].value_counts()\n",
        "county_freq.plot(kind='bar');\n",
        "plt.xlabel('County Code');\n",
        "plt.ylabel('Count');"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAABP8AAAE+CAYAAADyErXtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XtcVVX+//E3CISICIaMYYIJaqGk\nFhqWipm3qZlmprz0ncbRRp3GabTb18ppTM34qjnlZKiTXdTGZgzvNmoXb2WjWOYFxCuYiiKiIooo\nKHB+f/jjDCBnA/vAAbav5+Ph4yF7s85nnb3X2nvtD2vv7Waz2WwCAAAAAAAAYDnutV0BAAAAAAAA\nADWD5B8AAAAAAABgUST/AAAAAAAAAIsi+QcAAAAAAABYFMk/AAAAAAAAwKJI/gEAAAAAAAAW5VGT\nH378+HHNmDFDjzzyiAYMGKCzZ89q7ty5KigokIeHh8aMGSN/f39t2bJFa9eulZubm/r06aPevXvX\nZLUAAAAAAACAm0KNJf/y8vI0f/58dejQwb5s8eLFeuihh3T//ffr888/17///W8NHDhQS5cu1dSp\nU+Xh4aHx48era9eu8vX1Nfz89PT0cpcHBwc7XFcRs2WJWTPliFn3yhHTWjHrU12JWffKEdNaMetT\nXYlZ98oR01ox61Ndb5aY9amuxKx75YhprZgVlQsODi53eY3d9uvp6anx48crICDAvmzkyJGKjo6W\nJPn5+SknJ0cpKSkKCwuTj4+PvLy81K5dOx04cKCmqgUAAAAAAADcNGos+degQQN5eXmVWubt7S13\nd3cVFRXpiy++UPfu3ZWdnS0/Pz/77/j5+Sk7O7umqgUAAAAAAADcNNxsNputJgPEx8fLz89PAwYM\nkCQVFRXp3XffVXBwsAYNGqRvv/1WKSkpGj58uKTrtwYHBgaqT58+NVktAAAAAAAAwPJq9IUf5Zkz\nZ45uu+02DRo0SJIUEBBQaqZfVlaW2rRpU+Hn8My/+hOzPtX1ZolZn+pKzLpXjpjWilmf6krMuleO\nmNaKWZ/qSsy6V46Yda8cMa0Vsz7VlZi1V87lz/wrz5YtW+Th4aHBgwfbl7Vp00apqanKzc1VXl6e\nDh48qLvuusuV1QIAAAAAAAAsqcZm/h05ckQff/yxzpw5owYNGighIUEXLlyQl5eXJk2aJEm6/fbb\nNXLkSD355JOKjY2Vm5ubBg4cKB8fn5qqFgAAAAAAAHDTqLHkX+vWre1JvopER0fb3wIMAAAAAAAA\noHq49LZfAAAAAAAAAK5D8g8AAAAAAACwKJJ/AAAAAAAAgEWR/AMAAAAAAAAsqsZe+OFKhaMetf8/\nrcy6Bu+vdm1lAAAAAAAAgDqCmX8AAAAAAACARZH8AwAAAAAAACyK5B8AAAAAAABgUST/AAAAAAAA\nAIsi+QcAAAAAAABYFMk/AAAAAAAAwKJI/gEAAAAAAAAWRfIPAAAAAAAAsCiSfwAAAAAAAIBFkfwD\nAAAAAAAALIrkHwAAAAAAAGBRJP8AAAAAAAAAiyL5BwAAAAAAAFgUyT8AAAAAAADAokj+AQAAAAAA\nABZF8g8AAAAAAACwKJJ/AAAAAAAAgEWR/AMAAAAAAAAsiuQfAAAAAAAAYFEk/wAAAAAAAACLIvkH\nAAAAAAAAWBTJPwAAAAAAAMCiSP4BAAAAAAAAFuVR2xWoTb/45ECZJf/9edWTd7q2MgAAAAAAAEA1\nY+YfAAAAAAAAYFEk/wAAAAAAAACLIvkHAAAAAAAAWBTJPwAAAAAAAMCiSP4BAAAAAAAAFkXyDwAA\nAAAAALAokn8AAAAAAACARXnU5IcfP35cM2bM0COPPKIBAwbo7NmziouLU1FRkfz9/TVmzBh5enpq\ny5YtWrt2rdzc3NSnTx/17t27JqsFAAAAAAAA3BRqbOZfXl6e5s+frw4dOtiXxcfHq3///nr99dfV\nvHlzbdq0SXl5eVq6dKkmTJigSZMmac2aNbp06VJNVQsAAAAAAAC4adRY8s/T01Pjx49XQECAfVly\ncrKioqIkSVFRUUpMTFRKSorCwsLk4+MjLy8vtWvXTgcOHKipagEAAAAAAAA3DTebzWaryQDx8fHy\n8/PTgAEDNHLkSH3wwQeSpIyMDMXFxWnAgAFKSUnR8OHDJUmLFy9WYGCg+vTpU+kYaY9EOVzXcs0O\nh+u6zNjocN3347j1GAAAAAAAAPVbjT7zryalp6dX6+85Uy44ONh0HLNl61PM+lTXmyVmfaorMete\nOWJaK2Z9qisx6145YlorZn2qKzHrXjli1r1yxLRWzPpUV2LWXrng4OByl7v0bb/e3t66evWqJCkr\nK0sBAQEKCAhQdna2/XeKlwMAAAAAAABwjkuTf5GRkUpISJAkJSQkqFOnTmrTpo1SU1OVm5urvLw8\nHTx4UHfddZcrqwUAAAAAAABYUo3d9nvkyBF9/PHHOnPmjBo0aKCEhASNHTtWs2fP1vr16xUYGKiY\nmBh5eHjoySefVGxsrNzc3DRw4ED5+PjUVLUAAAAAAACAm0aNJf9at26tSZMm3bB8woQJNyyLjo5W\ndHR0TVUFAAAAAAAAuCm59LZfAAAAAAAAAK5D8g8AAAAAAACwKJJ/AAAAAAAAgEWR/AMAAAAAAAAs\niuQfAAAAAAAAYFEk/wAAAAAAAACLIvkHAAAAAAAAWBTJPwAAAAAAAMCiSP4BAAAAAAAAFkXyDwAA\nAAAAALAokn8AAAAAAACARZH8AwAAAAAAACyK5B8AAAAAAABgUST/AAAAAAAAAIsi+QcAAAAAAABY\nFMk/AAAAAAAAwKJI/gEAAAAAAAAWRfIPAAAAAAAAsCiSfwAAAAAAAIBFkfwDAAAAAAAALIrkHwAA\nAAAAAGBRJP8AAAAAAAAAiyL5BwAAAAAAAFgUyT8AAAAAAADAokj+AQAAAAAAABZF8g8AAAAAAACw\nKJJ/AAAAAAAAgEWR/AMAAAAAAAAsiuQfAAAAAAAAYFEk/wAAAAAAAACLIvkHAAAAAAAAWBTJPwAA\nAAAAAMCiSP4BAAAAAAAAFkXyDwAAAAAAALAokn8AAAAAAACARXm4MlheXp7i4uKUm5ura9euaeDA\ngfL399cHH3wgNzc3hYSEaNSoUa6sEgAAAAAAAGBZLk3+bd68WcHBwfr1r3+trKwsvf766woICNDw\n4cMVHh6ud955R7t27VLnzp1dWS0AAAAAAADAklx622/jxo2Vk5MjScrNzZWvr68yMzMVHh4uSbr3\n3nuVlJTkyioBAAAAAAAAluXS5N8DDzygs2fPasyYMZo4caKGDh2qRo0a2dc3adJE58+fd2WVAAAA\nAAAAAMtys9lsNlcF++abb7R//349/fTTOnr0qP7617/Kx8dHb775piQpMTFRmzZt0rPPPlulz017\nJMrhupZrdjhc12XGRofrvh/Xu0p1AAAAAAAAAOoalz7z7+DBg+rYsaMkqVWrVrp69aoKCwvt67Oy\nshQQEFCpz0pPT6/W33OmXHBwsOk4ZsvWp5j1qa43S8z6VFdi1r1yxLRWzPpUV2LWvXLEtFbM+lRX\nYta9csSse+WIaa2Y9amuxKy9csHBweUud+ltv82bN1dKSook6cyZM2rYsKFatGihAwcOSJK+++47\nderUyZVVAgAAAAAAACzLpTP/+vbtqzlz5mjixIkqKirSqFGj5O/vr3nz5slmsyk8PFx33323K6sE\nAAAAAAAAWJZLk3/e3t564YUXblj++uuvu7IaAAAAAAAAwE3Bpbf9AgAAAAAAAHAdkn8AAAAAAACA\nRZH8AwAAAAAAACyK5B8AAAAAAABgUST/AAAAAAAAAIsi+QcAAAAAAABYFMk/AAAAAAAAwKJI/gEA\nAAAAAAAWRfIPAAAAAAAAsCiSfwAAAAAAAIBFkfwDAAAAAAAALIrkHwAAAAAAAGBRJP8AAAAAAAAA\niyL5BwAAAAAAAFgUyT8AAAAAAADAokj+AQAAAAAAABZF8g8AAAAAAACwKJJ/AAAAAAAAgEWR/AMA\nAAAAAAAsiuQfAAAAAAAAYFEk/wAAAAAAAACLIvkHAAAAAAAAWBTJPwAAAAAAAMCiSP4BAAAAAAAA\nFkXyDwAAAAAAALCoSiX/Tp48ecOyQ4cOVXtlAAAAAAAAAFQfw+Rfbm6uMjIyNGfOHGVmZur06dM6\nffq00tPTNXv2bFfVEQAAAAAAAIAJHkYrDx06pDVr1ujo0aOaPHmyfbm7u7s6duxY45UDAAAAAAAA\nYJ5h8q9z587q3LmzvvzyS/Xr189VdQIAAAAAAABQDQyTf8W6du2qtWvX6tKlS7LZbPblQ4YMqbGK\nAQAAAAAAAHBOpV74MXXqVB09elRubm5yd3e3/wMAAAAAAABQd1Vq5p+3t7f++Mc/1nRdAAAAAAAA\nAFSjSk3fa9OmjU6ePFnTdQEAAAAAAABQjSo182/Pnj1as2aN/Pz8St3uO3fu3BqrWF332afZJX7K\nLrXu50P8XVsZAAAAAAAAoByVSv699NJLNV0PAAAAAAAAANWsUsm/pKSkcpf37t27WisDAAAAAAAA\noPpUKvl34MAB+/8LCgp0+PBhtWvXzlTyb8uWLVq9erXc3d01ZMgQhYSEKC4uTkVFRfL399eYMWPk\n6elZ5c8FAAAAAAAAUFqlkn9l3/Sbn5+vOXPmVDlYTk6Oli5dqmnTpikvL0/x8fFKSEhQ//791a1b\nN/3zn//Upk2b1K9fvyp/NgAAAAAAAIDSKvW237JuueUWZWRkVLlcUlKSIiMj1bBhQwUEBOjpp59W\ncnKyoqKiJElRUVFKTEw0UyUAAAAAAAAAZVRq5t9rr70mNzc3+89ZWVkKCQmpcrDMzEzl5+dr+vTp\nys3N1aBBg5Sfn2+/zdfPz0/Z2dkVfAoAAAAAAACAynCz2Wy2in5p3759pX728fFRaGhoqYRgZaxc\nuVIHDhzQuHHjdObMGU2ePFlXr17Vhx9+KEnKyMhQXFyc3njjjSp9btojUQ7XtVyzw+G6LjM2Olz3\n/Tjj5xm+N3Ofw3VPPx9hWBYAAAAAAABwhUrN/IuIiND+/fuVmpoqSWrbtm2VE3+S1KRJE7Vr104N\nGjRQ8+bN1bBhQzVo0EBXr16Vl5eXsrKyFBAQUKnPSk9Pr9bfq65yVS0bHBxsKpbZcrURsz7V9WaJ\nWZ/qSsy6V46Y1opZn+pKzLpXjpjWilmf6krMuleOmHWvHDGtFbM+1ZWYtVcuODi43OWVeubfp59+\nqkWLFun8+fPKysrS/PnztWLFiipXsmPHjtq7d6+KioqUk5OjvLw8RUZGKiEhQZKUkJCgTp06Vflz\nAQAAAAAAANyoUjP/kpOTNWXKFLm7X88VFhYWauLEifrVr35VpWBNmzZVdHS0Xn31VUnS7373O4WF\nhSkuLk7r169XYGCgYmJiqvgVAAAAAAAAAJSnUsk/m81mT/xJUoMGDUzd9itJffv2Vd++fUstmzBh\ngqnPAgAAAAAAAOBYpZJ/d9xxh6ZPn67IyEhJUmJiosLCwmq0YgAAAAAAAACcU2HyLzMzU8OHD9e2\nbdt0+PBhubm5KSIiQo8++qgr6gcAAAAAAADAJMMXfiQlJWnChAnKy8vTAw88oOHDh+vBBx/Ul19+\nqSNHjriqjgAAAAAAAABMMEz+LV26VK+++qp8fHzsy0JCQvTyyy9r8eLFNV45AAAAAAAAAOYZJv9s\nNptCQkJuWN6yZUtdu3atxioFAAAAAAAAwHmGyb+8vDyH63Jycqq9MgAAAAAAAACqj2Hyr2XLlvry\nyy9vWL5q1SqFh4fXWKUAAAAAAAAAOM/wbb9Dhw7VjBkz9M033ygsLExFRUU6ePCgGjZsqPHjx7uq\njgAAAAAAAABMMEz++fv7KzY2VklJSUpLS5O7u7u6deumiIgIV9UPAAAAAAAAgEmGyb9ikZGRioyM\nrOm6AAAAAAAAAKhGhs/8AwAAAAAAAFB/kfwDAAAAAAAALIrkHwAAAAAAAGBRlXrmH6rPrFmzHK4b\nO3asC2sCAAAAAAAAq2PmHwAAAAAAAGBRJP8AAAAAAAAAiyL5BwAAAAAAAFgUyT8AAAAAAADAokj+\nAQAAAAAAABZF8g8AAAAAAACwKJJ/AAAAAAAAgEWR/AMAAAAAAAAsiuQfAAAAAAAAYFEk/wAAAAAA\nAACLIvkHAAAAAAAAWBTJPwAAAAAAAMCiPGq7AqicoJTxpX4uSJGCSvycGT7VtRUCAAAAAABAncfM\nPwAAAAAAAMCiSP4BAAAAAAAAFkXyDwAAAAAAALAokn8AAAAAAACARZH8AwAAAAAAACyK5B8AAAAA\nAABgUST/AAAAAAAAAIsi+QcAAAAAAABYFMk/AAAAAAAAwKI8aiPo1atX9eKLL+rxxx9Xhw4dFBcX\np6KiIvn7+2vMmDHy9PSsjWoBAAAAAAAAllIrM/+WLVsmX19fSVJ8fLz69++v119/Xc2bN9emTZtq\no0oAAAAAAACA5bg8+Xfy5EmdOHFCnTt3liQlJycrKipKkhQVFaXExERXVwkAAAAAAACwJDebzWZz\nZcCpU6dqxIgR2rx5s4KCgrRo0SJ98MEHkqSMjAzFxcXpjTfeqNJnpj0S5XBdyzU7HK7rMmOjw3Xf\nj+ttGPO9mfscrnv6+QiH61555RWH66ZNm+ZwXcE3wwzr49FzoeF6AAAAAAAA3Hxc+sy/r7/+Wm3b\ntlVQUJDTn5Wenl6tv1dd5WoqZkVbrCoxg4ODTdXR1eWIWffKEdNaMetTXYlZ98oR01ox61NdiVn3\nyhHTWjHrU11vlpj1qa7ErHvliGmtmBWVCw4OLne5S5N/O3fuVGZmpnbu3Klz587J09NT3t7eunr1\nqry8vJSVlaWAgABXVgkAAAAAAACwLJcm/55//nn7/+Pj4xUUFKSDBw8qISFBPXv2VEJCgjp16uTK\nKgEAAAAAAACW5dLkX3kGDx6suLg4rV+/XoGBgYqJiantKgEAAAAAAACWUGvJv8GDB9v/P2HChNqq\nBgAAAAAAAGBZ7rVdAQAAAAAAAAA1g+QfAAAAAAAAYFEk/wAAAAAAAACLIvkHAAAAAAAAWBTJPwAA\nAAAAAMCiSP4BAAAAAAAAFkXyDwAAAAAAALAokn8AAAAAAACARZH8AwAAAAAAACyK5B8AAAAAAABg\nUST/AAAAAAAAAIsi+QcAAAAAAABYlEdtVwA179PkoaUXJP/3v0Pa/8O1lQEAAAAAAIDLMPMPAAAA\nAAAAsCiSfwAAAAAAAIBFkfwDAAAAAAAALIrkHwAAAAAAAGBRJP8AAAAAAAAAiyL5BwAAAAAAAFgU\nyT8AAAAAAADAokj+AQAAAAAAABZF8g8AAAAAAACwKJJ/AAAAAAAAgEWR/AMAAAAAAAAsiuQfAAAA\nAAAAYFEk/wAAAAAAAACLIvkHAAAAAAAAWBTJPwAAAAAAAMCiPGq7Aqi7gncnlV6wO0nBJX5M7xTp\n0voAAAAAAACgapj5BwAAAAAAAFgUyT8AAAAAAADAokj+AQAAAAAAABbFM/9QIwpHPWr/f1qZdQ3e\nX+3aygAAAAAAANykmPkHAAAAAAAAWBTJPwAAAAAAAMCiXH7b76JFi7R//34VFRXpl7/8pcLCwhQX\nF6eioiL5+/trzJgx8vT0dHW1AAAAAAAAAMtxafJv7969SktLU2xsrHJycvTSSy8pMjJS/fv3V7du\n3fTPf/5TmzZtUr9+/VxZLQAAAAAAAMCSXHrbb0REhJ5//nlJUqNGjZSfn6/k5GRFRUVJkqKiopSY\nmOjKKgEAAAAAAACW5dLkn7u7u7y9vSVJGzduVOfOnZWfn2+/zdfPz0/Z2dmurBIAAAAAAABgWW42\nm83m6qDff/+9VqxYob/85S8aO3asPvjgA0lSRkaG4uLi9MYbb1Tp89IeiXK4ruWaHQ7XdZmx0XEd\nx/U2jPnezH0O1z39fITDda+88orDddOmTXO4ruCbYYb18ei50OG6mV895HDd8303OP7QtV8YxtTD\n/R2uMrtPAAAAAAAAUH1c/sKP3bt3a/ny5Xr11Vfl4+Mjb29vXb16VV5eXsrKylJAQEClPic9Pb1a\nf6+6ytVUzKBaiBlcCzFvqENwsKk4ZsvdLDHrU12JWffKEdNaMetTXYlZ98oR01ox61NdiVn3yhGz\n7pUjprVi1qe6ErP2ygUHl5/Jceltv5cvX9aiRYv0yiuvyNfXV5IUGRmphIQESVJCQoI6derkyioB\nAAAAAAAAluXSmX9bt25VTk6OZs6caV/2zDPP6O9//7vWr1+vwMBAxcTEuLJKAAAAAAAAgGW5NPnX\np08f9enT54blEyZMcGU1AAAAAAAAgJuCS2/7BQAAAAAAAOA6JP8AAAAAAAAAiyL5BwAAAAAAAFgU\nyT8AAAAAAADAokj+AQAAAAAAABZF8g8AAAAAAACwKJJ/AAAAAAAAgEWR/AMAAAAAAAAsyqO2KwCU\nVDjq0VI/p5X4f4P3VxuW/cUnB0r8dKDUulVP3umw3GefZpdZ8t+ffz7E3zDmrFmzHK4bO3asYVkA\nAAAAAICaxsw/AAAAAAAAwKJI/gEAAAAAAAAWxW2/QC0IShlf6ueCFCno//8/M3yqw3KfJg8tvSC5\n9I9D2v/DYdng3Un//WF3koJLrEvvFOmwnNGt2JLx7dgly1alHAAAAAAAqB7M/AMAAAAAAAAsipl/\nAOqc0i9vkUq+wMXo5S0AAAAAAKA0Zv4BAAAAAAAAFsXMPwCW8tmn2SV+yi617udD/B2WmzVrlsN1\nY8eOdbZaAAAAAADUCmb+AQAAAAAAABbFzD8AcFLJtzeXfHOzZP7tzUZvbgYAAAAAoLKY+QcAAAAA\nAABYFDP/AKCeCd6dVHrB7iQFl/gxvVOkw7KFox61/z+tzLoG76+uVLmqli399ubSb3I2entz6ec3\nSiWf4Wj0/EYAAAAAwH8x8w8AAAAAAACwKJJ/AAAAAAAAgEVx2y8AwHJmzZrlcN3YsWMdriv58hap\n9AtcjF7eIpV5gUty6XW8wAUAAABAbWHmHwAAAAAAAGBRzPwDAKAWGb3ApbIvb5HMv8DF/MtbpKq8\nwAUAAABA7WDmHwAAAAAAAGBRzPwDAAA17rNPs0v8lF1q3c+H+DssZ/b5jVLpZziWfH6jZPwMx1LP\nb5RKPcPR6PmNRrM4JeOZnAAAAEBNYeYfAAAAAAAAYFHM/AMAAKhlRs9wrOzzG8uWq6hs6Wc4Vv75\njaVncUolZ3IazeKUeBM3AABAbWDmHwAAAAAAAGBRJP8AAAAAAAAAi+K2XwAAAKAMoxe4GL28xZlb\nsc3e/l36Fm7J/G3c1nwZj1Rmf1bhZTy1f0u+VHJ/mr0lX6qZ/Vkbt+Sb7ZsAcDNj5h8AAAAAAABg\nUcz8AwAAAABYXskZmTfjrFyp8jM5rTort2xZV8zKBeoCZv4BAAAAAAAAFsXMPwAAAAAAgBpQ157H\nKVVhJqeLnsfp6lm5UtmZnJWflVtf1Znk34IFC3T48GG5ublp+PDhCg8Pr+0qAQAAAAAAAPVanUj+\n7du3TxkZGYqNjdWJEyc0d+5cxcbG1na1AAAAAAAAAEnGMznr2vM4S6oTz/xLSkpSly5dJEm33367\ncnNzdfny5VquFQAAAAAAAFC/1YnkX3Z2tvz8/Ow/+/n5KTu7bDYVAAAAAAAAQFW42Ww2W21X4r33\n3tM999xjn/03YcIEjR49WsHBwRWUBAAAAAAAAOBInZj5FxAQUGqm3/nz5xUQEFCLNQIAAAAAAADq\nvzqR/OvYsaMSEhIkSUeOHFFAQIAaNmxYy7UCAAAAAAAA6rc6cduvJH3yySfav3+/3NzcNGLECLVq\n1aq2qwQAAAAAAADUa3Um+QcAAAAAAACgetWJ234BAAAAAAAAVD+SfwAAAAAAAIBFkfwDAAAAAAAA\nLIrkXx2Qm5tb21XA/3fx4sXargJuUjx+FfUNx8vy2Ww2Xbx4URcuXKjtqgBVtnfvXpfEycvLU0ZG\nhjIyMpSXl+eSmEBJhYWFOnPmjAoLC2u7KqgF5Y27z507Vws1AVynwaRJkybVdiXMKioq0tatW/Xv\nf/9b69ev17fffquDBw/Kzc1Nt912m8Nyu3btsq+/dOmSFi1apNWrV+vIkSMKCwvTLbfc4rDs5cuX\ntXbtWqWnpys0NFSff/65vvzyS504cUKhoaHy9PSs8veYOnWqevXqZfg7hw4d0qZNm7R161bt3LlT\naWlp8vX1lZ+fn8MyBQUF2rBhg1atWqW1a9dq06ZN2r17t65cuaLQ0FC5ubk5LLt792599dVX+vrr\nr7V9+3YdPnxY7u7uCgoKMqxnTk6O1q1bp8zMTIWGhmrFihVatWqVjhw5ojvuuMNw25Y1efLkCreL\nZH5/7tq1S5999pnuvfde7d27V5MmTdJ//vMfrV69WrfddpthGzL7PS9fvqzExEQFBwcrNzdX8fHx\nWrt2bY1uH7NttqCgQNu2bVNubq6aNWum//znP/rqq690+vRphYaGqkGDBg5jXr16VRs2bNCePXvU\npEkTNW7c2L5u2bJlioiIcFh23rx58vf3V9OmTSv8biWZPR5UZ5stac2aNWrbtm256/bs2aPp06dr\n27ZtatmypWbMmKElS5Zow4YNCgkJUbNmzRx+rtl+bXa7OrN9Lly4oKVLl2rr1q3y9vYudfz48MMP\ndc8995Rbzmw/ceaYV57K9jFXb9uCggL9+OOPatq0qQoKCrR+/Xpt3ry5Un2zPJX5njt37tS0adO0\nY8cO3XHHHXrjjTe0adOmSh0vze5PZ46XZtuekU8++UR33323w/Xp6emaPXu2Vq1apSVLligxMVGr\nVq3SwYMH1a5dOzVs2LDccmbbjzPjEbN9xZnxkxGjbetMv165cqWCgoIcbntHqvNcXayi9mO27LBh\nw3T+/Hm1bt1a3t7eVfrcr7/+WseOHbP/O3r0qBYtWiQfHx8dO3ZMrVq1qtLnVeZYkpqaqrffflvr\n1q1TUlKSvv/+e61evVrbt29XWFiYmjRpUm45Z/p0YWGhEhISdOHCBQUFBemHH37Qt99+q6ysLLVs\n2dJhG3JmHGS2f5o9HjjTT3788UcFBARIkq5du6Z169Zp/fr1ysjI0B133OHwezpzDDI7ZjO7fSRp\n/vz56ty5syQpMTFRU6ZM0e5cETOOAAAgAElEQVTdu7Vy5Uq1aNFCzZs3d1jW7HWR2e9ptr07c4x2\n5tqmOq+N33nnHUVHRxv+jjPHg++++06xsbFasWKF0tLSFBkZaa/jjBkzKjXuK6ui47vZ66LqPOdW\ndkzr6jbrzPm2uscklWl7ZrePEWfGB0bXm+Wp18m/999/XzabTffff78KCgrUtGlTRUZGavPmzTp0\n6JA6duxYbrl58+bZG//cuXN122236Wc/+5ny8/O1evVqde/e3WHMt99+W02bNlVWVpaWLFmiRo0a\nqUuXLjp//ry++uor3X///eWW++KLL5SamnrDv5SUFCUlJemnP/2pw5iLFy/W/v371aZNG2VkZKhR\no0YKCAjQ8uXLlZub63CHz5o1S+7u7urVq5d69uypLl26qEWLFkpOTtbWrVvVtWvXcst98MEHSklJ\nUfv27XXXXXepVatW8vLy0ubNm5WamupwuxZvn4CAAKWlpWnjxo3y9vZWv379lJ+fr1WrVqlHjx7l\nlhsyZIg2bNigzz//XGvWrNHatWt18uRJbdy4UevWrdPDDz/sMKbZ/Tlr1iw9/fTT8vb21uzZszVu\n3Dg99thj6tmzp95//3099NBD1f49p02bpmbNmql169aaO3euAgMD1bt3bxUUFNTY9jHbZt955x2d\nPXtWe/bs0eHDh3Xs2DF16tRJaWlp2rx5s7p16+Yw5syZM+Xm5qaGDRvqX//6lxo2bKjQ0FBJ0pIl\nSwxPPosXL1ZWVpa++eYb+fv7GybCSjJ7PDC7Lyvyr3/9y+H3nD17tl555RV16NBBb7zxhl588UX9\n5je/0b333quPPvpIDz74oMPPNduvzW5XZ7bPX//6V7Vp00ahoaH67LPPdPbsWfsA57PPPnO4fcz2\nE7PbRnKuj7l6286cOVN5eXmKiIjQhx9+qPPnzysyMlKnTp3Sxo0bDfum2e8ZFxenl19+WXfccYfe\nfvttPffccxoyZIgeeOABzZs3T3369HEY0+z+NFtOMt/28vPzVVhYWO6/lStXGvbNt99+W6NGjdJj\njz2m++67T+fOndNf/vIXNWzYUAsXLnQY05n2Y+bYLpnvK86Mn8xuW2f69dy5c3Xw4EGlpqYqODhY\njRo1cvi7JZlte860H7NlExMT1b9/fy1YsEDJyclq1KiRmjZtKnf3im/ueeutt5SWlqYmTZooNzdX\nly9fVmpqqoKCgpSbm6v27ds7LGv2WDJz5kyNHj1agwYN0oMPPqgHH3xQDz/8sEJDQzV//nyH39Ns\nn5auH7+ysrJ04MABbdu2TcePH1erVq108OBBff/99+rSpUu55ZwZB5ntn2aPB870k1mzZtm330cf\nfaTc3Fx17dpVaWlp2rJli+67775q/Y6S+TGb2e0jSatXr7Z/zzlz5mjcuHH6xS9+oQceeMBw3O/M\ndZHZ72m2vTt7jWtmTOJMO3jmmWe0bt06rV271n4cSU1NrfA44szxYPbs2Xr99df1+OOPKycnRx9/\n/LG6desmT09Pff311zUyPjB7XWR2fzozpnV1m3VmrOdMezfb9sxuH2fajxGj683yeJiKUkecOnVK\nv//97yVJ4eHhmjJligYOHKiOHTvqz3/+c6U+48KFC/rlL38pSbr99tu1bds2w9/Py8vTr371K0nS\n888/r9/85jeSpI4dO2ry5MkOy/373/9WZGSk/S9rJVU03fzAgQMqztF269ZNsbGxevLJJ9W7d2/9\n+c9/1s9+9rNyy50/f17PPfdcqWXNmzdXRESEJk6c6DDe8ePH9frrr9+wPCYmRq+99pphXa9evaqB\nAwfKZrPpueee07hx4yRd3z8JCQkOy40fP16rVq1S//797Rn3V199VbGxsYbxyqrK/iwoKLDPBmjU\nqJH9LxS+vr4V3oJp9nteuXLFPrjIzs7Ws88+K0kKCwvTli1bHJZzZvuYbbMXL17UxIkTVVhYqDFj\nxiguLk7u7u7q1q2bYfuRrt/KXhynX79+mjFjhoqKiip1cPL19dXo0aOVnp6utWvXav78+QoPD1do\naKiaNGnicLBt9nhgdl9K0siRI8tdbrPZdOXKFYflPDw8FBAQoICAADVq1Mg+AGjWrFmFF2xm+7XZ\n7erM9ikoKFD//v0lSffdd5/i4uK0dOlSDRw40LCc2X5idttIzvUxV2/brKwsvfjii5KkEydO2Ptx\ndHR0jX1PT09PBQYGKjAwUE2bNrXPCPL395eXl5dhWbP702w5yXzbGz58+A0zStzc3GSz2Sq8jbeg\noEA/+clPJF1vd8ePH5ckderUSfHx8Q7LmW0/Zo/tknN9pVhVx09mt60zdW3atKnGjx+vpKQkzZ8/\nX1evXlX79u3VqlUr+fn5KTw8vNxyZtueM+3HmbIRERGKiIhQamqqNmzYoHnz5qlhw4by8/PT+PHj\nHZZ76623tGzZMh09elTDhg1Ts2bNtHv3bg0aNMgwnmT+WGKz2dSiRYsblrdu3VpFRUUOy5nt09L1\nNlTcVsaMGaN3331XktSnT58aGweZ7Z9mjwfV0ael0ueUzp07y2iOiDPHILNjNrPbp6xGjRrZj9f+\n/v7y8HB8SezMdZHZ7+lMey9W1WO02TGJM+3gl7/8pb777jsNGTLEfjyuzHHEme3j7u4uX19fSdeP\nAU2aNFFsbKxefvllw3LOHKOduS4qVpX96cyY1tVt1pmxXklVbe9m257Z7eNM+zF7vVmeep38Kyoq\n0p49exQWFqadO3fap+zu3LnTsFxOTo79dzw8PHTs2DGFhoYqMzOzwueOFBYWKiMjQxcvXtSlS5d0\n6NAhtW3bVidPnlRBQYHDcuPGjdP8+fP11FNP3TD9ed++fYYxr127pvT0dAUHB+vAgQP2wdGJEycM\nk1Tu7u7avn277r33XvtJ7dq1a0pISDCcgm2z2XTkyBG1bt261PLiKa1Gip+f0axZMz311FP25UeP\nHjVMcnbq1Ent27fXihUrtGXLFg0bNqzSt+iV3J+enp6V3p+PPvqoXnrpJd19993y9fXVjBkz1LZt\nWyUnJ6t3796mvuexY8cMv2fz5s21YMEC9ejRQ+3bt9fWrVsVERGh3bt3y9/f32E5Z7aPozZ74sQJ\nwzZ77do15eXlydvbW4MHD7YnpbKzs3Xt2jXDmEVFRfY25O3trXHjxmnGjBnKysoyjFlScHCwRo4c\nqYKCAu3bt0+pqak6deqUw4Gd2eOB2TYrSQ8++KCaNWumfv363bDOaMDTqFEj/etf/1JOTo6aN2+u\nefPmqVOnTjp06JDDW56Kubu7KyEhQVFRUaX69fbt2yt1a0VVt6vZti5JDRo0UEJCgu677z65u7vr\nT3/6k+bMmaP33nvP8GRltp+YPeZJzvWxYq7atr6+vlq7dq26d++uu+++WykpKQoPD1dycnKFiTiz\n37NJkyZavXq1Hn30UU2ZMkXS9WfjfPbZZ7r11lsNy5bdn9u2bdNdd91V4f402w4k821v6NChunjx\nop544okb1lV0EdOyZUv97W9/U3h4uPbs2WOfNTV37lzdfvvthmUlc+3HzLFdcnwcqaivmD3fSua3\nraN+XdljniRFRkYqMjJSWVlZ2rVrl3bs2KHs7Gy99NJL5f6+2TbrTPtxpmyxsLAwhYWFSbqeDDp/\n/rzh73t5eel//ud/lJ6erg8//FBt27at9PNnzR5L2rRpo+nTp6tLly72R9dkZ2dr+/btuuuuuxyW\nc9Sn//73v1d48VM8nrl06ZJyc3OVmZmpoKAg5eTk6OrVqxWWMzMOcqZ/SlU/Hjhz/svLy9OJEyck\nSX5+fvbtc/nyZcN+7cx3NDtmK1bV7SNJaWlpevvttyVJGRkZ2rZtm7p166bPPvtMPj4+Dss5ui46\ncOBAhW3e7Pc0ew5z5hhtdjxs9tpYkvr27av77rtPixYt0ubNm/XrX/+6UscRZ44Hd955p6ZNm6YX\nXnhBXl5e6tKlizw9PTVlyhTl5OQ4LOfMMdrsdZHZ/enMmNbVbdaZsZ4z7d1s2zO7fZxpP2avN8vj\nZqvHT5lPT0/XokWLlJGRoZCQEP32t7+Vh4eHPv/8c/Xo0cPhfddz5syRdH3nubu7q2fPnurQoYPe\neustdenSRT179nQYc9++ffrHP/4hPz8//fa3v9VHH32k48ePy9/fX6NGjTK85zo/P1+enp43zOwp\n74RS0qFDh/T+++8rOztbQUFBGj16tG6//XbNnTtXvXv3Vrt27cotd+7cOX366afat2+f8vPzJUne\n3t6KjIzU4MGDHXaoY8eOacGCBTpz5oz9FpmcnBy1aNFCw4YNM7yIOXTokNasWaPnn39e0vUTwuef\nf64NGzboT3/6k+H3lK6fYDMyMrRw4UKlp6dr1qxZOnfunOHFZfH+LNazZ0+FhIRo3rx5evjhhw2f\nL3fp0iUlJibqzJkzstls8vf31913313hs0TKfk/p+jMklixZoj/84Q/2wXdZhYWF2rBhg3744Qdl\nZmZKup7179mzpx5++OEKL9yl639xKLl9KrJ//359/PHHpdrs0aNH1bRpU8M2u2PHDq1bt04TJkyw\nL9u1a5fef/99/f73v1enTp0cxjx+/Ljmz5+vl19+2f4MooKCAq1evVpffPGF3nvvPYdl3377bb3w\nwgvlrsvNzXV421bZ48HQoUN16623Kj4+XpGRkQ4vKsrbl9u2bdPy5cs1evRowzZrs9m0atUqDRgw\n4IZnLX300Uf63e9+V265vLw8bd68WX5+frr//vu1ZcsWHTx4UM2bN1efPn0Mn9tU3K/3799vP7EV\n9+vHHntMgYGB5ZYrb7vu3btXHTp0cBhLMm7rI0aM0J133mlY18WLF2vUqFH2tp2UlKQLFy7oiy++\nsCeRyjLbT4y2jdExr6xTp05pwYIFSkxM1CeffFLhbMzytu3FixcNn8kqlb9tExIStHz5cv3hD39w\n2PYuX76s5cuX64cffrBfhAYGBqpDhw4aPHhwhXFLfs+FCxfq5MmT9pkwjly9elU7duwodfvOkSNH\ntH//fvXt29fw2FVyf5Y81nbq1Endu3d3eHwv2w4KCgpUVFSk3r176+c//7lhzLNnz+rTTz8t1fYk\nacuWLVq5cqXeeusth2W/+eYbde3a9YZ+WNFfsG02m77//nv7Maj4GHns2DGFhIQ4HFQWtx+bzXbD\n7xid/8o7ticnJ6tly5YVjkfM9pXyzreS9OWXX2rAgAGG51vpxm1b3E+WL1+uxx57zLCu5Y1lBg0a\nVO4dFcUmT55cpVlPxco7BhW32ZiYGMNjiZnv6KhssRUrVthn1pS1cePGCv9gWVnffPONdu7cecMM\nsoqkp6fr448/to9LjM7V0vWx9N69e5WdnS1JCggIUMeOHSvVZkeOHGnv0xcvXlRiYqLh+USStm/f\nrvnz56tx48Z66qmnNH/+fEnXx4AjRoxQVFRUueXKGwft3r1b7733np5++mnDcVDZ/vnhhx8qLS1N\nfn5+evrppx1+17Lnk8LCQmVlZalp06aGzxgs20+Kiop06dIlPfjgg3riiScM22zZi8eHH35YXbp0\n0ZQpU9S3b1+Hz8By5hhUPGY7ffq0WrZsWeoarnv37goODq7U9pEqd76Vbpxs0bx5czVt2lTffvut\noqKiHI69jh07poULFyozM/OG66Lhw4eXO5O17Pcsb2x6zz33OJx9XHb8VNwO9u3bp/Xr1zts72WP\n0d27d9dtt92mhQsXVnhNZHY8XNwOfH19NXz48FLXxhWNE0tKTk5WfHy8Ll68qJkzZxr+rjPHg+JY\nERERpc65ly9f1tatWw0fY1LyGF2yb65atcrw+F7edVFRUZFWrlxpeF1U3jm3OGfRo0cPw9v5i1Vl\nrCfd2GZ/85vfyGazaf369YqJiXHYN8vbJ9L1cVd8fLzD2OWdbxs3bqx7771XP/3pTw3Hes7kAEra\nt2+f4uPjdf78eb3zzjuGv1vesSs9PV379++39zdHyp7ji6/DKhofFF9v/vSnP7U/x7B4zGh0vVme\nej3z7/Tp0/L399dLL72kvXv32p+tk5eXp/DwcIcbPyoqSgsWLFB+fr46d+5sT9S8+OKLmjx5smHy\nLyIiQlOnTrX/XHJAUNEg09FDJytKiMXGxiomJkaPP/54qVlBo0ePNix3/PhxeXh4KC4uTnv37tWc\nOXPk5eWlPXv2KCoqyuGDN1977TXFxMToT3/6k/2i18/Pr1IPkf/Pf/5jP2kkJibq73//u/z9/ZWf\nn2/4Zsjvvvuu1D559tln7X8liIuLM9yu0dHR2rFjh37/+9/bv2dxO6jorz++vr6Gz6JwpG3btjcM\narp27aquXbsatoPExEQdPXpU48ePL1XXr776SiEhIQ73yddff13q5+joaC1btsy+PCYmxmFdMzMz\n7dOwU1JS1KNHD2VmZupnP/uZTp065XBw5unpqaysLE2cOFFDhw7Vhx9+qNOnT6tRo0YV/hX5/Pnz\nys7O1tSpU+1lz58/L29vb40dO9awrKPEn3T9eRKOtu2JEyd09OhR5efny8PDw/5X3MGDBxvuk0uX\nLtkHc2XbT/GFiSPDhw9XTEyM8vPzbxg0Gh2Ivb29NWDAAPvPPXr0qPSzBVNTU7V37157XxkxYoT9\n9nWj73nPPfeUakc2m03Lly/X448/LslxG8rOztbhw4c1YsQIe7zKtHXpvye14qn3JWMaDa727Nmj\nNWvWKDAwUM8//7zeffddFRYWVthP9u7dq/bt25c62S9btkzt2rXTnj17DPtJ8cxsSTpz5ozS0tIU\nFBSkMWPGaNSoUYYXeb169dK8efPKPQaNGDHCYX1PnTpVar+U3D5paWkOzw2jR49WTEyMJk2aVOFM\n0bJ27typhQsXKjAwUMOGDdO5c+dks9n0zDPPGNbVy8vrhmNl69atKzx/SdIPP/yglStXVrnNfvzx\nx3rqqafUr1+/UueTjRs3KiwszHCfBAYG6plnnrlheY8ePbRx40aH5b777jstXrxYCxcuvKGuycnJ\nhsk/Nze3cgfixbf1O9K9e3f98Y9/LHf7GJ3/rly5osuXL8vHx0eFhYW6ePGi/YHZly5dMox57Ngx\nHTx4UEFBQRo2bJhmzZqloqIi7d6923B8UDybsXjQefbsWXubPXPmjGHMXbt26cCBA+rZs+cN/cTo\neHnrrbfqj3/8Y7nrKjoGHT58WB999NEN46eKNGjQQP369XP4V3ZHMXfu3Gk/Nw8bNsx+7MrPz9eI\nESMqjOto/JmYmOgw+Vddib/i+EZjYEeCg4P1yiuvKCsrS5LxuVr6723KVVHcZqdPn27ftkVFRfbj\nrJHZs2erV69e9nbw5ptv6uLFi2rcuLHhH3eKn4V34cIFe/tp3769Zs+eXeEfhRISEuzXC4mJicrI\nyFCzZs104cIFXb582WG5ksnskse9CxcuGJ6LVq9ebe8nxeV+8pOfKCkpSdHR0YbHy9TU1FLbp1jJ\na5zyXLlyRXfccYf93Hfq1Cm1aNFCly9fNpxBJV2feXfy5EkFBgbqscce09SpU+19JTw83GGCoVev\nXnr22Wft57CyfczoRQ+O2pzRc8Gk6+Og8+fPKyAgQEOHDrXHzMzMVHZ2tmHyLzg4uNwZxhWNTVev\nXm0/hxXvz4CAAGVnZ2vUqFEO4zVs2NA+linbfoxmuUrXx8NHjx7VlClTSm3bisbDOTk5unDhgjIz\nM7VixQr97//+b6XO8dL1FxYVX+O2b99eEyZM0NGjRw3rKf33GtfLy6tK55Ni5T3P1MfHx3BsWnKc\nWF7fNDJhwoQbrhfc3d312GOPGSZ9HJ37inMWjpJ/O3futF8bnzt3TsePH5e7u3uFYz3p+nsKitts\nYmKiXnvtNfv3vOuuuxz2zeLrk9GjR5cay1Q07vr22291yy23lJq1u3z5cvn7+2vbtm2GY/fo6OhS\nY9qSbbaimX979uzRggUL5Ofnp6FDh+rKlSs6d+6cnnvuOY0aNcrhM28PHz6s++67zz4OSkpK0rJl\ny/T444/r0KFDDvNPxWP97du3Syo95q/oDpqkpCR9/fXX2rVrl/26OisrSw0bNqyw7ZVVr5N/8fHx\n9vvzlyxZookTJ+onP/mJcnJyNG3aNIcNe+XKlXrzzTfl4+OjjRs36o033tCrr75qOOW72BdffOFw\nXfGgp7q1bt1a0dHRmjVrlm699Vb16tVL7dq1qzAZZ3b7FMebM2eOPV5lZ5IUP+NIun7RXRwzOztb\nM2bMcDj4cLRPKmPJkiWG3/Pee++t1OdUhdl2YHafLF26VI0bN9Y999xjvyWnoKDA/hcSI47KVnSh\ntnTpUk2cOFE5OTmaNGmSXnvtNYWGhurMmTN69913y33+SXWUdbRtbTab4bY1268raj9GJ0izfdMZ\nZr/nsmXLbmgH165dq7ANOXO8LK/tFcc0mlq/bNkyTZgwQWfPntX06dP10ksvKTQ0VNnZ2Zo+fXqV\n+0nxbDMjZo9dkvk25Gj7VFRfZ9qd2W3rDLNtyJl9YvYY7Ux7N8uZPm12X1ZnHytusxWp6Pzn6Fzt\nzLirTZs2pvqK2Zhlt+u4cePUqlWrSu0Ts+e+2uLq+jqzbcPCwsptBxUl8MyWk8wfv1xdzuh7VtRP\nnBl/m92fzrQDs8qOaSdOnKiQkBCnxrSS8bGkNtqB2W3rzHnT0XimImbPJ85wZtuaHbeZbT/OXNtU\n97V8Re3A6FqhIs6Mg5YuXarXXntNly5dqlK/Njt2N3tNJF3fhyXrWpXr6rLqdfLP7AsbzD7oU3Lu\nxR3OMPMw57Lbp/iNWJV5oYXZh0eXVJWH6TqzT5x5cYdZZtuB2X3izEO5HZWt6AGsHh4e8vf3l7+/\nf5VfSuFMWbPb1mwbcrb9VEdfqQqz39NsG3Kmb5qNWfblEsXtp6KXS5ht62VV5dglmW9DztTXbLsz\nu22d4UwbKlbVfeLq44gzzMZ0Zl/WRh8z20+cHXeZ6StmYzrzYpzaGl+a5er6OrNtJfPHTFePh2ur\nnKuvM8zuT2fbgRllx7QhISGSanZMW5Kr2oHZbevsebM62p4rrv1KclUfc/X1ZlmuuJZ35hrXmXFQ\n2ZcvVrZfmx0HOfM9nXlR5A2fVaXfrmPMvrDB7IM+Jede3FEdqvIw57Lb569//WulX2hhJp5k/mG6\nzuwTZ17cYZbZdmB2nzjzUG6zZZ15KYUzZc1uW7NtqLraT1X7illmv6fZduBM3zQb0+zLJZzpJ2aP\nXZL5NuRMfYtVtd058+IOs8y2IWf2iauPI84wG9OZfVkbfcxsP6mucVdV+orZmM7sk9oeX1aVq+tb\nXccus+dqV42HXV3Ome/pzHWG2f1ZG+ew2hjT1kY7MLttq+u86Uzbc8W1X230MVdfb0quv5Z3Zlzh\nzPHAbL82W19nvqczx6Cy6vULPyTzL2ww+6BPyfyLO8xy5mHOZraPM/HMPkxXcm6fmG0HzjDbDqqj\nrmYfyl2Vss68lMKZspL5bWu2DZndJ9X5oPWqcKavFKtKG6qOeFWJ6czLJczEk5w7dkmu7dfOtLvq\n2rZVZaYNObtPXH0ccYaZmM7sy9roY5L5fmJ2XzrTV8zEdHa7unp86SxX1teZbWu2HdTGeNjV5STX\nX2dI5vdnbZzDamNMWxvtwJlta/a8WRttz6za6mOuvt6srWv5YlUZVzjTZp3t12bqa7ZcddVVskDy\nDwAAAAAAAED5qnaTMAAAAAAAAIB6g+QfAAAAAAAAYFH1+oUfAAAAuP7g8H/84x9KS0uzP/9l0KBB\nuvvuu6s1zs6dO9W2bVv7W/0qkp6erkWLFuns2bO65ZZb5OXlpSeffLJKz4SbMGGCnnjiCbVv395s\ntQEAAG5qJP8AAADqMZvNphkzZqhnz54aO3asJOn48eOaMmWKpkyZoubNm1dbrDVr1ig4OLhSyb+r\nV6/q//7v/zR8+HBFRUVJuv5A8KlTp+qdd96p0hsSAQAAYB7JPwAAgHosKSlJbm5uGjBggH1ZSEiI\nZs6cKV9fXxUVFWnBggU6cuSI3Nzc1L59ez3xxBNKTk7W4sWLNWXKFEnS7NmzdeeddyoyMlJvvvmm\nOnbsqJSUFF25ckWvvPKKduzYoQMHDmjWrFnq3r27fvzxRz3zzDOSpK1btyohIUEvvPCCvQ5btmxR\neHi4PfEnSe3bt9fMmTPl4+OjvLw8zZs3T+fOnVNBQYFiYmLUr18/5efn629/+5suXryo2267Tdeu\nXbOXX7dunbZt26bCwkK1aNFCI0eOrLG3VAMAAFgFz/wDAACox06cOKGwsLAblhfPztu6dasyMzM1\nZcoUTZ48WYmJidq3b1+Fn9mrVy9NnjxZrVq10tatW9WvXz/5+/tr7Nix6t27txITE5WXl2eP8dBD\nD5X6jLS0NMN6rVu3Tj4+Ppo8ebImTpyoVatW6fTp09qyZYu8vLwUGxurJ598UmlpaZKklJQUfffd\nd5o8ebJiY2Pl4+OjDRs2VH2DAQAA3GSY+QcAAFCPubu7q6ioyOH6lJQURUZGys3NTW5ubrrzzjuV\nkpJSbmKuWOPGjdWyZUtJUmBgoC5dulRqvbe3t6KiopSQkKDo6GidPHlSkZGRVarX4cOH1atXL0mS\nl5eXWrdurR9//FHHjx/XnXfeKUkKCAhQcHCwpOu3DGdkZGjy5MmSpPz8fHl4MJQFAACoCCMmAACA\neiwkJKTcGXDHjx9XUFBQuWWKE4ElFRQU2P/foMH/a9+OXRrZwjCMP+46StTGRkQkiBgDERHESmwU\nLawtTGFnqX+AINgH1FoRC1NaWdikTyUWmiKFSiwCohCDVXQUZ261ubt4F9nm3mt4fs3AmTN8H6d8\n53zfP627sLBAPp+nvb2d6elpvn37daAkmUxydnb24btKpUIymfxQP47j5vPndz8CxCAImJqaYnV1\n9dPeJEmS9DfHfiVJkr6wTCZDIpHg5OSkuVatVsnlctTrdUZHRymVSsRxzPv7O+VymVQqRSKRoF6v\nE8cxYRhyc3Pzaa22tuMQtMwAAAGgSURBVLZmSDg0NMTr6yuFQoHZ2dkPe2dmZqhWqxSLxeZauVxm\nZ2eHRqNBKpXi8vISgJeXF25vbxkeHmZwcJCrqysAarUad3d3AKTTaS4uLpqjxoVCoblPkiRJv9cW\n//jNKkmSpC+p0WiQz+e5vr6mp6eHjo4OlpeXGRkZIYoijo6OqFQqRFHE5OQkS0tLRFHE9vY2tVqN\n/v5+giAgk8kwPj7O1tYWe3t7ABwfHxNFEdlsloODA0qlEuvr66TTaQqFAufn52xubv5jX09PTxwe\nHnJ/f09XVxfd3d2srKwwMDBAGIbs7+/z+PjI29sb8/PzzM3N8fz8zO7uLmEY0tfXx8PDA9lslrGx\nMU5PTykWiwRBQG9vL2tra3R2dv6bRy1JkvTlGP5JkiTpj8VxTC6XY3FxkYmJif+6HUmSJP2GY7+S\nJEn6I5VKhY2NDZLJpMGfJEnS/5w3/yRJkiRJkqQW5c0/SZIkSZIkqUUZ/kmSJEmSJEktyvBPkiRJ\nkiRJalGGf5IkSZIkSVKLMvyTJEmSJEmSWpThnyRJkiRJktSi/gKCBihOwz1DqQAAAABJRU5ErkJg\ngg==\n",
            "text/plain": [
              "<matplotlib.figure.Figure at 0x7fc811421990>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "VxtwlODcdJZe"
      },
      "cell_type": "markdown",
      "source": [
        "This is worrisome.  If we attempted to fit this model, the `weight_county` vector would likely end up memorizing the results for counties which had only a few training samples. Not good.\n",
        "\n",
        "GLMM's offer a happy middle to the above two GLMs.  I.e., we might consider fitting,\n",
        "\n",
        "    radon ~ 1 + floor + (1 | county)\n",
        "\n",
        "This model is exactly the same as the first, except we've introduced the random-effect, `(1 | county)`.  Adding this random effect has the effect of allowing per-county random fluctuations in radon  In words, `(1 | county)` means that across samples and within a county, the random fluctuation in observed radon will be the same draw from a random Normal . Furthermore, since the covariance is shared among a group (i.e., a \"random effect\"), the counties with more observations provide a hint at the variance of counties with few observations (since the variance is estimated from the whole group)."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "QjvAR2-ZYgxP"
      },
      "cell_type": "markdown",
      "source": [
        "## Experiment"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "bioH0_7ZfC4Z"
      },
      "cell_type": "markdown",
      "source": [
        "We'll now try to fit the `radon ~ 1 + floor + (1 | county)` GLMM using variational inference in TensorFlow. Our only remaining trick is to use stochastic gradient descent. (For brevity, we omit additional details.)"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "17Lf6EEfGMSJ"
      },
      "cell_type": "markdown",
      "source": [
        "### Imports"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "YKuwO2YmH9gG",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import tensorflow_probability as tfp\n",
        "\n",
        "tfd = tfp.distributions\n",
        "tfb = tfp.bijectors"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "J1qaYmpKYdX_"
      },
      "cell_type": "markdown",
      "source": [
        "The following code allows `Session` customization.  Feel free to play with the `session_options` arguments and see how computational performance changes."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "hX3rIfVZCAIk",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def session_options(enable_gpu_ram_resizing=True,\n",
        "                    enable_xla=False):\n",
        "  \"\"\"Convenience function which sets common `tf.Session` options.\"\"\"\n",
        "  config = tf.ConfigProto()\n",
        "  config.log_device_placement = True\n",
        "  if enable_gpu_ram_resizing:\n",
        "    # `allow_growth=True` makes it possible to connect multiple colabs to your\n",
        "    # GPU. Otherwise the colab malloc's all GPU ram.\n",
        "    config.gpu_options.allow_growth = True\n",
        "  if enable_xla:\n",
        "    # Enable on XLA. https://www.tensorflow.org/performance/xla/.\n",
        "    config.graph_options.optimizer_options.global_jit_level = (\n",
        "        tf.OptimizerOptions.ON_1)\n",
        "  return config\n",
        "\n",
        "def reset_sess(config=None):\n",
        "  \"\"\"Convenience function to create the TF graph and session, or reset them.\"\"\"\n",
        "  if config is None:\n",
        "    config = session_options()\n",
        "  tf.reset_default_graph()\n",
        "  global sess\n",
        "  try:\n",
        "    sess.close()\n",
        "  except:\n",
        "    pass\n",
        "  sess = tf.InteractiveSession(config=config)\n",
        "\n",
        "reset_sess()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "VkJaseCAGes0"
      },
      "cell_type": "markdown",
      "source": [
        "### Setup"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "1dwMTHinKHe5"
      },
      "cell_type": "markdown",
      "source": [
        "For transparency, we'll record all knob settings here."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "QFy2NhmkBfGL",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "hparams = tf.contrib.training.HParams(\n",
        "    train_test_split = 0.8,\n",
        "    train_batch_size = 100,\n",
        "    test_batch_size = 10,\n",
        "\n",
        "    dtype=np.float32,\n",
        "    init_raw_scale=np.log(np.expm1(1.)),  # approx= 0.5413\n",
        "    scale_diag_offset=1e-3,\n",
        "\n",
        "    surrogate_posterior_rank=1,\n",
        "\n",
        "    num_monte_carlo_draws=2,\n",
        "\n",
        "    train_iterations = int(3e3),\n",
        "\n",
        "    learning_rate_start = 1e-2,\n",
        "    learning_rate_num_epochs_per_decay = 10,\n",
        "    learning_rate_decay_factor = 0.99,\n",
        ")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "iJz56O8kGPsM"
      },
      "cell_type": "markdown",
      "source": [
        "### Data Munging"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "InypzQj8KViU"
      },
      "cell_type": "markdown",
      "source": [
        "We'll use `tf.data.Dataset` to feed the data into the TensorFlow graph. For a nice tutorial on different data loading patterns, see [\"How to use Dataset in TensorFlow\"](https://towardsdatascience.com/how-to-use-dataset-in-tensorflow-c758ef9e4428)."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "zA_Lt3BO9vZa",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def dataset(df, batch_size):\n",
        "  with tf.name_scope(name='dataset'):\n",
        "    feature_cols=['county', 'log_uranium_ppm']\n",
        "    #by Shihao Shen for testing\n",
        "    #feature_cols=['county', 'log_u2']\n",
        "    label_cols=['log_radon']\n",
        "    dataset = tf.data.Dataset.from_tensor_slices((\n",
        "        df[feature_cols].values.astype(np.int32),\n",
        "        df[label_cols].values.astype(np.float32),\n",
        "    )).repeat().batch(batch_size)\n",
        "    iter_ = tf.data.Iterator.from_structure(\n",
        "        dataset.output_types, dataset.output_shapes)\n",
        "    init_op = iter_.make_initializer(dataset)\n",
        "    features, labels = iter_.get_next()\n",
        "    features.set_shape([batch_size, 2])\n",
        "    labels.set_shape([batch_size, 1])\n",
        "    return init_op, features, labels\n",
        "\n",
        "\n",
        "DataStats = collections.namedtuple(\n",
        "    'DataStats',\n",
        "    [\n",
        "        'num_train',\n",
        "        'num_test',\n",
        "        'num_unique_county',\n",
        "        'num_unique_floor',\n",
        "    ])\n",
        "\n",
        "\n",
        "def split_dataset(df, train_test_split, train_batch_size, test_batch_size):\n",
        "  \"\"\"Creates train/test split data.\"\"\"\n",
        "  with tf.name_scope(name='split_dataset'):\n",
        "    train_df = df.sample(frac=train_test_split, random_state=42)\n",
        "    test_df = df.drop(train_df.index)\n",
        "    data_stats = DataStats(\n",
        "        num_train=len(train_df),\n",
        "        num_test=len(test_df),\n",
        "        num_unique_county=1 + df['county'].max(),\n",
        "        #num_unique_floor=1 + df['floor3'].max(),\n",
        "        num_unique_floor=3,\n",
        "    )\n",
        "    train_init_op, train_features, train_labels = dataset(\n",
        "        train_df,\n",
        "        train_batch_size)\n",
        "    test_init_op, test_features, test_labels = dataset(\n",
        "        test_df,\n",
        "        test_batch_size)\n",
        "    return [\n",
        "        tf.group([train_init_op, test_init_op]),  # dataset_init\n",
        "        train_features,\n",
        "        train_labels,\n",
        "        test_features,\n",
        "        test_labels,\n",
        "        data_stats,\n",
        "    ]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "231v5VR9FgRm",
        "colab_type": "code",
        "outputId": "62f3a203-c1e0-4ff8-fced-fc1fb6754b92",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 85
        }
      },
      "cell_type": "code",
      "source": [
        "#test paddings\n",
        "t = tf.constant([[1, 2, 3], [4, 5, 6]])\n",
        "paddings = tf.constant([[1, 1,], [2, 2]])\n",
        "tf.pad(t, paddings, \"CONSTANT\").eval()\n",
        "#test gather\n"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[0, 0, 0, 0, 0, 0, 0],\n",
              "       [0, 0, 1, 2, 3, 0, 0],\n",
              "       [0, 0, 4, 5, 6, 0, 0],\n",
              "       [0, 0, 0, 0, 0, 0, 0]], dtype=int32)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 17
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "ZWvgUUAbGgkc"
      },
      "cell_type": "markdown",
      "source": [
        "### Specify Model"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "rcl2AuF1GePi",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "def make_prior():\n",
        "  def _fn():\n",
        "    dims = data_stats.num_unique_county\n",
        "    prior_raw_scale = tf.get_variable(\n",
        "        name='prior_raw_scale',\n",
        "        initializer=np.array(hparams.init_raw_scale, hparams.dtype))\n",
        "    scale = tf.nn.softplus(prior_raw_scale) + hparams.scale_diag_offset\n",
        "    return tfd.Independent(\n",
        "        tfd.Normal(loc=np.zeros(dims, hparams.dtype), scale=scale),\n",
        "        reinterpreted_batch_ndims=1,\n",
        "        name='prior')\n",
        "  return tf.make_template('make_prior', _fn)()\n",
        "\n",
        "\n",
        "def make_likelihood(x, u):\n",
        "  def _fn():\n",
        "    intercept = tf.get_variable(\n",
        "      name='intercept',\n",
        "      initializer=np.zeros(1, hparams.dtype))\n",
        "    \n",
        "    weights_floor = tf.pad(\n",
        "        tf.get_variable(\n",
        "            name='weights_floor',\n",
        "            initializer=np.zeros(data_stats.num_unique_floor - 1,\n",
        "                                 hparams.dtype)),\n",
        "        paddings=[[1, 0]])\n",
        "\n",
        "#    weights_floor2 = tf.pad(\n",
        "#        tf.get_variable(\n",
        "#            name='weights_floor2',\n",
        "#            initializer=np.zeros(data_stats.num_unique_floor2 - 1,\n",
        "#                                 hparams.dtype)),\n",
        "#        paddings=[[1, 0]])\n",
        "    \n",
        "   \n",
        "    random_effect = tf.transpose(tf.gather(tf.transpose(u), x[:, 0]))\n",
        "    #fixed_effect = intercept + tf.gather(weights_floor, x[:, 1]) + tf.gather(weights_u,x[:,2])\n",
        "    fixed_effect = intercept + tf.gather(weights_floor, x[:, 1])\n",
        "    #fixed_effect = intercept + tf.gather(weights_floor, x[:, 1]) + tf.gather(weights_floor2, x[:, 5])\n",
        "    predicted_linear_response = fixed_effect + random_effect\n",
        "  \n",
        "    likelihood_raw_scale = tf.get_variable(\n",
        "        name='likelihood_raw_scale',\n",
        "        initializer=np.array(hparams.init_raw_scale, hparams.dtype))\n",
        "    \n",
        "    scale = tf.nn.softplus(likelihood_raw_scale) + hparams.scale_diag_offset\n",
        "    \n",
        "    return tfd.Independent(\n",
        "        tfd.Normal(loc=predicted_linear_response[..., tf.newaxis],\n",
        "                   scale=scale),\n",
        "        reinterpreted_batch_ndims=2,\n",
        "        name='likelihood')\n",
        "\n",
        "  return tf.make_template('make_likelihood', _fn)()\n",
        "\n",
        "\n",
        "def make_surrogate_posterior():\n",
        "  dims = data_stats.num_unique_county\n",
        "  def _tril():\n",
        "    loc = tf.get_variable(\n",
        "        name='surrogate_loc',\n",
        "        shape=[dims],\n",
        "        initializer=tf.zeros_initializer())\n",
        "    raw_scale_tril = tf.get_variable(\n",
        "        name='surrogate_raw_scale_tril',\n",
        "        initializer=np.zeros(dims * (dims + 1) // 2, hparams.dtype))\n",
        "    scale_tril = tfd.fill_triangular(raw_scale_tril)\n",
        "    new_diag = hparams.scale_diag_offset + tf.nn.softplus(\n",
        "        hparams.init_raw_scale + tf.diag_part(scale_tril))\n",
        "    scale_tril = tf.linalg.set_diag(scale_tril, new_diag)\n",
        "    return tfd.MultivariateNormalTriL(\n",
        "        loc=loc,\n",
        "        scale_tril=scale_tril,\n",
        "        name='surrogate_posterior')\n",
        "  def _lowrank():\n",
        "    rank = 1\n",
        "    loc = tf.get_variable(\n",
        "        name='surrogate_loc',\n",
        "        shape=[dims],\n",
        "        initializer=tf.zeros_initializer())\n",
        "    raw_scale = tf.get_variable(\n",
        "        name='surrogate_raw_scale',\n",
        "        initializer=np.zeros(dims * (1 + rank), hparams.dtype))\n",
        "    return tfd.MultivariateNormalDiagPlusLowRank(\n",
        "        loc=loc,\n",
        "        scale_diag=tf.nn.softplus(hparams.init_raw_scale + raw_scale[:dims]),\n",
        "        scale_perturb_factor=tf.reshape(raw_scale[dims:], [dims, rank]),\n",
        "        name='surrogate_posterior')\n",
        "  return tf.make_template(\n",
        "      'make_surrogate_posterior',\n",
        "      _tril if hparams.surrogate_posterior_rank >= dims else _lowrank)()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "8cd1whNpMPwL"
      },
      "cell_type": "markdown",
      "source": [
        "### Main"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "85FS-MhiFHjF",
        "outputId": "82861324-9df6-42a1-846d-2efb9b4b74ac",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 258
        }
      },
      "cell_type": "code",
      "source": [
        "reset_sess()\n",
        "\n",
        "[\n",
        "  data_init,\n",
        "  train_features,\n",
        "  train_labels,\n",
        "  test_features,\n",
        "  test_labels,\n",
        "  data_stats,\n",
        "] = split_dataset(\n",
        "    df,\n",
        "    hparams.train_test_split,\n",
        "    hparams.train_batch_size,\n",
        "    hparams.test_batch_size,\n",
        ")\n",
        "\n",
        "prior = make_prior()\n",
        "print(prior)\n",
        "\n",
        "surrogate_posterior = make_surrogate_posterior()\n",
        "print(surrogate_posterior)\n",
        "\n",
        "minibatch_correction_factor = data_stats.num_train / hparams.train_batch_size\n",
        "\n",
        "def unnormalized_approx_posterior_log_prob(u):\n",
        "  likelihood = make_likelihood(train_features, u)\n",
        "  print(likelihood)\n",
        "  return (likelihood.log_prob(train_labels) * minibatch_correction_factor\n",
        "          + prior.log_prob(u))\n",
        "\n",
        "elbo_loss = tfp.vi.monte_carlo_csiszar_f_divergence(\n",
        "    f=tfp.vi.kl_reverse,  # same as: Evidence Lower BOund\n",
        "    p_log_prob=unnormalized_approx_posterior_log_prob,\n",
        "    q=surrogate_posterior,\n",
        "    num_draws=hparams.num_monte_carlo_draws,\n",
        "    name='elbo_loss')\n",
        "print(elbo_loss)\n",
        "\n",
        "global_step = tf.train.get_or_create_global_step()\n",
        "learning_rate = tf.train.exponential_decay(\n",
        "    learning_rate=hparams.learning_rate_start,\n",
        "    global_step=global_step,\n",
        "    decay_steps=minibatch_correction_factor * hparams.learning_rate_num_epochs_per_decay,\n",
        "    decay_rate=hparams.learning_rate_decay_factor,\n",
        "    staircase=True)\n",
        "opt = tf.train.AdamOptimizer(learning_rate=learning_rate)\n",
        "train = opt.minimize(elbo_loss, global_step=global_step)\n",
        "\n",
        "pprint(tf.trainable_variables())\n",
        "\n",
        "var_init = tf.global_variables_initializer()"
      ],
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "tfp.distributions.Independent(\"make_prior/prior/\", batch_shape=(), event_shape=(85,), dtype=float32)\n",
            "tfp.distributions.MultivariateNormalDiagPlusLowRank(\"make_surrogate_posterior/surrogate_posterior/\", batch_shape=(), event_shape=(85,), dtype=float32)\n",
            "tfp.distributions.Independent(\"elbo_loss/expectation/make_likelihood/likelihood/\", batch_shape=(2,), event_shape=(100, 1), dtype=float32)\n",
            "WARNING:tensorflow:From /usr/local/lib/python2.7/dist-packages/tensorflow_probability/python/monte_carlo.py:161: calling reduce_mean (from tensorflow.python.ops.math_ops) with keep_dims is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "keep_dims is deprecated, use keepdims instead\n",
            "Tensor(\"elbo_loss/expectation/Mean:0\", shape=(), dtype=float32)\n",
            "[<tf.Variable 'make_prior/prior_raw_scale:0' shape=() dtype=float32_ref>,\n",
            " <tf.Variable 'make_surrogate_posterior/surrogate_loc:0' shape=(85,) dtype=float32_ref>,\n",
            " <tf.Variable 'make_surrogate_posterior/surrogate_raw_scale:0' shape=(170,) dtype=float32_ref>,\n",
            " <tf.Variable 'make_likelihood/intercept:0' shape=(1,) dtype=float32_ref>,\n",
            " <tf.Variable 'make_likelihood/weights_floor:0' shape=(2,) dtype=float32_ref>,\n",
            " <tf.Variable 'make_likelihood/likelihood_raw_scale:0' shape=() dtype=float32_ref>]\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "e9IaaqZPK9V8",
        "outputId": "c30c5c8d-c148-4e12-a6d1-b7cc8f24047d",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 2071
        }
      },
      "cell_type": "code",
      "source": [
        "loss_ = np.zeros(hparams.train_iterations)\n",
        "\n",
        "sess.run([var_init, data_init])\n",
        "\n",
        "for iter_ in range(hparams.train_iterations):\n",
        "  [\n",
        "      _,\n",
        "      train_features_,\n",
        "      train_labels_,\n",
        "      test_features_,\n",
        "      test_labels_,\n",
        "      loss_[iter_],\n",
        "  ] = sess.run([\n",
        "      train,\n",
        "      train_features,\n",
        "      train_labels,\n",
        "      test_features,\n",
        "      test_labels,\n",
        "      elbo_loss,\n",
        "  ])\n",
        "  if iter_ % 200 == 0 or iter_ == hparams.train_iterations - 1:\n",
        "    print(\"iter:{:>4}  loss:{:.3f}\".format(\n",
        "        iter_, loss_[iter_]))"
      ],
      "execution_count": 20,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "iter:   0  loss:1645.816\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "error",
          "ename": "InvalidArgumentError",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m\u001b[0m",
            "\u001b[0;31mInvalidArgumentError\u001b[0mTraceback (most recent call last)",
            "\u001b[0;32m<ipython-input-20-1f45783de46f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     17\u001b[0m       \u001b[0mtest_features\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m       \u001b[0mtest_labels\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m       \u001b[0melbo_loss\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m   ])\n\u001b[1;32m     21\u001b[0m   \u001b[0;32mif\u001b[0m \u001b[0miter_\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;36m200\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m0\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0miter_\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mhparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrain_iterations\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python2.7/dist-packages/tensorflow/python/client/session.pyc\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m    927\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    928\u001b[0m       result = self._run(None, fetches, feed_dict, options_ptr,\n\u001b[0;32m--> 929\u001b[0;31m                          run_metadata_ptr)\n\u001b[0m\u001b[1;32m    930\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0mrun_metadata\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    931\u001b[0m         \u001b[0mproto_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf_session\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTF_GetBuffer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrun_metadata_ptr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python2.7/dist-packages/tensorflow/python/client/session.pyc\u001b[0m in \u001b[0;36m_run\u001b[0;34m(self, handle, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m   1150\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mfinal_fetches\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mfinal_targets\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mhandle\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mfeed_dict_tensor\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1151\u001b[0m       results = self._do_run(handle, final_targets, final_fetches,\n\u001b[0;32m-> 1152\u001b[0;31m                              feed_dict_tensor, options, run_metadata)\n\u001b[0m\u001b[1;32m   1153\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1154\u001b[0m       \u001b[0mresults\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python2.7/dist-packages/tensorflow/python/client/session.pyc\u001b[0m in \u001b[0;36m_do_run\u001b[0;34m(self, handle, target_list, fetch_list, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m   1326\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mhandle\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1327\u001b[0m       return self._do_call(_run_fn, feeds, fetches, targets, options,\n\u001b[0;32m-> 1328\u001b[0;31m                            run_metadata)\n\u001b[0m\u001b[1;32m   1329\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1330\u001b[0m       \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_do_call\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_prun_fn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhandle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfeeds\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfetches\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python2.7/dist-packages/tensorflow/python/client/session.pyc\u001b[0m in \u001b[0;36m_do_call\u001b[0;34m(self, fn, *args)\u001b[0m\n\u001b[1;32m   1346\u001b[0m           \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1347\u001b[0m       \u001b[0mmessage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0merror_interpolation\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpolate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmessage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_graph\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1348\u001b[0;31m       \u001b[0;32mraise\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnode_def\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mop\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmessage\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1349\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1350\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0m_extend_graph\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mInvalidArgumentError\u001b[0m: indices[80] = -2 is not in [0, 3)\n\t [[node elbo_loss/expectation/make_likelihood/GatherV2_1 (defined at <ipython-input-18-3d0114a40e06>:38)  = GatherV2[Taxis=DT_INT32, Tindices=DT_INT32, Tparams=DT_FLOAT, _device=\"/job:localhost/replica:0/task:0/device:CPU:0\"](elbo_loss/expectation/make_likelihood/Pad, elbo_loss/expectation/make_likelihood/strided_slice_1, elbo_loss/expectation/make_likelihood/GatherV2_1/axis)]]\n\nCaused by op u'elbo_loss/expectation/make_likelihood/GatherV2_1', defined at:\n  File \"/usr/lib/python2.7/runpy.py\", line 174, in _run_module_as_main\n    \"__main__\", fname, loader, pkg_name)\n  File \"/usr/lib/python2.7/runpy.py\", line 72, in _run_code\n    exec code in run_globals\n  File \"/usr/local/lib/python2.7/dist-packages/ipykernel_launcher.py\", line 16, in <module>\n    app.launch_new_instance()\n  File \"/usr/local/lib/python2.7/dist-packages/traitlets/config/application.py\", line 658, in launch_instance\n    app.start()\n  File \"/usr/local/lib/python2.7/dist-packages/ipykernel/kernelapp.py\", line 477, in start\n    ioloop.IOLoop.instance().start()\n  File \"/usr/local/lib/python2.7/dist-packages/tornado/ioloop.py\", line 832, in start\n    self._run_callback(self._callbacks.popleft())\n  File \"/usr/local/lib/python2.7/dist-packages/tornado/ioloop.py\", line 605, in _run_callback\n    ret = callback()\n  File \"/usr/local/lib/python2.7/dist-packages/tornado/stack_context.py\", line 277, in null_wrapper\n    return fn(*args, **kwargs)\n  File \"/usr/local/lib/python2.7/dist-packages/zmq/eventloop/zmqstream.py\", line 536, in <lambda>\n    self.io_loop.add_callback(lambda : self._handle_events(self.socket, 0))\n  File \"/usr/local/lib/python2.7/dist-packages/zmq/eventloop/zmqstream.py\", line 450, in _handle_events\n    self._handle_recv()\n  File \"/usr/local/lib/python2.7/dist-packages/zmq/eventloop/zmqstream.py\", line 480, in _handle_recv\n    self._run_callback(callback, msg)\n  File \"/usr/local/lib/python2.7/dist-packages/zmq/eventloop/zmqstream.py\", line 432, in _run_callback\n    callback(*args, **kwargs)\n  File \"/usr/local/lib/python2.7/dist-packages/tornado/stack_context.py\", line 277, in null_wrapper\n    return fn(*args, **kwargs)\n  File \"/usr/local/lib/python2.7/dist-packages/ipykernel/kernelbase.py\", line 283, in dispatcher\n    return self.dispatch_shell(stream, msg)\n  File \"/usr/local/lib/python2.7/dist-packages/ipykernel/kernelbase.py\", line 235, in dispatch_shell\n    handler(stream, idents, msg)\n  File \"/usr/local/lib/python2.7/dist-packages/ipykernel/kernelbase.py\", line 399, in execute_request\n    user_expressions, allow_stdin)\n  File \"/usr/local/lib/python2.7/dist-packages/ipykernel/ipkernel.py\", line 196, in do_execute\n    res = shell.run_cell(code, store_history=store_history, silent=silent)\n  File \"/usr/local/lib/python2.7/dist-packages/ipykernel/zmqshell.py\", line 533, in run_cell\n    return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n  File \"/usr/local/lib/python2.7/dist-packages/IPython/core/interactiveshell.py\", line 2718, in run_cell\n    interactivity=interactivity, compiler=compiler, result=result)\n  File \"/usr/local/lib/python2.7/dist-packages/IPython/core/interactiveshell.py\", line 2822, in run_ast_nodes\n    if self.run_code(code, result):\n  File \"/usr/local/lib/python2.7/dist-packages/IPython/core/interactiveshell.py\", line 2882, in run_code\n    exec(code_obj, self.user_global_ns, self.user_ns)\n  File \"<ipython-input-19-40b328b5c697>\", line 36, in <module>\n    name='elbo_loss')\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow_probability/python/vi/csiszar_divergence.py\", line 903, in monte_carlo_csiszar_f_divergence\n    use_reparametrization=use_reparametrization)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow_probability/python/monte_carlo.py\", line 161, in expectation\n    return tf.reduce_mean(f(samples), axis=axis, keep_dims=keep_dims)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow_probability/python/vi/csiszar_divergence.py\", line 900, in <lambda>\n    f=lambda q_samples: f(p_log_prob(q_samples) - q.log_prob(q_samples)),\n  File \"<ipython-input-19-40b328b5c697>\", line 26, in unnormalized_approx_posterior_log_prob\n    likelihood = make_likelihood(train_features, u)\n  File \"<ipython-input-18-3d0114a40e06>\", line 54, in make_likelihood\n    return tf.make_template('make_likelihood', _fn)()\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/ops/template.py\", line 368, in __call__\n    return self._call_func(args, kwargs)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/ops/template.py\", line 311, in _call_func\n    result = self._func(*args, **kwargs)\n  File \"<ipython-input-18-3d0114a40e06>\", line 38, in _fn\n    fixed_effect = intercept + tf.gather(weights_floor, x[:, 1])\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/ops/array_ops.py\", line 2675, in gather\n    return gen_array_ops.gather_v2(params, indices, axis, name=name)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/ops/gen_array_ops.py\", line 3332, in gather_v2\n    \"GatherV2\", params=params, indices=indices, axis=axis, name=name)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/framework/op_def_library.py\", line 787, in _apply_op_helper\n    op_def=op_def)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/util/deprecation.py\", line 488, in new_func\n    return func(*args, **kwargs)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/framework/ops.py\", line 3274, in create_op\n    op_def=op_def)\n  File \"/usr/local/lib/python2.7/dist-packages/tensorflow/python/framework/ops.py\", line 1770, in __init__\n    self._traceback = tf_stack.extract_stack()\n\nInvalidArgumentError (see above for traceback): indices[80] = -2 is not in [0, 3)\n\t [[node elbo_loss/expectation/make_likelihood/GatherV2_1 (defined at <ipython-input-18-3d0114a40e06>:38)  = GatherV2[Taxis=DT_INT32, Tindices=DT_INT32, Tparams=DT_FLOAT, _device=\"/job:localhost/replica:0/task:0/device:CPU:0\"](elbo_loss/expectation/make_likelihood/Pad, elbo_loss/expectation/make_likelihood/strided_slice_1, elbo_loss/expectation/make_likelihood/GatherV2_1/axis)]]\n"
          ]
        }
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "Sp-9-fdSK36r",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "train"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "6H7RVU-SXexC"
      },
      "cell_type": "markdown",
      "source": [
        "### Results"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "OmiHUfBDUZg4",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "surrogate_posterior_cov = surrogate_posterior.covariance()\n",
        "surrogate_posterior_mean = surrogate_posterior.mean()\n",
        "surrogate_posterior_scale = tf.cholesky(surrogate_posterior_cov )\n",
        "\n",
        "with tf.variable_scope('make_likelihood', reuse=True):\n",
        "  intercept = tf.get_variable(name='intercept')\n",
        "  weights_floor = tf.get_variable(name='weights_floor')\n",
        "\n",
        "[\n",
        "    surrogate_posterior_mean_,\n",
        "    surrogate_posterior_scale_,\n",
        "    surrogate_posterior_cov_,\n",
        "    intercept_,\n",
        "    weights_floor_,\n",
        "    prior_scale_,\n",
        "] = sess.run([\n",
        "    surrogate_posterior_mean,\n",
        "    surrogate_posterior_scale,\n",
        "    surrogate_posterior_cov,\n",
        "    intercept,\n",
        "    weights_floor,\n",
        "    prior.distribution.scale,\n",
        "])\n",
        "\n",
        "\n",
        "print('    intercept: ', intercept_)\n",
        "print('weights_floor: ', weights_floor_)\n",
        "print('  prior_scale: ', prior_scale_)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "aK9CaGtSflo4",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "print(len(surrogate_posterior_mean_))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "eG8Z-HVwiz63"
      },
      "cell_type": "markdown",
      "source": [
        "We will now plot the training-set loss as a function of training iteration. Note that although the `learning_rate` decays to `0`, there's still noise in the loss owing to the learned parameters being evaluated over different mini-batches."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "qa63HEgPu2Zn",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "plt.plot(loss_, 'b-');\n",
        "plt.xlabel('iter');\n",
        "plt.ylabel('loss');"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "WJHWXM33i_3o"
      },
      "cell_type": "markdown",
      "source": [
        "We now plot a nonparametric estimate of the density of the diagonal of the surrogate posterior covariance. Roughly speaking, this shows us how much intrinsic variability there is in each county's `log(radon)` readings."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "Pa47BzL0u2KA",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "surrogate_posterior_diag_scale_ = np.diag(surrogate_posterior_scale_)\n",
        "with warnings.catch_warnings():\n",
        "  warnings.simplefilter(\"ignore\")\n",
        "  sns.kdeplot(surrogate_posterior_diag_scale_, shade=True, color=\"r\");\n",
        "  plt.xlabel('posterior std. deviation');\n",
        "  plt.ylabel('density');"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "L9WYd1WxjQ2u"
      },
      "cell_type": "markdown",
      "source": [
        "We now show a heat map of the covariance matrix. This shows us that just looking at the diagonal was a perfectly reasonable thing to do; the matrix is clearly [diagonally dominant](https://en.wikipedia.org/wiki/Diagonally_dominant_matrix)."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "-etGwnAHu2Br",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "sns.heatmap(surrogate_posterior_cov_)\n",
        "plt.title('Surrogate Posterior Covariance')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "unU7DMKVjel7"
      },
      "cell_type": "markdown",
      "source": [
        "We now conjecture that the county-wise standard deviations might be log-correlated with the number of observations for that county."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "PqpaImuYu14c",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "cnt = np.histogram(\n",
        "    df['county'],\n",
        "    np.arange(data_stats.num_unique_county+1))[0]\n",
        "\n",
        "y = np.stack([cnt, surrogate_posterior_diag_scale_]).T\n",
        "y = y[y[:,0].argsort()]  # sort by zero-th col\n",
        "sns.regplot(x=np.log(y[:,0]), y=y[:, 1], color='g')\n",
        "plt.xlabel('county log-count');\n",
        "plt.ylabel('posterior std. deviation');"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "JYJL76jkcFjl"
      },
      "cell_type": "markdown",
      "source": [
        "Indeed--we do appear to have learned a log-linear relationship between std. deviation and county frequency. This is neat because we never expressed this idea in the model; it is simply apparently true."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "AxhvAIeFs2cL"
      },
      "cell_type": "markdown",
      "source": [
        "\n",
        "## Comparing to `lme4` in R"
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "hvRDd7T-s7qY",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "%%shell\n",
        "exit  # Trick to make this block not execute.\n",
        "\n",
        "radon = read.csv('srrs2.dat', header = TRUE)\n",
        "radon = radon[radon$state=='MN',]\n",
        "radon$radon = ifelse(radon$activity==0., 0.1, radon$activity)\n",
        "radon$log_radon = log(radon$radon)\n",
        "\n",
        "# install.packages('lme4')\n",
        "library(lme4)\n",
        "fit <- lmer(log_radon ~ 1 + floor + (1 | county), data=radon)\n",
        "fit\n",
        "\n",
        "# Linear mixed model fit by REML ['lmerMod']\n",
        "# Formula: log_radon ~ 1 + floor + (1 | county)\n",
        "#    Data: radon\n",
        "# REML criterion at convergence: 2171.305\n",
        "# Random effects:\n",
        "#  Groups   Name        Std.Dev.\n",
        "#  county   (Intercept) 0.3282\n",
        "#  Residual             0.7556\n",
        "# Number of obs: 919, groups:  county, 85\n",
        "# Fixed Effects:\n",
        "# (Intercept)        floor\n",
        "#       1.462       -0.693"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "n5IqiHERv91u"
      },
      "cell_type": "markdown",
      "source": [
        "The following table summarizes the results."
      ]
    },
    {
      "metadata": {
        "colab_type": "code",
        "id": "D0sUh3NNuqlw",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "print(pd.DataFrame(data=dict(intercept=[1.462, intercept_[0]],\n",
        "                             floor=[-0.693, weights_floor_[0]],\n",
        "                             scale=[0.3282, prior_scale_]),\n",
        "                   index=['lme4', 'vi']))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "nVjHJxVdwBXb"
      },
      "cell_type": "markdown",
      "source": [
        "This table indicates the VI results are within ~10% of `lme4`'s.  This is somewhat surprising since:\n",
        "- `lme4` is based on [Laplace's method](https://www.jstatsoft.org/article/view/v067i01/) (not VI),\n",
        "- we used mini-batch SGD,\n",
        "- no effort was made in this colab to actually converge,\n",
        "- minimal effort was made to tune hyperparameters,\n",
        "- no effort was taken regularize or preprocess the data (eg, center features, etc.)."
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "ApP0PtwYN_ah"
      },
      "cell_type": "markdown",
      "source": [
        "## Conclusion"
      ]
    },
    {
      "metadata": {
        "colab_type": "text",
        "id": "eIFHW00tOJwo"
      },
      "cell_type": "markdown",
      "source": [
        "In this colab we described Generalized Linear Mixed-effects Models and showed how to use variational inference to fit them using TensorFlow. Although the toy problem only had a few 100 training samples, the techniques used here are identical to what's needed at scale."
      ]
    }
  ]
}