{ "cells": [ { "cell_type": "markdown", "id": "ad947718", "metadata": {}, "source": [ "## 2. Model building\n", "\n", "
\n", "\n", "A **Bayesian statistical model** is specified by the **likelihood model** (also called observational model) and **prior model**. This can be expressed as the joint probability density of the data, $x$, and model parameters, $\\theta$, such that\n", "\n", "
\n", "\n", "$$\n", "p(x, \\theta) = p(x | \\theta) p(\\theta).\n", "$$\n", "\n", "
\n", "\n", "We perform inference by conditioning on the data and calculating expectation values of the resulting posterior distribution. Recalling **Bayes' theorem**\n", "\n", "
\n", "\n", "$$\n", "p(\\theta | x) = \\frac{p(x | \\theta)p(\\theta)}{p(x)} \\propto p(x, \\theta),\n", "$$\n", "\n", "
\n", "\n", "we note that the results of our inference are *completely dependent on the assumed model*. Bayes' theorem is not some magical guarantee, and if our model does not represent the actual *data generating process* our inference will yield results that are meaningless in practice. Even worse, if we are not careful, a misspecified model may produce results that \"look good\" but are far from the truth - leading us to draw the wrong conclusions from our analysis. \n", "\n", "
\n", "\n", "Our workflow begins with building our first sensible model as a starting point, that we can later check and iterate upon to make sure our inferences are useful. The first model doesn't need to be perfect, and in fact should be as simple as possible while still including the scientific phenomenon of interest. This process is what we will focus on here.\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "be4cd990", "metadata": {}, "source": [ "### Going from a scientific question to a statistical model\n", "\n", "
\n", "\n", "So, model building is very important, but where to start? It can be helpful to break down the problem into the following components:\n", "\n", "* Scientific question\n", "* Data\n", "* Parameter(s) of interest\n", "* Key nuisance parameter(s)\n", "* Relevant domain expertise\n", "\n", "The model-building problem then comes down to **defining a connection between your parameters and data** (the likelihood/observational model) and **incorporating relevant domain expertise** (the prior model). \n", "\n", "To get a feel for how this can work, consider the following examples:\n", "\n", "
\n", "\n", "\n", "\n", "
\n", "\n", "**Exercise 1 (10 points):** Fill out this table for some more complicated examples such as:\n", "\n", "* What is the value of the Hubble constant?\n", "* What is the mass of the Higgs boson?\n", "* An example from your own field/interests\n", "Below are some relevant plots of data used to make such measurements.\n", "\n", "The original Hubble diagram from [this paper](https://www.pnas.org/content/15/3/168):\n", "\n", "\n", "\n", "First presentation of ATLAS Higgs to $\\gamma\\gamma$ results on 4th July 2012:\n", "\n", "" ] }, { "cell_type": "markdown", "id": "00d845ba", "metadata": {}, "source": [ "### The impossible model\n", "\n", "
\n", "\n", "When starting the model-building process, it can also be helpful to think about the ideal model that you would fit if you had no restrictions on model complexity or computational resources. The purpose of this is to consider each step of the *data generating process* in detail, from start to finish. We are forced to think about complexities of the physical system and measurement process that may or may not be important for our later inference. The model itself may be impossible to actually write down or implement, but in doing this, we give can identify important steps and give context to the simpler models we will later go on to create. \n", "\n", "Taking the first example in the table above, consider what other effects may be important...\n", "\n", "
\n", "\n", "\n", "\n", "
\n", "\n", "**Exercise 2 (5 points):** Can you expand the complexity for one of the other models that we discussed? List possible factors." ] }, { "cell_type": "markdown", "id": "b9011f2c", "metadata": {}, "source": [ "### Probabilistic graphical models\n", "\n", "
\n", "\n", "For both of the above steps, it can be useful to sketch out connections between data and parameters using directed acyclic graphs (DAGs), often also called probabilistic graphical models (PGMs). These can be thought of as a visualisation of the data generating process, with parameters in unshaded circles, data or observations in shaded circles, and the arrows representing the connections between them. \n", "\n", "Here is a trivial example for that of data sampled from a normal distribution, such as what we looked at when introducing Stan in the *Are you ready?* section of the [introduction notebook](introduction.ipynb). The normal distribution is defined by two parameters, $\\mu$ and $\\sigma$. $N$ independent samples from this distribution are then taken, shown by $x_i$. The box around $x_i$ shows that there are a number of repeated samples, that we can think of as observations or measurements.\n", "\n", "
\n", "\n", "\n", "\n", "
\n", "\n", "**Exercise 3 (5 points):** Now it's your turn! Sketch graphical models for the first two examples in the above table (measurement of gravitational acceleration and radioactivity)." ] }, { "cell_type": "markdown", "id": "51d17676", "metadata": {}, "source": [ "### An example: the period-luminosity relation of Cepheid variable stars\n", "\n", "
\n", "\n", "Cepheid variables are a type of star that pulsates in observed flux, such that the apparent magnitude of the star varies in time with a well-defined period and amplitude. This behaviour is interesting in its own right, but also because it has been observed that there is a strong correlation between this period and the intrinsic luminosity of these stars. This **period-luminosity** relation has the form\n", "$$\n", "M = \\alpha + \\beta \\log_{10} \\bigg( \\frac{P}{\\mathrm{days}} \\bigg), \n", "$$\n", "where $M$ is the absolute magnitude, related to the apparent magnitude, $m$, by the luminosity distance, $d_L(z)$\n", "$$\n", "m = M + 5 \\log_{10} \\bigg( \\frac{d_L(z)}{\\mathrm{Mpc}} \\bigg) + 25.\n", "$$\n", "\n", "
\n", "\n", "\n", "\n", "
\n", "\n", "In astronomy, it is common to work with magnitudes instead of flux and luminosity, so the above expression is basically the equivalent of $F = L / 4 \\pi d_L^2(z)$. The magnitude is related to the logarithm of flux/luminosity and it is important to remember that magnitude *decreases* with increasing brightness, although we typically invert the y-axis when plotting magnitude for easier interpretation. So, **intrinsically brighter Cepheids have longer periods**.\n", "\n", "
\n", "\n", "We know this because we can measure the distance to nearby Cepheids using *parallax*. By determining the period-luminosity relation, we can then measure the distance to Cepheids that are further away, extending the [cosmic distance ladder](https://en.wikipedia.org/wiki/Cosmic_distance_ladder). Measuring distances to far away objects allows us to compare their luminosity distance and redshift, giving us an estimate of the [Hubble constant](https://en.wikipedia.org/wiki/Hubble%27s_law). \n", "\n", "
\n", "\n", "In the coming notebooks, we will explore a series of scientific questions related to Cepheid variable observations:\n", "\n", "1. What is the nature of the period-luminosity relation?\n", "\n", "2. Is the relationship universal? (i.e. the same in different galaxies)\n", "\n", "3. What is the value of the Hubble constant?" ] }, { "cell_type": "markdown", "id": "aabba55c", "metadata": {}, "source": [ "### Generative modelling\n", "\n", "
\n", "\n", "Using what we know about Cepheid variables and the period-luminosity relation, we can start to build a simple initial model for what we think could be the data-generating process. \n", "\n", "
\n", "\n", "Towards an initial model\n", "\n", "Measuring the period of a Cepheid is relatively straightforward from studying its lightcurve, and so to start with we can assume that the period is a known quantity. However, the observed characteristic apparent magnitude, $\\hat{m}$, of a Cepheid can be more tricky to reconstruct exactly, so let's say that it is reported with some fixed uncertainty, $\\sigma_m$, that is the same for all Cepheids. We can assume that \n", "$$\n", "\\hat{m} \\sim \\mathrm{Normal}(m, \\sigma_m),\n", "$$ \n", "where $m$ is the true apparent magnitude that is not observed. Knowing that our measured Cepheids are going to be pretty nearby ($z < 0.01$), we can approximate the luminosity distance as \n", "$$\n", "d_L(z) = \\frac{cz}{H_0},\n", "$$\n", "where $c$ is the speed of light, and $H_0$ is the Hubble constant. Again, for simplicity, let's assume that for now we just look at $N_c$ Cepheids in a single galaxy with a known redshift, $z$.\n", "\n", "
\n", "\n", "**Exercise 4 (5 points):** With these simplifying assumptions in mind, sketch the corresponding graphical model for our initial Cepheid variable model. The parameters and known quantities are summarised below, just arrange them into a graph.\n", "\n", "
\n", "\n", "\n", "\n", "
\n", "\n", "**Exercise 5 (10 points):** Using the parameter and fixed quantity values given below, as well as your graphical model, write a simulation for the observed apparent magnitudes, $m_i$. Plot your results in $\\log_{10}(P/\\mathrm{days})$-$m$ space, comparing with the true period-luminosity relation. \n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "47a38052", "metadata": { "ExecuteTime": { "end_time": "2021-09-11T11:36:20.799040Z", "start_time": "2021-09-11T11:36:20.543182Z" } }, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "id": "98edc7e0", "metadata": { "ExecuteTime": { "end_time": "2021-09-11T11:36:20.803888Z", "start_time": "2021-09-11T11:36:20.800912Z" } }, "outputs": [], "source": [ "# known quantities\n", "c = 3.0e5 # km/s\n", "H0 = 70 # km/s/Mpc\n", "z = 0.005\n", "\n", "Nc = 10\n", "P = [71.5, 45.0, 19.9, 30.6, 8.8, 75.2, 21.5, 22.0, 34.5, 87.4] # days\n", "sigma_m = 0.5\n", "\n", "# high-level parameter values\n", "alpha = -2\n", "beta = -4" ] }, { "cell_type": "markdown", "id": "2ced8589", "metadata": {}, "source": [ "Recall: \n", "$$\n", "M = \\alpha + \\beta \\log_{10} \\bigg( \\frac{P}{\\mathrm{days}} \\bigg), \n", "$$\n", "\n", "$$\n", "m = M + 5 \\log_{10} \\bigg( \\frac{d_L(z)}{\\mathrm{Mpc}} \\bigg) + 25,\n", "$$\n", "\n", "$$\n", "d_L(z) = \\frac{cz}{H_0}.\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "id": "0e0a6e34", "metadata": { "ExecuteTime": { "end_time": "2021-09-11T11:36:20.810637Z", "start_time": "2021-09-11T11:36:20.807982Z" } }, "outputs": [], "source": [ "# define some useful functions\n", "# dL: luminosity distance\n", "def luminosity_distance(z):\n", " # to be completed...\n", " pass\n", "\n", "# M: absolute magnitude\n", "def absolute_magnitude(alpha, beta, P):\n", " # to be completed...\n", " pass\n", "\n", "# m: apparent magnitude\n", "def apparent_magnitude(M, z):\n", " # to be completed...\n", " pass" ] }, { "cell_type": "code", "execution_count": 4, "id": "0668a9e4", "metadata": { "ExecuteTime": { "end_time": "2021-09-11T11:36:20.813982Z", "start_time": "2021-09-11T11:36:20.811762Z" } }, "outputs": [], "source": [ "# simulate the observed apparent magnitudes\n", "# M_true = ...\n", "# m_true = ...\n", "# m_obs = ...\n", "# M_obs = ...\n", "pass" ] }, { "cell_type": "code", "execution_count": 5, "id": "f6a040de", "metadata": { "ExecuteTime": { "end_time": "2021-09-11T11:36:20.963115Z", "start_time": "2021-09-11T11:36:20.815378Z" } }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "No handles with labels found to put in legend.\n" ] }, { "data": { "text/plain": [ "Text(0, 0.5, 'm')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAASGklEQVR4nO3dfZBdd13H8feHJiEjVHDSyGC3kCgpEEEprLUOKFWqpp0xHXmQdqwCdojCVB1hcOrAFCz/iB11Rq0P4amIQikozKrBqFAeRFq6tVBIa3EJD90idom1M4qhLXz9456Q63b3l7vbnL03m/drJjPn4XfP/d7fbPaz5/zu+Z1UFZIkLedh4y5AkjTZDApJUpNBIUlqMigkSU0GhSSpyaCQJDX1GhRJdiW5I8lcksuX2P/wJO/q9t+YZFuf9UiSVq63oEhyCnA1cD6wE7g4yc5FzS4F7qmqJwC/B7yhr3okSavT5xnF2cBcVR2sqvuAa4ELF7W5EHhbt/we4DlJ0mNNkqQV2tDjsU8H7hxanwd+cLk2VfVAknuBLcBXhxsl2QPsAXjEIx7xjCc96Ul91SxJ69LNN9/81arauprX9hkUx01V7QX2AkxPT9fs7OyYK5KkE0uSL672tX1eeroLOGNofarbtmSbJBuARwGHeqxJkrRCfQbFTcCOJNuTbAIuAmYWtZkBXtQtPx/4YDlLoSRNlN4uPXVjDpcB+4FTgLdU1YEkVwKzVTUDvBl4e5I54D8ZhIkkaYL0OkZRVfuAfYu2XTG0fBh4QZ81SNJ6c//99zM/P8/hw4cftG/z5s1MTU2xcePG4/Z+J8RgtiTpqPn5eU499VS2bdvG8B0FVcWhQ4eYn59n+/btx+39nMJDkk4whw8fZsuWLSy+7SwJW7ZsWfJM46EwKCTpBLTcvcl93LNsUEiSmgwKSVKTQSFJJ6Dlbjnr41Y0g0KSTjCbN2/m0KFDDwqFI9962rx583F9P78eK0knmKmpKebn51lYWHjQviP3URxPBoUknWA2btx4XO+TOBYvPUmSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVJTr0GRZFeSO5LMJbl8if2vSHJbkluTfCDJ4/usR5K0cr0FRZJTgKuB84GdwMVJdi5qdgswXVXfB7wH+O2+6pEkrU6fZxRnA3NVdbCq7gOuBS4cblBV11fV17rVG4CpHuuRJK1Cn0FxOnDn0Pp8t205lwLvX2pHkj1JZpPMLiwsHMcSJUnHMhGD2UkuAaaBq5baX1V7q2q6qqa3bt26tsVJ0kluQ4/Hvgs4Y2h9qtv2/yQ5D3g18Oyq+nqP9UiSVqHPM4qbgB1JtifZBFwEzAw3SHIW8KfA7qq6u8daJEmr1FtQVNUDwGXAfuB24LqqOpDkyiS7u2ZXAY8E3p3kk0lmljmcJGlM+rz0RFXtA/Yt2nbF0PJ5fb6/JOmhm4jBbEnS5DIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTb0GRZJdSe5IMpfk8ka75yWpJNN91iNJWrnegiLJKcDVwPnATuDiJDuXaHcq8KvAjX3VIklavT7PKM4G5qrqYFXdB1wLXLhEu9cDbwAO91iLJGmV+gyK04E7h9bnu23fkuTpwBlV9betAyXZk2Q2yezCwsLxr1SStKyxDWYneRjwu8Arj9W2qvZW1XRVTW/durX/4iRJ39JnUNwFnDG0PtVtO+JU4CnAh5J8ATgHmHFAW5ImS59BcROwI8n2JJuAi4CZIzur6t6qOq2qtlXVNuAGYHdVzfZYkyRphXoLiqp6ALgM2A/cDlxXVQeSXJlkd1/vK0k6vjb0efCq2gfsW7TtimXanttnLZKk1fHObElSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkppGmGe+eOvdq4PHdawJUVX1fj7VJkibAqM+j+AvgVcCngW/2V44kadKMGhQLVTVz7GaSpPVm1KB4bZI3AR8Avn5kY1X9VS9VSZImxqhB8RLgScBGjl56KsCgkKR1btSg+IGqemKvlUiSJtKoX4/95yQ7e61EkjSRRj2jOAf4ZJLPMxij8OuxknSSGDUodvVahSRpYo0UFFX1xb4LkSRNJqfwkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1NRrUCTZleSOJHNJLl+mzc8kuS3JgSTv6LMeSdLKjTp77IolOQW4GvhxYB64KclMVd021GYH8BvAM6vqniTf2Vc9kqTV6fOM4mxgrqoOVtV9wLXAhYvavBS4uqruAaiqu3usR5K0Cn0GxenAnUPr8922YWcCZyb5WJIbkiz53Iske5LMJpldWFjoqVxJ0lLGPZi9AdgBnAtcDLwxyaMXN6qqvVU1XVXTW7duXdsKJekk12dQ3AWcMbQ+1W0bNg/MVNX9VfV54LMMgkOSNCH6DIqbgB1JtifZBFwEzCxq8z4GZxMkOY3BpaiDPdYkSVqh3oKiqh4ALgP2A7cD11XVgSRXJtndNdsPHEpyG3A98KqqOtRXTZKklUtVjbuGFZmenq7Z2dlxlyFJJ5QkN1fV9GpeO+7BbEnShDMoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTb0GRZJdSe5IMpfk8iX2Py7J9UluSXJrkgv6rEeStHK9BUWSU4CrgfOBncDFSXYuavYa4LqqOgu4CPijvuqRJK1On2cUZwNzVXWwqu4DrgUuXNSmgG/vlh8FfLnHeiRJq9BnUJwO3Dm0Pt9tG/Y64JIk88A+4JeXOlCSPUlmk8wuLCz0UaskaRnjHsy+GLimqqaAC4C3J3lQTVW1t6qmq2p669ata16kJJ3M+gyKu4Azhtanum3DLgWuA6iqjwObgdN6rEmStEJ9BsVNwI4k25NsYjBYPbOozZeA5wAkeTKDoPDakiRNkN6CoqoeAC4D9gO3M/h204EkVybZ3TV7JfDSJJ8C3gm8uKqqr5okSSu3oc+DV9U+BoPUw9uuGFq+DXhmnzVIkh6acQ9mS5ImnEEhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU0GhSSpyaCQJDUZFJKkJoNCktRkUEiSmgwKSVKTQSFJajIoJElNBoUkqcmgkCQ1GRSSpCaDQpLUZFBIkpoMCklSk0EhSWoyKCRJTQaFJKnJoJAkNRkUkqQmg0KS1GRQSJKaDApJUpNBIUlqMigkSU29BUWStyS5O8lnltmfJL+fZC7JrUme3lctkqTV6/OM4hpgV2P/+cCO7t8e4I97rEWStEq9BUVVfQT4z0aTC4E/q4EbgEcneWxf9UiSVmfDGN/7dODOofX5btu/L26YZA+Dsw6Ary93OeskdBrw1XEXMSHsi6Psi6Psi6OeuNoXjjMoRlZVe4G9AElmq2p6zCVNBPviKPviKPviKPviqCSzq33tOL/1dBdwxtD6VLdNkjRBxhkUM8DPd99+Oge4t6oedNlJkjRevV16SvJO4FzgtCTzwGuBjQBV9SfAPuACYA74GvCSEQ+997gXe+KyL46yL46yL46yL45adV+kqo5nIZKkdcY7syVJTQaFJKlpYoMiya4kd3RTfFy+xP6HJ3lXt//GJNvGUOaaGKEvXpHktm4qlA8kefw46lwLx+qLoXbPS1JJ1u1XI0fpiyQ/0/1sHEjyjrWuca2M8H/kcUmuT3JL9//kgnHU2bfepk6qqon7B5wCfA74bmAT8Clg56I2Lwf+pFu+CHjXuOseY1/8KPBt3fLLTua+6NqdCnwEuAGYHnfdY/y52AHcAnxHt/6d4657jH2xF3hZt7wT+MK46+6pL34EeDrwmWX2XwC8HwhwDnDjKMed1DOKs4G5qjpYVfcB1zKY8mPYhcDbuuX3AM9JkjWsca0csy+q6vqq+lq3egODe1LWo1F+LgBeD7wBOLyWxa2xUfripcDVVXUPQFXdvcY1rpVR+qKAb++WHwV8eQ3rWzPV09RJkxoUy03vsWSbqnoAuBfYsibVra1R+mLYpQz+YliPjtkX3an0GVX1t2tZ2BiM8nNxJnBmko8luSFJa5LOE9koffE64JLuq/r7gF9em9Imzkp/nwAnyBQeGk2SS4Bp4NnjrmUckjwM+F3gxWMuZVJsYHD56VwGZ5kfSfLUqvqvcRY1JhcD11TV7yT5IeDtSZ5SVd8cd2Engkk9oxhleo9vtUmygcHp5KE1qW5tjTTVSZLzgFcDu6vq62tU21o7Vl+cCjwF+FCSLzC4BjuzTge0R/m5mAdmqur+qvo88FkGwbHejNIXlwLXAVTVx4HNDCYMPNmsauqkSQ2Km4AdSbYn2cRgsHpmUZsZ4EXd8vOBD1Y3WrPOHLMvkpwF/CmDkFiv16HhGH1RVfdW1WlVta2qtjEYr9ldVaueDG2CjfJ/5H0MziZIchqDS1EH17DGtTJKX3wJeA5AkiczCIqFNa1yMqxq6qSJvPRUVQ8kuQzYz+AbDW+pqgNJrgRmq2oGeDOD08c5BoM3F42v4v6M2BdXAY8E3t2N53+pqnaPreiejNgXJ4UR+2I/8BNJbgO+AbyqqtbdWfeIffFK4I1Jfo3BwPaL1+Mfln1NneQUHpKkpkm99CRJmhAGhSSpyaCQJDUZFJKkJoNCktRkUGhdSPLfD+G1l3WzaVZ3v8GR7cvOtJnksUn+pls+N8m9ST6Z5PYkr110/JuTPHzRthcn+cPV1rzEZ3hqkmuO1/GkYQaFBB8DzgO+uGj7+QzuZN4B7AH+eGjfK4A3Dq1/tKqexmAKlUuOhEqS7cBdfd8tX1WfBqaSPK7P99HJyaDQutKdBVyV5DNJPp3khd32hyX5oyT/muQfkuxL8nyAqrqlqr6wxOFaM20+D/i7xS+oqv8Bbgae0G3adaRdkpck+WySTwDPHKr5pzJ4psotSf4xyWO6ev8tydah+ueSbE3ygu7zfSrJR4be/q9ZpzeearwMCq03zwWeBnw/g7OEq7pf7s8FtjF4FsHPAT80wrGWnGmzO0u4Z6mzhCRbGMwxdaDbtAv4u66G32QQEM/q6jjin4BzquosBlNk/3o3Wd2fAz/btTkP+FRVLQBXAD9ZVd8PDN+BPwv88AifS1oRg0LrzbOAd1bVN6rqP4APAz/QbX93VX2zqr4CXP8Q3uOxPHieoB9Ocgvw98BvdVNIbAKmquog8IPAh6pqoXtmwruGXjsF7E/yaeBVwPd2298C/Hy3/AvAW7vljwHXJHkpgykrjrgb+K6H8LmkJRkU0vKWm2nzfxlMKjfso1V1VlU9o5tTBwZ/3f/TCO/zB8AfVtVTgV88cuyquhP4jyQ/xuDhPO/vtv8S8Jqutpu7sxi61/3vyj6idGwGhdabjwIvTHJKd33/R4BPMPgr/Hndtf7H0M2qegzLzbT5WQaXsY5lF0cfInUj8OwkW5JsBF4w1O5RHJ3q+UX8f29icAnq3VX1DYAk31NVN1bVFQzObI6E2ZnAks9Klh4Kg0LrzXuBWxk8N/mDDK73fwX4SwZjDLcx+MX7LwyeikiSX+lm2pwCbk3ypu5Y+xhMyz3H4BtOL4dvDVh/LsmRAevlnMvg0hddwLwO+DiD0Lp9qN3rGMz8ezPw1UXHmGEwM/Bbh7Zd1Q3Ufwb45+6zwuDZ6ev9yX4aA2eP1UkjySOr6r+7SzWfAJ7ZhchqjvXTwDOq6jXL7J8C3lhV56++YsjgoUu/V1XNQeruPo0PA8/qHg0sHTcT+TwKqSd/k+TRwCbg9asNCYCqeu/Q2MBS++cZ3IexakkuB17G0W8+tTwOuNyQUB88o5AkNTlGIUlqMigkSU0GhSSpyaCQJDUZFJKkpv8DUmoVL6Uas5wAAAAASUVORK5CYII=\n", "text/plain": [ "