{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# MLE and MAP Estimation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this short tutorial we review how to do Maximum Likelihood (MLE) and Maximum a Posteriori (MAP) estimation in Pyro." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "import torch\n", "from torch.distributions import constraints\n", "import pyro\n", "import pyro.distributions as dist\n", "from pyro.infer import SVI, Trace_ELBO\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider the simple \"fair coin\" example covered in a [previous tutorial](http://pyro.ai/examples/svi_part_i.html#A-simple-example)." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "data = torch.zeros(10)\n", "data[0:6] = 1.0\n", "\n", "def original_model(data):\n", " f = pyro.sample(\"latent_fairness\", dist.Beta(10.0, 10.0))\n", " with pyro.plate(\"data\", data.size(0)):\n", " pyro.sample(\"obs\", dist.Bernoulli(f), obs=data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To facilitate comparison between different inference techniques, we construct a training helper:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "def train(model, guide, lr=0.005, n_steps=201):\n", " pyro.clear_param_store()\n", " adam_params = {\"lr\": lr}\n", " adam = pyro.optim.Adam(adam_params)\n", " svi = SVI(model, guide, adam, loss=Trace_ELBO())\n", "\n", " for step in range(n_steps):\n", " loss = svi.step(data)\n", " if step % 50 == 0:\n", " print('[iter {}] loss: {:.4f}'.format(step, loss))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## MLE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our model has a single latent variable `latent_fairness`. To do Maximum Likelihood Estimation we simply \"demote\" our latent variable `latent_fairness` to a Pyro parameter." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "def model_mle(data):\n", " # note that we need to include the interval constraint; \n", " # in original_model() this constraint appears implicitly in \n", " # the support of the Beta distribution.\n", " f = pyro.param(\"latent_fairness\", torch.tensor(0.5), \n", " constraint=constraints.unit_interval)\n", " with pyro.plate(\"data\", data.size(0)):\n", " pyro.sample(\"obs\", dist.Bernoulli(f), obs=data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can render our model as shown below." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "%3\n", "\n", "\n", "cluster_data\n", "\n", "data\n", "\n", "\n", "\n", "latent_fairness\n", "\n", "latent_fairness\n", "\n", "\n", "\n", "obs\n", "\n", "obs\n", "\n", "\n", "\n", "latent_fairness->obs\n", "\n", "\n", "\n", "\n", "\n", "distribution_description_node\n", "obs ~ Bernoulli\n", "latent_fairness ∈ Interval(lower_bound=0.0, upper_bound=1.0)\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pyro.render_model(model_mle, model_args=(data,), render_distributions=True, render_params=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we no longer have any latent variables, our guide can be empty:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [], "source": [ "def guide_mle(data):\n", " pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see what result we get." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[iter 0] loss: 6.9315\n", "[iter 50] loss: 6.7693\n", "[iter 100] loss: 6.7333\n", "[iter 150] loss: 6.7302\n", "[iter 200] loss: 6.7301\n" ] } ], "source": [ "train(model_mle, guide_mle)" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our MLE estimate of the latent fairness is 0.600\n" ] } ], "source": [ "mle_estimate = pyro.param(\"latent_fairness\").item()\n", "print(\"Our MLE estimate of the latent fairness is {:.3f}\".format(mle_estimate))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also compare our MLE estimate with the analytical MLE estimate which is given as: $\\frac{\\#Heads}{\\#Heads + \\#Tails}$. As we encode `Heads` as 1 and `Tails` as 0, we can directly find the analytical MLE as `data.sum()/data.size(0)` or `data.mean()`." ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The analytical MLE estimate of the latent fairness is 0.600\n" ] } ], "source": [ "print(\"The analytical MLE estimate of the latent fairness is {:.3f}\".format(\n", " data.mean()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Thus with MLE we get a point estimate of `latent_fairness` which matches the analytical MLE estimate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You may be wondering how to interpret the loss numbers in our experiment above. The loss is equivalent to the negative log likelihood (NLL) of observing the data under the Bernoulli likelihood. Thus, the above procedure was equivalent to minimizing the NLL. We confirm the same below." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The negative log likelihood given latent fairness = 0.600 is 6.7301 which matches the loss obtained via our training procedure.\n" ] } ], "source": [ "nll = -dist.Bernoulli(mle_estimate).log_prob(data).sum()\n", "print(f\"The negative log likelihood given latent fairness = {mle_estimate:0.3f} is {nll:0.4f} which matches the loss obtained via our training procedure.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## MAP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With Maximum a Posteriori estimation, we also get a point estimate of our latent variables. The difference to MLE is that these estimates will be regularized by the prior. We can understand the difference between the model we use for MLE and MAP via the rendering below, where we can see `latent_fairness` is a `pyro.sample` in original model." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "%3\n", "\n", "\n", "cluster_data\n", "\n", "data\n", "\n", "\n", "\n", "latent_fairness\n", "\n", "latent_fairness\n", "\n", "\n", "\n", "obs\n", "\n", "obs\n", "\n", "\n", "\n", "latent_fairness->obs\n", "\n", "\n", "\n", "\n", "\n", "distribution_description_node\n", "latent_fairness ~ Beta\n", "obs ~ Bernoulli\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pyro.render_model(original_model, model_args=(data,), render_distributions=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To do MAP in Pyro we use a [Delta distribution](http://docs.pyro.ai/en/stable/distributions.html#pyro.distributions.Delta) for the guide. Recall that the `Delta` distribution puts all its probability mass at a single value. The `Delta` distribution will be parameterized by a learnable parameter. " ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [], "source": [ "def guide_map(data):\n", " f_map = pyro.param(\"f_map\", torch.tensor(0.5),\n", " constraint=constraints.unit_interval)\n", " pyro.sample(\"latent_fairness\", dist.Delta(f_map))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see how this result differs from MLE." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[iter 0] loss: 5.6719\n", "[iter 50] loss: 5.6007\n", "[iter 100] loss: 5.6004\n", "[iter 150] loss: 5.6004\n", "[iter 200] loss: 5.6004\n" ] } ], "source": [ "train(original_model, guide_map)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our MAP estimate of the latent fairness is 0.536\n" ] } ], "source": [ "map_estimate = pyro.param(\"f_map\").item()\n", "print(\"Our MAP estimate of the latent fairness is {:.3f}\".format(map_estimate))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To understand what's going on note that the prior mean of the `latent_fairness` in our model is 0.5, since that is the mean of `Beta(10.0, 10.0)`. The MLE estimate (which ignores the prior) gives us a result that is entirely determined by the raw counts (6 heads and 4 tails). In contrast the MAP estimate is regularized towards the prior mean, which is why the MAP estimate is somewhere between 0.5 and 0.6. We can also understand these from the plot below. Infact, we can also analytically calculate the MAP estimate given the `Beta` prior and `Bernoulli` likelihood.\n", "\n", "Our `Beta` prior is parameterised by $\\alpha_{Heads}$ (= 10 in our example) and $\\alpha_{Tails}$ (= 10 in our example). The closed form expression for MAP estimate is:\n", "$\\frac{\\alpha_{Heads} + ~\\#Heads}{\\alpha_{Heads} + ~\\#Heads +~ \\alpha_{Tails} + ~\\#Tails}$ = $\\frac{10 + 6}{10 + 6 + 10 + 4}$ = $\\frac{16}{30} = 0.5\\bar{3}$" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgcAAAEGCAYAAADxFTYDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABIpUlEQVR4nO3dd3yUVfY/8M+ZyaSR3kNICIFUAgQJQXqAqIAiLiDCCoplsStr2a+iK4u6rvsTXCsIuEjRFRBBAREEIRRBIECAkEIJJQnppJI6mfv7YyaYxCQzIc/MMzM579drXply597zUDJn7nOfe0gIAcYYY4yxRgq5A2CMMcaYeeHkgDHGGGPNcHLAGGOMsWY4OWCMMcZYM5wcMMYYY6wZG7kG9vLyEsHBwXINz5jZycjIAACEh4fLHAkzZ8ePHy8SQnhL0I+PjY3NFwCiwV8UuxoNgBS1Wv34oEGDClprIFtyEBwcjKSkJLmGZ6xNWWVZAIBA18D2Gx46pP05bJgk48bHxwMAEhMTJemPWSciuiJFPzY2Nl/4+flFent7lygUCr6mvQvRaDRUWFgYlZeX9wWAe1trI1tywJi5mr15NgAgcU5i+w3nz9f+5A9zZpmiOTHomhQKhfD29i7Ly8uLbqsNJweMtfDGqDfkGfcNecZlXZaCE4OuS/d33+bpJE4OGGshISRBnnET5BmXMcZa4kUojLWQWZKJzJJMk4+bnJyM5ORkk4/LGGMtcXLAWAuP/vAoHv3hUZOPO2/ePMybN8/k4zImF6VSOSgiIiIqNDS074QJE0IqKipa/UwaOHBghNRjx8XFhQcHB0dHREREhYSE9F20aJGXvve89dZbPm3F2NL48eNDUlNTbQHgueeeC/Dz8+vv6Og4sGmb6upquvvuu0OCgoKi+/fvH5GRkWGrr9/7778/2MPDY0BoaGjfps/n5+crhw0bFtqzZ8/oYcOGhRYWFioB4JtvvnGdN29ed0NiboqTA8ZaWBi/EAvjF8odBmNWz87OTpOenp56/vz5syqVSixevLjZJZr19fUAgJMnT6Yb2mfjewyxZs2azPT09NTDhw+nv/XWWz1qamqovfbLli3zrays1Pu5mZSUZN/Q0EBRUVF1AHDfffeVHjlyJK1lu48++sjL1dVVffXq1ZRnn302/8UXX+yhr+9HH320aMuWLedbPr9gwQL/+Pj4iitXrqTEx8dXvPnmm34A8MADD5Tt3LnTzdCkppHeNQdEZA9gPwA7XfuNQogFLdrMAfA+gBzdU58KIb7oSCCMmYvRwaPlDoExk3pl46nAc3kVjlL2GebnXPX+tAFZhrYfMWJE5enTpx22bdvmvGDBgu6urq4NmZmZ9pcvX05xdHQcWFVVdVKj0eCpp57qsWfPHlciEq+88kruX/7yl5LW3tORWMvLy5UODg4aGxsbAQCbNm1yeeutt7rX1dVRz549a9etW3f5k08+8SooKFCNHj06zN3dXX3kyJFzDz74YNCpU6e61dTUKCZNmlTyn//85xoArFq1ynPSpEmljf2PGzfuRmvjbtu2ze0f//jHNQB45JFHSv7v//4vSKPRQKFo+3N8woQJla3NMOzYscNt3759GQDwxBNPFI8ePTocQI5CocCwYcMq1q9f7/r444+XGPpnYsiCxFoAY4UQlUSkAnCQiH4SQvzWot16IcSzhg7MmDlp0AjsSs1H6rUyFFRr1xsM7N4Xk2MCYK9SyhwdY9atvr4eO3fudLnzzjvLASA1NdXx5MmTZyMiIuqatluzZo3bmTNnHNLS0s7m5ubaxMXFRd55552Vrb1n4cKFPt99952HSqUSs2fPLho3blzlt99+6zZq1KjKhISEGwDw0EMPhdja2mquXr1q//bbb1+1sbFBbm6uzbvvvuu/f//+cy4uLprXX3/d7+233/ZdtGhR7tKlS3337dt3zt/fXw0AH3zwQY6vr2+DWq3GsGHDwo8cOeIwZMiQ6iNHjjg99NBD1/Udd35+vm2vXr3qAEClUsHJyakhPz/fprH/jiguLrbp2bNnPQAEBgbWFxcX3/x8j42NvXHgwAEnSZMDIYQAUKl7qNLd+PIXZhXq1Bp8fzIHn++/iMxCbXKfZ/cqIAC/w+/h/+3IwKMjemHW7T3h6qCSOVrGjKMj3/ClVFtbq4iIiIgCgCFDhlS88MILRbt373bq37//jZaJAQAcOHDAefr06ddtbGwQGBioHjJkSOXBgwcdXV1dNS3fk5+frzp+/Hh6amqq3YIFC7p/9NFHfpMnTy4ZM2bMzW/xa9asyRw1alTVtWvXbIYOHRoxefLk8uPHjztcvHjRPi4uLgIA6uvradCgQZUtYwGA1atXe6xatcpLrVZTYWGh6tSpU/ZDhgypLiwsVPn5+Rl+fkNiCoUCRL+fIfHz81Pn5eXpXc/QlEGXMhKREsBxAH0AfCaEONJKs6lENArAOQB/FUL84R8bEc0FMBcAgoKCOhInY5IrqKjBQ/89ivS8CvTt7oJP/zwQE6L9cSTHHUII2KgjsCTxIt7fmYGVBy9h9aNxiA5w/b2Dd9+VNJ53Je6PMXPXuOag5fOOjo6ajvbV8j1LlizJAYABAwbUfv/995fae2/37t3V0dHRVfv37+/m6OioGTFiRPnWrVvbfU96errtp59+6nv8+PE0b2/vhqlTpwbX1NQoGo+rurpa7zl+X1/fukuXLtn27t27vr6+HpWVlUpfX98OzxoAgKenp/rKlSuqnj171l+5ckXl4eFxs5/q6mqyt7fv0J+pQQsUhBANQogYAD0AxBFRy12VtgIIFkL0B7ALwOo2+lkuhIgVQsR6e3d6a3DGblnW9Src//lhXL1ehc9n3YZtz43APf27Q6kgDAschuFBwzEkxBOrH43D1mdHwF6lxMzlv+FIZvHvnQwbJtnWydruhmGYhP0xZm1GjRpVsXHjRg+1Wo1r167ZHD161GnkyJGtns/viIqKCsXZs2cdw8PDa+Pj428kJSU5paSk2AFAeXm54vTp03YA0K1bt4aysjIFAJSUlCgdHBw0Hh4eDVlZWTaJiYk3vzmEhobWpKWl2ekb9+677y5duXKlJwB8+eWX7kOHDq1QKBS4dOmSaujQoWEdOYa77rqrdNmyZZ4AsGzZMs/x48eXNr6WkZFh37dv3+qO9Neh1YtCiFIAewGMb/F8sRCiVvfwCwCDOtIvY6Z0Pr8C0z4/hNKqenz1+BCMj/ZvNgWXUpCClILf1zP16+GKb58cCh8XOzy08ij2puvqlBw69Ht9BQkcOnQIhyTsjzFrM3v27NK+fftWR0ZG9o2Pjw9buHBhdlBQ0C190wa0aw4iIiKiBgwYEDljxoyikSNHVnXv3l29bNmyyzNmzAgJCwuLio2NjThz5ow9ADz88MNF48ePDxsyZEjY0KFDq6Ojo6t69+4dPX369JCmpx4mTJhQumfPHufGx08++WQPX1/f/jU1NQpfX9/+L774YncAeOGFF4pKSkpsgoKCoj/55BO/RYsWZQNAVlaWSqlUtnr6ftKkSb1GjBgRcenSJTtfX9/+//nPf7wAYOHChbl79+516dmzZ3RiYqLLwoULcxvfs3//fuf77ruvrCN/NqRdUtBOAyJvAPVCiFIicgDwM4B/CyG2NWnjL4TI1d3/E4D/E0Lc3l6/sbGxggsvMVO7fqMOEz7aD40A1j4Whwg/lz+0iV8VD+CPtRWKK2vx8JdHkZFXgU1PDUe/WZO1L0hUW4ELLzFDENFxIURsZ/s5derU5QEDBhRJERNrrrKykoYPHx5+/PjxdBubjm9E/O6773r37Nmz7sEHH+zQB3prsrKybKZPnx5y+PDhcy1fO3XqlNeAAQOCW3ufIVH7A1itW3egALBBCLGNiN4CkCSE2ALgeSK6F4AawHUAc27xOBgzGiEEXvn2FEpu1GPzM8NaTQwA4P073m/1eU8nO6x9dAgmfnwAz31zAj9+8hm62fIO5Iyx5pycnMSbb7557dKlS7ahoaF/WFipz/z58wuliiUzM9N28eLFHV5wasjVCqcBDGzl+Teb3H8NwGsdHZwxU1p96DJ+SS/APyZFoW931zbbDQ4Y3OZr7t1s8Z8HYvDnFb/hzTQ1Fk/v22ZbxljXNXXq1HK5YwCA0aNHV93K+3iHRNYlpF4rx7vb0zEuwgcPDwtut21yXjKS85LbfP32EE88OzYUpRu+w28frpI0TsYYMwc8J8qsXn2DBi+sOwk3RxXev39As8WHrZm3Yx6AP645aOr5sX2Q8cw2VCWpkffIDPi52ksYMWOMyYuTA2b1/nfkKs4XVGLFQ7Hw6KZ/H5APx3+ot42NUoHe3t1wKrsMi3/OwPv3D+h0nB9+qH9cxhgzBT6twKxaeU09Ptx9DkNDPJEQ6WPQe2L8YhDjF6O3nb1KCT8Xe2w8kY2z1zq9qBgxMTGIidE/LmPWgogGTZ48uVfj4/r6eri7uw8YM2ZMHwD4+OOPPR966KE/7JgXEBDQLywsLCoiIiIqIiIias6cOYGmjLsr4JkDZtU+23sBpdX1eP3uSL2nExodyzkGoP2FiY0C3B3g6qDCu9vT8NVjQwweozW7d+8GACQkJNxyH4xZEgcHB01GRoZDZWUlOTk5ic2bN7v4+voatO1w0xoHTHo8c8CsVtb1Knz562X8aWBA822P9Xhl1yt4ZdcrBrW1URCeHxuKXy8UIzGjc1cfvfPOO3jnnXc61QdjliYhIaHs22+/dQOAb775xmPq1Kl6CxYx4+OZA2a13t+ZAQUBr9wV3qH3fTrx0w61n3V7T6w5fBn/3J6GkaFesFFyzs0sT1xc3B/+o0yZMuX6q6++WlhRUaEYN25caMvXZ82aVfT8888X5+bm2kyePLl309eOHj2aYci4s2fPvr5gwQL/Bx54oDQtLc3xscceKz506JCTvveNHj06rLG08cyZM4sWLFhQYMh4zDCcHDCrdKGgAltOXcMzY3rD39WhQ++N9mlZOqR9tjYKvDohAk9+dQI/nsnF5JiADr2fsa5syJAh1dnZ2XYrVqzwSEhIMHjxDp9WMC5ODphVWr4/E/YqBR4bEdLh9x7K0tY3GBZoeBGkO6P80Nu7Gz7fl4l7B3Tv1NoDxuTQ3jd9Z2dnTXuv+/v7qw2dKWjN+PHjSxcsWBD4888/ZxQUFPDnkhng+U9mdfLLa7D5ZA6mxwYadOliS/N/mY/5v8zv0HsUCsITo3ojLbccBy/wdvWMdcRTTz1V9PLLL1+Li4vrUOVAZjycoTGr8+Wvl9GgEXj8FmYNAGDZPcsMbNi83eSB3bHo5wws25eJkaEdL0m+bJmB4zJmZXr37l3/xhtvtLpmYOPGjZ47d+50a3x86NChNKD5moPIyMiqzZs3XzZBqF2G3qqMxsJVGZkxVNTUY9i/9mB0uDc+/fNtJh//830X8d5P6dj23IgOXSHBmKG4KiOTSntVGfm0ArMq3xy9iopaNZ4Y1Vt/4zbsu7wP+y7v099w61btrYk/DwmCk50Nlu/P7PC4W7duxdYW/THGmBz4tAKzGvUNGqw8eBnDenuiX49b/9a+IHEBgPZrKwAAFi/W/pw06eZTLvYq/HlIEP578BJeuSscgR6OBo+7WNffpCb9McaYHHjmgFmNX9LykVdeg8dG9NLfuB0rJ6/Eyskr9TfcuFF7a2HOsGAIIbD+WIdLqDPGmFng5IBZjf8dzYK/qz3iww2rodCWEPcQhLgbsJjRy0t7a6G7mwPGhPtgQ1IW6hs0nYqFMcbkwMkBswpZ16tw4HwhpscGQqno3B4DuzN3Y3fmbv0NV63S3loxMy4IBRW12JPOm7YxxiwPrzlgVmFDUhYIwPTBnS/O9s5+bX2DhBA9BZAaE4M5c/7wUny4N3xd7LDu6FXc1dev0zExxpgp8cwBs3jqBg3WH8vC6DBvBLh1bKvk1qz901qs/dPaTvVho1TggdhAJJ4rRE6pYfu6rF27FmvXdm5cxiyJvpLNjRISEnoPGDAgoulzL774YncfH5/+ERERUaGhoX2//vprvnZYQnqTAyKyJ6KjRHSKiM4S0cJW2tgR0XoiukBER4go2CjRMtaKPekFKKioxcy4P5R9vyWBroEIdO38DETjLIahCxMDAwMRGMhl6VnX0bRkMwC0VrK5qKhImZKS0q2iokKZmprabMvTJ598Mj89PT11/fr1F5999tnghoYGU4Zv1QyZOagFMFYIMQBADIDxRHR7izaPASgRQvQB8B8A/5Y0Ssbase5YFnyc7TA2onMLERvtuLADOy7s6HQ/PdwdMSrUG98mZUFtwMLE9evXY/369Z0elzFLoq9k81dffeWWkJBQ+qc//en6mjVrPFrr47bbbqtRKpXIy8vjU+US0ZscCK1K3UOV7tZyW8XJAFbr7m8EMI648gwzgdyyaiRmFOD+2B6SlUp+7+B7eO/ge5L0NTMuCLllNdh/vlBv26VLl2Lp0qWSjMtYh8XFheu9vfmmb7P2H3/sCQDIzbX5Q1sDzZ49+/r69evdq6qqKC0tzXHo0KE3mr6+YcMGj1mzZl1/+OGHr2/atKnV5GDPnj3dFAqF4CqN0jEoyyIiJYDjAPoA+EwIcaRFkwAAWQAghFATURkATwBFLfqZC2AuAAQFSTMFzLq2LcnXoBHA/YOkm45fN22dZH2Ni/SBu6MKm09ew9gIX/1vYKyLaa9kc1ZWls2VK1fs77zzzkqFQgEbGxtx7Ngx+8GDB9cAwOeff+67YcMGz27dujWsWbMms7HWAus8g5IDIUQDgBgicgOwmYiihRApHR1MCLEcwHJAW1uho+9nrKXvk69hQKAbgr26Sdann5N0VxeolArc3d8fG49no7JWDSc7nvVkZqqjJZebtvf3V3f4/U20VbJ5zZo1HuXl5crAwMB+AFBZWalcs2aN5+DBg3MA7ZqDt956K/9Wx2Vt61CaJYQoBbAXwPgWL+UACAQAIrIB4AqgWIL4GGvTufwKpOWW476Y7pL2uzVjK7ZmSFfj4L6YANTUa/Dz2TzJ+mTMmrRVsnnjxo0emzdvPp+Tk3MmJyfnzJEjR1K///57d7ni7EoMuVrBWzdjACJyAHAHgPQWzbYAeFh3fxqAPUKuco+sy/j+ZA6UCsI9/aVNDhYfXozFhxdL1t+gnu7o4e6A75OvSdYnY9aktZLNGRkZtjk5ObZjx469uQYhIiKiztnZuWHPnj3STRWyVhkyx+kPYLVu3YECwAYhxDYiegtAkhBiC4D/AlhLRBcAXAcww2gRMwZAoxH4IfkahvfxgreznaR9b5z+x3oJrTc0rB0RYXJMdyxNvIjCito2491oYH+MWYuqqqqTLZ+75557Ku65554KACgoKDjd8vXU1NQ0AGiaNDDp6U0OhBCnAQxs5fk3m9yvAXC/tKEx1rbjV0uQU1qNl+4Mk7xvL8c/1ktovaGB7aA9tfDZ3ovYdvoaHhneemEorw70xxhjxsRLO5lF+v5kDuxVCtxphK2JN6Vtwqa0TfobtlNboaVQX2dE+bu0e2ph1apVWGVgf4wxZkycHDCLU6fW4Mczubgjys8oq/8/PvIxPj7ysf6GHUgOAOC+gd1xKqsUl4panw3l5IAxZi44OWAW59cLRSitqsfkAdIuRGz0w4wf8MOMH/Q3TEzU3gx074AAAMC2U7wwkTFm3jg5YBbnp5RcONvZYGSYcc7Ru9q7wtVe+houfq72GNTTHT+l8CWNjDHzxskBsyj1DRr8nJqPcZE+sLNRGmWM9SnrsT7FgBoHixZpbx0wIdoPqbnluNzGqQXGGDMHnBwwi/JbZjFKq+oxoZ+/0cZYmrQUS5MMqHGwbZv21gGNcfPsAWOAUqkc1FhyecKECSEVFRWtfiYNHDgworXnDRUQENBv0KBBzeo9NI7bmX6tGScHzKL8lJIHR1slRod5G22M7Q9ux/YHtxul7wA3Bwzo4YodKbl/HHf7dmzfbpxxGTNHdnZ2mvT09NTz58+fValUYvHixc3+Y9fXa6s3nzx5suXGe21qfE9LN27cUF64cEEFACdOnLC/9ai7Bk4OmMVo0AjsTMnD2Agf2KuMc0oBABxVjnBUORqt/wn9/HEquwzZJVXNx3V0hKOj8cZlzJyNGDGi8sKFC3bbtm1zHjRoUPjYsWP7hIaGRgOAo6PjQADQaDR44okneoSGhvYNCwuLWrFihTsAtPaelu67776bJZ/XrFnTrDS0Wq3GE0880SM6OjoyLCws6v333/cCgLKyMsXQoUPDoqKiIsPCwqK++uorN0C7e2NISEjfGTNm9OzTp0/f4cOHh1ZWVlpVJWJODpjFOHrpOopv1GFCtPFOKQDAV6e/wlenvzJa/xOitXsz7GhxamHJkiVYsmSJ0cZlrD1xK+LCPz6iLcFcq66luBVx4UuOLfEAgIraCkXcirjwFSe0H8bFVcXKuBVx4auTV7sBQG5Frk3cirjw/535nysAXC272qFrjOvr67Fz506Xfv36VQNAamqq45IlS65evny5WYG/NWvWuJ05c8YhLS3t7C+//HLuzTff7HHlyhVVe+9pNHPmzJKtW7e6A8DOnTvdpkyZUtr42ocffujl6urakJKSknbq1Km01atXe6enp9s6Ojpqfvzxxwupqalp+/btOzd//vweGo0GAHD16lX7559/vuDChQtnXV1dG9asWWNVNR84OWAWY0dKLuxVCsSHG++UAgB8ceILfHHiC6P139OzG6L8Xf6w7mDDhg3YsGGD0cZlzNzU1tYqIiIiovr16xfVo0ePuhdeeKEIAPr3738jIiKirmX7AwcOOE+fPv26jY0NAgMD1UOGDKk8ePCgY3vvaeTj49Pg6uqqXr58uXufPn2qnZycNI2v7d6922XDhg2eERERUQMHDowsKSmxSU1NtddoNDRv3rweYWFhUWPGjAkrKCiwzc7OtgGAgICA2mHDhlUDwMCBA6suX74s7T7uMuP6scwiaDQCP6XkIT7MB92MXPZ41+xdRu0fACb288Oin88hr6wGfq58+pPJ7+hffi+5bGdjJ5o+drZz1jR97Ono2dD0sb+zv7rp4yDXILUhYzauOWj5vKOjo6a19u0x5D3Tpk0r+dvf/tZzyZIll5o+L4SgxYsXX506dWp50+c//vhjz+LiYpszZ86k2dnZiYCAgH7V1dUKALC1tb1ZXFCpVIrG562FVR0Ms14ns0pQUFGLCf2k3y65JZVSBZVSZdQxGq9aaG1hImOsdaNGjarYuHGjh1qtxrVr12yOHj3qNHLkSIOvC37wwQdLnnnmmbwpU6Y0SwLuuOOOsqVLl3rX1tYSAJw+fdquvLxcUVZWpvTy8qq3s7MTW7dudb527Zqt1MdkrnjmgFmEn8/mQ6UkjInwMfpYq5JXAQDmxMwx2hi9vZ3Qx8cJu9LyMaeNQkyMseZmz55deujQIafIyMi+RCQWLlyYHRQUpD59+g/FG1vl7u6u+ec///mH64j/+te/Fl2+fNmuX79+kUII8vDwqN++ffvFxx9//PqECRP6hIWFRfXv37+qV69eNZIflJkiIYT+VkYQGxsrkpKSZBmbWZ6xixIR4O6AtY8NMfpY8aviAQCJcxL1NNS268gWyk2991M6vjiQieN/vwOuDirE6/pLvMX+WNdARMeFELGd7efUqVOXBwwYUCRFTMwynTp1ymvAgAHBrb3GMwfM7F0oqERm0Q3MGR5skvH0JgU3GxrYrg13RPni830XkZhRgMkxAZwUMMbMBq85YGZvV2o+ACAh0lfmSKQ1MNANXk52+Fl3fIwxZi44OWBmb3daPqIDXNDdzcEk4604vgIrjq/Q3/AWais0pVAQEiJ9sC+jEHVqDRYtWoRFneiPsQ7SaDQaq9q4hxlO93ff5hUenBwws1ZYUYsTV0tMOmuw/ux6rD9rQOGlw4e1t064I8oXlbVq/JZZjG3btmFbB2s1MNYJKYWFha6cIHQ9Go2GCgsLXQG0umEUwGsOmJnbk54PIbQfoqay+6HdhjX87rtOjzW8jxccVMqbp04YMxW1Wv14Xl7eF3l5edHgL4pdjQZAilqtfrytBnqTAyIKBLAGgC8AAWC5EOKjFm3iAfwAoHFjiU1CiLduLWbGfrcrNR8Bbg6I8neROxSjsFcpMSrMC7vT8mFV26sxszdo0KACAPfKHQczT4bMHKgBvCSEOEFEzgCOE9EuIUTLXa0OCCHukT5E1lVV1alx4HwRZsYFgch0M59LjmnrGzw9+On2G772mvbnv/7VqfHuiPLDzrP58K5VG333R8YYM4Te30RCiFwAubr7FUSUBiAAwB+2vGRMSgfPF6FWrTH5VQpbz20FYEBy0Mn1Bo3GRvhAQUClWgEvEy26ZIyx9nToawoRBQMYCOBIKy8PJaJTAK4BeFkIcbbz4bGubE96AZztbBDXy8Ok4/704E8mHc+jmy0G9XRH9UP/xLbnRpp0bMYYa43Bi1CIyAnAdwDmCSHKW7x8AkBPIcQAAJ8A+L6NPuYSURIRJRUWFt5iyKwr0GgE9qQXYFSYN2xtrH+t1NgIX6TklCOvrMvszsoYM2MG/dYlIhW0icHXQohNLV8XQpQLISp197cDUBGRVyvtlgshYoUQsd7exi27yyzb2WvlKKioNUkthZY++u0jfPTbR/obSmhshA9Kf/0Gz/3tdZOOyxhjrdGbHJB2Jdh/AaQJIT5oo42frh2IKE7Xb7GUgbKuZU96AYiA+HDTJ5G/XPoFv1z6xaRjhvk6QeSk4OD+RJOOyxhjrTFkzcFwALMBnCGiZN1z8wEEAYAQ4nMA0wA8RURqANUAZgi5Kjoxq7AnPR8xuu2FTW3LzC0mH5OI4OaoQmFFLWrqG2CvUpo8BsYYa2TI1QoHAbR7HZkQ4lMAn0oVFOvaCipqcCq7DC/dESZ3KCbl7miL/PIa/JZZjPhw059OYYyxRta/0otZnMQM7WLVsZHyfEAuOrQIiw6ZvsaBi4MNFETYm15g8rEZY6wpTg6Y2dmTVgA/F3vZdkU8nH0Yh7MN2MPA01N7k4i3lxd8vb3wS3oB+KwcY0xOvB0bMyu16gYcOF+Ie2MCTLorYlPfTTewZoIEtRWad/cdvj5yBa9vTsH5gkqE+TpL2j9jjBmKZw6YWTl2qQQ36howToZLGM3BWN1x/5LGpxYYY/LhmQNmVvakF8DWRoFhfaSbru+o9w6+BwB4dcSr7TeUqLbC791p+4v0n4i9GQV4Kr63JP0yxlhHcXLAzEpiRgFuD/GEo618/zST85INa1gs7VYeh3W1Gia++jCW7c9EWXU9XB1Uko7BGGOG4OSAmY0rxTeQWXQDs4f2lDWOddPWGdZw+XKjjD8mwgdLEi/i1wtFmNjP3yhjMMZYe3jNATMbjZcwjuni1/gPDHSDi70NX9LIGJMNJwfMbOzNKEAvr24I9uomaxxv73sbb+97W3/DuXO1N4nZKBUYFeaNxHOF0Gj4kkbGmOnxaQVmFqrrGnD4YjH+PCRI7lCQUZxhWMNz5yQdt0ePHjfvjwn3wbbTuUjNLUd0gKuk4zDGmD6cHDCz8FtmMWrVGrM4pfDVlK/kGfer38cdrSs4tTe9gJMDxpjJ8WkFZhb2ZhTAQaVEXC8PuUMxC15OdhjQwxV7M3jdAWPM9Dg5YLITQiAxoxDD+3iaRTXCN/e+iTf3vmnycefNm4d58+bdfBwf7oPkrFKU3KgzeSyMsa6NkwMmu8yiG7h6vQqjzeCUAgBklWchqzzL5OMmJycjOTn55uP4cG9oBLD/fKHJY2GMdW285oDJrvGSvfgwb5kj0fpy8pdyhwAA6N/DDR7dbLE3vQCTYwLkDocx1oXwzAGTXWJGIfr4OCHQw1HuUMyKUkEYHeaN/eeL+JJGxphJcXLAZHWjVo2jl65jTLh5zBoAwGu7X8Nru1+TOwwA2lML12/U4XROmdyhMMa6ED6twGR1+GIx6ho0iDeT9QYAUFxtYM2EsDBJxw1rpb+Rod4g0taciAl0k3Q8xhhrCycHTFZ7MwrQzVaJ2GB3uUO5afkkA2smSFxbYXkr/Xl0s0VMoBv2ZhRiXoK0yQhjjLVF72kFIgokor1ElEpEZ4nohVbaEBF9TEQXiOg0Ed1mnHCZNWm8hHFYHy/Y2ch/CaO5ig/zwensUhRX1sodCmOsizBkzYEawEtCiCgAtwN4hoiiWrSZACBUd5sLYKmkUTKrdLGwEjml1Yg3o/UGAPDyzy/j5Z9f1t9Q4toKc+fOxdxW+osP94YQwIHzRZKNxRhj7dGbHAghcoUQJ3T3KwCkAWh5XdVkAGuE1m8A3IiIa82ydjVWYTSn9QYAUF1fjer6av0NPT21N4mcO3cO51qp19AvwBWe3WyRyLslMsZMpENrDogoGMBAAEdavBQAoOmuMdm653I7ExyzbnszChDm64QANwe5Q2nms7s/M6zhv/5l3EB0FLpLGvdmFKBBI6BUkEnGZYx1XQZfykhETgC+AzBPCFF+K4MR0VwiSiKipMJC3vWtK7tRq8axSyVmN2tgrkaHe6Okqh6ns0vlDoUx1gUYlBwQkQraxOBrIcSmVprkAAhs8riH7rlmhBDLhRCxQohYb2/zOs/MTOtQ4yWMZrIrYlPzdszDvB3z9DecOlV7M4FRod5Q0O+nYhhjzJgMuVqBAPwXQJoQ4oM2mm0B8JDuqoXbAZQJIfiUAmvT75cwWnAVxuJi7U0iMTExiImJafU1d90ljbzugDFmCoasORgOYDaAM0SUrHtuPoAgABBCfA5gO4CJAC4AqALwiOSRMqshhMC+jEIM7+MFWxvz26Tzw/EfyjPuh+2PGx/ug//sPoeiylp4OdmZJijGWJekNzkQQhwE0O4KKCGEAPCMVEEx63a+QHsJ47Nj+8gdikUZE+6DD3adw/5zhZhyWw+5w2GMWTHz+9rGrN7NKoxmtr9Bo2d+fAbP/Gj6XHfWrFmYNWtWm6/37e4CLyc77OV1B4wxI+Ptk5nJ7c0oQISfM/xdzesSxkYOKnniys7Obvd1hYIQH+6NXan5UDdoYKPk3J4xZhz824WZVHlNPZIul2BMhPlewrjozkVYdOciucNo1ZhwH5RV1yM5q1TuUBhjVoyTA2ZSv54vglojzPISRkswItQLSgVhL1+1wBgzIk4OmEklZhTC2d4Gt/U0nyqMLc3dOhdzt0pXM0FKrg4qDApy5/0OGGNGxWsOmMkIIbA3owCjQr2hMuPz5Z4OBtZLGDpU0nGHGthffIQ3/t+ODOSX18DXxV7SGBhjDODkgJlQam45CipqzfYqhUb/SjCwZoLEtRX+ZWB/Y8J98P92ZGBfRiGmDw7U/wbGGOsg8/36xqxO41T4aDNPDsxdhJ8z/Fzsed0BY8xoODlgJrMnvQDRAS7wcTbvqfBHfngEj/xgwCafEtdWmDp1KqYa0B8RYUyENw6cL0KdWiPZ+Iwx1oiTA2YSJTfqcPJqCcZG+Modil6BLoEIdDFgun7oUEnXHRQXF6PYwFoNY8J9UFmrRtLl65KNzxhjjXjNATOJfecKoRHAWDPe36DRW2PeMqzhyy8bN5B2DO/jBVulAr+kF2BYHy/Z4mCMWSeeOWAm8Ut6AbycbNE/wFXuUKxCNzsb3N7b8+ZW1IwxJiVODpjRqRs02JdRgPhwHygU7dbwMguzNs3CrE1t1zi4KT5ee5PJ2HBvZBbdwKWiG7LFwBizTpwcMKM7cbUU5TVqjLOAUwoAEO4ZjnDPcJOPO27cOIwbN87g9o3rN/bw7AFjTGK85oAZ3S/p+VApCSNCLePc+N9H/12ecf/esXGDPB3Rx8cJe9ML8NiIXkaKijHWFfHMATO6vekFiOvlAWd7ldyhWJ1xET44cqkYlbVquUNhjFkRTg6YUWVdr8K5/EqMCbeMUwoAMGPjDMzYOMPk406YMAETJkzo0HvGRPigvkHg4HmutcAYkw6fVmBG1biL37hI89/foFGMX4ws41ZXV3f4PYN6usPF3ga/pBVgfLS/EaJijHVFnBwwo9qTXoBeXt3Qy6ub3KEY7NURr8odgsFUSgVGhXljb0YhNBphEVeDMMbMH59WYEZzo1aNQxeKLeYqBUuVEOmLospanMoulTsUxpiV0JscENFKIiogopQ2Xo8nojIiStbd3pQ+TGaJDpwvRF2DBglRlnNKAQCmbpiKqRukq5lgbPHh3lAqCLvT8uUOhTFmJQw5rbAKwKcA1rTT5oAQ4h5JImJWY1dqAVwdVIjt6S53KB0ytIeB9RLukfaf/D232J+boy0GB7tjd2oBXrkrQtKYGGNdk97kQAixn4iCTRALsyINGoE96fkYG+EDG6Vlnb16eZiBNRMkrq3wcif6S4j0xTs/puFqcRWCPB0ljIox1hVJ9Vt7KBGdIqKfiKhvW42IaC4RJRFRUmEhX3plzY5fKUFJVT0SLOgqBUt2h+7UDZ9aYIxJQYrk4ASAnkKIAQA+AfB9Ww2FEMuFELFCiFhvb28JhmbmaneadlfEUWGWsStiU/d+cy/u/eZe/Q0lrq0QHx+P+Fvsr6dnN4T6OHFywBiTRKcvZRRClDe5v52IlhCRlxCiqLN9M8u1OzUft4d4WuSuiON6GVjfYM4co8bRUQlRvli+PxNlVfVwdbS8P3fGmPnodHJARH4A8oUQgojioJ2NKO50ZMxiXSysRGbRDcwZHix3KLfkhdtfMKyhuSUHkb5YmngRiecKMDkmQO5wGGMWTG9yQETfAIgH4EVE2QAWAFABgBDicwDTADxFRGoA1QBmCCGE0SJmZm93qnZq25J2RbwlRbrJMS/zOHUSE+gGLydb7E7j5IAx1jmGXK0wU8/rn0J7qSNjAIBdqfmI9HdBgJuD3KHckglfa+sb/PTgT+03nDZN+zMx0bgBGUipIIyN8MFPZ/JQp9bA1sayrhJhjJkP3j6ZSaqgogbHr5bghXGhcodyyyaFTZJl3OnTp3e6jzuj/LAhKRuHM4sxOowX/TLGbg0nB0xSu1LzIQQwPtpP7lBu2dODn5Zn3Kc7P+6IUC842iqxIyWPkwPG2C3jeUcmqR0peQj2dES4r7PcoVicqqoqVFVVdaoPe5USYyJ8sCs1Dw0aXvrDGLs1nBwwyZRV1+PwxWLcFe0HIsutDpiwJgEJaxJMPu7EiRMxceLETvczvq8fiirrcOJqiQRRMca6Ij6twCSzJz0fao3AXX0t95QCADzQ9wG5Q+iU+HBv2CoV2JGSh8HBHnKHwxizQJwcMMnsSMmDr4sdYnq4yR1Kp/xl0F/kDqFTnO1VGBHqhR0peXjj7kiLnsVhjMmDTyswSVTVqbHvXCHu6usHhYI/jOQ2vq8fckqrcfZauf7GjDHWAicHTBL7zxWipl6D8RZ+SgEA4lfFI35VvNxhdEpClC8UpJ3NYYyxjuLTCkwSO8/mw81Rhbheln+Oe07MHAMbGtjO0HEl7M+jmy2G9PLEzrN5ePmucMn6ZYx1DZwcsE6rVTdgd2o+xkf7wUZp+ZNR1pAcANq9JhZsOYvz+RUI5UtLGWMdYPm/yZns9mUUoqJWjbv7+8sdiiTqG+pR31Cvv2FR0e/1FSRQVFSEIgn7mxDtByJg6+lcyfpkjHUNnBywTtt2OhfujioM72MeBYg66461d+COtXfobzht2u/1FSQwbdo0TJOwPx8Xewzp5YFtp6+Ba6ExxjqCTyuwTqmua8DutHxMjukOlRWcUgCAx2973LCGL71k3EAkcE//7njj+xSk5VYgqruL3OEwxiwEJwesU/ZmFKCqrgH39O8udyiSmdV/lmENJ8lToKkjJujWHWw7fY2TA8aYwazjqx6TzbbT1+DlZIshVnCVQqOq+ipU1RtQ4yAjQ3szY55OdhjW2xPbTufyqQXGmME4OWC37EatGnvSCzAh2t8qrlJoNPHriZj4tQE1Dp54Qnszc/f098fV61U4k1MmdyiMMQvBpxXYLdudlo+aeg3usZKrFBo9FfuUPOM+ZZxx7+rrh9c3p2Db6Vz0t/CtrRljpsHJAbtl207nwtfFzuqK+zwQLU/hpQceMM64bo62GBnqhR9P5+K1CRFca4Exppf1zAUzkyqrqse+jEJM7OdvdbUUymrKUFZj+in4rKwsZGVlGaXvSQO6I6e0GsevcBlnxph+epMDIlpJRAVElNLG60REHxPRBSI6TUS3SR8mMzc/nslFXYMGUwb2kDsUyU1eNxmT1002+bizZ8/G7NmzjdL3nX394KBSYtPJHKP0zxizLobMHKwCML6d1ycACNXd5gJY2vmwmLnbfDIbfXycEB1gfZfHPT/keTw/5Hm5w5CUk50N7urri22nrqFW3SB3OIwxM6c3ORBC7AdwvZ0mkwGsEVq/AXAjIutaocaauVpchWOXSzDltgCrPH89JXIKpkROkTsMyU25rQfKa9TYk1YgdyiMMTMnxZqDAABNT5Rm6577AyKaS0RJRJRUWFgowdBMDptP5oAIuC+m1b9mi1dUVYSiKulqHJiL4X284ONsx6cWGGN6mXRBohBiuRAiVggR6+3tbcqhmUSEENh0MhtDQzzR3c1B7nCMYtqGaZi2QboaB+ZCqSDcNzAAe9MLcP1GndzhMMbMmBSXMuYACGzyuIfuOWaFTlwtwZXiKjw7po/coRjNS0MNrJkgcW2Fl0xQq+FPAwOwfH8mtp66hoeHBRt9PMaYZZIiOdgC4FkiWgdgCIAyIQTXiLVSm07kwF6lwIR+1rusZFK4gTUTJK6tMMkEtRoi/V0Q6e+CTSdzODlgjLXJkEsZvwFwGEA4EWUT0WNE9CQRPalrsh1AJoALAFYAeNpo0TJZ1dQ3YNvpXNzV1w9Odta7f1ZeZR7yKvP0N5S4tkJGRgYyTFCrYcrAAJzKKsWFgkqjj8UYs0x6f8MLIWbqeV0AeEayiJjZ2nk2D2XV9Zg2yPr2NmhqxsYZAIDEOYntN2ysq5Cop52BntD1lyhRf225b2AA/r0jHeuPXcXrd0cZdSzGmGWy3q9/THL/O3IVQR6OGN7bS+5QjOrVEa8a1vDdd40biJF4O9vhjihfbDyejZfvCoedjVLukBhjZoa3T2YGuVhYiSOXrmNGXKDVbZfc0vg+4zG+T3v7fukMG6a9WaCZcUEoqarHzrP5cofCGDNDnBwwg6w7ehU2CsL9gwL1N7ZwWWVZyCozoMbBoUPamwUa0ccLgR4O+ObIVblDYYyZIU4OmF616gZsPJ6NO/v6wtvZTu5wjG725tmYvdmAGgfz52tvFkihIMwYHITDmcXILOSFiYyx5njNAdNr59l8lFTVY2ZckNyhmMQbo96QZ9w3TDvu/bE98J9d57D+WBZemxhp0rEZY+aNkwOm1zdHriLQw8HqFyI2SghJkGfcBNOO6+Nsj4RIX3x7PBsv3hnGCxMZYzfxaQXWrgsFlTicWYwZg4OsfiFio8ySTGSWZJp83OTkZCQnJ5t0zJlDgnD9Rh12pBiwrwNjrMvgmQPWrpW/XoKtjQIzBlv/QsRGj/7wKAAD9jmQ2Lx587TjGnmfg6ZG9vFCiHc3rDx4CfcO6G6VVTYZYx3HyQFrU8mNOmw6kY0pAwPg6WT9CxEbLYxfKHcIJqNQEB4Z3gt//z4Fx6+UIDbYQ+6QGGNmgE8rsDb97+hV1NRr8MjwXnKHYlKjg0djdPBoucMwmam3BcDVQYWVv16SOxTGmJng5IC1qk6twZrDlzEy1Avhfs5yh2NSGUUZyCgyfo0Dc+Foa4OZcUHYkZKHrOtVcofDGDMDnBywVm0/k4v88lo8OqJrzRoAwBPbnsAT256QOwyTenhYTxARVh+6LHcojDEzwGsO2B8IIbDy10sI8e6G0aHecodjcu+OM7BmgsS1Fd6VsVaDv6sDJvbzx/pjWZh3R5hVV91kjOnHvwHYH/yWeR2ns8vwzn3RXebyxaaGBRpYL0HiugrDZK7T8NiIXth66hq+OXIVfxkVImssjDF58WkF9gcf/3IePs52Vl+auS0pBSlIKUjR31Di2gqHDh3CIRlrNcQEumF4H08s25+JmvoG2eJgjMmPkwPWzNFL13E4sxhPjO4Ne1XX3DHv2e3P4tntz+pvKHFthfnz52O+zLUanh8biqLKWnxzlAsyMdaV8WkF1swne87Dy8kWf+4idRRa8/4d7xvWcNky4wYigyEhnhjSywOf77uImXFBXTZBZKyr45kDdtPxKyU4cL4Ic0eFwMG2634oDA4YjMEBg/U3DA/X3qzMC+NCkV9ei2+TDChbzRizSpwcsJs+2XMeHt1s8eCQnnKHIqvkvGQk5yXrb7h1q/ZmZYb29kRsT3csSbyIWjWvPWCsKzIoOSCi8USUQUQXiOjVVl6fQ0SFRJSsuz0ufajMmE5cLUFiRiEeH9kL3br4ZWzzdszDvB3z9DdcvFh7szJEhOfHhSK3rAYbjvHsAWNdkd5PASJSAvgMwB0AsgEcI6ItQojUFk3XCyEMWMXFzI0QAv/8MQ1eTnZ4aGiw3OHI7sPxH8oz7ofyjNuakaFeiOvlgQ93n8d9AwPgbK+SOyTGmAkZMnMQB+CCECJTCFEHYB2AycYNi5nSTyl5OH6lBC/dyZvfAECMXwxi/GJMP25MDGJiTD9ua4gIb9wdieIbdViaeFHucBhjJmZIchAAoOncYrbuuZamEtFpItpIRK3W9yWiuUSURERJhYWFtxAuk1qtugHv/ZSOcF9nTI/tOmWZ23Ms5xiO5Rwz+bi7d+/G7t27TT5uW/r3cMN9Md3xxcFLyC7hmguMdSVSLUjcCiBYCNEfwC4Aq1trJIRYLoSIFULEent3vW15zdGaQ1dw9XoVXr87EsouuBtia17Z9Qpe2fWKycd955138M4775h83Pa8Mj4CBOD9nV2nEBVjzLDkIAdA06+UPXTP3SSEKBZC1OoefgFgkDThMWMquVGHT/acx+gwb4wK42St0acTP8WnEz+VOwyzEODmgMdG9MIPydeQnFUqdziMMRMxJDk4BiCUiHoRkS2AGQC2NG1ARP5NHt4LIE26EJmxvLs9DTfqGjB/YqTcoZiVaJ9oRPtEyx2G2Xgqvje8nOzwxvdnoG7QyB0OY8wE9CYHQgg1gGcB7IT2Q3+DEOIsEb1FRPfqmj1PRGeJ6BSA5wHMMVbATBoHzxfh2+PZmDsqBOF+znKHY1YOZR3CoSz5ahyYG2d7Ff5xbxRScsrx34OX5A6HMWYCBi1NF0JsB7C9xXNvNrn/GoDXpA2NGUt1XQNe23wavby64YVxoXKHY3bm/6Ktb5A4J1HeQMzI3f388X3kNXyw6xzu6uuHYK9ucofEGDMivm6tC/pgVwayrldj3dzbee/8Viy7x8CaCRLXVlhmxrUaiAjv3BeNOz7Yh/mbz+Drx4eAiBewMmatePvkLiY5qxT/PXgJM+OCcHuIp9zhmKVwr3CEexlQM0Hi2grh4eEIN+NaDX6u9nh1YgQOXSzGOt45kTGrxslBF1JWXY/nvjkBPxd7vDohQu5wzNa+y/uw7/I+/Q0lrq2wdetWbDXzWg0zBwdhaIgnFm49i3P5FXKHwxgzEhJCyDJwbGysSEpKkmXsrkgIgSe/Oo5f0gqw/omhGNTTXe6QzFb8qngABqw5iNe2Q6KedoaOq+svUaL+jKWgvAYTPz4AN0dbbHl2OBxt+eykKRHRcSFErNxxMOvG/6u7iC9/vYydZ/Px+sRITgz0WDl5pWENN240biBmysfFHh8+MBCzVx7BG9+nYPH9A3j9AWNWhk8rdAEnrpbgXz+lISHSB4+P7CV3OGYvxD0EIe4h+ht6eWlvXdCIUC88NzYUm07kYD2vP2DM6nByYOUuF93AX1Ynwc/VHov4G55Bdmfuxu5MA2ocrFqlvXVRL4wLxYg+Xnjj+xQcOM+1UhizJpwcWLGiylo8/OVRaITA6kfi4OZoK3dIFuGd/e/gnf0G1Djo4smBUkFYMus29PFxwpNrjyMlp0zukBhjEuE1B1bqRq0aj646hvzyGvzvL7cjxNtJ7pAsxto/rZVn3LXyjNsZLvYqrH40DlOWHMIjq45h01PDEOjhKHdYjLFO4pkDK1ReU49HVh1DSk4ZPp15G24L4gWIHRHoGohAV9OXrw4MDERgoOWVzfZ1scfqRwejTq3BzBW/4UrxDblDYox1EicHVqa4shZ/XvEbTlwpwYczBiIhylfukCzOjgs7sOPCDpOPu379eqxfv97k40qhj48z1j4Whxu1akz7/DDScsvlDokx1gmcHFiRnNJq3L/sMM7nV2LFQ7G4d0B3uUOySO8dfA/vHXzP5OMuXboUS5cuNfm4Uunfww0bnhgKJREeWHYYx69clzskxtgt4uTASuw/V4hJnxxEYXkt1j42BGMifOQOyWKtm7YO66atkzsMixTq64xvnxwKj262mLn8CNYevgy5NlpjjN06Tg4sXING4INd5/Dwl0fh5WSLzc8MR1wvD7nDsmh+Tn7wc/KTOwyLFejhiM1PD8ewPp74+w9n8fy6ZFTWquUOizHWAXy1ggVLzyvHG5tTkHSlBNMG9cDbk6PhYMtVFjtra4a2vsGk8EkyR2K53LvZYuXDg7F030Us/jkDp7NLsfDevogP5xktxiwBJwcWqKKmHh/uPo9Vhy7D1UGFD6YPwJTbesgdltVYfHgxAE4OOkuhIDwzpg8GB3vg1U2nMefLYxjf1w9/nxSFADcHucNjjLWDCy9ZkLKqeqw+fBlf/noJpdX1mBkXhL/dFc6bG0msqKoIAODlqGdr5CJtO6m2UC7S9edlhVsy16ob8MWBS/hkz3loBPBAbCDmjgrhPRFuARdeYqbAyYGZE0IgNbcc3x3PwfpjV3GjrgEJkT54flwo+vdwkzs8xjoku6QKn+29gI3Hs6ERwD39/fFAbCBuD/GEQsFbexuCkwNmCpwcmCGNRpsQ7DtXiC3J15CRXwGVkjAh2h9PxfdGpL+L3CFatU1pmwAAUyKntN+wcevkOXMkGXeVrr85EvVnznLLqrFi/yVsSMpCZa0a/q72uHdAd4yJ8MFtQe6wteG10m3h5ICZgkHJARGNB/ARACWAL4QQ77V43Q7AGgCDABQDeEAIcbm9Pjk5+F15TT3O5pQjJacMp7JLcfhiMYpv1AEAYgLdMPW2ANzTvzvcu/HpA1OIXxUPAEick6inobYdEvW0M3RcXX+JEvVnCarrGrArLR+bT2Rj//kiNGgEHG2VuD3EEwN6uKF/D1f0DXCBt5MdFw3T4eSAmYLeBYlEpATwGYA7AGQDOEZEW4QQqU2aPQagRAjRh4hmAPg3gAeMEbAlaNAIVNWpUV3fgKraBpTX1KOsuh6lVfUorKhFQUUtCsprcOV6Fa4U30BRZd3N93Z3tcfIUC+MDPXGyFAv+LjYy3gkXdMPM34wrGEX+hA3FgdbJe4d0B33DuiO8pp6HL5YjAPnC3H4YjH2ZhSg8buLs70Ngj27oaenI/xc7OHjYgcfZ3u4Oarg6qCCi4MKTnY2sFcp4WirhErJMw+MdYYhVyvEAbgghMgEACJaB2AygKbJwWQA/9Dd3wjgUyIiYYRzFvvOFeKdban6GxqgveCEEL+/LrRtG5/TCAGNRvtTrRHQaATqGzSoa9CgvkGgQdP+YauUBG8nOwR6OCIh0hc9Pbsh0t8Z0QGu8HKyk+TY2K1ztXeVO4QuycVehbv6+uGuvto9Jipr1TibU4az18pxufgGLhdX4UxOGXan5aOmXtNuXwoCVEoFbG0UUCkVUCoINgqCgghEuPmTABARCAB0j6F7rjVSzV08MDgQj48Mkag3xqRnSHIQACCryeNsAEPaaiOEUBNRGQBPAEVNGxHRXABzASAoKOiWAnays0Gor3QVBqm9/+4tfllof5EASiIQERQE2ChJ94un8RcRwVap/fZib6tEN1slXOy132xcHVTwdraDm4OKF1+ZsfUp2voGD0TrmfxatEj78+WXjRxR1+RkZ4MhIZ4YEuLZ7HkhBCpq1Sgor0VZdT3Kq7Uzc1V1DdoZu7oG1DVoUKfWoFatQYNGm8SrGzTQiOZJvmiR+GsHaD0e0e7XiY7hLwHM3Jl0nwMhxHIAywHtmoNb6WNQT3cM6jlI0rgYa2ppkra+gd7kYNs27U9ODkyKiLQJt71K7lAYs1qGJAc5AJrWke2he661NtlEZAPAFdqFiYxZnO0Pbpdn3O3yjMsYYy0ZsmrnGIBQIupFRLYAZgDY0qLNFgAP6+5PA7DHGOsNGDMFR5UjHFWm35zH0dERjo68KRBjTH56Zw50awieBbAT2ksZVwohzhLRWwCShBBbAPwXwFoiugDgOrQJBGMW6avTXwEAZvWfZdJxlyxZAgB4+umnTTouY4y1ZNCaAyHEdgDbWzz3ZpP7NQDulzY0xuTxxYkvAJg+OdiwYQMATg4YY/LjwkuMtbBr9i65Q2CMMVlxcsBYCyolr4JnjHVtvI0YYy2sSl6FVcmr5A6DMcZkw8kBYy1wcsAY6+pkq8pIRIUArtzi273QYvfFLoCPuWvgY+4aOnPMPYUQ3lIGw1hLsiUHnUFESV2tKhkfc9fAx9w1dMVjZpaFTyswxhhjrBlODhhjjDHWjKUmB8vlDkAGfMxdAx9z19AVj5lZEItcc8AYY4wx47HUmQPGGGOMGQknB4wxxhhrxqyTAyIaT0QZRHSBiF5t5XU7Ilqve/0IEQXLEKakDDjmF4kolYhOE9EvRNRTjjilpO+Ym7SbSkSCiCz+EjBDjpmIpuv+rs8S0f9MHaPUDPi3HUREe4nopO7f90Q54pQKEa0kogIiSmnjdSKij3V/HqeJ6DZTx8hYm4QQZnmDtjz0RQAhAGwBnAIQ1aLN0wA+192fAWC93HGb4JjHAHDU3X+qKxyzrp0zgP0AfgMQK3fcJvh7DgVwEoC77rGP3HGb4JiXA3hKdz8KwGW54+7kMY8CcBuAlDZenwjgJwAE4HYAR+SOmW98a7yZ88xBHIALQohMIUQdgHUAJrdoMxnAat39jQDGERGZMEap6T1mIcReIUSV7uFvAHqYOEapGfL3DABvA/g3gBpTBmckhhzzXwB8JoQoAQAhRIGJY5SaIccsALjo7rsCuGbC+CQnhNgP4Ho7TSYDWCO0fgPgRkT+pomOsfaZc3IQACCryeNs3XOtthFCqAGUAfA0SXTGYcgxN/UYtN88LJneY9ZNtwYKIX40ZWBGZMjfcxiAMCL6lYh+I6LxJovOOAw55n8AmEVE2QC2A3jONKHJpqP/3xkzGS7ZbKGIaBaAWACj5Y7FmIhIAeADAHNkDsXUbKA9tRAP7ezQfiLqJ4QolTMoI5sJYJUQYjERDQWwloiihRAauQNjrKsx55mDHACBTR730D3XahsisoF2KrLYJNEZhyHHDCJKAPA6gHuFELUmis1Y9B2zM4BoAIlEdBnac7NbLHxRoiF/z9kAtggh6oUQlwCcgzZZsFSGHPNjADYAgBDiMAB7aAsUWSuD/r8zJgdzTg6OAQglol5EZAvtgsMtLdpsAfCw7v40AHuEEJa8q5PeYyaigQCWQZsYWPp5aEDPMQshyoQQXkKIYCFEMLTrLO4VQiTJE64kDPm3/T20swYgIi9oTzNkmjBGqRlyzFcBjAMAIoqENjkoNGmUprUFwEO6qxZuB1AmhMiVOyjGADM+rSCEUBPRswB2QrvSeaUQ4iwRvQUgSQixBcB/oZ16vADtwp8Z8kXceQYe8/sAnAB8q1t7eVUIca9sQXeSgcdsVQw85p0A7iSiVAANAF4RQljsrJiBx/wSgBVE9FdoFyfOseRkn4i+gTbB89Kto1gAQAUAQojPoV1XMRHABQBVAB6RJ1LG/oi3T2aMMcZYM+Z8WoExxhhjMuDkgDHGGGPNcHLAGGOMsWY4OWCMMcZYM5wcMMYYY6wZTg6Y7IiosgNt44loWCfGciOip9t5vYGIkpvcgttpe+hW42CMMXPGlzIy2RFRpRDCycC2/wBQKYRYdItjBQPYJoSI7mwsbbzfRlfngzHGLBbPHDCzRESTiOgIEZ0kot1E5Kv7YH8SwF913+pHEpE3EX1HRMd0t+G69/+DiFYSUSIRZRLR87qu3wPQW/f+9w2Iw4mIfiGiE0R0hogmN3mtUvcznogOENEWAKm6x4lEtJGI0ono68ZqoUQ0iIj2EdFxItrZWIWPiJ4nolQiOk1E63TPjW4yg3GSiJyl+xNmjLG28cwBk11r39aJyB1AqRBCENHjACKFEC+1nDkgov8BWCKEOEhEQQB2CiEide3uBDAG2voMGQD8oK16197MQQOAM7qHlwDcD8BRCFGu28b4NwChurgqhRBORBQP4EcA0UKIS7rHPwDoC23Z4V8BvALgCIB9ACYLIQqJ6AEAdwkhHiWiawB6CSFqichNCFFKRFsBvCeE+JWInADU8KwEY8wUzHb7ZNbl9QCwXvfN2hbaD+rWJACI0n0xBwAX3QcpAPyoK0xVS0QFAHwNGLdaCBHT+ICIVADeJaJRADTQJhe+APJavO+orkBS08fZuj6SAQQDKIW2iNQuXbxKAI176Z8G8DURfQ9tXQVAm1R8QERfA9jU2B9jjBkbn1Zg5uoTAJ8KIfoBeALaIjytUQC4XQgRo7sFCCEaFzg2rVjZgFtLhh8E4A1gkC5pyG8jlhstHrc2NgE42yTWfkKIO3Vt7gbwGYDbABzTrV14D8DjABwA/EpEEbcQP2OMdRgnB8xcueL38rUPN3m+AtrTBI1+BvBc4wMiitHTb8v3GxJHgRCinojGAOjZgfe2lAHAm4iGAtpZCSLqS0QKAIFCiL0A/k83phMR9RZCnBFC/BvaqoacHDDGTIKTA2YOHIkou8ntRQD/gLby5HEARU3abgXwp8YFiQCeBxCrW8iXCu2CxTbpKhv+SkQphixIBPC1rv8zAB4CkN7xw7s5dh20pcX/TUSnACQDGAbt6YWvdGOcBPCxEKIUwDxdnKcB1AP46VbHZoyxjuAFiYwxxhhrhmcOGGOMMdYMJweMMcYYa4aTA8YYY4w1w8kBY4wxxprh5IAxxhhjzXBywBhjjLFmODlgjDHGWDP/H2gJF/UDzMkYAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = torch.linspace(0.0, 1.0, 100)\n", "plt.plot(x, dist.Beta(10, 10).log_prob(x).exp(), label='Prior~Beta(10, 10)')\n", "plt.xlabel(\"Latent Fairness\")\n", "plt.axvline(mle_estimate, color='k', linestyle='--', label='MLE')\n", "plt.axvline(map_estimate, color='r', linestyle='-.', label='MAP')\n", "plt.axvline(0.5, color='g', linestyle=':', label='Prior Mean')\n", "plt.legend(bbox_to_anchor=(1.44,1), borderaxespad=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Doing the same thing with AutoGuides\n", "\n", "In the above we defined guides by hand. \n", "It's often much easier to rely on Pyro's [AutoGuide machinery](https://docs.pyro.ai/en/stable/infer.autoguide.html?highlight=autoguide). \n", "Let's see how we can do MLE and MAP inference using AutoGuides.\n", "\n", "To do MLE estimation we first use [`mask(False)`](https://docs.pyro.ai/en/stable/poutine.html?highlight=mask#pyro.poutine.handlers.mask) to instruct Pyro to ignore the `log_prob` of the latent variable `latent_fairness` in the model. \n", "(Note we need to do this for every latent variable.)\n", "This way the only non-zero `log_prob` in the model will be from the Bernoulli likelihood and ELBO maximization will be equivalent to likelihood maximization." ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "def masked_model(data):\n", " f = pyro.sample(\"latent_fairness\", \n", " dist.Beta(10.0, 10.0).mask(False))\n", " with pyro.plate(\"data\", data.size(0)):\n", " pyro.sample(\"obs\", dist.Bernoulli(f), obs=data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we define an [`AutoDelta`](https://docs.pyro.ai/en/stable/infer.autoguide.html?highlight=autodelta#autodelta) guide, which learns a point estimate for each latent variable (i.e. we do not learn any uncertainty):" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[iter 0] loss: 7.0436\n", "[iter 50] loss: 6.8213\n", "[iter 100] loss: 6.7467\n", "[iter 150] loss: 6.7319\n", "[iter 200] loss: 6.7302\n", "Our MLE estimate of the latent fairness is 0.598\n" ] } ], "source": [ "autoguide_mle = pyro.infer.autoguide.AutoDelta(masked_model)\n", "train(masked_model, autoguide_mle)\n", "print(\"Our MLE estimate of the latent fairness is {:.3f}\".format(\n", " autoguide_mle.median(data)[\"latent_fairness\"].item()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To do MAP inference we again use an `AutoDelta` guide but this time on the original model in which `latent_fairness` is a latent variable:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[iter 0] loss: 5.6889\n", "[iter 50] loss: 5.6005\n", "[iter 100] loss: 5.6004\n", "[iter 150] loss: 5.6004\n", "[iter 200] loss: 5.6004\n", "Our MAP estimate of the latent fairness is 0.536\n" ] } ], "source": [ "autoguide_map = pyro.infer.autoguide.AutoDelta(original_model)\n", "train(original_model, autoguide_map)\n", "print(\"Our MAP estimate of the latent fairness is {:.3f}\".format(\n", " autoguide_map.median(data)[\"latent_fairness\"].item()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also quickly verify that had we chosen a uniform prior in our original model, our MAP estimate would be the same as the MLE estimate." ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "def original_model_uniform_prior(data):\n", " f = pyro.sample(\"latent_fairness\", dist.Uniform(low=0.0, high=1.0))\n", " with pyro.plate(\"data\", data.size(0)):\n", " pyro.sample(\"obs\", dist.Bernoulli(f), obs=data)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[iter 0] loss: 6.7490\n", "[iter 50] loss: 6.7302\n", "[iter 100] loss: 6.7301\n", "[iter 150] loss: 6.7301\n", "[iter 200] loss: 6.7301\n", "Our MAP estimate of the latent fairness under the Uniform prior is 0.600 matching the MLE estimate\n" ] } ], "source": [ "autoguide_map_uniform_prior = pyro.infer.autoguide.AutoDelta(original_model_uniform_prior)\n", "train(original_model_uniform_prior, autoguide_map_uniform_prior)\n", "print(\"Our MAP estimate of the latent fairness under the Uniform prior is {:.3f} matching the MLE estimate\".format(\n", " autoguide_map_uniform_prior.median(data)[\"latent_fairness\"].item()))" ] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:root] *", "language": "python", "name": "conda-root-py" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 4 }