Skip to content

Instantly share code, notes, and snippets.

@iamirmasoud
Created July 31, 2022 20:26
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save iamirmasoud/2b83de47ee93b9a27a7694495182164d to your computer and use it in GitHub Desktop.
Save iamirmasoud/2b83de47ee93b9a27a7694495182164d to your computer and use it in GitHub Desktop.
Importance Sampling in Python
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Importance Sampling\n",
"---\n",
"Using sampling to approximate a distribution\n",
"\n",
"$$E[f(x)] = \\int f(x)p(x) dx \\approx \\frac{1}{n}\\sum_{i} f(x_i)$$\n",
"where $ x \\sim p(x)$\n",
"\n",
"$$E[f(x)] = \\int f(x)p(x) dx = \\int f(x)\\frac{p(x)}{q(x)}q(x) dx \\approx \\frac{1}{n} \\sum_{i} f(x_i)\\frac{p(x_i)}{q(x_i)}$$\n",
"\n",
"where $ x \\sim q(x)$\n",
"\n",
"Idea of importance sampling: draw the sample from a proposal distribution and re-weight the integral using importance weights so that the correct distribution is targeted\n",
"\n",
"$$Var(X) = E[X^2] - E[X]^2$$\n",
"\n",
"**Reference**\n",
"\n",
"- [1](https://www.youtube.com/watch?v=3Mw6ivkDVZc)\n",
"- [2](https://astrostatistics.psu.edu/su14/lectures/cisewski_is.pdf)"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jeremy.zhang/anaconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
" return f(*args, **kwds)\n",
"/Users/jeremy.zhang/anaconda3/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88\n",
" return f(*args, **kwds)\n"
]
}
],
"source": [
"import numpy as np\n",
"import scipy.stats as stats\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns"
]
},
{
"cell_type": "code",
"execution_count": 161,
"metadata": {},
"outputs": [],
"source": [
"def f_x(x):\n",
" return 1/(1 + np.exp(-x))\n",
"\n",
"def distribution(mu=0, sigma=1):\n",
" # return probability given a value\n",
" distribution = stats.norm(mu, sigma)\n",
" return distribution"
]
},
{
"cell_type": "code",
"execution_count": 121,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x115b111d0>"
]
},
"execution_count": 121,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAESCAYAAAABl4lHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VOX9/vH3h5CQBEgCJCwh7IRFFhHCJlqXulC02ir2ay1aFUQtVr+1topff9paW6ut3dTWilJUVrUu4FKXinUHEgz7FtaEAFlJCEnI9vz+SKQxBJlgZs4kc7+uK5eZmSeZ20Nm7jnnOYs55xAREfFFG68DiIhIy6HSEBERn6k0RETEZyoNERHxmUpDRER8ptIQERGfqTRERMRnKg0REfGZZ6VhZnPNLMfM1h/ncTOzv5hZhpmtNbPRgc4oIiJf1tbD554HPAY8e5zHvwUk132NB/5W99+vFB8f7/r27ds8CUVEQkRaWlqecy7hROM8Kw3n3Adm1vcrhlwKPOtqz3PymZnFmVkP59y+r/q9ffv2JTU1tRmTioi0fma225dxwTyn0RPIrHc7q+4+ERHxSDCXhjVyX6NnVzSzmWaWamapubm5fo4lIhK6grk0soBe9W4nAdmNDXTOPemcS3HOpSQknHCTnIiInKRgLo2lwDV1e1FNAIpONJ8hIiL+5dlEuJktAs4G4s0sC7gPCAdwzj0BvAFMATKAUuA6b5KKiMgXvNx76vsneNwBs5r7eYuLi8nJyaGysrK5f3Wr0b59e5KSkmjTJphXREXEC14epxFwxcXFHDhwgJ49exIVFYVZY3Ptoa2mpoa9e/eSl5dH165dvY4jIkEmpEojJyeHnj17Eh0d7XWUoNWmTRu6devG7t27VRoiQaqquob8wxUcKC4np/gIOYeOcKC4nG8O7crIpDi/PndIlUZlZSVRUVFexwh64eHhVFVVeR1DJCSVV1aTfbCMfUXl7Csq50BxOfvrf19cTn7JEWoaOQAhoWM7lUZz0yapE9MyEvEP5xwHSyvJKiwjq7CUrMIy9h4sI/tgGdlFZWQfLKfgcMUxPxcbFU6P2Ei6xURySo8YusW0IyEmkm4d29E1JpKuHdsR36EdEW39Pw8ZcqUhIuJPR6qqySwoI7OglN35h9ldUEpmQSmZBbVFcbii+kvjO7RrS2JcJIlxUYxMiiMxtvb7HrFRR4siKiLMo/+bY6k0RESaqKKqhj0FpezKO8yu/MPszKv92pV3mH3F5bh6m46iI8Lo3Tma3l2iOX1gF5I6RZPUKaruK5rYqHDv/kdOgkpDROQ4isoq2Z5bQkZOCdtzS9ieU8L23MPsKSilut6kQlx0OH27tGd8/y706RJN787Rdf9tT3yHiFa1yVel0cLMmTOHBx98kD179jB9+nReeeUVPvnkEwYMGODTz0+dOpXTTz+d22+/3c9JRVqOkiNVbDtwiK0HDrFlfwnbcg6xZf8hcg4dOTomIqwNfeOjGdqjIxeN6EH/hPb0jW9Pvy7t6dQ+wsP0gaXSaEE2b97MzTffzAsvvMCECRN4+OGHmTJlis+FAXDfffdx1llnMX36dGJjY/2YViT41NQ4MgtL2bSvmI37DrFpXzGb9hWTVVh2dExkeBuSu3bkzOQEkrt1YGBCBwZ07UCvTlG0DdMBryqNFmTp0qUMHz6c7373u5SWljJv3jyWLVvWpN8xYsQI+vfvz/z585k1q9kPuBcJGlXVNWzPPcy6vUWsr/vatK/46ER0G4N+8e0Z1SuOK8f2YnD3GAZ160BSp2jC2rSezUnNTaXRQgwaNIht27YB/90ltnPnzkyaNOlL41544QWmTZvG1q1b6dOnDwC33XYbr732Gp988gndunXjkksuYdGiRSoNaTVqahw78w+zJvMgazIPsm5vERv3FVNeWQNAVHgYwxJjmDomiaE9YhjaI4ZB3ToG1V5JLYVKo4X46KOPOPPMM7nmmmuYPn06s2fPZu/evcdMsE2dOpWHHnqIBx54gDlz5vD73/+eRYsW8fHHH9OtWzcAxo0bxwMPPEBZWZkOdpQWqeBwBZ/vKSQ98yDpdUVRXF57QGr7iDCGJcZy1bg+jEiKYUTPWPrFd9DaQzMJ+dL45bINbMwuDuhznpIYw33fHtakn4mJiWHHjh1MmjSJ7t27U1hYSI8ePY4ZZ2b85je/4aKLLmLAgAH8+te/5r333iM5OfnomMTERCorK8nOzm7SfIiIF5xz7Mg7TNquQlJ3F5C6u5AduYeB2k1Mg7vHcNHIREb1imVUr04M7KqC8KeQL42WYv369VRVVTFq1CgAysrKjq45NHTBBRcwduxY7rnnHpYtW8bYsWO/9PgXaxdlZWWN/biIp2pqHJv3H+KzHfms2JnPql2FR4+SjosOZ0zvTlw+OokxfToxMimW6Ai9jQVSyC/tpn7i90p6ejp9+vQhLq72vDLx8fEUFhY2Ova9995jzZo1OOcaLZaCggIAdJVDCQY1NY5N+4v5bEcBn+3IZ+XOAorKai9d0KtzFOcM7kpK306k9OnEgIQOtNFahKdCvjRaivT09KNrGQCnnXYa8+bNO2bcmjVruOyyy3j00Ud5/fXXmT17Nm+99daXxqxfv57ExMTjrqmI+FtWYSkfZ+TxUUY+n2TkkV+3JtGnSzQXDuvGhP5dGN+/Cz3jNOcWbFQaLUR6ejrnnXfe0dsXXnghd955J/n5+XTp0gWA3bt3M2XKFG6//Xauv/56xo0bx8iRI3n//fc5++yzj/7shx9+yOTJkwP9vyAhrLSiik8y8vnP1lw+yshjZ17tnETXju04a1ACkwbGM3FAFxJVEkFPpdECOOdYu3Ytd9xxx9H7RowYwbhx41i8eDGzZs2ioKCAyZMnc/HFF3PvvfcCMHz4cK644gpmz57Np59+CkB5eTkvv/zyMWsfIs3JOcfOvMMs35LL+1tyWLGjgIrqGqIjwpjYvwtXT+jDGcnxJHft0KpOsREKVBotgJlRXHzsHl733Xcft912GzfddBOdO3dm06ZNx4xZsmTJl24//fTTjB8/ngkTJvgtr4SmquoaVu4q4N2NOfx78wF255cCMLBrB66Z2IdzhtTOTbRrq2MjWjKVRgs2efJkZs2aRVZW1tED+U4kPDycRx991M/JJFSUHKnig625vLPxAO9tzqGorJKItm2YNKALM87ox9mDu9Krs66U2ZqoNFq4W2+9tUnjZ86c6ackEiqKyip5d+MB3li3jw+35VFRXUNcdDjfHNqVC07pxpnJCbRvp7eW1kr/siJyQkVllbxztChyqax2JMZGMm1CHy4Y1o2UPp10Mr8QodIQkUaVVVTzzqYDvPr5Xj6oK4qecVH8cGJfpozswaikOB0zEYJUGiJyVFV1DR9l5PFqejZvbdhPaUU1PWIjufb0vlw0MpFTk2K1t1OIC7nScM7pj/4EXP1rVUpI2JBdxItpWSxbk01eSQUxkW255NRELh3Vk/H9OmuNQo4KqdIIDw+nrKyM6GjtzfFVKisrads2pP40QtLB0gpeTc/m+dRMNmQXExHWhnOHdOU7p/XknCEJ2jVWGhVS7wxdu3Zl79699OzZk6ioKK1xNKKmpoYDBw7oqn6tVE2N48OMPJ5PzeSdDQeoqK5hWGIMv7xkGJecmhhSly2VkxNSpRETEwNAdnY2lZWVHqcJXu3btyc+Pt7rGNKM8kqO8EJqFgtX7iazoIxO0eFcNb43V6QkMSxRHxDEdyFVGlBbHF+Uh0hr5pxj5c4CFqzYw5vr91FZ7RjfrzM/u3AIFw7rps1PclJCrjREWrvSiipeWr2XZz/dxdYDJXSMbMsPxvdh2oTeDOza0et40sKpNERaieyDZTz76W4WrdxDUVklw3vG8PDlI/n2qYm6FrY0G5WGSAu3ek8hcz/ayZvr9+Oc48Jh3bn+jH6k9OmknT2k2XlaGmY2GfgzEAY85Zz7bYPH+wBzgQSgAJjmnMsKeFCRIFNT4/j35hye+M920nYX0jGyLdPP6MfVE/roBIHiV56VhpmFAY8D5wNZwCozW+qc21hv2O+BZ51zz5jZucCDwNWBTysSHCqqali6Jpu//2c723JKSOoUxS++fQpXpPTSSQIlILz8KxsHZDjndgCY2WLgUqB+aZwC/KTu++XAKwFNKBIkDh+pYvGqTJ7+cAfZReUM6d6RP185iotG9NCJAiWgvCyNnkBmvdtZwPgGY9YAl1O7Ceu7QEcz6+Kcyw9MRBFvHT5SxbOf7ubJD7ZTWFrJuH6d+fVlIzh7UILmK8QTXpZGY3/xDU96dAfwmJldC3wA7AWqjvlFZjOBmQC9e/du3pQiHig5UsWzn+5izgc7KCyt5OzBCfz43GTG9OnkdTQJcV6WRhbQq97tJCC7/gDnXDZwGYCZdQAud84VNfxFzrkngScBUlJSdLY9abFKjlTxzCe7mPPhDg6WVnLO4ARuO28Qo3rFeR1NBPC2NFYByWbWj9o1iCuBq+oPMLN4oMA5VwPMpnZPKpFWp7yymvmf7ebx5RkUqiwkiHlWGs65KjO7BXiL2l1u5zrnNpjZ/UCqc24pcDbwoJk5ajdPzfIqr4g/VFXX8NLqvfzp3a1kF5VzZnI8P71gsMpCgpa1tmsnpKSkuNTUVK9jiHwl5xxvbdjP797awvbcw5zaK447LxzM6QN1okjxhpmlOedSTjROO3aLBNiqXQU88Pom1mQeZEBCe56YNoYLh3XT3lDSIqg0RAIks6CUB9/cxBvr9tM9JpKHLx/JZaN76jgLaVFUGiJ+dqi8kseWZ/CPj3YR1sb4yXmDuOEb/YiO0MtPWh791Yr4SXWNY/GqPfzh7a3kH67g8tFJ/OzCwXSPjfQ6mshJU2mI+EHa7kLufXU9G7KLGde3M/OuO4URSbpCnrR8Kg2RZpRfcoSH/rWZ51Oz6B4TyWNXncZFI3pokltaDZWGSDOornEsXLmH3/1rM6UV1dx4Vn9uPTdZZ56VVkd/0SJfU3rmQe55ZR3r9xYzsX8XfvWdYbqsqrRaKg2Rk3T4SBWPvL2Vf3yyk4QO7fjL90/j2yO1KUpaN5WGyElYviWHe15ez96DZUyb0JufTx5CTGS417FE/E6lIdIE+SVH+NVrG3klPZsBCe154aaJjO3b2etYIgGj0hDxgXOOV9L3cv+yjZQcqeK2bybzo3MG0K5tmNfRRAJKpSFyArmHjnD3y+t4Z+MBTusdx0OXj2RQN010S2hSaYh8hdfX7uOeV9ZxuKKau6cMYfoZ/Qlro4luCV0qDZFGFB6u4N6lG1i2JpuRSbE8csWpJGvtQkSlIdLQuxsPMPvldRwsreCn5w/i5rMH6Ey0InVUGiJ1yiqqeeD1jSxYsYch3Tsy77qxDEvU+aJE6lNpiAAbsou4bXE6GTklzPxGf356wSDtGSXSCJWGhLSaGsfcj3fy8L+2EBcdzvzp4zkjWZdcFTkelYaErJxD5dzxwlo+2JrLeUO78fDUkXRuH+F1LJGgptKQkPTB1lx+siSdkiNVPPCd4fxgfG+dM0rEByoNCSnVNY4/v7uVR5dnMKhrRxbPnKBdaUWaQKUhISPnUDm3LUrn0x35XDEmifsvHU5UhCa7RZpCpSEh4dPt+dy6+HMOlVfy8NSRfC+ll9eRRFoklYa0ajU1jr/9ZzuPvL2FvvHteW76OIZ0j/E6lkiLpdKQVquorJLbl6Tz7805fPvURB68bAQddPlVka9FryBplbYdOMTM59LILCjll5cM45qJfbR3lEgzUGlIq/PWhv3cviSdqIgwFt4wgXH9dJEkkeai0pBWo6bG8cd3t/LoexmcmhTLE1ePoUdslNexRFoVlYa0CkVllfxkSTrvbc7hijFJ/Oo7w4kM1+60Is1NpSEt3s68w0x/ZhV78kv51XeGM01Hd4v4jUpDWrRPt+dz0/w02hgsmDGe8f27eB1JpFXz9MoyZjbZzLaYWYaZ3dXI473NbLmZfW5ma81sihc5JTgtWbWHq59eQULHdrwya5IKQyQAPFvTMLMw4HHgfCALWGVmS51zG+sNuwd43jn3NzM7BXgD6BvwsBJUqmscv31zE3M+3MmZyfE8/oPRxESGex1LJCR4uXlqHJDhnNsBYGaLgUuB+qXhgC8O340FsgOaUIJOyZEqblv0Of/enMMPJ/bh/118ii7FKhJAXpZGTyCz3u0sYHyDMb8A3jazHwPtgfMCE02C0f6icq79x0q25ZRw/6XDuGZiX68jiYQcLz+iNbZ7i2tw+/vAPOdcEjAFeM7MjslsZjPNLNXMUnNzc/0QVby2Zf8hvvvXj8ksKGXutWNVGCIe8bI0soD6pxpN4tjNT9OB5wGcc58CkcAx1+J0zj3pnEtxzqUkJCT4Ka545dPt+Ux94hOqaxzP3zSRswbp31jEK16Wxiog2cz6mVkEcCWwtMGYPcA3AcxsKLWloVWJELJ0TTY/nLuSbjGRvPSj0xmWGOt1JJGQ5tmchnOuysxuAd4CwoC5zrkNZnY/kOqcWwr8FJhjZj+hdtPVtc65hpuwpBVyzjHnwx385o3NjOvXmTlXpxAbrT2kRLzm6cF9zrk3qN2Ntv5999b7fiMwKdC5xFs1NY77X9vIvE92cdHIHjxyxak6JYhIkNAR4RJUKqpq+OkLa1i2JpsZZ/Tj7ilDadNGpwQRCRYqDQkapRVV3Dx/Nf/Zmstd3xrCTWcN8DqSiDSg0pCgUFRayfXPrOLzPYX89rIRXDmut9eRRKQRKg3xXE5xOVc/vZKdeYd5/KrRfGtED68jichxqDTEU7vzDzPt6RXkl1Qw99qxnJF8zGE4IhJEVBrimS37DzHt6RVUVdew8IYJjOoV53UkETkBlYZ4Yv3eIq5+egURbdvw/I0TSe7W0etIIuIDlYYE3Od7Crlm7kpiIsNZeMN4+nRp73UkEfGRziktAbVyZwHTnlpB5/YRLLlxggpDpIXRmoYEzMcZecx4JpUecZEsnDGB7rGRXkcSkSbSmoYExPItOVw3bxW9O0ezZOZEFYZIC6U1DfG7dzYe4EcL0hjUrSPPTR9P5/YRXkcSkZOk0hC/ereuME5JjOXZ68cRG6Uz1Yq0ZNo8JX7z700HuHlBGkN7xKgwRFoJlYb4xfLNOdw8fzVDusfw3PXjVRgirYRKQ5rd+1tyuPG5NAZ178D86eN18SSRVqRJpWFmW83sTjPr7q9A0rJ9sDWXmc+lMbCrCkOkNWrqmkYl8CCwx8xeMbOLzUxrKwLAR9vyuOHZVAYkdGDBjPHERWsvKZHWpklv+M65YcDpwDPAOcCrQKaZ/drMdMWcELZyZwEznl1Fv/j2LJgxnk7arVakVWryWoJz7jPn3A1AD2AGsBOYDWw1s/fM7Coza9fMOSWIrc06yPXzVpEYF6XjMERauZPetOScK3XO/cM5dwYwBFgMnA08B2Sb2R/NTJdfa+U27y/mmrkriYsOZ8GM8SR01OcFkdbsa81HmFmYmX0X+APwP4ADlgOfAT8GNpnZpV87pQSlHbklTHtqJe3atmHhjAn0iI3yOpKI+NlJlYaZDTGz3wF7gX8CKcDvgUHOufOccxdRu/axBXi4ucJK8MgqLGXaUytwzrFgxgR6d4n2OpKIBECTTiNiZtcD04EJdXe9CzwJvOqcq6o/1jmXYWZ/AZ5qjqASPHKKy/nBUysoOVLFopkTGNi1g9eRRCRAmnruqaeA/cBvgTnOuV0nGL+R2jkOaSUOllYw7ekV5B06wnMzxjMsMdbrSCISQE0tjcuBpc65al8GO+dWAiubnEqCUmlFFdfNW8Wu/FLmXTeW0b07eR1JRAKsSaXhnHvZX0EkuFVU1XDT/NWsyTzI36aN4fQB8V5HEhEP6NTockI1NY47XljDB1tzeejyEVw4TGeREQlVOgWIfCXnHL9ctoGla7K5c/IQ/mesDr0RCWUqDflKf/73Np75dDc3nNmPm87q73UcEfGYSkOO67lPd/Gnd7dx+egk7p4yFDPzOpKIeEylIY16Y90+7l26gfOGduWhy0eoMEQE8Lg0zGyymW0xswwzu6uRx/9oZul1X1vN7KAXOUPNyp0F/O+SdEb37sRjV42mbZg+W4hILc/2njKzMOBx4HwgC1hlZkudcxu/GOOc+0m98T8GTgt40BCz7cAhZjyziqROUTx1TQqR4WFeRxKRIOLlR8hxQIZzbodzroLas+R+1ckNvw8sCkiyEHWguJxr/7GKiLZhPHPdOF0TQ0SO4WVp9AQy693OqrvvGGbWB+gHvBeAXCHpUHklP5y7koOlFcy7biy9OusEhCJyLC9Lo7GZVXecsVcCLx7v9CVmNtPMUs0sNTc3t9kChorao73TyMgp4a/TxjC8p84nJSKN87I0soBe9W4nAdnHGXslX7Fpyjn3pHMuxTmXkpCQ0IwRWz/nHHf+cy0fZ+Tz28tHctYgLT8ROT4vS2MVkGxm/cwsgtpiWNpwkJkNBjoBnwY4X0h45O2tvPz5Xu64YBBTxyR5HUdEgpxnpVF3/Y1bgLeATcDzzrkNZna/mV1Sb+j3gcXOueNtupKT9HxqJo8tz+DKsb2Ydc5Ar+OISAvg6QkLnXNvAG80uO/eBrd/EchMoeKTjDzufmkdZwyM51ffGa6D90TEJzpqKwRl5Bzixvlp9Itvz1+njSZcB++JiI/0bhFicg8d4dp/rKJd2zDmXjuWmMhwryOJSAui0ggh5ZXV3PBsKnklR3jqhyk6FkNEmkwXYQoRNTWO259PZ03WQf72gzGM6hXndSQRaYG0phEifvf2Ft5Yt5+7vzWUycN15T0ROTkqjRDwz7Qs/vb+dr4/rjczzuzndRwRacFUGq1c2u4CZr+0jtMHdOH+S4dp11oR+VpUGq1YVmEpM59NIzEukr/+QLvWisjXp3eRVqrkSBUznkmlorqGp344lrhoneZcRL4+7T3VClXXOG5b9DnbckqYd91YBnbt4HUkEWkltKbRCj38r838e3MO9337FM5M1llrRaT5qDRamRdSM/n7Bzu4ekIfrpnY1+s4ItLKqDRakbTdhfzfy+uZNLAL9377FK/jiEgrpNJoJfYVlXHjc2n0iIvk8au0p5SI+IcmwluB8spqbnwujbKKKhbeMF57SomI36g0WjjnHLNfWsfarCKevHoMg7p19DqSiLRi2obRwj314U5e/nwvt58/iAuG6ZxSIuJfKo0W7D9bc3nwzU1MGdGdH5+ry7WKiP+pNFqonXmH+fHC1Qzq1pHfTT1V55QSkYBQabRAh8orueHZVMLaGHOuSaF9O01NiUhg6N2mhampcfz0+TXszDvMc9PH6ep7IhJQWtNoYf76fgZvbzzA3VOGcvqAeK/jiEiIUWm0IMs35/DIO1v5zqhErp/U1+s4IhKCVBotxK68w9y6+HOGdo/hwctGauJbRDyh0mgBDh+p4sbn0ghrY/z96jFERYR5HUlEQpRKI8g55/j5i2vZlnOIR79/mia+RcRTKo0g9/cPdvD6un38fPIQXRtDRDyn0ghiH23L4+F/beaiET248Rv9vY4jIqLSCFZ7D5Zx6+LPGZDQgYenauJbRIKDSiMIHamq5kfz06ioquGJq8foiG8RCRp6NwpCv1y2kTVZRTwxbQwDEjp4HUdE5CitaQSZF1IzWbhiDzee1Z/Jw3WqcxEJLp6WhplNNrMtZpZhZncdZ8z3zGyjmW0ws4WBzhhI6/cWcc8r65nYvws/u2Cw13FERI7h2eYpMwsDHgfOB7KAVWa21Dm3sd6YZGA2MMk5V2hmXb1J639FpZXcvCCNTtERPHrVabTVNb5FJAh5+c40Dshwzu1wzlUAi4FLG4y5AXjcOVcI4JzLCXDGgKipcfzvks/ZX1TOX6eNJr5DO68jiYg0ysvS6Alk1rudVXdffYOAQWb2sZl9ZmaTA5YugB5bnsHyLbnce/EpjO7dyes4IiLH5eXeU40deOAa3G4LJANnA0nAh2Y23Dl38Eu/yGwmMBOgd+/ezZ/Ujz7clssf3609c+20CX28jiMi8pW8XNPIAnrVu50EZDcy5lXnXKVzbiewhdoS+RLn3JPOuRTnXEpCQss51ca+ojJuW5xOctcO/OayETqAT0SCnpelsQpINrN+ZhYBXAksbTDmFeAcADOLp3Zz1Y6ApvSTiqoaZi1YzZHKav42bQzRETpkRkSCn2el4ZyrAm4B3gI2Ac875zaY2f1mdkndsLeAfDPbCCwHfuacy/cmcfP67ZubWb3nIA9NHakD+ESkxfD0461z7g3gjQb33VvvewfcXvfVary+dh9zP97Jtaf35eKRiV7HERHxmQ4GCLDtuSX8/MU1nNY7jrunDPU6johIk6g0Aqi0ooofzV9Nu/AwHr9qNBFttfhFpGXR7GuAOOe455X1bM05xLPXjyMxLsrrSCIiTaaPugHyfGomL63ey63nJusKfCLSYqk0AmBjdjH3vrqBMwbGc+s3jznMRESkxVBp+Nmh8kpmLVxNbFQ4f7pyFGFtdACfiLRcmtPwI+ccd720jj0FpSycMV4nIhSRFk9rGn40/7PdvL52H3dcMJjx/bt4HUdE5GtTafjJ2qyD/Oq1TZwzOIEbv9Hf6zgiIs1CpeEHRaWV/GjBauI7RPCH742ijeYxRKSV0JxGM3POcceLa9hfVM7zN02kU/sIryOJiDQbrWk0s7kf7+KdjQe461tDdEElEWl1VBrNKD3zIL99cxPnn9KN6Wf08zqOiEizU2k0k6LSSmYtWE3XjpH8fuqpuqCSiLRKmtNoBl/MY+QcKueFm04nNjrc60giIn6hNY1m8N95jKGM6hXndRwREb9RaXxN9ecxrp/U1+s4IiJ+pdL4GjSPISKhRnMaJ8k5x880jyEiIUZrGifpHx/v4u2NB7hz8hDNY4hIyFBpnIQ1mQd58M1NnDdUx2OISGhRaTRRUVkltyyqm8e4YqTmMUQkpGhOowmcc8x+aS37Dpaz5MaJxEXrvFIiElq0ptEE8z/bzRvr9vOzCwczpo/OKyUioUel4aP1e4uOXh/jhjN1fQwRCU0qDR8cKq/kloWr6dw+gkd0fQwRCWGa0zgB5xx3v7yezMIyFt0wgc66PoaIhDCtaZzA4lWZLFuTze3nD2KdPS/LAAAHxUlEQVRcv85exxER8ZRK4yts3l/ML5Zu4MzkeG4+a4DXcUREPKfSOI7SiipmLVhNTFS4rvMtIlJHcxrHce+rG9iRd5gF08eT0LGd13FERIKC1jQa8c+0LF5My+LH5yZz+sB4r+OIiAQNT0vDzCab2RYzyzCzuxp5/FozyzWz9LqvGf7OlJFTwv97dT3j+3Xmtm8m+/vpRERaFM82T5lZGPA4cD6QBawys6XOuY0Nhi5xzt0SiEzlldXcsnA1keFh/PnK0wjTPIaIyJd4uaYxDshwzu1wzlUAi4FLPczD/a9tZPP+QzzyvVPpHhvpZRQRkaDkZWn0BDLr3c6qu6+hy81srZm9aGa9/BVm2ZpsFq7Yw41n9eecwV399TQiIi2al6XR2LYf1+D2MqCvc24k8C7wTKO/yGymmaWaWWpubu5JhekUHcH5p3TjjgsGn9TPi4iEAi9LIwuov+aQBGTXH+Ccy3fOHam7OQcY09gvcs496ZxLcc6lJCQknFSYM5LjmXNNCuFh2qFMROR4vHyHXAUkm1k/M4sArgSW1h9gZj3q3bwE2BTAfCIi0oBne08556rM7BbgLSAMmOuc22Bm9wOpzrmlwK1mdglQBRQA13qVV0REwJxrOI3QsqWkpLjU1FSvY4iItChmluacSznROG3AFxERn6k0RETEZyoNERHxmUpDRER8ptIQERGftbq9p8wsF9h9kj8eD+Q1Y5zmolxNo1xNF6zZlKtpvk6uPs65Ex4d3epK4+sws1RfdjkLNOVqGuVqumDNplxNE4hc2jwlIiI+U2mIiIjPVBpf9qTXAY5DuZpGuZouWLMpV9P4PZfmNERExGda0xAREZ+FZGmY2WQz22JmGWZ2VyOPtzOzJXWPrzCzvkGS61ozyzWz9LqvGQHKNdfMcsxs/XEeNzP7S13utWY2OkhynW1mRfWW170ByNTLzJab2SYz22BmtzUyJuDLy8dcXiyvSDNbaWZr6nL9spExAX89+pjLk9dj3XOHmdnnZvZaI4/5d3k550Lqi9rTsG8H+gMRwBrglAZjfgQ8Uff9lcCSIMl1LfCYB8vsG8BoYP1xHp8CvEnt1RgnACuCJNfZwGsBXlY9gNF133cEtjby7xjw5eVjLi+WlwEd6r4PB1YAExqM8eL16EsuT16Pdc99O7CwsX8vfy+vUFzTGAdkOOd2OOcqgMXApQ3GXMp/Ly37IvBNM2vs8rSBzuUJ59wH1F7P5HguBZ51tT4D4hpcQMurXAHnnNvnnFtd9/0hai8c1rPBsIAvLx9zBVzdMiipuxle99VwojXgr0cfc3nCzJKAi4CnjjPEr8srFEujJ5BZ73YWx754jo5xzlUBRUCXIMgFcHndJo0XzaxXI497wdfsXphYt4nhTTMbFsgnrtsscBq1n1Lr83R5fUUu8GB51W1qSQdygHecc8ddXgF8PfqSC7x5Pf4J+DlQc5zH/bq8QrE0Gmvchp8gfBnT3Hx5zmVAX+fcSOBd/vtpwmteLC9frKb21AinAo8CrwTqic2sA/BP4H+dc8UNH27kRwKyvE6Qy5Pl5Zyrds6NApKAcWY2vMEQT5aXD7kC/no0s4uBHOdc2lcNa+S+ZlteoVgaWUD9TwRJQPbxxphZWyAW/28GOWEu51y+c+5I3c05wBg/Z/KVL8s04JxzxV9sYnDOvQGEm1m8v5/XzMKpfWNe4Jx7qZEhniyvE+XyannVe/6DwPvA5AYPefF6PGEuj16Pk4BLzGwXtZuwzzWz+Q3G+HV5hWJprAKSzayfmUVQO1G0tMGYpcAP676fCrzn6maVvMzVYLv3JdRulw4GS4Fr6vYKmgAUOef2eR3KzLp/sS3XzMZR+/ee7+fnNOBpYJNz7g/HGRbw5eVLLo+WV4KZxdV9HwWcB2xuMCzgr0dfcnnxenTOzXbOJTnn+lL7HvGec25ag2F+XV5tm+sXtRTOuSozuwV4i9o9luY65zaY2f1AqnNuKbUvrufMLIPahr4ySHLdamaXAFV1ua71dy4AM1tE7Z418WaWBdxH7cQgzrkngDeo3SMoAygFrguSXFOBm82sCigDrgxA+U8CrgbW1W0PB7gb6F0vlxfLy5dcXiyvHsAzZhZGbUk975x7zevXo4+5PHk9NiaQy0tHhIuIiM9CcfOUiIicJJWGiIj4TKUhIiI+U2mIiIjPVBoiIuIzlYaIiPhMpSEiIj5TaYiIiM9UGiIi4jOVhoifmFlbM/vYzErMbEiDx2aamas7/YNIi6HTiIj4kZn1AdKB3dRe+a287joVq4A04GznXLWXGUWaQmsaIn7knNsNTAdOBX5fd8bUxUA58AMVhrQ0WtMQCQAz+ytwM/AJcDpw+XGutSES1FQaIgFgZpHAemAAMMc5N9PjSCInRZunRAJjJHXXrgCG111RTaTFUWmI+JmZxVA7j5EH/B8wEfilp6FETpI+7Yj439+BPsD5zrn3zGwUcJeZveucW+5xNpEm0ZyGiB+Z2XTgKeA3zrn/q7svjtrdcMOBkc45v16HW6Q5qTRE/KTugL40agviLOdcVb3HJgIfAG865y7xKKJIk6k0RETEZ5oIFxERn6k0RETEZyoNERHxmUpDRER8ptIQERGfqTRERMRnKg0REfGZSkNERHym0hAREZ+pNERExGf/H13riWwhJrHSAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=[6, 4])\n",
"x = np.linspace(0, 4, 50) # x ranges from 0 to 4\n",
"y = [f_x(i) for i in x]\n",
"\n",
"plt.plot(x, y, label=\"$f(x)$\")\n",
"\n",
"plt.xlabel(\"x\", size=18)\n",
"plt.ylabel(\"y\", size=18)\n",
"plt.legend(prop={\"size\": 14})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sampling"
]
},
{
"cell_type": "code",
"execution_count": 169,
"metadata": {},
"outputs": [],
"source": [
"# pre-setting\n",
"n = 1000\n",
"\n",
"mu_target = 3.5\n",
"sigma_target = 1\n",
"mu_appro = 3\n",
"sigma_appro = 1\n",
"\n",
"p_x = distribution(mu_target, sigma_target)\n",
"q_x = distribution(mu_appro, sigma_appro)"
]
},
{
"cell_type": "code",
"execution_count": 170,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jeremy.zhang/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
" warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n",
"/Users/jeremy.zhang/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
" warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x117979f98>"
]
},
"execution_count": 170,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAEKCAYAAAAo4eD5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl4VdXV+PHvvjfzPCeEjEASAgECRIIgkwIOWCapWqVqbevcvq3Fqq3i0Fp9tfJaf1XrSKuVQWtVRBxAQOYhzGMgCSEJCRnJPCf798cNGCCQm3DPJcP6PE+e5J6zz147wcR199lnbaW1RgghhBBCGMt0uQcghBBCCNEbSNIlhBBCCGEHknQJIYQQQtiBJF1CCCGEEHYgSZcQQgghhB1I0iWEEEIIYQeSdAkh2qWUuksppVt9VCmlMpVSnyqlblZKmVq1jWppc1cH+p+olHq6dT8dGFNUq2OZSql/W9tHZ8fVme9RCCEk6RJCdMSPgSuBG4AngTpgMfCtUsq1pU1eS5svO9DvROApOvY36cuWOHkduKajJtL2uDrzPQohejmHyz0AIUS3sltrndbq9QdKqY+Bj4EXgV9preuALUYNQCnlCDRqrQuBQqPiXIzR36MQomeSmS4hxCXRWn8CfA78Uinl1tatN6XUFUqplUqpYqVUtVIqQyn1esu5p7HMJgE0nL6F2XLudF8PKKVeVErlYpld82nr9mKreL9USqUppWqVUjuVUpPOOb9WKbW2jesylVL/7MC47jrn+rlKqT0tcYuUUh8opfq0EePfSqlblVKHWm7Vpiilrjqn3QV/ZkKI7klmuoQQtrACmAkkAVmtTyilPIBvgG3AXUAFEAWMaWnyDhAG/By4Cmhqo/8/AtuBewAzUHuRsUwARrZcUwc8CnyllBqmtU7twPdkzbjOUErdA7wJLAUeB0KBvwDJSqkRWuvKVs3HAXFYbtHWAn8CliulorTWpVb8zIQQ3ZAkXUIIWzidaPXhnKQLGAj4Ar/XWu9tdfyfAFrrHKVUTsuxrVrrxjb6zwdm6VabxSqlLjSWYGCs1jqrpd13wHHgCeCn1n5DVo7r9FjMWBKntVrrW1sdPwysB+4GXm11iReQqLU+1dLuJJak8gZgEe38zIQQ3ZPcXhRC2MLpDEi3ce4oUAq82XL7LbwT/X/WOuFqx5bTCReA1rqCHxbdGyUOCAI+bH1Qa70BS8I34Zz2m08nXC32tXyOaPlsi5+ZEKKLkaRLCGELp5OC854k1FqXAZOAXOB1IEsptV8pdVMH+u/IE4r5FzjWtwN9dJRfy+e2xnmy1fnTSlq/aFmYD+DS8toWPzMhRBcjSZcQwhamYVmbtKOtk1rr3Vrrm7AkH1cC6cBHSqkEK/u3dpYLLLcX2zp2otXrWsCpjXbnJkfWOp1EhbRxLgQo7miHNviZCSG6GEm6hBCXRCk1G5gO/ENrXX2xtlrrRq31FiwLyE1AfMup0zM9rm1e2DGjW9+OU0p5YkkKN7dqcxyIVUo5tWo3HvA8py9rx5WKZTbt1tYHlVJjgEjg+458A61d5GcmhOhmZCG9EKIjEpVSAVhmiSKAG7EUTF2J5Ym98yilbsTy1OFnwDHAHfg1lifyTidCB1s+/04p9RXQpLVO6eQY87EUa32aH55edMey0P20JS1jeq+lREQ08DBQdk5fVo1La92klJqPZQ3Wv4F/Y7md+RyW9VkLO/INWPkzE0J0M5J0CSE64uOWz7VAAbATy+zOfy6y0P0oUINlpqYPlsRhOzBFa3366cDlWNYuPQDMx7Iw/4KPJ7bje2AtlnINYVgSp+u11kdON9Bar1FK3QfMA24CdgFzgU/O6cvqcWmt31JKVQOPYKlbVomllMbvzykXYQ1rfmZCiG5GWf9AkBBCCCGE6CxZ0yWEEEIIYQeSdAkhhBBC2IEkXUIIIYQQdiBJlxBCCCGEHXS5pxcDAgJ0VFTU5R6GEEIIIUS7duzYUaS1DrSmbZdLuqKiokhJ6Wx5HiGEEEII+1FKHbe2rdxeFEIIIYSwA0m6hBBCCCHsQJIuIYQQQgg76HJruoQQQoieoKGhgZycHGpray/3UIQNuLi4EBYWhqOjY6f7kKRLCCGEMEBOTg6enp5ERUWhVGe3EhVdgdaa4uJicnJyiI6O7nQ/cntRCCGEMEBtbS3+/v6ScPUASin8/f0vedZSki4hhBDCIJJw9Ry2+LeUpEsIIYQQwg4k6RJCCCGEsANZSC+EEO1YtDXL6ra3JUcYOBIhRHcmSZcQQghhBx1J3q3RmQT/6aefxsPDg3nz5jFmzBg2bdrUZrvS0lIWLVrEAw88cMG+Tl+fmZnJjTfeyP79+60aQ1t9X2wstlJTU8N1113H6tWrMZvNbbapr69n8uTJrF69GgcH26dIcntRCCGE6IUuluSUlpby+uuvt3lOa01zc3Onk6S2+jY64QJ47733mD179gUTLgAnJyeuueYali5dasgYJOkSQggherDnnnuOuLg4Jk+eTGpq6pnjHh4eAFRVVTFt2jSGDRtGQkICS5cu5bHHHiM9PZ3ExEQeeeQRMjMziY+P54EHHmDEiBFkZ2efuR6gsbGRO++8k6FDhzJnzhyqq6vJzMwkISHhTJu//vWvPP300+f13XosAAsWLCAhIYGEhAReeeUVgDPxf/nLXzJ48GCmTp1KTU1Nm9/vrbfeyi233EJycjKRkZF8+eWXAHz44YfMmDHjTLtJkyaxcuVKAJ544gl+/etfAzBz5kw+/PDDzv/AL0JuLwohhBA91I4dO1iyZAm7du2isbGRESNGMHLkyLPafP3114SGhp5JTsrKykhOTmb//v3s3r0bsCQ9qampLFy4sM0ZsNTUVN59913Gjh3L3Xffzeuvv86cOXPaHNMLL7xwVt/njnfhwoVs3boVrTXJyclMmDABX19fjh49yuLFi3n77be5+eab+eSTT5g7d+55fezZs4eZM2eydOlSNmzYwMMPP8yUKVPIyMggKirqTLtnnnmG+fPnU1BQwK5du1i2bBkACQkJbN++3bofcAfJTJcQQgjRQ61fv55Zs2bh5uaGl5cX06dPP6/NkCFDWLVqFY8++ijr16/H29u7zb4iIyMZPXp0m+fCw8MZO3YsAHPnzmXDhg2dGu+GDRuYNWsW7u7ueHh4MHv2bNavXw9AdHQ0iYmJAIwcOZLMzMzzrq+pqaGoqIinnnoKgEGDBnHq1CmKiorw8fE5q+348ePRWrNgwQKWLFly5raj2WzGycmJioqKTn0PFyNJlxBCCNGDtVfUMzY2lh07djBkyBAef/xxnn322Tbbubu7Wx1DKYWDgwPNzc1njllTzV1rfcFzzs7OZ742m800Njae12b//v3ExMTg4uICwM6dOxk2bBiurq7nxd+3bx95eXk4Ozvj6el51rm6urozfdiSVUmXUuo6pVSqUipNKfXYRdrNUUpppVRSq2OPt1yXqpS61haDFkIIIUT7xo8fz6effkpNTQ0VFRV88cUX57XJzc3Fzc2NuXPnMm/ePHbu3Imnp2eHZnqysrLYvHkzAIsXL+aqq64iODiYgoICiouLqaurY/ny5QAX7Xv8+PF89tlnVFdXU1VVxaeffsq4ceOsHseePXvIysqitraWqqoqnnrqKX7729/i6+tLU1PTmcQrLy+P22+/nc8//xx3d3e++eabM30UFxcTGBh4SRtbX0i7a7qUUmbgNWAKkANsV0ot01ofPKedJ/BrYGurY4OAW4HBQCiwSikVq7Vust23IIQQQnR9l6OG24gRI7jllltITEwkMjKyzQRm3759PPLII5hMJhwdHXnjjTfw9/dn7NixJCQkcP311/Pggw9eNE58fDz/+te/uPfee4mJieH+++/H0dGR+fPnk5ycTHR0NAMHDgQ4r++XXnrprPHeddddjBo1CoBf/OIXDB8+vM1biW3Zs2cPt99+OxMnTqS8vJw//OEPZ257Tp06lQ0bNjBmzBhmz57Nyy+/THx8PE8++SSPPvoo115rmRdas2YNN9xwg1XxOkpdbCoPQCl1JfC01vraltePA2itnz+n3SvAKmAeME9rnXJuW6XUNy19bb5QvKSkJJ2SktL570gIIWysdX2l+sZmiirrKKioQ2tNiLcLQZ4umE0d25dNiqj2fIcOHSI+Pv5yD6NXGT9+PG+//TZxcXHnndu1axcLFizggw8+uGgfs2fP5vnnn2+zj7b+TZVSO7TWSec1boM1Ty/2BbJbvc4Bks8JOBwI11ovV0rNO+faLedc2/fcAEqpe4B7ACIi5A+REMI+rC1WWVnXyMa0IvbmlFJa3cC5b1UdTIpgLxeG9PXmyv7+OJpluawQl0N6ejoxMTFtnhs+fDiTJk2iqanposVRZ86c2WbCZQvWJF1tvX078zdHKWUC/g+4q6PXnjmg9VvAW2CZ6bJiTEIIYbhTVfWsTyskJfMUTc2auBBPRkT6EuTpQqCnMwrIK6slr6yGrJJqvj5wks0ZxUwZFExiuA+mdhYwCyFs68SJExc9f/fdd1/0vJOTE3fccYcth3QWa5KuHCC81eswILfVa08gAVjb8vRCCLBMKTXdimuFEKLLadaajWlFfHsgH4DECB/GxwQS6Ol8XttgLxcSwy2PomcUVvLV/pP8Z0cOG9OK+PHIcEK8bf8ElBCie7Im6doOxCilooETWBbG33b6pNa6DAg4/VoptZYf1nTVAIuUUguwLKSPAbbZbvhCiG4rZaH1bZN+Ztw4zlFZ18h/dmRzJL+SwaFe3Dg0FG9X655i6hfowf0T+7PvRBkr9uXx5rp05o6OpH+gR/sXCyF6vHaTLq11o1LqIeAbwAy8p7U+oJR6FkjRWi+7yLUHlFIfAQeBRuBBeXJRCNFVpRdW8lFKNjX1TcxIDGVUlF+7NY7OZVKKYWE+RPq58c9NmfxzYyZzRoYxLNyn/YuFED2aVdsAaa1XACvOOTb/Am0nnvP6OeC5To5PCCHsYlfWKT7ZmYO/hzM/GxN9ybcFfdycuHd8fz7YcpylKdmU1TQwPjbQRqMVQnRH8oiNEKLX25xexMc7cogKcOeBCf1ttg7L1cnMz8ZGMaSvN18fOMmWjGKb9CuE6J5kw2shRK+ltWZNaiGrDuUT38eLW68It3m5B0eziVuuCKehqZkv9uTi6+ZIXIiXTWMIIboHmekSQvRKWmu+PZjPqkP5DA/34bZREYbV1zIpxS1XhNPH24XF27PJK6sxJI4QomuTpEsI0XVpDU0NhnT9xvfpfH+kkFHRftw0MqzDFeU7ytnBzE+vjMLFwcT7m4+TX97+5r9CiJ5Fbi8KIbqG+irI2wu5O6Ei15JsNdVbzq1/GfqOgNAREJEMkWPB1HZFaWv8e8txXvw6lcRwH6YPC7VbEVNvV0fuuDKKt9ZlcM/7KXx83xicHOS9b6/RkTIp1uhEKZWnn34aDw8P5s2bx5gxY9i0aVOb7UpLS1m0aBEPPPDABfs6fX1mZiY33ngj+/fvt2oMbfV9sbHYSk1NDddddx2rV6++aEX6yZMns3r1ahwcbJ8iSdIlhLi8ynPh8HIoPAy6GdwDIWQoOLiA2cmSXFUWQE4KpLY8RO3qCxFXQngyuHh3KNyybGee3OrFNQODmRgXZPeq8aE+rvw4KYwPt2bx129T+cMNsjefuDwuluSUlpby+uuvt5l0aa3RWnc6SWqrb6MTLoD33nuP2bNnXzDhAktF+muuuYalS5dy++2323wM8hZLCHF56GZIXwMbXobSLIieCOPmwcQ/wNBbYNAMiLseYqbC8Lkw6Q9w7fMw4i5LYpa6Ar57BnZ+AFVFVoVcd9KJh7d5cUVAA6/dPsLwW4oXMjjUm7mjI3hrXQbrjhReljGI3uO5554jLi6OyZMnk5qaeua4h4elaG9VVRXTpk1j2LBhJCQksHTpUh577DHS09NJTEzkkUceITMzk/j4eB544AFGjBhBdnb2mesBGhsbufPOOxk6dChz5syhurqazMxMEhISzrT561//ytNPP31e363HArBgwQISEhJISEjglVdeATgT/5e//CWDBw9m6tSp1NS0vTby6NGjTJw4kaSkJH7/+98zYMAAAD788ENmzJhxpt2kSZNYuXIlAE888QS//vWvAZg5cyYffvhh53/gFyEzXUII+yvLgS1vQPFRCE6AobeCsxVV2x1dITTR8lFVCMc3wfGNkLcbIsdAzLUX7OdwmZkHtngR49XIO2PLcHHs/O1JW3hi2iC2HSvh4Y/28PVvxhHgcf4WQ0Jcqh07drBkyRJ27dpFY2MjI0aMYOTIkWe1+frrrwkNDeXLL78EoKysjOTkZPbv38/u3bsBS9KTmprKwoULef3118+Lk5qayrvvvsvYsWO5++67ef3115kzZ06bY3rhhRfO6vvc8S5cuJCtW7eitSY5OZkJEybg6+vL0aNHWbx4MW+//TY333wzn3zyCXPnzj3r+qamJu644w5ee+01RowYwa9+9SsGDx5MfX09GRkZREVFnWn7zDPPMH/+fAoKCti1axfLlllqvSckJLB9+3brf8gdIEmXEMK22lu3UpkPm/6fZb3W0Fsttwg7c4vPPdAyG9ZvIhz52pJ85WyD2OshejyoHybyC2pN/HyjD+4OmvfGluHlqDsez8ZcHM28+pPhTP/7RuZ9vIf37rwC02WaeRM91/r165k1axZubm4ATJ8+/bw2Q4YMYd68eTz66KPceOONjBs3jlOnTp3XLjIyktGjR7cZJzw8nLFjxwIwd+5cXn311QsmXRezYcMGZs2ahbu7OwCzZ89m/fr1TJ8+nejoaBITEwEYOXIkmZmZ513/2WefMWjQIEaMGAFAfHw8Pj4+FBUV4eNz9q4Q48ePR2vNggULWLt27ZnbjmazGScnJyoqKvD09Ozw93AxcntRCGE/NacsM1wA434HEaM7l3C15uJtuR054VHw7QcHP4MN/2eZTQNqGuGXG70pqTPx7tgy+rg1X+I3YTsDQ7x4clo8a1ML+eemzMs9HNFDtbeVVWxsLDt27GDIkCE8/vjjPPvss222O50IWRNDKYWDgwPNzT/8vtXWtv/ErtYXfkPk7PzDbLDZbKaxsfG8Nrt27TqTmAHs2bOHYcOG4erqel78ffv2kZeXh7Oz83nJVV1dHS4utt+sXpIuIYR91FVaEq7GWki+DzyCO9zF1mMlF/4odGRr4ByOht0EtaWwYQH64DJ+v82FvacceDW5jATf8/9IX25zR0dyzcAgXvzmMJlFVZd7OKKHGT9+PJ9++ik1NTVUVFTwxRdfnNcmNzcXNzc35s6dy7x589i5cyeenp5UVFRYHScrK4vNmzcDsHjxYq666iqCg4MpKCiguLiYuro6li9fDnDRvsePH89nn31GdXU1VVVVfPrpp4wbN87qcfj7+3P48GEAtm7dyvvvv8/QoUPx9fWlqanpTOKVl5fH7bffzueff467uzvffPPNmT6Ki4sJDAzE0dG6je47Qm4vCiGM11AL2/5hmekafT94hxkTRylKvAfDoJFw+AtUxmp+03yIqbE/ZUpoqDExL5FSiudmDWHKgu957L97WfSL0XKbsafqRImHSzVixAhuueUWEhMTiYyMbDOB2bdvH4888ggmkwlHR0feeOMN/P39GTt2LAkJCVx//fU8+OCDF40THx/Pv/71L+69915iYmK4//77cXR0ZP78+SQnJxMdHc3AgQMBzuv7pZdeOmu8d911F6NGjQLgF7/4BcOHD2/zVmJbfvrTnzJt2jSGDBnCDTfcgL+//5mF9FOnTmXDhg2MGTOG2bNn8/LLLxMfH8+TTz7Jo48+yrXXXgvAmjVruOGGG6yK11HqYlN5l0NSUpJOSUm53MMQQnTWuWu6tIYdCyF/PyT9AoIHdbrrrcdKrGqXHO3H8mxnFm3P5jXXN/FpKkFFj4eB0yxlKACSfsairVnt9tU/62Orx5ce8WOr296WHHHW6yXbsnjsv/t4blYCtydHWt2P6LoOHTpEfLyUBLlcsrOzmTNnDlu3bgUstx4XLFjABx98cNHrZs+ezfPPP09cXNx559r6N1VK7dBaJ1kzJpnpEkIYK2c7nNwL8dMvKeHqiAOlDjyS4sUg31jcxjyCOrIcjn0P+Qdg2E/Av79dxtERt1wRzhd7c3l+xWEmxQUR6uMKYFViCOcncUL0dnv27GHo0KFnXg8fPpxJkybR1NR00eKoM2fObDPhsgVZ0yWEME51CRz4BPz6W54ytIPyBjP3bPLG26mZN64sx9nZBYbMgdEPAho2/x32/9dSAb8LUUrxwuyhNDVr/vDpvosuKBZCtO/GG2/k7bffPuvY3Xff3W5x1DvuuMOwMVmVdCmlrlNKpSql0pRSj7Vx/j6l1D6l1G6l1Aal1KCW41FKqZqW47uVUv+w9TcghOiidDPsWWT5OvG2s0o4GKVJw/9lhFJYa+LNK8sIcmn1pGJADIz/PURdBZnr4I0xhBQZXwW7I8L93Hj0ujjWphby+e7cyz0cIYSNtftXUCllBl4DrgcGAT85nVS1skhrPURrnQi8CCxodS5da53Y8nGfrQYuhOjijn0PxWkweDa4+dsl5OITgRysdOf5kRUM82vjSUUHZ0i4Ca58CJSZq7ffy9jdj+BaW2CX8VnjjiujSAz34c9fHqKsxpjNvoUQl4c1bz1HAWla6wytdT2wBJjRuoHWurzVS3dA5sWF6M0qC+Hwl5Zq82Gj7BJy6ykPvsj3Z0rgKW6KbKcekP8AuH8Tewc8QFj+am5cN524zH+jmi9/SQmTSfHnmQmUVNWx4NvU9i8QXZrcJu45bPFvac1C+r5AdqvXOUDyuY2UUg8CDwNOwNWtTkUrpXYB5cATWuv1nR+uEKJbOPS5ZaPqITdfevFTK+TWOvJGZh/6u9VwZ1gBW4+1/8cxvakAYu4nM3QaSQf/wshD/0v0iWVsH/wExT5D273eSAl9vbnjyije35zJ/RMG0NfX9bKOR3SOi4sLxcXF+Pv7t1ugVHRtWmuKi4svuWCqNUlXW/+lnPcXTWv9GvCaUuo24AngTiAPiNBaFyulRgKfKaUGnzMzhlLqHuAegIgIeQJHiG6tMNVSHmLgj8DFy/BwtU2KBel9cVCah/ufwNHUsXejle4RrE16g/CTKxl56H+ZunkuaeFz2BP7P9Q7eRs06vY9PDWW5Xvz+HzPCe6b0B+T/E+72wkLCyMnJ4fCQtnUvCdwcXEhLOzSagxak3TlAOGtXocBF1vhuQR4A0BrXQfUtXy9QymVDsQCZxXi0lq/BbwFljpd1g5eCNHFNDVatuFx84foCXYJ+V52CDm1zjw2IIcAp07eHlSK7D5TyQscy9CjrxF7fBHhJ1exe+BvLXXGDEh4rCkFcfXAID5KyWZ7ZgnJ0fZZFydsx9HRkejo6Ms9DNGFWLOmazsQo5SKVko5AbcCy1o3UErFtHo5DTjacjywZSE+Sql+QAyQYYuBCyG6oJ3/goo8S00us/FlANcVe/F9sTez+xST6H3pJSAaHdzZGf97vh6zlAr3SEbvm8+gzH/iWptvg9F23LAwb/oFuPPNgZNU1l3+9WZCiEvTbtKltW4EHgK+AQ4BH2mtDyilnlVKnd6u/CGl1AGl1G4s67rubDk+HtirlNoD/Ae4T2ttXUlpIUT3UlMKa56z1OQKMX5NVF6tI+9khTDQo5qb+hTZtO9SrzhWjv4XW4Y8i0tdMUPS3yIsfw2qucmmcdqjlGL6sFDqG5tZdfDyJH5CCNux6q2o1noFsOKcY/Nbff0/F7juE+CTSxmgEKKbWPeSpRjqyJ8Zvni+oVnxt2OWdVy/is7FbEQ4ZSIjbBYO9RVE5n9L36L1+FQcIaPvDKpdQwwI2LYgLxeS+/mzJb2Y5H5+9PGWRfVCdFdSkV4IcenK82Db25YiqEZtZt3K4hOBHKt24f6ovM6v47JSk4MrGX1nkBpxC45NVQzOeIfQwg2WtV52MnlgMK5OZr7cmyclCIToxiTpEkJcuo1/g+ZGGP+I4aF2lbnzZYEf1wae4gqfSsPjnVbqGcfe/vdzyiue8ILVxGR/jKmpzi6xXZ3MTI4PJqOoigO55e1fIITokiTpEkJcmoqTsGMhDLsV/Ix9Uqui0cQ/MvsQ7lLL3DD7V5FvcnAlLWw2x4On4luRyuBj7+FcZ59lqldE+RHs5cxX+/NoaGpu/wIhRJcjSZcQ4tJsfBWaGmDc7wwP9W5WCBVNZh6KzsOpg/W4bEYpTgaM5nDk7Tg2VpKQ8Q6eVZmGhzWbFNOGhHKquoGNabZ9cEAIYR+SdAkhOq+yAFLeg6E3g39/Q0Mty3Zm8ykv5vQpIsrNPrf1Lqbcox8H+v2CBgcP4o4vtkviNSDIg/g+Xqw9UiglJITohiTpEkJ03qZXoanO8LVc+TUmntzpyQD3GmaEFBsaqyPqnHw5GH0HdU4+xB1fTFDxdsNjXj84hMamZlYflhISQnQ3knQJITqnshC2vwtDfmzoLJfW8OgOT+qaFQ9GGVQe4hI0OnhwKOqn1Dn5MDHlAYKKtxkaL8DTmSui/Nh2rISiiss/4yeEsJ4kXUKIztn8/6ChBsbNMzTMx8ddWHvSmceGVBLq0mBorM6yJF53UOkWxsQdD+FTftjQeFcPDMLBZOKbgycNjSOEsC3j9+kQQvQ8VcWw7R1IuAkCYw0LU1ireG6PB6MC6rmjfw3bM23Xd/+sj23XGZYthFaPepvrNt7C+J3/wzdXLqbO2c+mMU7zdHFkXGwA3x0qIKu4igh/d0PiCCFsS2a6hBAdt/nv0FANE35vaJhndntS06T4y4gKTF3stmJbap0DWDfyVVzrirhq9zxUs3Ezc1cNCMDT2YGvDpyUgqlCdBOSdAkhOqa6BLa9BYNnQWCcYWFW5TqxPMeFX8VXMcDLvnseXooS78FsTXia4JLtjDj8kmFxnB3MXB0fxPHiag7lScFUIboDub0ohOiYza9BfaWhs1wVDYond3kS59XIvXHVhsUxSmbfH+Fbfpj4zPcp8RrEsbCZZ53vyK3N9IgfX/BcUqQfm9KK+fZgPs/MSMDcHaYDhejFZKZLCGG96hLY+iYMmgFB8YaF+et+d07WmHh+ZDlO3fTUmDR8AAAgAElEQVSv1O6433LSbxRJB5/HrSbXkBhmk2LyoGAKKur4Yo8xMYQQtiMzXUKI9qUstHxOXQH1FRAQ98MxG9tb4sD76a7cOaCGEf7dtwCoNjmwZeifmLZ+Fsn7nmLNFW+Bsv1M1OBQL/p4u/DKqiNMG9oHR3M3zVKF6AXkt1MIYZ36aji2DkKGgleoISGaNczf7Ym/czMPD64yJIY9VbuGsmvg7+hTvIUB2bZ9WvI0k1JMjg8ms7iaT3bkGBJDCGEbMtMlhLBO5jporIXYaw0L8Z9MF3aXOLLginK8HLvfE3ltrdXSGsrcoxl56H9xrC+j3snH5nEHhniSGO7Dq98dZdaIvjg7mG0eQwhx6aya6VJKXaeUSlVKpSmlHmvj/H1KqX1Kqd1KqQ1KqUGtzj3ecl2qUsq4v9ZCCOM01MCx7yF4CHj1NSREWb3ihX0eJPnXMyui1pAYl4VSZIT+CIB+uV9YsjCbh1DMmxpHblktS7Zl27x/IYRttJt0KaXMwGvA9cAg4Cetk6oWi7TWQ7TWicCLwIKWawcBtwKDgeuA11v6E0J0J8fWWRKv2KmGhVhwwJ3SesUzwyuNWPp0WdU7+ZAVPBnvqmMElO0zJMbYAf4kR/vx9zVp1NR3nxIbQvQm1sx0jQLStNYZWut6YAkwo3UDrXXrIjHuwOm3cjOAJVrrOq31MSCtpT8hRHdRW94yyzUYvMMNCXGw1IEP0l2Z27+GwT7dd/H8xRT4jqTSJZTw/O8wGVA0VSnFvGvjKKyo4/3NmTbvXwhx6axJuvoCreerc1qOnUUp9aBSKh3LTNevO3jtPUqpFKVUSmFhobVjF0LYw7Y3LdXnY64zpHut4endHvg4aX7XAxbPX5BSZIVMwamxgpCizYaEuCLKjwmxgfzj+3QqarvmPpVC9GbWJF1tTfSftyhBa/2a1ro/8CjwRAevfUtrnaS1TgoMDLRiSEIIu6irsBRDDRoEPsbMcn2b68S2IidmhxRw+EQxW4+VXPCju6twj6TEK57Qoo04NlQYEuN3U2M5Vd3Awo2ZhvQvhOg8a5KuHKD1X9sw4GJV+JYAp8svd/RaIURXsu1tqDkFscbMcjU0wwv7POjrUsfVAaWGxOhqsoKvQdFEWMEaQ/ofGubD1EHBvL0ug9LqekNiCCE6x5qkazsQo5SKVko5YVkYv6x1A6VUTKuX04CjLV8vA25VSjkrpaKBGGDbpQ9bCGG42nLY9CoMmAI+EYaEWJThyrFKB+aGFWDuYYvnL6TOyY98v1EElu7GreakITEenhpLZX0jb63LMKR/IUTntJt0aa0bgYeAb4BDwEda6wNKqWeVUtNbmj2klDqglNoNPAzc2XLtAeAj4CDwNfCg1loeqxGiO9j6D8ss16THDem+vEHxykF3rgysZ7hXD17L1YYTgeNpNLsSkb/SkP4Hhnjxo6GhLNyYSVFlnSExhBAdZ1VxVK31CmDFOcfmt/r6fy5y7XPAc50doBDiMqgphU1/h7gboO9IyNvbqW4utg5rUU4gp+pNTA840eNKRLSnyexCbuA4Ik9+i2fVcSrcI20e4zeTY1i+N5c31qbz5I3nVvkRQlwOsg2QEOJ8m1+DujKY9AdDui+sc2BFgS/j/Mro59Y7Z2LyfUdS7+BO38J1hvTfL9CDm0aE8cGW4+SV1RgSQwjRMZJ0CSHOVl0CW96AQTMgZIghIT7OszylfGvf3lsiRpscyfMfg3fVMTyqswyJ8etrYtBa89qaNEP6F0J0jCRdQoizbfwb1FfCRGPWcp2odWJdsRdTA0sJcOqZhVCtVeA3kgazG30LjJntCvdz4+akcJZuzybnVLUhMYQQ1pMNr4UQP6gshG1vQcJNEBRvSIiPcwNwMmlmhBQb0n930mxyIi9gDBH5q/CozqHSLazTfS3a2vZsWZivG80afrNkN7NHhHFbsjFPogoh2iczXUKIH3z/AjTWwcTz9rW3iePVzmw+5cUNQSV4O8qDzAD5vkmW2S6D1nZ5uzoyKsqPnVmnKJYnGYW4rGSmS4jeLGXhD19XnISU9yBiDGRusHzY2Md5AbiZm7gxuPtXl7eVZrMTef6jiShYjXv1Carcztsp7ZJNiAtke2YJqw8X8KtrYtq/QAhhCJnpEkJYHFoGZifDqs+nV7mwvdSTaUEleDg0GxKju8r3u4JGswuhRRsN6d/LxZHR/fzZnV1KWkGlITGEEO2TpEsIAYWpUHDQUn3e2cOQEB/lBuBhbuKG4FOG9N+dNZudyfdNwrfiMM51xswCjo8NxNFs4m/fHW2/sRDCEJJ0CdHb6WY4+Dm4+kH0eENCHKl0YXe5B9NDinEzyyxXW/L9RqGVmT7FWwzp38PZgSv7+7N8by6pJ43ZbFsIcXGSdAnR22Vvg4pciP8RmB0NCfHpyQA8zY1cGyizXBfS4OhBkfdQAkp349BoTHmHcTEBeDg58MqqI4b0L4S4OEm6hOjNGusgdQX4RkGfRENCZFY7s7PMg+uDT+Fi1obE6ClO+o/GrBsJKkkxpH83Jwfuviqar/af5EBumSExhBAXJkmXEL1Z+ndQVw6DZmLUBoifnfTH1dQks1xWqHEJ5JTHAEJKtqOajSkce/dV0Xi5OPB/K2W2Swh7k6RLiN6q7ASkr4HQ4ZaZLgPk1Tqy5ZQnUwNL5YlFK50MuBLHpioCyjq3yXh7vF0duWd8P1YdKmB3dqkhMYQQbZOkS4jeavWfAA0DbzQsxLJ8fxyU5gapy2W1crcoqlz60KdoC2hjbsfeNTYaXzdHFshslxB2JUmXEL3RiZ2wZzFETwA3f0NC5FWb+L7Ym6sDyvCR6vPWU4o8/9G41hfhU2lMeQcPZwfum9CfdUcKScmUhFgIe5GkS4jeRmv49glwC4ABkw0L8/YRN9Dwo2DZY7GjSrwHUe/gSUjxNsNi3HFlFAEezrz8rcx2CWEvViVdSqnrlFKpSqk0pdR5m7IppR5WSh1USu1VSn2nlIpsda5JKbW75WOZLQcvhOiEw1/C8Y0w6XFwdDUkxKk6xeJjroz1KyfQ2ZgF4T2ZVmby/ZLwrsrAqyLdkBiuTmYenNSfzRnFbEwrMiSGEOJs7SZdSikz8BpwPTAI+IlSatA5zXYBSVrrocB/gBdbnavRWie2fEy30biFEJ3RWA8rn4SAOBhxl2FhFmW4UtOkmB4is1ydVeA7kmblQNzxDw2LcVtyBKHeLrz0TSraoPVjQogfWDPTNQpI01pnaK3rgSXAjNYNtNZrtNanq/ltAcJsO0whhE2kvAslGTD1z2A2Zr/7+mb4V7or44LrCHetNyRGb9Do4EaR9xCiT3yBU70xNbWcHcz8+poYdmeX8t2hAkNiCCF+YE3S1RfIbvU6p+XYhfwc+KrVaxelVIpSaotSamZbFyil7mlpk1JYWGjFkIQQHVZdAmtfgH6TIGaKYWGWZ7tQUGvmFzE1hsXoLU76j8KhuZb+2f8xLMZNI8OI8nfjr9+m0twss11CGMmat7ptVUxs8zdTKTUXSAImtDocobXOVUr1A1YrpfZprc9apKC1fgt4CyApKUl+64W4FCkL2z5+4FOoLYOwK2DHPw0JrTW8c9SVGK9GxgfXsy3TkDC9Ro1LMCf9k4nNWsLh6DvRJtvPTjqaTfx2Siz/s2Q3X+7L40fDQm0eQwhhYc1MVw4Q3up1GJB7biOl1GTgj8B0rXXd6eNa69yWzxnAWmD4JYxXCNEZlYWQuQHCk8HLuP+pbi505GCpIz+PqTaqwH2vkxp5O+61JwnL/86wGD8aGkpcsCf/t/IIjU1SxFYIo1iTdG0HYpRS0UopJ+BW4KynEJVSw4E3sSRcBa2O+yqlnFu+DgDGAgdtNXghhJUOfwEmM8TdYGiYd4+64e/czMyIWkPj9Ca5QeOpcAtnYOa/DYthMikenhpLRlEV/911wrA4QvR27c5Va60blVIPAd8AZuA9rfUBpdSzQIrWehnwEuABfKwsb2+zWp5UjAfeVEo1Y0nwXtBaS9IlhD0Vp8HJvZaEy8XrkrvbeqztYpq5tY58lxfEnD5F7MmSgpu2opWZI5G3MfLQ/+JXdoAS78GX1N+irVltx9GaMF9X/rLiELX1TTiYTdyWHHFJsYQQZ7OqTpfWeoXWOlZr3V9r/VzLsfktCRda68la6+BzS0NorTdprYdorYe1fH7XuG9FCHEe3QwHPwcXH+g30dBQKwr8cFDNTJGNrW0uve9MGszuxBk426WUYkp8MKXVDWw/Lv+GQhhBKtIL0ZOd2AFl2Zb9Fc1OhoWpajSxrtibsX7lsuWPARodPcgIm0lE3te41Br3hPeAIA+i/N1Ze7iA+kZZ2yWErUnSJURP1VQPh5eDdzj0HWFoqO+LvalrNnGdzHIZJjXyNky6iZispYbFUEoxdVAwFXWNbMmQwrZC2JokXUL0VBlrLSUiBs0EZdyvutbwbaEPA9xr6Ode1/4FolMq3SM4ETSBmOyPMTUZ93OOCnAnNtiD748UUlHbYFgcIXojSbqE6IlqyyBtFYQMBf/+hobaV+FGXp0z18osl+FSI2/Hpb6EyLyv2m98CabEh1DT0MS7G44ZGkeI3kaSLiF6otSvoLkJ4n9keKhvC33xdGhktG+F4bF6u3z/ZEo9BljKRxi4V2JfX1cGh3rxzvpjnKqSrZyEsBVJuoToaSryIHsrRF0F7oGGhiqqdyCl1IOrA8pwMslmEoZTitSoufhWpBJUkmJoqMnxwVTVN/KP79PbbyyEsIokXUL0NIeXg4MzxEw1PNSqQh8ApgTIrUV7yQydRq2jD3HHPzQ0TrCXCzMT+/KvzZkUlEuxWyFsQZIuIXqS45sg/wAMuAac3A0N1dCs+K7Ih5HelQQ6NxoaS/ygyexCevgc+uavwb06x9BYv5kcQ2OT5u9r0gyNI0RvYfvdU4UQl4fWsHI+uHhD9IT221+iLac8KW90YGqQzHLZ25GIW4g/9k9ijy9mV/wjhsXZmFbM8AhfPtySRbCnC77ubdd6k8r1QlhHZrqE6CkOLYOc7RB7vaGFUE/7ttCHPs71DPGsNjyWOFuNawhZIVPon/MpDo3G/vyvHhiEUrD6cEH7jYUQFyVJlxA9QVMDfPcsBA6EsCsMD3es2pkjVW5MCTyFSRkeTrQhNfJ2nBoriD7xuaFxvF0dSY72Y2fWKQorpA6bEJdCki4heoJdH1g2tr7mKTCZDQ/3TYEvzqZmJvqXGR5LtK3YdxhF3kOJy/zQssemgSbEBeFoNrHqUL6hcYTo6STpEqK7a6iF71+C8GSIu97wcJWNJjaWeHGVXznuDrI/3+WUGnU7XtXHCS3cYGgcD2cHxgzwZ9+JMvLKagyNJURPJkmXEN3djn9CRS5M+iMo4+/1rS32pl6bmCoV6C+7rJApVDsHGV4+AmDcgEBcHE2sPCizXUJ0liRdQnRn9dWw/mWIGgf9jH9isVlbKtDHuVcT5Sbrey43bXLkSOSt9CnahHfFEUNjuTqZGR8TyOGTFWSVyMMTQnSGVUmXUuo6pVSqUipNKfVYG+cfVkodVErtVUp9p5SKbHXuTqXU0ZaPO205eCF6lZSF5398ei9UFUDo8B+OGWhdvhP5dU5cK2Uiuoy08JtpMLsyKOOfhse6sr8/7s4OrDx40vBYQvRE7SZdSikz8BpwPTAI+IlSatA5zXYBSVrrocB/gBdbrvUDngKSgVHAU0opX9sNX4herLEW0r+DwDjw62eXkB+ku+Lt0Eiyj+yz2FXUO3mTHn4TkXlf4VaTZ2gsZwczE2MDSS+sIr2w0tBYQvRE1sx0jQLStNYZWut6YAkwo3UDrfUarfXp+eYtQFjL19cCK7XWJVrrU8BK4DrbDF2IXu7Yeqivgtgb7BIuu8rE6jwnrgkoxUEWJnQph6PuAGBg5vuGxxoV7Ye3qyMrD+ajDdx0W4ieyJo/nX2B7Favc1qOXcjPga86cq1S6h6lVIpSKqWwsNCKIQnRyzVUQ8ZqCBoMvpHtt7eBf6e7YlIwObDULvGE9apd+3C8z/X0z/4Ep3pjy3g4mk1Migsiq6Sa1HyZ8RSiI6zZBqitx6HafHujlJoLJAGnV/Rada3W+i3gLYCkpCR56yREezK+h4Yau5SIAKhtgo8yXZkSWoe/k+yzaE/9sz62qt3BfncRnfsFMVlLOTDgHkPHNDLSl3VHC1l5MJ+4YE9DYwnRk1gz05UDhLd6HQbknttIKTUZ+CMwXWtd15FrhRAdUF8Fx9ZCyFDwDmu3uS0sz3bhVL2JO/pLjaauqswzltzAq4g7/iHmplpDY5lNiqsHBpFXVsuB3HJDYwnRk1iTdG0HYpRS0UopJ+BWYFnrBkqp4cCbWBKu1ht0fQNMVUr5tiygn9pyTAjRWemrobHesseinXyQ7soAz0auDGywW0zRcQej78alvoToHGO3BgJIDPchwMOZ7w7n09wsNyiEsEa7SZfWuhF4CEuydAj4SGt9QCn1rFJqekuzlwAP4GOl1G6l1LKWa0uAP2FJ3LYDz7YcE0J0Rl0FZK6H0ETw6mOXkHtKHNhzypGf9q+xR+1VcQkK/JIo8h7KoGMLUc3GJsgmpbgmPoj88jqW7zP2qUkhegqrnkHSWq/QWsdqrftrrZ9rOTZfa306uZqstQ7WWie2fExvde17WusBLR/GFhESoqdL/86yuXWs/R4C/iDdFXeHZmZHGnvLStiAUhzo/0s8ak4Qlful4eGG9PUmyNOZV1YdoUlmu4Rolzz4LUR3UVsGmRshLAk8gu0S8lSd4otsF2ZF1OLpKP9T7Q5OBE2gxHMgg9PfQekmQ2OZlGJyfDAZhVV8vvuEobGE6Akk6RKiu0hbBboJYq61W8iPMl2pa1aygL47UYoDA+7Bq/o4EXlfGx5uUKgXg/p48bfvjtLQJBugC3ExknQJ0R2UZkPWJggbBe4BdgnZpOHfGa4kB9QT623sjImwrezgayj1GMDg9LdBG5sImZTi4SmxHC+u5r87cwyNJUR3Z02dLiHE5bbuJcvnWPvNcn2X50R2lZnHh8h2L93BufW8inyGMiDnv1yx/0+UeJ+9c1t6xI9tGvua+CCGhXnz6ndpzBoehpNsWSBEm+Q3Q4iuriQDdn8IEWPA1X5bl75zxI2+bk1MDa1rv7Hocoq9BlHj5E9o0XoweLsepRS/nRLLidIaPkrJbv8CIXopSbqE6Oq+fxFMDjBgst1C7i1xYFuREz8bUC37LHZXykRu4FW41+bjW5FqeLgJsYGMjPTl76vTqG2Q29FCtEX+nArRlRUegb1L4YpfgIu33cK+e9QND4dmbomWMhHdWZH3EGqc/AgrWGuX2a7fTYnlZHkti7dlGRpLiO5Kki4hurK1z4ODK1z1W7uFzKs28WWOM7dGS5mIbk+ZOBE0Ebe6AvzKDxgebsyAAEb38+O1NenU1MtslxDnkqRLiK7q5H448F8Yfb/dnlgE+Fe6K80a7hxQbbeYwjjFXoOpdg4irOB7w59kBHh4ShxFlXV8sCXT8FhCdDfy9KIQXdXa58HZG8Y8ZLeQVY2KRRmuXB9WR7i71FzqEZQiJ2gisdkfEVC6lyLfRJuHWLT17NuJA4I8eGXVURzNJpwdzGeO35YcYfPYQnQnMtMlRFd0YiccXm5JuOz4xOJ/Ml0obzDx8xiZ5epJTnnGUekSSljh96jmRsPjTYkPprq+iY1pxYbHEqI7kZkuIbqiNX+xJFvJ99kt5KaMEl4/2I8Y9xoaygvYWm630MJoSpETPJGBxxcRWLqLNH5iaLhwPzfiQzxZf7SQ0f38cHOS/9UIATLTJUTXk7UV0lbC2N+Ai5fdwm4s8SK/3omZITI70ROVufen3C2CvoXrMTcaP5M5ZXAI9Y3NrDtSaHgsIboLSbqE6GrW/Bncg2DUL+0WsknDp3n+RLjWMtJbKtD3SEqRHXwNTo2VxGe+b3i4EC8XhoX7sCm9mPKaBsPjCdEdSNIlRFdybJ3lY9zD4ORut7ArcpzJrXNmdp9ilLJbWGFnlW7hlHgOJD5jIS51RYbHmxwfTLPWrE4tMDyWEN2BJF1CdBVaw+rnwDMURv7MbmGbNfz9kDt9XepI9qmwW1xxeWQHX4O5uZ6EtH8YHsvP3YkrovxIySyhuFK2kxLCqqRLKXWdUipVKZWmlHqsjfPjlVI7lVKNSqk555xrUkrtbvlYZquBC9HjpK6A7C0w4RFwdLFb2G9znUgtd2BWSDEmmeXq8Wqd/UkLn8OA7P/gVZlheLxJA4MwmxSrDuUbHkuIrq7dpEspZQZeA64HBgE/UUoNOqdZFnAXsKiNLmq01oktH9MvcbxC9ExNDbByPgTEwfA77BZWa/h/h9yJdG9kjJ88rthb7BtwH00mF4Yd+ZvhsbxcHBnTP4C9OWXsP1FmeDwhujJrZrpGAWla6wytdT2wBJjRuoHWOlNrvReQaopCdMaOf0JxGkx5Fsz2e7x+7UknDpQ68uDAaswyy9Vr1Dn7c7Df3YTnryawZIfh8cbHBOLqZOaFrw6jDd4DUoiuzJqkqy+Q3ep1Tssxa7kopVKUUluUUjPbaqCUuqelTUphoTxeLHqZ2jJL9fmocRB7rd3CNmv46wF3+ro1MStSNrbubQ5H/5Rq5yBGHHrR8O2BXJ3MTIoLYkNaEeuOGr+AX4iuypqkq633vx15qxKhtU4CbgNeUUr1P68zrd/SWidprZMCAwM70LUQPcCGV6C6GKb+CXs+Orgs25kDpY7MG1yJozxS0+s0mV3ZHfdb/MsP0i/nc8PjJffzI8LPjedXHKKpWWa7RO9kzZ/aHCC81eswINfaAFrr3JbPGcBaYHgHxidEz1aWA1teh6G3QKj9fjXqmuCl/R4M8mlgRoQ8VdZbZYZOo9AnkWFHXsGxwdgnVx1MJn5/XRyHT1bw3505hsYSoquyJunaDsQopaKVUk7ArYBVTyEqpXyVUs4tXwcAY4GDnR2sED3Oqmcsq9mvftKuYT9Id+VEtZk/DKmUJxZ7M6VIGfQYLvWnSEh/0/Bw04b0YVi4Dy9/e4Sa+ibD4wnR1bS7Yldr3aiUegj4BjAD72mtDyilngVStNbLlFJXAJ8CvsCPlFLPaK0HA/HAm0qpZiwJ3gtaa0m6hADI3Aj7PoLxj4BPePvtbaSsXvH/DrkzLriOq4KlUnhv1D/r47NeF/okEnfsA+rN7tQ6B5w5nh7xY5vGVUrxxxviufnNzby38RgPThpg0/6F6OqsekxKa70CWHHOsfmtvt6O5bbjuddtAoZc4hiF6HmaGmHFI5ZNrd0DIWWh3UK/ftiN8gbFY0Oq7BZTdG3ZwVfjV36QyJPfkhp5m6GxRkX7MXVQMK+vSePHI8MI8rJfTTohLjdZPivE5bD9bSg4AINmgdnJbmFzq00sTHNjVmQtg30a7RZXdG2NDu6cCByPT2UavuWphsf747R4Gpo0L3x92PBYQnQlknQJYW8V+bDmL9D/Ggix70Twn/d4oIDfDZZZLnG2fP9RVDsHEXnya0zN9YbGivR35+fjovnvzhPsyjplaCwhuhJJuoSwt1VPQUMNXP+iXUtErDvpxIoTLjwUX0VfN6ljLM6mlZljoTfg3FBG34J1hsd7cNIAgjydeXrZAZqlhIToJSTpEsKeMjfAnsUw5lcQYL9FxLVNMH+3B9EejdwTW223uKJ7qXSLoMBnOCHFW3CtNXavRA9nBx69biB7csr4764ThsYSoquQpEsIe6mvhmW/At8oGD/PrqHfSnUjs9KBZ4dX4Gy2a2jRzWQHX0OT2ZnovBWGV6qfNbwvieE+/O/Xh6mskzWGouez3yZvQvR2a56Dkgy48wtwcrdb2KxKE68ddmdaWC3jpESEaEejgxtZwVPon7uMfjmfkRE+22Z9L9qadd6xK/v588b36dz3wQ5uGNIHgNuSI2wWU4iuRGa6hLCH7O2WyvNJd0P0eLuF1Rqe3u2Jg0nz5LBKu8UV3VuRzzDK3SIYnvoyLnXG7pUY7udGUqQvm9KLyCurMTSWEJebzHQJYYTWdbeaGmH9S+DsBf6xdq3J9WWOM6tPOvPHoRWEuMrieWElpTgWeiMJGe+QdOAvbBixwNBw1yWEcCivnM92neDeCedtzytEjyFJlxBGO/oNVObDqHvB0XaFILceK7no+fIGM384GM1Q3wZ+NkBmEETH1DoHsG/A/SQe+RvhJ1eSHTLFsFhuTg7cMKQPH+/IYXtmCXNHRxoWS4jLSW4vCmGkkgxIWwVhoyAo3q6hF2YHU9Vk5qWkchzkN110wqHoOynxiifpwHM41ZcZGisx3Id+ge58c+AkBRW1hsYS4nKRmS4hjFJfDTvfBzd/GGy7xcjW2F7qwaZTXtwcWkhpSTFbLz4pJkSbtMmRLUOe5bpNP2HE4RfZMvQ5w2IppZgxrC+vrj7Kn5cf4tWfDDcslhCXi7z/FcIIWsPepVBXDiPusOltxfZUNpp453gIka61zAgptltc0TOVeg3kYL+f0e/EMvoUbjA0VqCnMxNiA1m2J5c1qQWGxhLicpCkSwgjZG2Gk3tg4I3gY9/H39/PCaa80cz9UXk42K/gvejB9ve/j1KP/iTvewrHBmNvM06MDSQ22IPHP9lHWY2UOBE9iyRdQthawSE48CkExEG/iXYNvb3Ug++LvZkRUky0W51dY4ueq9nsxOahf8GlvoQrDvzF0FgOZhN//fEwCivr+NPyg4bGEsLeJOkSwpZqSmHJ7eDgDIm3g7Lfr1jZ/2/vvsPjqM7Fj3/fnS1a9S5ZsmXZci/YxsamxfQWAqYYMCWAgTjEBJJcbvKEJPcmIb+EC0+4IQWSAKGZFiAUQ8jFCSU2BvfeLdmWLEtW79JKW87vj1kZYQxeY61Wst/P88yzM7Mz8x49o919Z86Zc/wWj5bmUuj1MXtQdPtWUsefhpRxbBrxTQor32ZI5TtRjXXC4FTmn1nEK6vLeXdrdIcjUqovRYWuMV0AAB1VSURBVPSLICIXish2ESkWkR8e4v2ZIrJGRAIiMvug924SkZ3h6abeKrhS/U4oCK9+AxpLYepciEvus9DGwGNlubQHHdwxrEKfVlRRsXn4bdSlTOCkzf8v6p2m3nn2SMbkJvHDVzfS2N4V1VhK9ZXDPr0oIhbwMHAeUA6sFJGFxpie933LgJuB/zxo33Tgp8A0wACrw/s29E7xlepH3v8V7FwEFz8I0rcDHC6uT2ZlYxI35FdT4NUfKNU7ispe/sy68qyZTCx5lLNWzGNHwRwQu+FgScFVvRrb7bSrGS97eCk/XbiZ387RpxnVwBfJ9fB0oNgYs8sY0wW8CMzquYExZo8xZgNwcJfXFwD/NMbUhxOtfwIX9kK5lepfNr8OS35tP6k47dY+DV3b5eTJshzGJLZzcY72DaGiy+fJpCznHNJad5LdsDqqsSbkp3Dn2SN5Y10Fr6/dF9VYSvWFSJKufGBvj+Xy8LpIHM2+Sg0MFevg9fl2B6hf/fWBK/++EDLwyJ5BhBDmF1bi0KcVVR+oSp9OY2IRQ/cvwuuriWqsO84qYnphOj9+bSO7anT8UDWwRZJ0Hepr3ER4/Ij2FZF5IrJKRFbV1ET3A6xUr2rYA89dBfHpcM0CuwF9H3qrKp3NLQnMHVJFjkcfr1d9RIRdebMIOtwUlb+KhAJRC+W0HPz22sm4nA6+/fxafP5g1GIpFW2R9EhfDgzpsTwYqIjw+OXAmQft+8HBGxljHgUeBZg2bVqkCZ1SsdVeD8/OhmAX3PwWJOX2afjd7R5erMhiemoLZ2ZEt+8kpQ7mdyVSkj+LMWUvMKTqXYoLr+21Yz+/vOwz6y45IY8Fy0q5+cmVXDopD4DrZvRtH3hKHa1I7nStBEaKyDARcQNzgIURHv8d4HwRSRORNOD88DqlBjZ/B7wwBxrL4NoXIWt0n4bvCMDvd+eR7Awwb2hlX9ZoKnVAU9JI9qdPZ1D9cvKqF0c11thByZxWlMGyXXVs2qcXGWpgOmzSZYwJAN/GTpa2Ai8ZYzaLyL0icimAiJwkIuXAVcCfRWRzeN964BfYidtK4N7wOqUGrmAA/nYb7F0BVzwKQ0/p8yLctzGRfT4P8wsrSXIe/PyKUn2nLOdc2jw5nLzhJ3g79kc11gUTchmc5uVva8qpatZBsdXAE1FvPsaYt40xo4wxRcaYX4bX/bcxZmF4fqUxZrAxJsEYk2GMGd9j3yeMMSPC05PR+TOU6iOhILx+O2x7Cy66H8Zf1udFeL/SzTMl8VycXc8Jye19Hl+pnozDSfGQK7FCnZy2/gdIKHptC50OB9dNL8BlOViwrJSGNu0eRQ0s2oWiUpEyBt76Hmx8Gc75b5jxzT4vQq1P+P6qZMak+JmTrw+dqP7B58lkxYSfkd2wlkk7fh/VWKnxbm6YUUBTh5/5z63BH9Q7vWrg0KRLqUgYA+/8CNY8DV+5255iUIQfrEqm2S/8bnozboc+c6L6j9K8i9hRcA3jdj9JftUHUY1VkJHA5VPy+XhXHfe+qeMzqoFDky6lDscYePfnsOwRmHE7nP1fMSnGs7u8vLffwz0TWxmVoo/Nq/5nzZjvU588lpM3/piE9vKoxjqxII1vzhzOgmWlPLV0d1RjKdVbIukyQqlj26ovaGpoDGx7E0res8dTvOC+Pu38tFtxs8UvNyQyM6eTm0d09Hl8pSIRsjwsmfIgFy29hplrvsuiUxYQtLxRi/eDC8ewq7aNn7+1hbQEN7Mma9/bqn/TO11KfR5jYOtCO+Eaehpc/L/g6PuPTFcIvrsimXjL8OtpLdo9hOrX2uKHsHTyA6S27GDGxp/an6MosRzC76+dwvTCdO5+aT3vb6uOWiyleoPe6VLqUIyBLa/D7n/D0NNhwpV2e64Y+NWGRDY1unj01EayvdpoWPV/lVmns37UXUze8Vvqk8exbfjNUYsV57J4/KZpXPvYMr713GoW3DqDkwrToxZPqaOhd7qUOpgxsOU1O+EqnGknXDG6vfR2uYeniuO5ZUQ75+fp4/Fq4Ngy/FZKc89n8vbfkFvzUVRjJcW5eGrudPJSvNzy1Eo2lmvnqap/0qRLqZ6Mgc2vwu7FMGwmjL+8zxOu5bvrWb67nje2tXD3ikRGxHdwTlLZgfXdk1L9mgjLJv6CpqQiTl/3nyS37opquMxEDwtum0FynIvrHl/G2rKGqMZT6svQpEupbsbApr/BniUw7AwY1/cJV7eukPCbknwcAt8dvg+nflJVP1dU9vJnpsKKv7Mn9yLEhDh32U2M3vVUVMuQn+rlpdtPIT3Bzdf/soKVe/TiRPUv+lWuFIAJwaaXofRDGH4WjLssZgkXwDN7s9nTEccdhZVkeQIxK4dSR6vLncr2oXNwBVoZXfYiVjC6T9/mp3r567xTyE72cNMTK/i4pC6q8ZQ6Epp0KRUKwvoXoPQjKDoHxl4a04RrUU0q/6xN45KcOqamtsasHEr1ljZvPiWDryChYx+nrr8HMdHtZy43JY4X551MfqqXm59cwaLN0R0TUqlIiYni47xfxrRp08yqVatiXQx1vAh0wV/Og8p1MOoiGHl+TBOuJVUublqSyuSUNr5fVI5Du4dQx5CcuuUU7n+H7UOvZfXYe6L+WWvrDPD0x3uoaOzgF5dN4PoZQ6MaTx2fRGS1MWZaJNtqlxHq+OX3wUs32gnXuFl2tWIMFTdbzF+WQn5cJ3cOq9CESx1zqjJm0OHNZezupwlYCawf/Z2oxkvwOLnt9OEs2VnDj1/bxP4mH/9x3ihEO7tTMaJJlzo+dbXBC9faTylOvMru/DSGGjqF2z5Kwe2AH4zYR7yl/XGpY9Pa0XfjDLQzftfj+J3xbCn6RlTjuZ0O/vz1qfzk9U38/r1iKpt83HfFRFyWtq5RfU+TLnX88TXBc1dD+Qq4/E/gj+2wOq1+4dalqVS0W7wws4FAiz+m5VEqqkRYOf4nOIMdTN7xO4KWl+2FN0Q1pNNycN8VE8lNieOhf+2kpqWTR64/kQSP/gSqvqX/cer40l4PCy6Hqk0w+0kYf9kXj70YZW0BYe6HKaxvcPLwyc1MzQywvCVmxVGqb4iDZRN/gTPYwdSt9wNENfF6fnkZANlJcVw+JZ831u3j/N8s5sZThpIU5/rUttfNKIhaOZSK6P6qiFwoIttFpFhEfniI9z0i8tfw+8tFpDC8vlBEOkRkXXj6U+8WX6kj0FwBT10M1VthzvN2whVD7QGY+2EKa+pd/G5GMxfmd8a0PEr1JeNwsnTSA5TlnMfUrfczcecjUR2nsdtJhencMGMo1S0+/rx4F7Wt+rlTfeewSZeIWMDDwEXAOOBaERl30Ga3Ag3GmBHAb4D7e7xXYoyZHJ5u76VyK3VkarbD4+dB4164/mUYdUFMi9MegFuWprKq1sVD05u5eLB+8avjT8hys3TyA5TkX8bE4j9y4tYH7D7zomzMoGRuO304Pn+QP35Qwu7atqjHVAoiq16cDhQbY3YBiMiLwCxgS49tZgE/C8+/AvxB9PEQ1V+ULYfnrwanB+a+DYNOiGlx9rU7mPdRClsbnfxmejOXDNGESx0/ispe/sy62pSJxHdUMqb0WTIaN7A772sYh5OSgquiVo4h6fF864winv54D08s3c2VJw5m8pDUqMVTCiKrXswH9vZYLg+vO+Q2xpgA0ARkhN8bJiJrReTfIvKVQwUQkXkiskpEVtXU1BzRH6DUF9r6JjxzKcRnwK2LYp5wrahxcem76ZS2Wjx2ahOzCjThUgoRynLPZ2/2mWQ1bWBs6QKcgeh3DJyR6OH2M4oYkhbPS6v28t62Kvpb35Xq2BJJ0nWoO1YH/1d+3jaVQIExZgrwH8DzIpL8mQ2NedQYM80YMy0rKyuCIil1GMbAkgfhrzdAzng74UorjGlxnt8Vx/WLU0l2hXj97AbOyeuKWXmU6ndEqMiayc7Bs4nvqGTCrsdJbd4W9bDxbie3nFbIlCGp/GtrNXe/vJ6ugHbZoqIjkqSrHBjSY3kwUPF524iIE0gB6o0xncaYOgBjzGqgBBh1tIVW6gv5ffDaN+Hde2HCbLj575CQGbPiVLY7+MZHKfxoTTKnZnfx+tkNjEiO7jAoSg1U9Snj2DJsLhg4f9mNDK34e9RjOi0Hs6cO5tyx2by6Zh83PrGcxna9KFK9L5KkayUwUkSGiYgbmAMsPGibhcBN4fnZwHvGGCMiWeGG+IjIcGAksKt3iq7UITRXwNNfgw1/hbN/Alc+Di5vTIoSMrCgxMt5i9L5sNrNjya28MTpTaS4tfpCqS/S7h3E5uG3UZ88htPW/5AZG36CM9Ae1ZgiwtljcnjomsmsKW3kij9+RGmdNrBXveuwDemNMQER+TbwDmABTxhjNovIvcAqY8xC4C/AAhEpBuqxEzOAmcC9IhIAgsDtxpj6aPwhSrHzX/DaPPtO19S5EJ8Jq5+KSVE+rnZx/6ZE1tW7+Ep2F788sZmCRK2yUCpSflci705/ggnFf2JCyaNkNaxj6eT7aUgZH9W4l03JJy/Vy7wFq7j8kY947MapTB2aHtWY6vihA16rgS8YgA9+Zbfhyh4PVz8Nez6MSVHW1Dl5cHMiS6vd5HqDfH9CG1cU+I5oXN/lu/W6RCngwNOL2XUrOWXDPXg769g8/FY2F80jZLmjErO7c9TdtW3MfXIFFU0+HrxqEpdMyotKPDXw6YDX6vhRVwKvz4e9y+DEG+HC+8Ed36dJlzGwpMrNYzu9LKnykOkJ8V+TWrh+eAdx1ifbaTKl1JdTnXES/zjtFbsT1ZI/U7D/nyyf+HNq0yZHLeawzARem38a8xas4s4X1lJW3878M4t0sGx1VPROlxpYuofsMSF7sOptfweHZTeYHxzRhUZEIkmQukJCpWMQj++MZ0ezk6y4ILeM6ODGER0kOD/7udKkS6mjl9Kyk2GVb+P2N1GdNpXy7LPYPvymw+/4JQWCIV5du491exs5sSCNyybn4bQcOlyQOkDvdKljW2uV3VC+fhdkj4OJV4O37zo1bPZbLKpJZVFNGk0BJ0O9PuYXVnNqWgsuh2HT3sMfQyn15TQljWRD/O0Mrv6A3PoVZDRvxjic7Cy4BuPo/Z80p+XgqqmDSU9w8962ampafFw/Y2ivx1HHB0261MDha4Ytb8Duf9u9y0+6DgafxBE1mDoKFT43b1Wls7guGb9xcGJKKxdn1zM+qb2viqCUAkKWh7JBF1CTNoWh+99h2tb/YcTeV1g3+ntUZH2l178TRIRzx+aQmxzHK6vLefiDYqYPT+fEgrRejaOOfVq9qPq/UBDWvwD/+jm01cCQGTDma+BJjFrInlWB21u9LNyfzuqmRJxiOCOjia/mNJAfp/34KBVzxtAZl8GUrb8mqaOcqrSprB/93ai199rf7OPZZaW0+gL89NJxXDe9QNt5HeeOpHpRky7Vf4VCsPUNeP8+qN0Og6dD4emQGv22FB/vqmdNUyIL96ezvS2eRCvIBdkNXJDVQIpLOzZVqr+RUJCsxjXk1yzGHWijIWkU5Vln0O4d9KntemM8x/auAIt31rJ4Rw1fO2EQ910xkaQ411EfVw1M2qZLDWyhIGx/Gz64H6o2QtYYuPoZGHMJrHk6qqE7g/BGWRwPbRlGhc9DltvPzUOqOCujkTirf12gKKU+YRwW1eknUZsyidz6ZQyqXcbElsfCyddM2r291+VDvNvJUzefxJ8Wl/Dgoh1sKG/iD9dN4YTBOmC2+mKadKn+o6sd1j8PHz9sN5JPHw5XPAYTrrSfUIyixi7huRIvT5V4qfFZFHp93DVsHyentWBpzYFSA0bIclORNZOq9Onk1K1gUP0yJu56nMbEIiozT7X7eOmF6kCHQ5h/5gimF6Zz1wtrufKPH/Gdc0byzTOKcFmRDPaijkdavahir67E7jl+7bPQUQ/5U+HUu2DsJZ9Ntrq7jOglpa0WT+z08tIeLx1BYWZOJ98Y1Y6rrUobxyt1DLCCneTUryC3bgWuYBt1yePYOnwue3POPaqnHXt2GdHY3sWPX9/E3zdUMiY3iQdmn6B3vY4j2qZL9X9d7XYV4pqn7f62xILRF8Epd0DBKZ9/JdpLSdfqOieP7YjnnX0enAKzCnzcNqqdMSl2ey3tU0upY4uEAmQ2riezaRPJ7aW0e7IpHjKb4iGz8cVl9UqMLRXNLFy/jxZfgNNGZHL2mGxuOX1Yrxxb9V/apkv1T8EA7P4ANrwM296CrlZIKbAHpp7ydUjKjWr4jgC8vS+O50q8rKl3keIK8a3R7dw0ooMcr46LqNSxzDic1KRPZfmkX5JXvZiRZX/lhOJHmFDyKOU5Z7Er/zIqM0/BOL58g/hxeckMz0rgH5v282FxLatLG/AHQ9x0aiFxrug2kVADg97pUtHl74CS9+0ka/s/7OpDTwqMuxROuBqGngZrnolaeGNgQ4OTl/bEsbAsjpaAg8LEADeP6OCqQt8he44HvdOl1LGq59OLSW2ljCh7iWH73iTO30CHO53SQV+ldNCF1KVOBPnybbMqGjtYtGU/O6payU2OY/5ZRVxx4mASPXqv41ij1Ysqtup3QfG7dlVg3U4IdoEzDnLGQ+4kyB4LVvQer+5OtN4u9/D2vjj2tll4HIaLB/u4epiPGZn+w7bX0qRLqWPTobqMkJCfvJoPGbbvTfKrP8Ayfjo8mZRnn0V5ztlUp08laHm/VLzhWQk88H/bWFPWSILbYtaUfK6bXsCE/JSj/VNUP6FJl+pbTeX2ANN7ltivDXvs9fEZkDUWcidAxgiIwhAd3Wp9wtJqN0uq7KnKZ+EUw+k5XXw1v5ML8jtJcRtNppRSX8gK+kht2Yk72Mqgmg9xBTsIiovatMnsz5hBdfo06lPGRZyEXTejAGMM6/Y28tzyMt5cX0FnIMSonETOHJ3NmaOymFaYjtupTzwOVJp0qej5+BForoDGUntqKAVfo/2eKx7SiyBzpH03K6F3GqcezB+C4maLtfUue6pzUdxiJ3SJVpAJyW1MSW5jWmoLiU5tq6WUOnIlBVfhCHaSU7+SnLrl5NYtJ615G4IhJBaNSaOoS5lIQ9IomhOH05RYRKcn/TPHOXhg7KZ2P6+tLeefW6tYsbsef9CQ4LaYkJ/C+LwUxuUlM3ZQEgXp8drh6gChSZc6eh2NdlcOdTuhrhiqt8L+jXai1S0+w+4dPrUQMkdA0qCjagNxsBa/UNpqUdpmsafFYkezk+3NTkqaLfzGrh9Md4eYkuEnS5qYmNzO8HgfDu3qQSl1lA5VDenpaiCzYT0ZTRvIaNxIRtMm3IHWA+/7rXg6PJl0eLLp8GTic6czZlgBeNPs8WKnzf3U8Vo7A3xUXMuSnbVs3NfEtv3N+PyfXCimeF0MTvMyKCWO9AQ36Qke0hNcpCd4yEhwk5bgPvCa4LZ0OKIY6fWkS0QuBH4LWMDjxpj/Oeh9D/AMMBWoA64xxuwJv3cPcCsQBO4yxrzzRbE06YoCYyDoh4APOhqgvQ7a6+1G7e119tRWCy2VdlVh417obPpkf7EgfRjkToRAFyTn2cmWJ+mIihE00NwlNHY5aAi/HpjvdFDb6WB/h4PKDgf7Oywauz6dwOV5g4xOCTA6JcDYlACT0gMMTQgiom2wlFIxYAzuQDPezhq8vlr7tbMGb2ctzpDv09u64iGjCFKGQMpg+3s0ISs8ZfL6Tj8drjSqOoSq5k4a2rpoaO+isd1PU4ef9q4AbZ1Bgp/zm+10CAkeJ/FuC6/LYlROEileFynxLlK8LpK99muixyLOZW/jDW/rdVnEuS3inBYuSzR5O0K9mnSJiAXsAM4DyoGVwLXGmC09tpkPnGCMuV1E5gCXG2OuEZFxwAvAdCAP+BcwyhjzuYPXxTzpMgZMKDz1nA9PmE+/FwrYCU2wy34N9Zg/gvWmbAWYICb0yUQoAN3zJmgPjxMKQkImhPyYQFf4OAH7NWQfW0J+JNiFhAJIyI/DBL74T0bwOZNpcWfR5M6lwZ1LgzOHKtcQ9rvyqXYOoiNkEQiGCNSX4Q+BPyQEQuA34deQ4DcQCK/vCkHAdM/bryE+/4PswJDmMQzyBsn1hhgUHyTPG8Lf3kSOp4scjx+vpVWFSqkBwBhcgVY8/kbGpwbsi92ORvAkHvrCtge/lYDPnYbflUSXMwm/K9F+dSbR5Uygw5FIm3HTEXLRFrRoDTlpDThpDlg0B5y0BBw0+y2MM446n9DoC9IRMBgcBHEQQjAIofB89zI9vp9dluCyHAcmtyW4nA7au4I4HYJDBMthT07HJ/OWQ7BEGJ2bhMty4LQEd/cxnA7iXA48Tuszr56Dl52OA8d2hI9pWeHXHnEc/aRao7f76ZoOFBtjdoUP/iIwC9jSY5tZwM/C868AfxA7VZ4FvGiM6QR2i0hx+HgfR1K4qHlwjP0BOFQyFSMCBIyDABZdOAlg4cdpT6bHPBb+mir8xhne1oUf74H9utd/avvw+i6cNJJIg0mi3iTRSCL1JolmEgjhQDA4xZ4sAUsMTgnhlHIsMfayo3t993YGt4DX8el13cf49LIhwRki0QqS6AyS5AwemI+3QoeuFvT09ZlQSqmjJILflYTflcTyLsAaBonh9zKnQyY4Ql24Am04A224gm24Au2fWu50p+MKtJLUVoYr0Ior0PqpqsyIOYnol747+TLhZIzuZCwAJgAvOq/nleQrCYYMwVAo/Grw+UMEjTmwHAwZKpo68AcN/kCIrqA9RaMlkwhYIoiAdCeN8kn6KAIjshN5686v9H7wLymSpCsf2NtjuRyY8XnbGGMCItIEZITXLzto3/yDA4jIPGBeeLFVRLZHVHoVC5lAbawLob4UPXcDl567gekYOm8PhaeBZTsgd32pXY/k3A2N9KCRJF2Huv9wcM76edtEsi/GmEeBRyMoi4oxEVkV6W1U1b/ouRu49NwNTHreBq5onbtIHjUrB4b0WB4MVHzeNiLiBFKA+gj3VUoppZQ65kWSdK0ERorIMBFxA3OAhQdtsxC4KTw/G3jP2C30FwJzRMQjIsOAkcCK3im6UkoppdTAcdjqxXAbrW8D72B3GfGEMWaziNwLrDLGLAT+AiwIN5Svx07MCG/3Enaj+wBwxxc9uagGBK0GHrj03A1ceu4GJj1vA1dUzl2/6xxVKaWUUupYpIM9KaWUUkr1AU26lFJKKaX6gCZdKiIicqGIbBeRYhH5YazLoyIjIkNE5H0R2Soim0XkO7EukzoyImKJyFoReSvWZVGRE5FUEXlFRLaFP3+nxLpMKjIi8r3w9+UmEXlBROJ669iadKnDCg8F9TBwETAOuDY8xJPq/wLA3caYscDJwB167gac7wBbY10IdcR+C/yfMWYMMAk9hwOCiOQDdwHTjDETsB8gnNNbx9ekS0XiwFBQxpguoHsoKNXPGWMqjTFrwvMt2F/8nxkVQvVPIjIYuBh4PNZlUZETkWRgJvaT/RhjuowxjbEtlToCTsAb7nc0nl7sX1STLhWJQw0FpT/cA4yIFAJTgOWxLYk6Ag8BPwB0tPeBZThQAzwZrhp+XEQSYl0odXjGmH3Ar4EyoBJoMsYs6q3ja9KlIhHRcE6q/xKRROBvwHeNMc2xLo86PBH5GlBtjFkd67KoI+YETgT+aIyZArQB2hZ2ABCRNOyanGFAHpAgIjf01vE16VKR0OGcBjARcWEnXM8ZY16NdXlUxE4DLhWRPdhV+meLyLOxLZKKUDlQbozpvqv8CnYSpvq/c4HdxpgaY4wfeBU4tbcOrkmXikQkQ0GpfkhEBLtdyVZjzP/GujwqcsaYe4wxg40xhdifufeMMb12xa2ixxizH9grIqPDq87BHplF9X9lwMkiEh/+/jyHXnwI4rDDACn1eUNBxbhYKjKnAV8HNorIuvC6Hxlj3o5hmZQ6HtwJPBe+UN0FzI1xeVQEjDHLReQVYA32099r6cUhgXQYIKWUUkqpPqDVi0oppZRSfUCTLqWUUkqpPqBJl1JKKaVUH9CkSymllFKqD2jSpZRSSinVBzTpUkoppZTqA5p0KaWUUkr1gf8P7G8EABRWXxcAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=[10, 4])\n",
"\n",
"sns.distplot([np.random.normal(mu_target, sigma_target) for _ in range(3000)], label=\"distribution $p(x)$\")\n",
"sns.distplot([np.random.normal(mu_appro, sigma_appro) for _ in range(3000)], label=\"distribution $q(x)$\")\n",
"\n",
"plt.title(\"Distributions\", size=16)\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": 178,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"simulate value 0.9542816022260111\n"
]
}
],
"source": [
"# value\n",
"s = 0\n",
"for i in range(n):\n",
" # draw a sample\n",
" x_i = np.random.normal(mu_target, sigma_target)\n",
" s += f_x(x_i)\n",
"print(\"simulate value\", s/n)"
]
},
{
"cell_type": "code",
"execution_count": 172,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"average 0.9495227171370471 variance 0.3043862985463373\n"
]
}
],
"source": [
"# calculate value sampling from a different distribution\n",
"\n",
"value_list = []\n",
"for i in range(n):\n",
" # sample from different distribution\n",
" x_i = np.random.normal(mu_appro, sigma_appro)\n",
" value = f_x(x_i)*(p_x.pdf(x_i) / q_x.pdf(x_i))\n",
" \n",
" value_list.append(value)\n",
"\n",
"print(\"average {} variance {}\".format(np.mean(value_list), np.var(value_list)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Different $q(x)$"
]
},
{
"cell_type": "code",
"execution_count": 179,
"metadata": {},
"outputs": [],
"source": [
"# pre-setting\n",
"n = 5000\n",
"\n",
"mu_target = 3.5\n",
"sigma_target = 1\n",
"mu_appro = 1\n",
"sigma_appro = 1\n",
"\n",
"p_x = distribution(mu_target, sigma_target)\n",
"q_x = distribution(mu_appro, sigma_appro)"
]
},
{
"cell_type": "code",
"execution_count": 182,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/jeremy.zhang/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
" warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n",
"/Users/jeremy.zhang/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.\n",
" warnings.warn(\"The 'normed' kwarg is deprecated, and has been \"\n"
]
},
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x117564ac8>"
]
},
"execution_count": 182,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAEKCAYAAADQARsOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl4VdW9//H3OudknhMyQEIGIEAgzGEeFZzAKlLqgIjWsWLrbb1atb9q1dZOKvX2Vq0jWgeg1qIItl4QRRAISZinMIRMZCTzPJys3x8HaICQnCQn2TnJ9/U8eSBn7732Jwkk36y19lpKa40QQgghhOg8k9EBhBBCCCF6CymshBBCCCEcRAorIYQQQggHkcJKCCGEEMJBpLASQgghhHAQKayEEEIIIRxECishxHlKqbuUUrrZW5VSKl0ptVYpdbNSytTs3Oiz59zVjvbnKKWead5OOzJFN3stXSn1gb1tdDRXRz5GIUTfJoWVEKIlPwCmAvOBp4A6YBXwf0opj7Pn5J49Z0M72p0D/Ir2fe/ZcPY+ue24pr3m0HKujnyMQog+zGJ0ACFEj7RXa32i2fvvK6U+Bj4G/gj8RGtdB+zsqgBKKRegUWtdCBR21X1a09UfoxCi95EeKyGEXbTWnwCfAfcppTxbGiZTSk1USm1UShUppaqVUmlKqVfPHnsGW68QQMO54cazx861tVwp9UelVA62XjL/loYCm93vPqXUCaVUrVJqt1LqiouOf6OU+qaF69KVUu+2I9ddF12/VCm17+x9zyil3ldK9W/hHh8opW5VSh05O6yarJSacdF5l/2cCSGcj/RYCSHa4wtgIZAAZDY/oJTyBr4EdgF3ARVANDDt7ClvARHAPcAMwNpC+/8PSALuB8xAbStZZgMTzl5TBzwO/EspNUZrndqOj8meXOcppe4HXgfWAE8CA4DfApOVUuO11pXNTp8JDMM2nFoL/BpYr5SK1lqX2vE5E0I4GSmshBDtca6Y6s9FhRUwHAgAfq613t/s9XcBtNbZSqnss68laq0bW2g/H7hJN9vEVCl1uSyhwHStdebZ874CMoBfAnfY+wHZmetcFjO24ugbrfWtzV4/CmwF7gb+3OwSX2Cs1rrk7Hl52ArH+cBHtPE5E0I4HxkKFEK0x7kqp6Xd248DpcDrZ4fKBnag/U+1/TvD7zxXVAForSv4z0T3rjIMCAE+bP6i1nobtqJu9kXn7zhXVJ114OyfkWf/dMTnTAjRg0hhJYRoj3M/+C95Qk9rXQZcAeQArwKZSqmDSqnvt6P99jz5l3+Z18Lb0UZ7BZ79s6Wcec2On1Pc/J2zk+EB3M++74jPmRCiB5HCSgjRHguwzRVKaemg1nqv1vr72AqMqcBJ4O9KqXg727e3twpsQ4EtvXa62fu1gGsL511cANnrXKEU1sKxMKCovQ064HMmhOhBpLASQthFKbUIuAH4q9a6urVztdaNWuud2CZtm4C4s4fO9dh4tHhh+0xpPnSmlPLBVvjtaHZOBjBUKeXa7LxZgM9FbdmbKxVbr9itzV9USk0DooAt7fkAmmvlcyaEcCIyeV0I0ZKxSql+2Hp7IoHrsS0auhHbk3CXUEpdj+1pvk+BU4AX8DC2J93OFTuHz/7530qpfwFWrXVyBzPmY1uw9Bn+81SgF7bJ5eesPpvpnbPLK8QAjwBlF7VlVy6ttVUp9TS2OVEfAB9gG3p8Htt8qZXt+QDs/JwJIZyIFFZCiJZ8fPbPWqAA2I2tl+YfrUwuPw7UYOtx6Y+tOEgCrtJan3vqbj22uUTLgaexTYa/7GN/bdgCfINtqYMIbMXRdVrrY+dO0Fp/rZT6EfAo8H1gD7AU+OSituzOpbV+QylVDTyGbV2vSmzLUPz8oqUW7GHP50wI4USU/Q/gCCGEEEKI1sgcKyGEEEIIB5HCSgghhBDCQaSwEkIIIYRwECmshBBCCCEcxLCnAvv166ejo6ONur0QQgghhN1SUlLOaK2D2zrPsMIqOjqa5OSOLl8jhBBCCNF9lFIZ9pwnQ4FCCCGEEA4ihZUQQgghhINIYSWEEEII4SCypY0QQgjRCQ0NDWRnZ1NbW2t0FOEA7u7uRERE4OLi0qHrpbASQgghOiE7OxsfHx+io6NRqqNbX4qeQGtNUVER2dnZxMTEdKgNGQoUQgghOqG2tpagoCApqnoBpRRBQUGd6n2UwkoIIYToJCmqeo/Ofi2lsBJCCCGEcBAprIQQQgghHEQmrwshRC/xUWJmm+csmRzZDUmE6LuksBJCCCEcyJ4Ctz06Ugw/88wzeHt78+ijjzJt2jS2b9/e4nmlpaV89NFHLF++/LJtnbs+PT2d66+/noMHD9qVoaW2W8viKDU1NVx77bVs3rwZs9nc4jn19fXMmzePzZs3Y7E4thSSoUAhhBCiF2utkCktLeXVV19t8ZjWmqampg4XQi213dVFFcA777zDokWLLltUAbi6ujJ37lzWrFnj8PtLYSWEEEL0As8//zzDhg1j3rx5pKamnn/d29sbgKqqKhYsWMCYMWOIj49nzZo1PPHEE5w8eZKxY8fy2GOPkZ6eTlxcHMuXL2f8+PFkZWWdvx6gsbGRO++8k9GjR7N48WKqq6tJT08nPj7+/DkvvvgizzzzzCVtN88CsGLFCuLj44mPj+fll18GOH//++67j5EjR3L11VdTU1PT4sd76623cssttzB58mSioqLYsGEDAB9++CE33njj+fOuuOIKNm7cCMAvf/lLHn74YQAWLlzIhx9+2PFP+GXIUKAQQgjh5FJSUli9ejV79uyhsbGR8ePHM2HChAvO+fe//82AAQPOFyBlZWVMnjyZgwcPsnfvXsBW2KSmprJy5coWe7JSU1N5++23mT59OnfffTevvvoqixcvbjHT73//+wvavjjvypUrSUxMRGvN5MmTmT17NgEBARw/fpxVq1bx5ptvcvPNN/PJJ5+wdOnSS9rYt28fCxcuZM2aNWzbto1HHnmEq666irS0NKKjo8+f9+yzz/L0009TUFDAnj17WLduHQDx8fEkJSXZ9wluB+mxEkIIIZzc1q1buemmm/D09MTX15cbbrjhknNGjRrFpk2bePzxx9m6dSt+fn4tthUVFcWUKVNaPDZw4ECmT58OwNKlS9m2bVuH8m7bto2bbroJLy8vvL29WbRoEVu3bgUgJiaGsWPHAjBhwgTS09Mvub6mpoYzZ87wq1/9CoARI0ZQUlLCmTNn8Pf3v+DcWbNmobVmxYoVrF69+vwQodlsxtXVlYqKig59DJcjhZUQQgjRC7S1sOXQoUNJSUlh1KhRPPnkkzz33HMtnufl5WX3PZRSWCwWmpqazr9mz6rlWuvLHnNzczv/d7PZTGNj4yXnHDx4kNjYWNzd3QHYvXs3Y8aMwcPD45L7HzhwgNzcXNzc3PDx8bngWF1d3fk2HEUKKyGEEMLJzZo1i7Vr11JTU0NFRQWff/75Jefk5OTg6enJ0qVLefTRR9m9ezc+Pj7t6rHJzMxkx44dAKxatYoZM2YQGhpKQUEBRUVF1NXVsX79eoBW2541axaffvop1dXVVFVVsXbtWmbOnGl3jn379pGZmUltbS1VVVX86le/4mc/+xkBAQFYrdbzxVVubi633347n332GV5eXnz55Zfn2ygqKiI4OLjDmy1fjsyxEkIIIRzIiLXCxo8fzy233MLYsWOJiopqsUg5cOAAjz32GCaTCRcXF1577TWCgoKYPn068fHxXHfddTz00EOt3icuLo733nuPBx54gNjYWB588EFcXFx4+umnmTx5MjExMQwfPhzgkrZfeOGFC/LeddddTJo0CYB7772XcePGtTjs15J9+/Zx++23M2fOHMrLy/nFL35xfojy6quvZtu2bUybNo1Fixbx0ksvERcXx1NPPcXjjz/ONddcA8DXX3/N/Pnz7bpfe6jWuuO6UkJCgk5OTjbk3kII0RvJAqHGOHLkCHFxcUbH6FNmzZrFm2++ybBhwy45tmfPHlasWMH777/fahuLFi3id7/7XYtttPQ1VUqlaK0T2somQ4FCCCGEcConT54kNja2xWPjxo3jiiuuwGq1Xvb6+vp6Fi5c2GJR1VkyFCiEEH1IW71a0qMlnMHp06dbPX733Xe3etzV1ZVly5Y5MtJ50mMlhBBCCOEgUlgJIYQQQjiIFFZCCCGEEA4ihZUQQgghhINIYSWEEEII4SBSWAkhhBBCOIhdhZVS6lqlVKpS6oRS6olWzluslNJKqTYX0BJCCCGE6G3aLKyUUmbgFeA6YARwm1JqRAvn+QAPA4mODimEEEII4QzsWSB0EnBCa50GoJRaDdwIHL7ovF8DfwQedWhCIYQQwpkkr3Rsewk/bPclzzzzDN7e3jz66KNMmzaN7du3t3heaWkpH330EcuXL79sW+euT09P5/rrr+fgwYN2ZWip7dayOEpNTQ3XXnstmzdvxmw2t3hOfX098+bNY/PmzVgsjl0r3Z6hwHAgq9n72WdfO08pNQ4YqLVe31pDSqn7lVLJSqnkwsLCdocVQghhv/QzVfx1y0le+PIoL/1fKn/aeIxP956msq7R6GiiG7VWyJSWlvLqq6+2eExrTVNTU4cLoZba7uqiCuCdd95h0aJFly2qwLby+ty5c1mzZo3D729PYaVaeO38zs1KKRPwJ+C/22pIa/2G1jpBa50QHBxsf0ohhBB2q6m3snbPad7YmkZ5TQPRQV5EBHjQz9uV5PRiVmxMZfvJM1ibdNuNCafx/PPPM2zYMObNm0dqaur51729vQGoqqpiwYIFjBkzhvj4eNasWcMTTzzByZMnGTt2LI899hjp6enExcWxfPlyxo8fT1ZW1vnrARobG7nzzjsZPXo0ixcvprq6mvT0dOLj48+f8+KLL/LMM89c0nbzLAArVqwgPj6e+Ph4Xn75ZYDz97/vvvsYOXIkV199NTU1NS1+vMePH2fOnDkkJCTw85//nCFDhgDw4YcfcuONN54/74orrmDjxo0A/PKXv+Thhx8GYOHChXz44Ycd/4Rfhj39X9nAwGbvRwA5zd73AeKBb5RSAGHAOqXUDVrrZEcFFUII0baiyjre2JpGZW0jM4b0Y25cCG6W//zmXlBey/oDuazfn8u+rFLumhaDh+vlf7MXziElJYXVq1ezZ88eGhsbGT9+PBMmTLjgnH//+98MGDCADRs2AFBWVsbkyZM5ePAge/fuBWyFTWpqKitXrmyxJys1NZW3336b6dOnc/fdd/Pqq6+yePHiFjP9/ve/v6Dti/OuXLmSxMREtNZMnjyZ2bNnExAQwPHjx1m1ahVvvvkmN998M5988glLly694Hqr1cqyZct45ZVXGD9+PD/5yU8YOXIk9fX1pKWlER0dff7cZ599lqeffpqCggL27NnDunXrAIiPjycpKcn+T7Kd7OmxSgJilVIxSilX4FZg3bmDWusyrXU/rXW01joa2AlIUSWEEN2spt7K33ZkYG3SLJ8zhPmj+l9QVAGE+Lrzw2nR3JIwkJzSWt7+Lo1qGRp0elu3buWmm27C09MTX19fbrjhhkvOGTVqFJs2beLxxx9n69at+Pn5tdhWVFQUU6ZMafHYwIEDmT59OgBLly5l27ZtHcq7bds2brrpJry8vPD29mbRokVs3boVgJiYGMaOHQvAhAkTSE9Pv+T6Tz/9lBEjRjB+/HgA4uLiGD16NGfOnMHf3/+Cc2fNmoXWmhUrVrB69erzQ4RmsxlXV1cqKio69DFcTpuFlda6Efgx8CVwBPi71vqQUuo5pdSlXzkhhBDdrtHaxOqkTIqr6rl9chThAR6XPVcpxZiB/iydEklBeR1vf3dK5l31AmdHjS5r6NChpKSkMGrUKJ588kmee+65Fs/z8vKy+x5KKSwWC01NTedfq62tbTOr1pcfhnZzczv/d7PZTGPjpf829+zZc774Ati3bx9jxozBw8PjkvsfOHCA3Nxc3Nzc8PHxueBYXV0d7u7ubeZtD7vWsdJaf6G1Hqq1Hqy1fv7sa09rrde1cO4c6a0SQoju9ZsNRzheUMmNYwcQ0+/yPxibGxbmyx1ToyisqOOdbaeobbB2cUrRVWbNmsXatWupqamhoqKCzz///JJzcnJy8PT0ZOnSpTz66KPs3r0bHx+fdvXYZGZmsmPHDgBWrVrFjBkzCA0NpaCggKKiIurq6li/3vYcW2ttz5o1i08//ZTq6mqqqqpYu3YtM2fOtDtHUFAQR48eBSAxMZG//e1vjB49moCAAKxW6/niKjc3l9tvv53PPvsMLy8vvvzyy/NtFBUVERwcjIuLi933tYdjnzEUQgjR7Tbsz+Xd7enMGNKPhOjAdl0bG+LDsqnRvLv9FB8mZrB0ShSuFtmUo1M6sDxCZ40fP55bbrmFsWPHEhUV1WKRcuDAAR577DFMJhMuLi689tprBAUFMX36dOLj47nuuut46KGHWr1PXFwc7733Hg888ACxsbE8+OCDuLi48PTTTzN58mRiYmIYPnw4wCVtv/DCCxfkveuuu5g0aRIA9957L+PGjWtx2K8ld9xxBwsWLGDUqFHMnz+foKCg85PXr776arZt28a0adNYtGgRL730EnFxcTz11FM8/vjjXHPNNQB8/fXXzJ8/3677tYdqrTuuKyUkJOjkZOnYEqLD2lorx4Bv7qL7ldU0MG/FFsJ83bk5YSBmU+vDQZezO6OEf+zO5vvjI3jxB6PbHFYS/3HkyBHi4uKMjtFnZWVlsXjxYhITbeuT79mzhxUrVvD++++3et2iRYv43e9+x7Bhwy451tLXVCmVorVuc2cZ6bESQnQtexZLlCKww/7w76MUVdax8q6J7M8u63A746MCKKmu55Pd2UQGevJf82IdmFKIrrNv3z5Gjx59/v1x48ZxxRVXYLVaW10gdOHChS0WVZ0lhZUQfZkjVoiWosgwSenFfJSYyb0zYogP9+tUYQVw5fAQ/D1d+dOmY0QEePD9CREOSipE17n++uu5/vrrL3jt7rvvbvUaV1dXli1b1iV5pLASQthPa7DWgTKD2bETPkX71Dc28Yt/HiDc34OfXTXUIW0qpfjdolHkltXw+Cf76e/nzrQh/RzSthB9hRRWQohLWRugLNv2Vn4aKnKgthzqK6Dp7JNjJhdw9YRD/4TBV9reQkeBSSY+d4d3vjvF8YJK3rkrAS83x30rd7WYeG3pBH7w1+088EEKnzw4jaGhPm1fKIQApLASQgA01kLxKShOg+KTUJrxnwLKxQv8wiE4DFy9bW/aCvU10FAF1cWw6RnbW0A0zPkFjPqBFFhdqLCijr9sPsHc4SFcOTzU4e37ebiw8oeTWPjKd/xwZRJrl08jxNexa/30NlprmfDfS3T2oT4prIToi2pKIGMHHP4Uik7aeqbQoEzgGwHRsyAwBvwiwd0PWvuBkfBDqMiDk5th52uw9n7Y/meY9wzEXtVNH1DfsmLjMWobrPxiQdc9iRbu78HKuyZy8+s7uPu9JNbcP9WhPWO9ibu7O0VFRQQFBUlx5eS01hQVFXVq0VD5XyJEX1BTCunbIOM7SN8KeQcBDSYL+EfZCqDAwbYeJ4tbW61dyicMxi6B0bfC4bWw+Tfw4WKY+d/gG24r2IRDHM0rZ01SJsumRjM42LvtCzohPtyPvywZx73vJfOTVXt4444JWMzytbxYREQE2dnZFBYWGh1FOIC7uzsRER1/cEMKKyF6K60hOxmS37HNg2qsBYs7REyEOU9C9HTIP+zYSegmE8R/H4Z/D754FLa+BGGjYeztHSvYxAW01vxm/RF83F34aTcth3Dl8FCevTGepz49yLOfH+a5G0dKr8xFXFxciImJMTqG6CGksBKiNyrPgbfmwukU25yosUsgfjFEJFxY4Jw53vl7XW7Jhv5jYUS5bbhx+59h8oPg1rU9LL3d5qMFbDtxhqevH4G/p2u33feOKVFkF1fz+rdpDAz04P5Zg7vt3kI4GymshOhNdBOkbYHU9eARCAtWwOibwc2Ap7qUgkGzwTvYVnwlvw1TlssyDR3UYG3i+S+OMKifF3dMjer2+z9+7XCyS2r47RdHCff3ZMHo/t2eQQhnIIPlQvQWDTW2yeNHPoPgEbB8J0y8x5iiqrmQEbahwJJTsPdDW/En2u2DnRmkFVbxi/lxuBgwz8lkUrx08xgmRAXws7/vJSWjuNszCOEMpLASojewNkDSW7blEkbfCgl3g1cPWthxwFjbvKvcvXB0g9FpnE5pdT0vbzrO9CFBzI0LMSyHu4uZN5clEO7vwb3vJXPqTJVhWYToqWQoUAhnp5tgzwe29afG3QHhE4xO1LLBV0J1EZz8yvakYPh4oxM5jf/56jgVtQ38csEIwyeOB3q5smhcOK9tOcni17bzo9mDL1mGYcnkSIPSCWE86bESwplpDQf/CXn7YMTCnltUgW3OVfz3bcs7HPwH1HZuX7u+4mRhJe/vyOCWiQOJ6+9rdBwAgrzduGNKFGU1Dby/M4MGqwzvCnGOFFZCOLPsXZCxDQZdAYPmGJ2mbSYzjF1qG7rct8pWGIpWPb/hCO4uZh65apjRUS4QFeTFDxIGklVczccp2TTJ11IIQIYChXBedRVw+DMIHARx3zM6jf28gyHuBjj0CWTugKhpRifqET5KzLzktSO55Ww+WsB18WFsPJzf44bYRoX7URofxr8O5rHRy5VrRoYZHUkIw0mPlRDO6tCn0FgHo252vpXNo6dDv6G2Na6qzhidpkdqsDaxfn8OIT5uTBvcgx5EuMiMIf2YGB3AlmOFHM4pNzqOEIaTHishnNGJTZCTArHX2LaTcTbKBGNugy2/h4OfwOzHjE7U42w5VkhJdQP3zIjBbOq+Cest9Zy1RinF9aMHkFNay8cpWTzkO6SLkgnhHKSwEsLZ1FfD+kfAKwSGtLLJ8eVWRO8pPAJg6LW24cxjX8LQa4xO1GMUVdbx7bFCRkf4dfl+gI7gYjaxZHIkf9l8go8SM7l3ZgyervLjRfRNTjZ+IIRg+/9CaYZtRXWzk//wip4J3iHw7ydsw5oCrTXr9+diMinmxzvP6uYBnq7cOnEg+eW1PP3ZIaPjCGEYJ/+uLEQfU1sGO1+B4ddDUC8YcjFZYMRNsOt126rxM35qdCLDHc2rIDW/gvnxYfh6XLj9T3uH6bpbbKgPc4aF8I+UbK4eEcrVMpld9EHSYyWEM9n1hq24mtWL5iSFxMGw+fDtC1CRZ3QaQzWfsD61B09Yb80Vw4MZ0d+XX6w9QHFVvdFxhOh2UlgJ4SzqKmDHK7Z5SQPGGp3Gsa55Hqz18NVzRicx1LkJ6zeMGdCtE9YdyWIyseKWMZTVNPDUZweNjiNEt5OhQCGcRdLbUFMCs35udBLHS9sCUdNh70fgFwE+F80tSvihMbm60bkJ62Mi/BjkBBPWW7M7o5QrhoWwYX8uPm77GR3hf8k5PW1NLiEcRXqshHAG9VW2SeuD50JED962pjOGzAOLGxz9wugk3U5rzef7czCbFNc50YT11syMDSYiwIN1+3KoqbcaHUeIbiOFlRDOIOVdqD4Dsx83OknXcfWybdScfwBK0o1O0602Hs7nWH4lc+NCL5mw7qzMJsXCseHU1FvZdDTf6DhCdBsZChSiJ2q+BpVugm9fhMDBUHDY9tZbxcyGU9/Ckc9h6o9tGzf3cjX1Vp79/DChvm5MHRRkdByHGuDvwaSYQBLTipgYHUiYr7vRkYToctJjJURPV3AEaoptaz71dhY322ryxSeh8KjRabrFa9+c4HRpDTeMCXfaCeutuSouFDeLmfX7c9CyUbPoA6SwEqKnS98Kbn4QNsroJN0jaip4BELqBujlP4jTz1Tx1y1pLBw7gJh+XkbH6RKebhauGhFKWmEVh2QvQdEHSGElRE9WVWjruYmaCiaz0Wm6h8li296mLBsKeu8K3lprnvn8EK4WE7+YH2d0nC51bhjwiwO5NFibjI4jRJeSwkqInizjO9uGxZFTjU7SvcITwDPItodgL+21+jq1gG9SC/npvFhCevncI7NJMX9Uf0prGkhKLzY6jhBdSgorIXoqaz1kJULYaHD3MzpN9zKZbRtMl2X1ysn6DdYmfrPhCIOCvbhzWrTRcbrF4GAvYvp58U1qIfWN0mslei8prIToqU7vhoYaiJ5hdBJjREy0zbXqhb1WH+zMIK2wiv83Pw4Xc9/4NqyU4uoRoVTWNbIzrcjoOEJ0GVluQYieKmM7+ITZllnoi0xmiL0K9q+B4xth6NVGJ+qUcxsoV9c38tL/HWNIsDd5ZbU9fmNlR4oK8mJoqDdbjhVSUduAj3vvWLNLiOb6xq9KQjibijwoy4SBk/vEWk6XFTERPAJgy+97Ta/VV0cLqG2wMn9Uf1Qf/NpeFRdGTYOVd7alGx1FiC5hV2GllLpWKZWqlDqhlHqiheM/UkodUErtVUptU0qNcHxUIfqQ08mAggG9dPsae5kstrlWp1PgxFdGp+m0MxV1/1ks0693T1i/nPAAD0b09+WtrWmUVTcYHUcIh2uzsFJKmYFXgOuAEcBtLRROH2mtR2mtxwJ/BFY4PKkQfUVTE2QnQ/BwcPc1Oo3xBk4Cv4G9otfqq6P5mE2KuXEhRkcx1Ny4ECrqGnl/Z7rRUYRwOHt6rCYBJ7TWaVrremA1cGPzE7TWzVd98wKc+7ufEEbK2Aa1pRCRYHSSnsFkgRk/g+wkOLnZ6DQdlldey/7sMqYN7tfn5xb19/NgzrBgVn6XTm2DbNAsehd7Jq+HA1nN3s8GJl98klLqIeARwBW40iHphOiL9q2xbe3SV1Zat8e4pbD1Jfjm97aNmnvY3CR7JqB/dSQfV4uJmUP6dUOinu9Hswdz6xs7+Tg5izumRhsdRwiHsaewauk72CU9UlrrV4BXlFJLgF8Cd17SkFL3A/cDREZGti+pEL1F8w2WL2athwMfw4AxYHbtvkw9ncXN1mv1xaOQ9rWtuHIiOaU1HMop58rhIXi6ycPYAJNjAhkX6c8bW9O4bVIklj6y7ITo/ez5l5wNDGz2fgSQ08r5q4GFLR3QWr+htU7QWicEBwfbn1KIviLvAFjrIHyi0Ul6nvHLwGcAfPMHp5trtelIPh4uZmZIb9V5SimDl8NXAAAgAElEQVR+NHswWcU1bDiQa3QcIRzGnsIqCYhVSsUopVyBW4F1zU9QSsU2e3cBcNxxEYXoQ7KTbMsLBPXRtataY3GDmY9A1k5I32Z0GrtlFVdzNK+CmbH9cHfpI/s92umquFAGB3vx1y1paCcrloW4nDYLK611I/Bj4EvgCPB3rfUhpdRzSqkbzp72Y6XUIaXUXmzzrC4ZBhRCtKG2HApTbfvkKRkWadG4peAVDNuc58HjLccK8XAxM3VwkNFRehyTSfHA7MEcyS1n6/EzRscRwiHsGuzXWn8BfHHRa083+/t/OTiXEH1PTgqg5WnA1rh4wJTl8NWzkLMHBowzOlGrzlTUcSS3nDnDgnGzSG9Vc+cm/Ddam/B2s/D8hiNkl9ScP75ksszDFc5Jfi0WoqfITgK/SPAONTpJzzbxHnDzhW0vG52kTVtPnMFsUkwZJL1Vl2Mxm5gUE0hqfgVnKuuMjiNEp0lhJURPUJ5je4uQSettcvezFVeHP4MzJ4xOc1kVtQ3sySxhfGRAn1+3qi2TYwIxK8UO2ZxZ9AJSWAnRE2Qn2eZVhffsoa0eY8py22T273pur9WOk0VYmzQzYuVJwLb4uLswKsKP3RklsmCocHpSWAlhNN1k2wsvZAS4ehudxjl4h9gmsu9bDWWnjU5zibpGKztPFTFigC/9vN2MjuMUpg0Ooq6xid2ZJUZHEaJTpLASwmhnjkFdue1pQGG/aQ/bitIdrxid5BIpGSXUNjQxK1bW67NXRIAnkYGe7DhZRJMsvSCcmBRWQhgtO9n2tFvoSKOTOJeAKBi1GFLehepio9Ocp7Vm16liIgI8GBjoaXQcpzJtcBBFVfUcy68wOooQHSaFlRBGaqyDvP3QfyyYZYJzu03/KTRUwa43jE5yXmZxNQUVdUyKDjQ6itMZOcAPH3cLu071nEJZiPaSwkoII+Xts+0PKE8DdkzoCBh6HST+FeoqjU4DQFJ6MW4WE6Mi/IyO4nTMJsWEqABS8yrIKa1p+wIheiDZDVQII2Ung2cQBMQYnaRna23j6uChcOxfsPs9mPpQ92VqQU29lf3ZZYyPDJAFQTtoYnQgW1ILWZ2UxSNXDTU6jhDtJoWVEEapKYUzxyH2alDK6DTOKyAGAgfDlj+AyRXMl/m2lvDDLo+yN6uExibNxBgZBuyoAE9XYkO9WZOUycNXDsFiloEV4VzkX6wQRjmdjG0LGxkG7LQh86C27Ozn1Bhaa5LSSwj39yDc38OwHL3BpOgg8svr2Hy0wOgoQrSbFFZCGEFr2zBgQAx4yQKSnRY8HHz6w6ktts+tAbJLasgrr2WiTFrvtGFhPoT5uvPRrkyjowjRblJYCWGE8myozJMNlx1FKYiZDRW5UHTckAi70otxNZsYI5PWO81sUtw8cSBbjhWSVVxtdBwh2kXmWAlhhOwkMJmhv2xh4zDhE+Doekj7Bvp176Tn2gYr+7NLGRPhj5uLTFp3BHeLCTQ89dlBrh4R1uI5SyZHdnMqIdomPVZCdLcmK5zeDSHx4CoLSDqM2QWipkPBYajs3rk5+7JLabBqJsmkdYfx93RlWJgPKeklWJtkJXbhPKSwEqK7FR6F+koZBuwK0TNsPYGnvu22W55bab2/n7tMWnewSdGBVNQ1ciS33OgoQthNCishutvpZHDxgpA4o5P0Pm4+MGACZO+C+u6Zm3PgdBm5ZbZJ60qWzXCooWE++Hm4kJQuK7EL5yGFlRDdqbYM8g7AgHFgkimOXWLQbNtq9pk7uuV2q3Zl4mJWjB3o3y3360tMSpEQHcDxgkqKq+qNjiOEXaSwEqI7Hf4Mmhpl7aqu5BsOQbGQvtU2n60LVdY18tneHEaH++Muk9a7REJUICaF9FoJpyGFlRDdad9q8AoBf3maqUsNmg21pba9GLvQ5/tyqK63ykrrXcjPw4VhYb4kZ5TQ2NRkdBwh2iSFlRDdpSQDMr6zTVqXuThdK2QEePaDtC1deptVuzIZHubDwACZtN6VJkUHUlXXyJHcCqOjCNEmKayE6C77/277M1yeBuxyymRbMLQ0A0rSu+QWB0+XsT+7jNsmRcqk9S4WG+qNv6cLSadkOFD0fFJYCdEdtIZ9qyBqBnjKsFG3GDgJLO62BUO7wKpdmbhZTCwcG94l7Yv/MClFQlQAJwplErvo+aSwEqI7nE6B4pMw5lajk/QdFjeInAp5+6GmxKFNV9fbJq0vGN0fP08Xh7YtWjYhKhAFJMskdtHDSWElRHfY+yFYPGDEDUYn6VtiZtn+dPCCoev35VJZ18iSSfIQQnexTWL3ISVTVmIXPZsUVkJ0tYYaOPCJrahylw16u5VHAISNgqxEhy4Y+tGuTGJDvJkQFeCwNkXbEqICqaht5Fi+TGIXPZcUVkJ0tSProa4Mxt5udJK+KXoWNFTDgY8d0tzhnHL2ZpXKpHUDDAvzwcfdImtaiR5NCishutreD2zrVkXPNDpJ3xQ4CHz6w643bQ8RdNLqpExcLSYWjZdJ693NbFJMiAwgNa+CspoGo+MI0SIprIToSqVZtrWUxiwBk/x3M4RStl6r/AOQubNTTdXUW1m75zTz48Pw93R1UEDRHgnRgWggJUN6rUTPJN/phehK+1YBGsYuMTpJ3xY+3ja/bdcbnWpmw4FcKmobuU0mrRsm0MuVwcFeJGeU0CST2EUPJIWVEF2lqQn2fGB7Mi0gyug0fZvFDcbdAUfWQXluh5tZtSuTQcFeTJItbAw1MTqQ0uoGtp04Y3QUIS4hhZUQXSXjO9vK32OXGp1EAEy8x7Ypc8rKDl1+LL+ClIwSlsikdcON6O+Lp6uZ1UmZRkcR4hJSWAnRVfa8D26+EPc9o5MIsE1ij70akldCY/tX7161KxNXs4lF4yO6IJxoD4vZxPjIADYezudMZZ3RcYS4gBRWQnSF6mI49CmMvgVcPY1OI86ZdD9UFdiGBNuhtsHKP3ef5pr4MAK9ZNJ6T5AQFUCDVfNJSrbRUYS4gBRWQnSFvR+BtQ4Sfmh0EtHc4CshcDAkvt6uy/51MJeymgZumzSwi4KJ9grxdSchKoA1SVloByyjIYSjSGElhKNpDSnvQsQkCB1pdBrRnMkEk+6D7F2Qs8fuy1YlZhEd5MnUQUFdGE60162TIkk7U8WuU7L0gug5LEYHEKLXSd8GRcdh4WtGJxEtGXMbfPVr2PUWLHylzdNT8yrYlV7MtSPDWLUrqxsCCnstGNWfZ9cdYnVSFpOl6BU9hPRYCeFoKSttayaNvMnoJKIlHv4w5hbbFjfVbfd0vL8zHYtJyb6APZCHq5mF48LZcCCXkqr2P5AgRFewq7BSSl2rlEpVSp1QSj3RwvFHlFKHlVL7lVJfKaVk0R7RN1UWwuF1tpXWXTyMTiMuZ+J9tjlwu//W6mnltQ38c/dpRkf44+UmHfw90e1TIqlvbOIfMold9BBtfqdQSpmBV4CrgGwgSSm1Tmt9uNlpe4AErXW1UupB4I/ALV0RWIgebe+H0NQgk9Z7ouSL1q8KioXvXgY3H1Bnf8e86Ov2z5RsquutTBkkC4L2VMPDfJkYHcAHiRncMyMGk0nWGBPGsqfHahJwQmudprWuB1YDNzY/QWv9tda6+uy7OwFZ6EX0PU1WSH4boqZD8DCj04i2RM+AmhLIP9TiYa017+/MYMxAfyICZMmMnmzplCgyiqplJXbRI9hTWIUDzWdsZp997XLuAf7V0gGl1P1KqWSlVHJhYaH9KYVwBqn/gtJMmPyA0UmEPULjwd0f0re2eHj7ySJOFlaxbIrMbOjpro0PI8jLlfd3ZhgdRQi7CquW+lVbXDREKbUUSABeaOm41voNrXWC1johODjY/pRCOIPEv4LfQBi2wOgkwh4ms6138cwxqMi75PB729MJ9HJlwej+BoQT7eFmMXPLxIF8dSSfnNIao+OIPs6ewiobaL4qXgSQc/FJSql5wP8DbtBayx4Dom/JO2jr+Zh0H5hlkrPTiJxqK7DSt13wcnZJNZuO5HPLxIG4u5gNCifaY8nkSDS2rYeEMJI9hVUSEKuUilFKuQK3AhfsB6GUGge8jq2oKnB8TCF6uMS/gsUDxt1hdBLRHm7eMGA8ZCdBQ+35l1d+l45JKe6QYUCnERHgyZXDQli1K4v6xiaj44g+rM1frbXWjUqpHwNfAmbgHa31IaXUc0Cy1nodtqE/b+Djs7u+Z2qtb+jC3EL0HFVFtjWRxtwGnoGXPn0merbombbCKjsJeJCymgZW78rk+tH9GeAvS2Y4k2XTornznV1sOJDDTePkGSphDLvGLLTWXwBfXPTa083+Ps/BuYRwHrvfhcZamPwjo5OIjvCPtL2lbwWtWbUrk6p6K/fOHGR0MtGGjxIvHPbTWhPs48YLX6ZSXWfldulxFAaQldeF6IzGetvWKIPmQMhwo9OIjoqeBVUFNBzbyLvfpTN9SBDx4X5GpxLtpJRi2uAgckprSS+qbvsCIbqAzLIVoj0uHubL2gUVORD3PRkCdFKJp4pRTVGMtXhzat0fySv/GdfGh13SGyKcw7iBAfzfoXy+kzWthEGkx0qIjtJNkLYZfAZAsPRWOTNtMpMXMJFhVUlM9c4nNsTb6Eiig1wtJibFBHIkt5xM6bUSBpDCSoiOKjhiW/9o8JWgZBsNZ/e1ywxqtCs/9d6Ekq+nU5syKAil4N3t6UZHEX2QFFZCdNTJzbaVuweMMzqJcIBVBVGsU3NIKN+Ie50MIzkzPw8XRoX78ffkLCpqG4yOI/oYKayE6IiSdCg+aZu0bpIFJJ3d4QoPjlZ6khq9FHNTPbGZfzc6kuikGUOCqaxr5EOZKye6mRRWQnTEyc3g4mFbuVs4vX/m9sPP0kjU0DFkh8whNnMNZmtt2xeKHis8wIOZsf14a+spahusRscRfYgUVkK0V2U+5B2AqBlgcTM6jeikY5XuHKjw4nuhxbiYTRyNXoZ7fTExp9e1fbHo0ZbPGcKZyjo+Tsk2OoroQ6SwEqK9TmwCkwViZhudRDjA2rx++JgbuSq4BICCwASK/OKJO/UuSktPhzObMiiQcZH+vL7lJI1W2eZGdA8prIRoj6ozcDoFoqbZ9pkTTu1giYXdZd7MDy3B3axtLyrFoUH34FOdxcC8jcYGFJ2ilOKhOUPILqnh8/05RscRfYQUVkK0x8mvbEsrDL7S6CTCAf502Asvs5VrzvZWnZMdeiVlXtGMSHsHtDYonXCEK4eHMDzMh1e/PklTk3wtRdeTwkoIe5Wdtq20PnAKuMt2J84upcjCV7lufC+0CC/LRcNEysSRmLsJLD9C2JkdxgQUDmEyKR6cM5jjBZX83+E8o+OIPkAKKyHstf1/AS29Vb2A1vDCQW/6uTVxXUhJi+ekD1hAtVsII9Pe6uZ0wtEWjOrPoH5e/Gnjcem1El1OCish7FFZACnvQkQCeAYZnUZ00ncFLuwsdOWh4VX/mVt1kSazK0djlhFanERQyb5uTigcyWI28V/zYknNr2DDgVyj44heTgorIeyx7WWw1sPgq4xOIjrpXG/VAA8rSwbVtHruiYE/oM7Fj/iTb3ZTOtFVvjd6AENDvXl50zGs0mslupDF6ABC9HgVeZD8Noy+BbyDjU4jOmljriv7Slz444Ry3C5aNH9w5seXnF8QMJ6BBV8zKvV/ODDsv7oppXA0k0nxs3lDefDD3Xy29zSLxkcYHUn0UlJYCdGWbS+DtQFmPwZpW4xOIzrBquGlg94M8m5kUZR9K6vnB06k/5kdRBR+K4WVk/noou1smrSmv587v9lwhKo6K2aTYsnkSIPSid5KhgKFaE15LiS/A2Nvg8BBRqcRnbQ+y43Ucgs/G1mFxc7vflazO3lBkwmoOEZA2ZGuDSi6lEkp5sWFUlxVz+6Mlh9aEKKzpLASojXbVoC2wqzHjE4iOqmhCVYc8iLOr4EFEXXtujYvaDKNJjfiT77eRelEdxke5kNkoCebjuZT1ygr6wvHk8JKiMspO217EnDsEgiINjqN6KSP093JqLLw6MgqTKp9157rtRqY/xX+5aldE1B0C6UU18WHUVHbyLbjZ4yOI3ohKayEuJxv/2h7hGzmo0YnEZ1Ua4U/H/FifGADV/av71AbeUFTqLd4M+rEqw5OJ7pbVJAX8QN82Xr8DAXl9s21E8JeUlgJ0ZKik7D7fUi4GwKijE4jOiHxVDG/TdLk1Zi5vl8uu9KLSTz1nzd7Wc3uHI1exsD8zQSWHuzCxKI7XDMyDGuT5k+bjhkdRfQyUlgJ0ZKvfwsWN5glvVXOrqrRxD9z+zHap4qRPtWdautozDJqXQIYc+zPDkonjBLk7cbkQYGsScoiNa/C6DiiF5HCSoiL5R2Ag/+AKQ+Cd4jRaUQnrc0Lospq4vaIgk631Wjx4tDge+lftIOQol0OSCeMdOWwELzcLPzuX/K0p3AcKayEaC55JfzzfnDxAI8g2/vN34RTyaoy8a+CAGYFlRHt2b4nAS/neOQtVLmH2nqttKzg7cw83Sz85MohfJNaKBPZhcNIYSVEc8VpUHAYBs8FV0+j04hOeuGgNybglgGO+6HZZHbj4JAfEVy6jwGF3zqsXWGMZVOjiQjw4PkvjshWN8IhpLAS4hyt4cg6cPOF6JlGpxGdtK/YwrosdxaEFhPk2ujQttPCb6TCM5KxqS+jtKyF5MzcXcw8ds0wjuSWs3bPaaPjiF5ACishzjn8GZSkw7D5tonrwmlpDc/v96afWxM3htn/5J/d7Ztc2Dv0YfwrTxCT/ZnD2xfd63ujBzAmwo8Xv0ylpl4KZdE5slegEACN9bDpV+DTHwZOMjqN6KSNua7sOuPKb8aV42Fq6pJ7ZIVdzRn/0Yw+/hcy+19Lo0WGjp3Ruf0EJ8UE8ebWNB5evYcrhv3noRXZS1C0l/RYCQGQ9JattyruBlDy38KZNTTB7/d7M9inkVtjunDxR6XYPfwxPOsKGX7qva67j+gWMf28GNHfly2phZTVNBgdRzgx6bESoqbEtsr64CshJM7oNKKTVqV5kFZp4a1ppXZvtGyvwZkfX/JakW8cI0++SaPJjQYXH05G/sCxNxXdZv6o/ry86RhfHsrj5oSBRscRTkp+NRfi2xehphSu+rXRSUQnlTcoXj7sxZTgeuZ2cOua9soKmYvCSkTBN91yP9F1Ar1cmRHbj71ZpWQUVRkdRzgpKaxE31aYCol/hXFLISze6DSik/561JPiehP/b3Qlqp0bLXdUnVsg+QETCS7dg2dNbvfcVHSZOUND8HW38Pn+HJpknTLRAVJYib5La/j3E+DiBXN/ZXQa0UlZVSbeOu7JwshaRgU4dnmFtpwOmU2j2ZPovH/LoqFOztVi4rpR/ckprSUlvcToOMIJSWEl+q6jG+DkZrjiF+AdbHQa0Um/2eeNWcHj8ZXdfm+r2Z2s0Ln4VGcRnbOh2+8vHGt0uB/RQZ58eTiPkqruGVIWvYcUVqJvaqiBL5+E4DiYeI/RaUQnbc134cscd34cV0V/z65ZXqEthf5jqfQYwNjUFVgaZX6OM1NKccOYcGobrPz+X0eNjiOcjBRWom/a/r9QmgnX/QHMLkanER2UeKqY704W83iSJ6Fu9YxxOU3iqeIL3rqNUqSHXYtnXSHxJ17vvvuKLhHm5870wf1Yk5xFcno3/jsSTs+uwkopda1SKlUpdUIp9UQLx2cppXYrpRqVUosdH1MIBzi3kfKWP9re+o+17Q0omyw7tX8VBpJT68adEfm4mIyd31TlGcHJ8IUMS38f38o0Q7OIzrsyLoQBfu788tODNFiN6QkVzqfNwkopZQZeAa4DRgC3KaVGXHRaJnAX8JGjAwrhUFrDgX+AyQQjbzI6jeik4noLn+QEMc63kgn+PWP4be+wn9Jo8WLSwWdByw9jZ+ZmMfOrG0ZyNK+Cld+dMjqOcBL29FhNAk5ordO01vXAauDG5idordO11vsB+S4ierac3XAmFYZfD+5+RqcRnfRuVgiNWnHXwHyjo5xX5xbEnmH/TUjJbgZlf2p0HNFJV48IZe7wEP608ThZxdVGxxFOwJ7CKhzIavZ+9tnX2k0pdb9SKlkplVxYWNiRJoTouPoqOLQW/CMharrRaUQnbcpxJbHUl+/3LyLMvWdtQZIWsZD8gAmMS30Jt7oio+OITlBK8dzCeMwmxeOf7EfLchqiDfYUVi0ts9ehf1la6ze01gla64TgYHm8XXSzo+uhoRpG3Sz7ATq5qkbF03t8iHCv43uhPbBwUYqk+KexNFYz/ugLRqcRnRTu78GT84ez/WQRq3ZltX2B6NPs2SswG2i+aVIEkNM1cYToImlbIHMHDLoC/CKMTiM66aWDXuTWmHhuWK7D9wPsrOb7Ceb2m0ZMzgbqLL6U+cSef132E3Q+SyZFsmF/Lr/94gizhwUT7u9hdCTRQ9nzLSkJiFVKxSilXIFbgXVdG0sIB6qrhHU/Ac9+MOw6o9OITtpbbOHdEx7cPqiGod61RsdpVU6/GVS7BTMoZz1ma8/OKlqnlOIP3x9Nk9Y8+c8DMiQoLqvNwkpr3Qj8GPgSOAL8XWt9SCn1nFLqBgCl1ESlVDbwA+B1pdShrgwtRLts/rVtzaoxt4HZ1eg0ohNqrfBIki9hHk38fFTPeAqwNdpkIS38BlwaK4nM+z+j44hOGhjoyePXDufbY4V8sDPD6Diih7JnKBCt9RfAFxe99nSzvydhGyIUomfJ2AGJr8Ok+yFosNFpRCf98YA3aRUWPphZgq+Lc/QYVHmEk9NvOuFntlHiO5xSn6FGRxKdcMeUKL5OLeDXG44wMSaQ4WG+RkcSPUwPm50ghAPVV8NnD9meApwnmyw7ux0FLrxzwpNlg6uZEdqzngJsy+ngWVS7hRCTswGztcboOKITTCbFiz8Yg6+7Cw+v2kNNvdXoSKKHkcJK9F4bn7KtrH7jX8DVy+g0ohMqGxSPJfsS7d3IE6O6f5PlztImCyfDb8TSWEVMzgbbQrXCafXzdmPFzWM4ll/JbzYcNjqO6GHsGgoUwik035Im/zAkvWV7CrDopO1NOCWt4ak9PuRUm/h4TgmeTvpdq9qjP6dD5jCwYDMxpz/jVMRCoyMJO3yUmHnZYzNj+/FhYiZTBwdx/egB3ZhK9GTSYyV6n7pK2L8KfPrDsAVGpxGd9I8Md9ZmuvNfI6qY0K/R6DidktNvGuWeUSQc/i0+VTL52dldNSKUhKgAHvt4P0dyy42OI3oIKaxE76I17F9tWwh03B1gdtLuDQHAiXIzT+/xYUpwPT+O6wXbiSgTJyNuosnkwrR9j2Nqcq65YuJCFpOJV5eOx9fDwv3vJ1NaXW90JNEDSGElepesnZB/EIZ/D3yla96Z1Vrhx4l+eJg1/zOpHHNLe0A4oXoXXxLjnyWo7BBjjv2P0XFEJ4X4uPPXpRPIL6vjJ6v2YG2S+XN9nfw6L3qPqkLbXoD9hkLMLKPTiDYknipu9fjrGWEcLbOwcnopoR69a3/37LB5HIu8lbhT73HGfwxZYVcZHUl0wrjIAH69cCSPf3KAX68/zDM3jDQ6kjCQFFaid7A2wp4PwGSGMUtkL0Ant6nQj81n/HloeBVX9O+dwyu7435OYPlhpux/ilLvIVR4xxgdSXRA88ntM4b0493t6RRW1DF9SL/zry+ZHGlENGEQ+ekjeoetL0FpBoz6AXj4G51GdMKxSnfeyQpjjG8lj4zs+aurd1STyYWtY1/CanZj5p6fYWnsBXPI+rhr48MYOcCXLw7kcvB0mdFxhEGksBLOLysJtvwBwhNgwHij04hOKG0wsyItnCCXBh6Oyek186ouNjjzYwZnfsyAwq2c6r8Av8o05ibexeCMv1+wibNwLialuDlhIAMDPfl7chYZRb33FwNxeTIUKJxbdTH844fgFw7x3zc6jeiEhibFipPhVDWa+fXwDLwtTW3Ow+oNyr1jyAydR1T+RmoLNpMdOtfoSKITXMwm7pgSxevfnuS9HencM2OQ0ZFEN5MeK+G8mppg7Y+gMh9+8B64eBidSHSQ1rbJ6qlVniyPziXas87oSN0qL2gK+QHjCT/zHf1K9hodR3SSl5uFu6fH4O5iZuV3pziWX2F0JNGNpLASzmv7n+H4l3D18xAuQ4DO7NO8ILYW+3HzgEKmBvbBH0JKkdH/Osq8BhGTs57QokSjE4lO8vd05Z7pMZhNitvfSiT9jAwL9hVSWAnnlLEdvnoORiyESfcZnUZ0QmKJD6tzgpkeWMaisCKj4xhGKzPHBy6m1i2Imbt/SkCZ7EHn7IK83bh7egzWJs3tbyVyulQ24O4LpLASzqc0E9bcAQFRcMOfQfXSGc59QGqlB/97qj+xXjX8KCqvz38prWZ3UqOW0GDx4YqkB/CtkD0unV2orzt/u3sS5bUN3P7mTgrKa42OJLqYFFbCudRVwqolYG2A29aAu5/RiUQHna515Y8nIujn2sjPB2fjapIVqwHqXfz4atJbaJOFK5Puw6s6y+hIopP2Z5dx+6RIckpruf5/t/Hmt2l8lJh5/k30LlJYCefR1ARrH4CCQ7D4HQgeanQi0UEFNSZ+d3wgZqV5MjYLXxer0ZF6lEqvSDZPfANzUz1zd90rxVUvEBnkxbKpURRX1fPOd6eoqnPuDcXF5UlhJZxD8kr4cDEcXQ9xN0BZlu215m/CKZTVK+76zo/yRjNPDMki1E02Im5JmU8smye+jktjFVftvAvfyjSjI4lOGhTszdIpURRW1PH2tlNUSnHVK0lhJZzDiU1w8iuInAoxs41OIzqoskFx5zZ/TpRb+O/B2Qzy6lvLKrRXid9INk1eidJW5u28i4CyI0ZHEp00NNSHZVOjKaqq482taZTXyi8WvY0UVqLn2/WmradqwHjbljV9fYazk6pphHu+8+NAiYW/TCljjK9s4WKPMp9YNk15D6vZnbm77pGlGHqBISHe3DktmrLqBt78No2sYvm/0C5T3aIAABFFSURBVJsorY2ZMJqQkKCTk5MNubdwIns+gM8egtB4mPBD2ybLwunUWuH+7X5szXfl5Unl3BhZ1ydWVXck14YyhmWswr3uDOkDFrBzzG+NjiQ6KaOoivd2pOPr7sJ7d08irr+v0ZFEK5RSKVrrhLbOkx4r0TNpDd++aCuqBl0B4++UospJVTYo7tzqz9Z8V/4woYIbI2X4ryPqXfw4HHMX5d4xDMr5nDGpfwLdZHQs0QlRQV48MGuwbY/B13ewM63vruPWm0hhJXoeawOs+zFs/jWMuhmWrAGzi9GpRAeU1itu/9af5CIXXp5Uzs0xsoZPZ1jN7qRG3kZ+QAIj095hTspDuNaXGR1LdEKorzufLJ9GiI8by97exT93Zxsd6f+3d+fRcVX3Ace/v1k1o13IkhdJ2Ja3eAODMcY0AQKlZjksORDIAmnIciDN3pxSSJv2JOnCSU4LPQnJoQklJEBJKElpCRSTUprgQDAYY1vG8S7Ltqx1NJJGM5qZ9+sfb2wJr6o9mieNf59z3rnvvXka/3w1mvnNvffda06TJVZmYunbBz++0e0CvORu+MBDEAh7HZU5Bc9ujXPdmnI2x/x8afY+pmYP8NqunsObOUXiY/e0q/jdor+kvus1Vq/9INV9m72OypyGGVURnrpzFeefXc2Xf7qBv/3lFrKOzes2WVliZSYGVTeZevAi2PcG3PB9uOxeG6g+Sb3VE+CrW2bSmQpy95w2Lqga8Dqk4iLC9qYP8uLKHyGqXPnb25i353H378hMStWlIR79xApuv+hsHvrfndzxyOv0Dg57HZY5BZZYGe8dbIHHbnbHU01dDHe9Aud+yOuozCl6Zm+YW/6nmpDP4RsL9rDU7v4bN91VS3j+4idpr13J8pa/49J1d1GS6vI6LHOKgn4fX79+MX9z42LW7ujiqgd+zdod9vucbOyuQOOdgy3w8n3Q8gsIlcP7/wJWfBp8x8j3bQLQCS+Vhfs2lvHw9igraof51PQ9NqN6oahS17uOs9vXkPWF2D3tal5f/DVr8Z0kPnxh01HnNu3r4/NPrGdX9yB3XtLMl66YRyhgbSFeGutdgZZYmcJKDUDLv8Nbj8Oe37gJ1co7YeVnIFpz/J+zxGpC2x7387nXKtjSF+SP5yS4d+kA6/fYOKpCK0l1MqftF5QmD7C37jLWLbyXochUr8Myp2g44/Dsxv28vruX5imlfOOGxaxqrvU6rDOWJVZm4hjshm0vwNZfujOopxNQ0wznfhiW33HihOoQS6wmpIwDj2yP8O3NZUT9yrcviPP+ae64EBug7hF1mNb9KtM7f4OKn7fnfZbfN92C+uzO2slqa3s/L23toLUnwQ3nTueeq99DfUWJ12GdcSyxMoVzZNKjCoMd0L7JXTC5ZxegEK6E+kXQcAFUz7RuikluQ0+Ae98sZ3MsyGVTU9x3fj91kZF5lSyx8lZ77UpWbP4m07rWEi+dyfr5X2Zf3aX2dzdJfeC8GTz40na+//JOROC2lWdz56XN1JbZXdOFMtbEKlCIYMwZwMlCz043kWrfBIncgMuKBph7pTtzeuUMEBsjMNm1D/m4v6WUJ3eVUFfi8ODKPq6akbLP6wlmMNrIS8u/z/TOX7PsnW9xyZuf52DNcjbOuYuOmgsswZpkSoJ+vnzlfG5e3sj9L27j4Vd28fjvWvnIhU18bNVMGqqjXodocqzFypw6VXdqhDV/BfvXQ3rQnR39rHluy1T9IohUex2lyZPYsPC9d6I8sj2Ko3D7nCG+uHCQ8uCx30OsxcpbO5puPrwvTpq5e3/Goh3/TCTVRUf1MjY3f4oDtRfbl51J4sgB7ts7BvinX23j2Y0HUFVWL57K7RfN5MJZNYglzePCugLN+OnZBW//FN5+Enp2gC/oTpMw7VyYssAm9Cwy+xM+Ht4W5YldJSQywo1NSS6r2E9dOO11aOYERidWh/iyKZrbnmbhzocpTbbTH21iW9MH2TnjBoZDlR5EacbqWHcOAuyLDfHob3fzo7W7SaYdakpDLGuqYlljNTWloTE9hxkbS6xMfiX7YPPP4a0nYO+r7rmZ74Wlt0CqH4IRb+MzeaUKb3QHeWxnhP/YG0aBaxtS3LVgkAWVWWuNmgSOlVgd4nPSNLa/wJJt36UisRdH/MTK5tJduZhY+Vyc3ED3Ez2HKayTJUWPvLKbzfv7eKO1l52dgwBMqyxh0fQKFk6vpL48zEdWnl2IUIuWjbEyp8/JwvP3QNvr0L4RnDSU1cGCa2DGcrebz8lYUlVE2od8PNsW5sldEX4fD1AWcPho8xCfmJugsdQW/J1Mmlt/dtJrtsz6OJHkQep611MTb6Gm/x2yviB9pbPpK5vDgSl/QCIyrQDRmtMVCvhY1lTNsqZqehPDbNrXx+b9cX61pYMXt3RwVmmI1p4EVy6q55yGKgJ+6wIeL9ZiZd5NFdrfhk1Pu119/QcgGIXp50HjCqhstEGvRUQVdvT7eflgiOfawqzrdrsOmqNDXDElxqrqOCV+WybljKAOFYk91PRtoWpgG+G0u7jzQGQGXVXn0FV9Dp1V5xArn2dTN3jgZC1Wj7/Weszz/ck0LQfitOyPs7t7kHRWKQ8HuHB2Dauaa1k15yzm15fbuKwxsK5AM3aJHmh9FXb8CrY+B/F9IH73br7yqVC3CPzWuFkMEhloiQXZFAuwoSfIKx1BOpJ+ABZUprmmIcXVDSm6ujo9jtR4SpWS4W4ygTKm9K6nNraBaKoDgIyvhO7KRXRXLSFWPo/eivnES2dZsjUJXLN0Gr/e1snaHd2s3d7F7m53uanashArZ5/FeU3VLG2oZOH0CqIhe88/kiVW5miqMNgJne+4y8l0bIa2ddDR4j4ejELz+2H+1TDvj6C01ibmnGCON7YpmRVimQCxdIDa6io6hnx0JH10Jt3S3ffTkxIU95tpZSDDwvIESyoGWVKesMHo5iiHx1ipEk22UxvbwJTet6iNbaAqvhW/uq8ZR/wMhaeQKKknEa53y5J6MoGojdOawGKJYXZ0DrKjc4CdnQPEkxkAfAJz68pZPKOSpQ2VzJ9azuzaUqaUh8/olq28JlYishp4APADP1DVvz/i8TDwKHA+0A3coqq7T/SclljlUTYDQ72Q6HYn5hw4tB10E6mBgxDbC7FWyAyN/FykBqafC2evgqZVMON8CB4xm68lVp5Rdac4GEmO/KzbnySWDhzeetN+YukAQ47/qJ8PiFJb4lCX26aUOEyNZAmkYsyOJqkOZqxX15wy0SwlqW6iyYNEk+1EUx1EkwcJZQYOXzMcKKe3YgH90SYGog25spH+aCOZYJmH0ZtjiQ+l2RcbcrdetxxIZQ4/Hgr4mDOljFm1pcyqLaWpJkp9ZQn1FWHqy0uoigaLOvHKW2IlIn7g98AfAm3A68CHVLVl1DWfAZaq6p0icitwo6recqLn9TyxUgV1RpUOMGpfFTQLjuOW6riDuTWbK52R8vC5Udc7Ryw++64XmxznnEImCemkmwClR22puNtlN9QzUsZaYXjQ/ZljET+Ey90tUu0uHROpgbJ6t4svXDFpxkupggKOgpMrVcFByDiQygqpXJkctT/6fMrJHWch7QgZhUyubO1NklXJbZBVQXG/uflRpleG8QkEBHyiBHzgF/Dn9gOj9v3iJjV+HwRzpSpkVMg6kFbIOkJa3Zam+LDQl/YRGxb6hn30pd2yM+lj2Dn69xP2OVQFM1QHM+8qqwJZqoIZ3jczSl3EoTqk+I7x67U7+sx4CmQG3UQreZBosgPER1milcjwu193aX+UVKiaVLCSVKia4UNloIysP4zjC+H4gmR9YbKj9h1fKHccQMWPI35UAjgSQH0jx+9+zI/6AoePJ8v7ntdUlXgyQ0c8SdfgMN0DKUIBH7u7BtnbO0TWeXf+EPL7qKsIU1cepioaojISPGorDfspCfqJBP1EQm5ZEvQTDvjw+YSAT0ZKcUu/TyZEwpbPuwJXANtVdWfuif8VuB5oGXXN9cBf5/afAr4jIqJe9TOCO1boqTuOTp4OJVCTUUmlmyBFaiB6ljuxX7AUQtGRMlwxkkwFo5PqDWTPgJ/Va2rcpIlcEpVLnsaTXxQfETcZyiVIfnE7zNwkTpC4z024HDdBcnKJUj6EfUplyKEypFQGHaZFHBZUZpgSdqiLuC1Nh1qdWg92E/Gf+O68eO8g8d68hGbM/1smUEq8rJl4WTMw0p0YyAxSlthLeWIvZYOtRIa7CA/HCKdjVAzsJJAZIpBNEHBS4x6jkusSFwEk94kw8oXX/esXtjfexJsL7x73eCYqETmcEM094rGso/QNpelPpoknM8SH0vQnM9SWhejoT9HRn2RbRz+xhHv+dPkEN8FCRn5VjHzEfeXK+XzyvbNP+9/Jh7G0WN0ErFbVT+aObwMuVNXPjrpmU+6attzxjtw1XUc816eBT+cO5wNbTxJfLdB1kmtMflhdF5bVd+FYXReW1XdhWX0XznxVLT/ZRWNpsTrW1/Ijs7GxXIOqPgQ8NIZ/031SkXVjaXYzp8/qurCsvgvH6rqwrL4Ly+q7cERkTOOXxjJDWBvQOOq4Adh/vGtEJABUAjaQwxhjjDFnlLEkVq8Dc0VkloiEgFuBZ4645hngY7n9m4D/9nR8lTHGGGOMB07aFaiqGRH5LPBfuNMtPKyqm0Xk68A6VX0G+CHwYxHZjttSdWue4htzt6E5bVbXhWX1XThW14Vl9V1YVt+FM6a69myCUGOMMcaYYmOrMBpjjDHG5IklVsYYY4wxeTIpEisR+YqIqIjUeh1LMRORb4nIOyLytoj8XESqvI6p2IjIahHZKiLbReTPvY6nmIlIo4i8JCJbRGSziHzB65iKnYj4RWS9iPyn17EUOxGpEpGncu/ZW0TkIq9jKmYi8qXc+8gmEXlCREqOd+2ET6xEpBF3OZ1Wr2M5A6wBFqvqUtxljO7xOJ6iklse6rvAVcBC4EMistDbqIpaBvhTVX0PsBL4E6vvcfcFYIvXQZwhHgCeV9UFwDlYvY8bEZkBfB5YrqqLcW/kO+5NehM+sQL+EfgzJu06NJOHqr6gqofWHngVd84ykz+Hl4dS1WHg0PJQZhyo6gFVfTO334/7wTPD26iKl4g0ANcAP/A6lmInIhXA+3DvyEdVh1U15m1URS8ARHJzdUY5ej7PwyZ0YiUi1wH7VHWD17Gcge4AnvM6iCIzA9g76rgN+6AvCBGZCSwDXvM2kqJ2P+6X4BMvZmnyYTbQCfxLruv1ByJS6nVQxUpV9wHfxu05OwD0qeoLx7ve88RKRF7M9VkeuV0PfBX4mtcxFpOT1Peha76K243ymHeRFqUxLf1k8ktEyoB/A76oqnGv4ylGInIt0KGqb3gdyxkiAJwHfE9VlwGDgI3ZHCciUo3buzALmA6UishHj3f9WNYKHFeqesWxzovIEtz/xAZxl69uAN4UkRWq2l7AEIvK8er7EBH5GHAtcLnNnp93Y1keyuSRiARxk6rHVPVpr+MpYhcD14nI1UAJUCEiP1HV4374mNPSBrSp6qEW2KewxGo8XQHsUtVOABF5GlgF/ORYF3veYnU8qrpRVetUdaaqzsR9IZ1nSdX4EZHVwN3Adaqa8DqeIjSW5aFMnoj7jeyHwBZV/Qev4ylmqnqPqjbk3qtvxV3WzJKqcZL7HNwrIvNzpy4HWjwMqdi1AitFJJp7X7mcE9ws4HmLlZlQvgOEgTW5VsJXVfVOb0MqHsdbHsrjsIrZxcBtwEYReSt37l5V/aWHMRmTL58DHst9SdsJfNzjeIqWqr4mIk8Bb+IOk1nPCZa3sSVtjDHGGGPyZMJ2BRpjjDHGTDaWWBljjDHG5IklVsYYY4wxeWKJlTHGGGNMnlhiZYwxxhiTJ5ZYGWOMMcbkiSVWxhhjjDF58n9wUajkCy1G9gAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 720x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure(figsize=[10, 4])\n",
"\n",
"sns.distplot([np.random.normal(mu_target, sigma_target) for _ in range(3000)], label=\"distribution $p(x)$\")\n",
"sns.distplot([np.random.normal(mu_appro, sigma_appro) for _ in range(3000)], label=\"distribution $q(x)$\")\n",
"\n",
"plt.title(\"Distributions\", size=16)\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": 181,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"average 0.9959984807502844 variance 83.36158359644132\n"
]
}
],
"source": [
"# calculate value sampling from a different distribution\n",
"\n",
"value_list = []\n",
"# need larger steps\n",
"for i in range(n):\n",
" # sample from different distribution\n",
" x_i = np.random.normal(mu_appro, sigma_appro)\n",
" value = f_x(x_i)*(p_x.pdf(x_i) / q_x.pdf(x_i))\n",
" \n",
" value_list.append(value)\n",
"\n",
"print(\"average {} variance {}\".format(np.mean(value_list), np.var(value_list)))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment