diff --git a/Exercise sheet 3/exercise_sheet_03.ipynb b/Exercise sheet 3/exercise_sheet_03.ipynb index 0b9fa94..73dd500 100644 --- a/Exercise sheet 3/exercise_sheet_03.ipynb +++ b/Exercise sheet 3/exercise_sheet_03.ipynb @@ -231,7 +231,7 @@ "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEICAYAAABI7RO5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdY0lEQVR4nO3dfZzNdf7/8cfLGKELylWicpEbUYhBKEPC1ipMSt1ascWkJL7f32qlzVarbb9xayuFrK1QbRYp0m6bXbNqU7loulAhF635KmQzMYZmeP3+mJnzNcyM+ejMfM6ced5vt7lxznwuXuczc+Z5Xp+L98fcHRERkdKqEnYBIiJSsSg4REQkEAWHiIgEouAQEZFAFBwiIhJI1bALKA9169b1Jk2ahF2GiEiFsXbt2m/dvV5R34vr4DCza4BrLrjgAtasWRN2OSIiFYaZfVXc9+J6V5W7L3X31Fq1aoVdiohI3Ijr4BARkehTcIiISCBxfYxDyl9OTg4ZGRkcPHgw7FJEpBSqV69O48aNSUxMLPU8Cg6JqoyMDE4//XSaNGmCmYVdjoiUwN3Zs2cPGRkZNG3atNTzxfWuKjO7xsxmZWZmhl1KpXHw4EHq1Kmj0BCpAMyMOnXqBN5DENfBobOqwqHQEKk4Tub9GtfBISIi0adjHCfwyoavwy4BgJSWDcMu4aREe/sF3Q4PPPAAp512Gr/4xS+iVsPVV1/NSy+9BMBLL73EnXfeGZXljh8/njfeeIOrr76aKVOmRGWZP8Zvf/tbJk6cGHYZJVqzZg1z587lySefDDzv448/TmpqKjVr1gT+7+dau3btKFdZdtLS0pg6dSqvv/56ua5XHYdIQG+88Qa1a9dm7969TJ8+PWrLfeaZZ1i3bl1MhAbkBUcsyM3NLfZ7SUlJJxUakBccBw4ciDwu+LnKianjCGDS0OvKdX0PzVtUruuLFw8//DBz587l3HPPpV69enTs2BGAzZs3M3r0aHbv3k3NmjX5wx/+QKtWrRg+fDhnnHEGa9as4ZtvvuHRRx9l8ODBfP311wwZMoTvv/+e3NxcZsyYweWXX06TJk1Ys2YNEyZMYPPmzbRv354+ffrwzTffMHjwYAYMGADAzTffzJAhQ7j22msjtbk799xzD3/5y18wM371q19FpsnKyqJLly7ce++9DBkyJDLPBx98wLhx48jOzqZGjRo899xztGzZksOHD/PLX/6SN998EzNj5MiRjBkzhtWrVzN27FiysrI45ZRT+Pvf/07NmjWZMGECaWlpHDp0iNGjR3P77beTlpbGpEmTqFOnDhs2bKBHjx5Mnz6diRMnkp2dTfv27WnTpg0vvvgiAwcOZPv27Rw8eJCxY8eSmpoKwGmnncbYsWN5/fXXqVGjBq+99hoNGjRg586djBo1ii1btgAwY8YMunXrxgsvvMCTTz7JDz/8QJcuXZg+fToJCQmFfobPP/88y5Yt4+DBg2RlZbF06VLGjBnDJ598Qm5uLg888AADBgwo9Ik7KyuryGmK2k7uzo4dO+jVqxd169ZlxYoVkZ9r3bp1eeyxx3j22WcBGDFiBOPGjWPbtm1cddVVXHbZZbz77rs0atSI1157jRo1ahSqfcGCBTz44IMkJCRQq1YtVq5cybZt2xg6dChZWVkAPPXUU3Tr1o20tDR+/etf06BBA9LT00lJSeHiiy/miSeeIDs7m1dffZXmzZszfPhwqlevzvr169m5cyePPfYY/fv3L7Te4l7/+vXr+fnPf84PP/zAkSNHWLRoES1atPhxbzJ3j/uvjh07+sla9MWOyFebTl3L9evodVcUn332WaHHR7+GaHydyJo1a/yiiy7yrKwsz8zM9ObNm/uUKVPc3f2KK67wjRs3urv7e++957169XJ392HDhvngwYP98OHDvn79em/evLm7u0+dOtUnT57s7u65ubn+/fffu7v7+eef77t37/atW7d6mzZtIutOS0vzAQMGuLv73r17vUmTJp6Tk1OovoULF/qVV17pubm5/s033/i5557rO3bkva5TTz21yNeUmZkZWc5bb73lKSkp7u4+ffp0T0lJiXxvz549fujQIW/atKl/8MEHheZ95pln/De/+Y27ux88eNA7duzoW7Zs8RUrVvgpp5zimzdv9tzcXL/yyit9wYIFRdazZ88ed3c/cOCAt2nTxr/99lt3dwd8yZIl7u4+fvz4yHpuuOEG//3vfx/Zfnv37vXPPvvM+/fv7z/88IO7u99xxx0+Z86c417zc889540aNYqs89577/V58+a5u/t3333nLVq08P379/uKFSv8pz/9aYnTFLWdjv45Fih4XPA7tH//ft+3b5+3bt3a161b51u3bvWEhAT/8MMP3d39+uuvj6zvaBdddJFnZGRE6nB3z8rK8uzsbHd337hxoxf8TVqxYoXXqlXLd+zY4QcPHvRzzjnHJ02a5O7ujz/+uI8dO9bd835H+/Xr54cPH/aNGzd6o0aNPDs7u1Sv/6677vIXXnjB3d0PHTrkBw4cOK7mY9+37u7AGi/mb2pcdxxHD3IYDeoAYt/bb7/NoEGDIvutCz7t79+/n3fffZfrr78+Mu2hQ4ci/x84cCBVqlShdevW7Ny5E4BOnTpx6623kpOTw8CBA2nfvn2J605OTmb06NHs2rWLV155heuuu46qVQu/xd555x1uuukmEhISaNCgAcnJyaxevbpQV3KszMxMhg0bxqZNmzAzcnJyAFi+fDmjRo2KrOOss87ik08+oWHDhnTq1AmAM844A4C//e1vfPzxxyxcuDCyzE2bNlGtWjU6d+5Ms2bNALjpppt45513GDx48HF1PPnkkyxevBiA7du3s2nTJurUqUO1atUin347duzIW2+9BcA//vEP5s6dCxD59D1v3jzWrl0bqS87O5v69esX+br79OnDWWedFal/yZIlTJ06Fcg77fvf//53oemLm6ao7VSSd955h0GDBnHqqacCkJKSwttvv821115L06ZNI78HHTt2ZNu2bcfN3717d4YPH84NN9xASkoKkHdh7F133UV6ejoJCQls3LgxMn2nTp1o2DDv2F3z5s3p27cvABdffDErVqyITHfDDTdQpUoVWrRoQbNmzfjiiy9K9fq7du3Kww8/TEZGBikpKT++2yDOd1W5+1JgaVJS0siwa5HyU9TphUeOHKF27dqkp6cXOc8pp5wS+X/ehy3o0aMHK1euZNmyZQwdOpTx48dzyy23lLjuoUOH8uKLL/Lyyy9HdnUcrWDZQdx///306tWLxYsXs23bNnr27BlZ1rGvtajnCp6fNm0a/fr1K/R8WlracdMXNX9aWhrLly9n1apV1KxZk549e0bO/U9MTIzMk5CQUOIxCXdn2LBhPPLII4WeX7x4MQ8++CAAs2fPBoj84S6Yb9GiRbRs2bLQfAUhX9I0xW2TkmosztG/JwkJCWRnZx83zcyZM3n//fdZtmwZ7du3Jz09nWnTptGgQQM++ugjjhw5QvXq1YtcZpUqVSKPq1SpUmhbnujnVNzrv/DCC+nSpQvLli2jX79+zJ49myuuuKKkTXBCOjgucaVHjx4sXryY7Oxs9u3bx9KlS4G8T95NmzZlwYIFQN6b7KOPPipxWV999RX169dn5MiR3Hbbbaxbt67Q908//XT27dtX6Lnhw4fz+OOPA9CmTZsi65s/fz6HDx9m9+7drFy5ks6dO5dYR2ZmJo0aNQLy9v0X6Nu3LzNnzoz8cfnPf/5Dq1at2LFjB6tXrwZg37595Obm0q9fP2bMmBHpVjZu3BjZ3/7BBx+wdetWjhw5wvz587nsssuAvEAomD4zM5MzzzyTmjVr8sUXX/Dee++VWDNA7969mTFjBgCHDx/m+++/p3fv3ixcuJBdu3ZFav7qq68YNGgQ6enppKenk5SUdNyy+vXrx7Rp0yJ/1D/88MNST1PUdoKif36Q9zN69dVXOXDgAFlZWSxevJjLL7/8hK+3wObNm+nSpQsPPfQQdevWZfv27WRmZtKwYUOqVKnCvHnzOHz4cKmXV2DBggUcOXKEzZs3s2XLluMCorjXv2XLFpo1a8bdd9/Ntddey8cffxx43ceK645DwlfepxF36NCBIUOG0L59e84///xCb/gXX3yRO+64g8mTJ5OTk8ONN95Iu3btil1WWloaU6ZMITExkdNOOy2y26VAnTp16N69OxdddBFXXXUVU6ZMoUGDBlx44YUMHDiwyGUOGjSIVatW0a5dO8yMRx99lLPPPrvE13TPPfcwbNgwHnvssUKfFEeMGMHGjRtp27YtiYmJjBw5krvuuov58+czZsyYyMH05cuXM2LECLZt20aHDh1wd+rVq8err74KQNeuXZkwYQKffPIJPXr0YNCgQQCkpqbStm1bOnTowLPPPsvMmTNp27YtLVu25NJLLy2xZoAnnniC1NRU/vjHP5KQkMCMGTPo2rUrkydPpm/fvhw5coTExESefvppzj///BKXdf/99zNu3Djatm2Lu9OkSZPIKagFn7yLm6a47ZSamspVV11Fw4YNC+0S6tChA8OHD48E+ogRI7jkkkuK3C1VlPHjx7Np0ybcnd69e9OuXTvuvPNOrrvuOhYsWECvXr0KdVOl1bJlS5KTk9m5cyczZ84s1LWU9Prnz5/PCy+8QGJiImeffTaTJk0KvO5j2cm0zhVNUlKSn+yNnHQdRzCff/45F154YdhlhObAgQNcfPHFrFu3joowYkFY1wFEy6JFi1iyZAlz5swJu5QyNXz4cPr371/ksadoKOp9a2Zr3f349g/tqhKJmuXLl9OqVSvGjBlTIUKjoluyZAn33Xcft99+e9ilVDrqOE5AHUcwlb3jEKmI1HGIiEiZiuvg0LDqIiLRF9fB4RpWXUQk6uI6OEREJPoUHBJ3jh61Ni0t7bjB4Mra888/z44dOyKPR4wYwWeffRZ4OWHULlIaCg6JO9Ee7rwoJQ2rcWxwzJ49m9atW5dpPSLlScEhcefo4c7Hjx/P/v37GTx4MK1ateLmm2+ODMmwdu1akpOT6dixI/369ePrr/NOvU5PT+fSSy+lbdu2DBo0iO+++w6Anj17MnHiRJKTk3niiSeKnH/hwoWsWbOGm2++mfbt25OdnU3Pnj0pOB38r3/9Kx06dKBdu3b07t0byBvyo1u3blxyySV069aNDRs2hLDVREpPQ45I3Pnd737Hp59+Snp6OmlpaZF7Epxzzjl0796df/3rX3Tp0oUxY8bw2muvUa9ePebPn899993Hs88+yy233MK0adNITk5m0qRJPPjgg5Hxp/bu3cs///lPcnJySE5OLnL+p556iqlTpx435tLu3bsZOXIkK1eupGnTppExk1q1asXKlSupWrUqy5cvZ+LEiSxapJGYJXYpOCSqhg0bFhnSvDykpaWdcJrOnTvTuHFjANq3b8+2bduoXbs2n376KX369AHyBuFr2LAhmZmZ7N27l+TkZCDv9Rw9FHvBDZY2bNhQ5Pwlee+99+jRowdNmzYF/m947+KGTReJVQoOiXvHDoWdm5uLu9OmTRtWrVpVaNoTXfNTMDhdcfOXpLjhvYsbNl0kVik4JKrmzJkT+pAjxQ2XfbSWLVuye/duVq1aRdeuXcnJyWHjxo20adOGM888k7fffpvLL7+cefPmRbqP0s5f3Pq7du3K6NGj2bp1a2RX1VlnnVXssOkisSqugyPadwCUiuHo4c5r1KhBgwYNjpumWrVqLFy4kLvvvpvMzExyc3MZN24cbdq0Yc6cOYwaNYoDBw7QrFkznnvuuUDzDx8+nFGjRlGjRo1CHUm9evWYNWsWKSkpHDlyhPr16/PWW28VO2y6SKzSIIcnoEEOg9EghyIVjwY5FBGRMqXgEBGRQBQcEnWVYfenSLw4mfergkOiqnr16uzZs0fhIVIBuDt79uw57v7lJxLXZ1VJ+WvcuDEZGRns3r077FJEpBSqV68euUC2tBQcElWJiYmRK6NFJD5pV5WIiASi4BARkUAUHCIiEoiCQ0REAonr4DCza8xs1olGPBURkdKL6+Bw96XunlqrVq2wSxERiRtxHRwiIhJ9Cg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCievg0B0ARUSiL66DQ3cAFBGJvrgODhERiT4Fh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoFUDbsAKZ1XNnwddgkApLRsGHYJIhIydRwiIhKIOo4KaNLQ68p1fQ/NW1Su6xOR2KaOQ0REAlHHUQGpAxCRMKnjEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQmkwgWHmTUzsz+a2cKwaxERqYzKNTjM7Fkz22Vmnx7z/E/MbIOZfWlmE0pahrtvcffbyrZSEREpTnnfj+N54ClgbsETZpYAPA30ATKA1Wa2BEgAHjlm/lvdfVf5lCoiIkUp1+Bw95Vm1uSYpzsDX7r7FgAzexkY4O6PAP1Pdl1mlgqkApx33nknuxgRETlGLBzjaARsP+pxRv5zRTKzOmY2E7jEzO4tbjp3n+XuSe6eVK9evehVKyJSycXCrWOtiOe8uIndfQ8wquzKERGRksRCx5EBnHvU48bAjpBqERGRE4iF4FgNtDCzpmZWDbgRWBJyTSIiUozAwWFmp+afCRWYmf0JWAW0NLMMM7vN3XOBu4A3gc+BP7v7+pNZfhHru8bMZmVmZkZjcSIiQimOcZhZFfK6gJuBTsAh4BQz2w28Acxy902lWZm731TM82/kLyuq3H0psDQpKWlktJctIlJZlabjWAE0B+4Fznb3c929PnA58B7wOzP7WRnWKCIiMaQ0Z1Vd6e45xz7p7v8BFgGLzCwx6pWJiEhMOmHHURAaZjb52O8VHOsoKlhERCQ+BTk43sjMIscozKw+sDz6JUWPDo6LiERfkOC4HUg1s85m1gn4BzC1bMqKDndf6u6ptWrVCrsUEZG4UZqzquYC64APgdHAS0AuMNDdvyzb8kREJNaUpuOYkz/dreSFRhPgO+BnZja47EoTEZFYdMKOw93/Dvy94LGZVQVaA+2ASwHdUKkSeWXD12GXQErLhmGXIFKplWZXlbl7ZNDB/Cu9P87/mlfUNCIiEr9KdQGgmY0xs0I3tTCzamZ2hZnNAYaVTXk/js6qEhGJvtJcAPgT8o5v/MnMmgJ7gerk3aHvb8Dv3T29rAr8MTTkSNmZNPS6cl3fQ/MWlev6RKR4pTnGcRCYDkzPv0K8LpDt7nvLuDYREYlBgW7k5O45ZnaLu/9PWRUkFYM6AJHK62Tux/GlmT2ef4zjZ2amvyAiIpXIyQTHYvI6lR35/+paDhGRSuRkguNPwF+A7kBf4MyoViQiIjEtcHC4+xB3X+buG4D/Ap6OflnRodNxRUSir9TBYWZXmdn7ZrbBzP5sZl3dfScwogzr+1E0yKGISPQF6TimA/9N3jAjs4ApZnaTu2eVSWUiIhKTgpyOu9Pd/5X//+Vmtgp4n7xjHiIiUkkE6Ti2mdlkM6uW/zgH2FcGNYmISAwLEhwOpADbzewd4EsgzcxalEllIiISk0q9q8rdbwIws+rAReQNq94OmG1mzdz93LIpUUREYkmgIUcgMnbVmvwvERGpZE7mAsAKQ9dxiIhEX1wHh67jEBGJvrgODhERiT4Fh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQkk8JXjImF7ZcPXYZcAQErLhmGXIBKKuO44dOW4iEj0xXXH4e5LgaVJSUkjw65FysakodeV6/oemreoXNcnEoviuuMQEZHoi+uOQ+KfOgCR8qeOQ0REAlFwiIhIIAoOEREJRMEhIiKBKDhERCQQBYeIiASi4BARkUAUHCIiEoiCQ0REAonr4NAghyIi0RfXweHuS909tVatWmGXIiISN+I6OEREJPoUHCIiEoiCQ0REAlFwiIhIILofh8hJ0r3PpbJSxyEiIoGo4xCJAt37XCoTdRwiIhKIOg6RKFAHIJWJOg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISSFwHh27kJCISfXEdHLqRk4hI9MV1cIiISPQpOEREJBAFh4iIBKKxqkQqON0XRMqbOg4REQlEHYdIHNF9QaQ8qOMQEZFA1HGIxBF1AFIe1HGIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggunJcRKJCo/RWHuo4REQkEHUcIhJ1GqU3vqnjEBGRQNRxiEjUqQOIb+o4REQkEAWHiIgEouAQEZFAFBwiIhJIhQsOMxtoZn8ws9fMrG/Y9YiIVDblGhxm9qyZ7TKzT495/idmtsHMvjSzCSUtw91fdfeRwHBgSBmWKyIiRSjv03GfB54C5hY8YWYJwNNAHyADWG1mS4AE4JFj5r/V3Xfl//9X+fOJiEg5KtfgcPeVZtbkmKc7A1+6+xYAM3sZGODujwD9j12GmRnwO+Av7r6uuHWZWSqQCnDeeedF5wWISMyLhTGz4n28rFi4ALARsP2oxxlAlxKmHwNcCdQyswvcfWZRE7n7LGAWQFJSkkepVhGpQDT0SdmIheCwIp4r9g+9uz8JPFl25YiISEliITgygHOPetwY2BFSLSISRypLB1DeYuF03NVACzNrambVgBuBJSHXJCIixSjv03H/BKwCWppZhpnd5u65wF3Am8DnwJ/dfX2U1neNmc3KzMyMxuJERITyP6vqpmKefwN4owzWtxRYmpSUNDLayxYRqaxiYVeViIhUIAoOEREJRMEhIiKBxMLpuGXGzK4BrrngggvCLkVEKpFYuHodyu4K9rjuONx9qbun1qpVK+xSRETiRlx3HCIiYYvHYU/iuuMQEZHoU8chIlKG4nHYE3UcIiISSFwHh4YcERGJvrgODp1VJSISfXEdHCIiEn0KDhERCUTBISIigSg4REQkEAWHiIgEYu4edg1lzsx2A1+FXUcMqAt8G3YRMULbojBtj8K0PeB8d69X1DcqRXBIHjNb4+5JYdcRC7QtCtP2KEzbo2TaVSUiIoEoOEREJBAFR+UyK+wCYoi2RWHaHoVpe5RAxzhERCQQdRwiIhKIgkNERAJRcMQ5MzvXzFaY2edmtt7MxoZdUywwswQz+9DMXg+7ljCZWW0zW2hmX+T/jnQNu6Ywmdl/5b9PPjWzP5lZ9bBrikUKjviXC/w/d78QuBQYbWatQ64pFowFPg+7iBjwBPBXd28FtKMSbxMzawTcDSS5+0VAAnBjuFXFJgVHnHP3r919Xf7/95H3h6FRuFWFy8waAz8FZoddS5jM7AygB/BHAHf/wd33hlpU+KoCNcysKlAT2BFyPTFJwVGJmFkT4BLg/ZBLCdvjwD3AkZDrCFszYDfwXP5uu9lmdmrYRYXF3f8XmAr8G/gayHT3v4VbVWxScFQSZnYasAgY5+7fh11PWMysP7DL3deGXUsMqAp0AGa4+yVAFjAh3JLCY2ZnAgOApsA5wKlm9rNwq4pNCo5KwMwSyQuNF939lbDrCVl34Foz2wa8DFxhZi+EW1JoMoAMdy/oQBeSFySV1ZXAVnff7e45wCtAt5BrikkKjjhnZkbePuzP3f2xsOsJm7vf6+6N3b0JeQc+/+HulfJTpbt/A2w3s5b5T/UGPguxpLD9G7jUzGrmv296U4lPFihJ1bALkDLXHRgKfGJm6fnPTXT3N8IrSWLIGOBFM6sGbAF+HnI9oXH3981sIbCOvLMRP0RDjxRJQ46IiEgg2lUlIiKBKDhERCQQBYeIiASi4BARkUAUHCIiEoiCQ0REAlFwiIhIIAoOkRCY2ZVmNi/sOkROhoJDJBztyLsyWaTCUXCIhKMd8KGZnWJmz5vZb/PHRxKJeRqrSiQc7YBdwJvAbHevrCP0SgWksapEyln+MPffAl8Bt7v7qpBLEglEu6pEyl9rYDV5I7AeDrkWkcAUHCLlrx3wLnn3A3nOzBqEXI9IIAoOkfLXDvjU3TcCvwT+nL/7SqRC0DEOEREJRB2HiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigfx/vf8Mt7nTGkQAAAAASUVORK5CYII=\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEICAYAAABI7RO5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdWUlEQVR4nO3de3RU9dn28e9NCCcPoJxEUDmUBQICQgABJSAq1aJARNFlEaoQUUR43rdYxMqjVmsfYVkFBaRUBcRKAVEQWistKVpRQYwHVEAOlrwoIJUIIWAC9/tHknkIJCEbJ9mTyfVZKyuZPftw7z2Zuea3D79t7o6IiEhpVQm7ABERqVgUHCIiEoiCQ0REAlFwiIhIIAoOEREJpGrYBZSHevXqedOmTcMuQ0Skwvjggw++dff6RT0X18FhZtcC1/7kJz9h3bp1YZcjIlJhmNlXxT0X17uq3H2Zu6fWrl077FJEROJGXAeHiIhEn4JDREQCietjHFL+cnJyyMjI4NChQ2GXIiKlUKNGDZo0aUJiYmKpp1FwSFRlZGRwxhln0LRpU8ws7HJEpATuzt69e8nIyKBZs2alni6ud1WZ2bVmNiszMzPsUiqNQ4cOUbduXYWGSAVgZtStWzfwHoK4Dg6dVRUOhYZIxXEq79e4Dg4REYk+HeM4iVc2fh12CQCktGoUdgmnJNrbL+h2ePDBBzn99NP55S9/GbUarrnmGl566SUAXnrpJe66666ozHf8+PGsWLGCa665hsmTJ0dlnj/Gb3/7WyZOnBh2GSVat24dc+fOZerUqYGnffLJJ0lNTaVWrVrA/76uderUiXKVZSctLY0pU6bw+uuvl+ty1eIQCWjFihXUqVOHffv2MX369KjN99lnn2X9+vUxERqQFxyxIDc3t9jnkpKSTik0IC84Dh48GHlc8LrKyanFEcCkodeX6/Ienre4XJcXLx599FHmzp3LeeedR/369encuTMAW7ZsYfTo0ezZs4datWrxhz/8gdatWzN8+HDOPPNM1q1bxzfffMPjjz/O4MGD+frrrxkyZAjff/89ubm5zJgxg8suu4ymTZuybt06JkyYwJYtW+jYsSNXXnkl33zzDYMHD2bAgAEA3HLLLQwZMoTrrrsuUpu7c++99/KXv/wFM+PXv/51ZJysrCy6devGfffdx5AhQyLTvP/++4wbN47s7Gxq1qzJ888/T6tWrThy5Ai/+tWveOONNzAzRo4cyZgxY1i7di1jx44lKyuL6tWr8/e//51atWoxYcIE0tLSOHz4MKNHj+aOO+4gLS2NSZMmUbduXTZu3EivXr2YPn06EydOJDs7m44dO9K2bVvmz5/PwIED2bFjB4cOHWLs2LGkpqYCcPrppzN27Fhef/11atasyWuvvUbDhg3ZtWsXo0aNYuvWrQDMmDGDHj168OKLLzJ16lR++OEHunXrxvTp00lISCj0Gr7wwgssX76cQ4cOkZWVxbJlyxgzZgyffPIJubm5PPjggwwYMKDQN+6srKwixylqO7k7O3fupE+fPtSrV49Vq1ZFXtd69erxxBNP8NxzzwEwYsQIxo0bx/bt27n66qu59NJLeeedd2jcuDGvvfYaNWvWLFT7woULeeihh0hISKB27dqsXr2a7du3M3ToULKysgB4+umn6dGjB2lpafz3f/83DRs2JD09nZSUFC666CKeeuopsrOzefXVV2nRogXDhw+nRo0abNiwgV27dvHEE0/Qv3//Qsstbv03bNjAL37xC3744QeOHj3K4sWLadmy5Y97k7l73P907tzZT9XiL3ZGftp26V6uP8cuu6L47LPPCj0+dh2i8XMy69at83bt2nlWVpZnZmZ6ixYtfPLkye7ufvnll/umTZvc3f3dd9/1Pn36uLv7sGHDfPDgwX7kyBHfsGGDt2jRwt3dp0yZ4o888oi7u+fm5vr333/v7u4XXHCB79mzx7dt2+Zt27aNLDstLc0HDBjg7u779u3zpk2bek5OTqH6Fi1a5FdccYXn5ub6N9984+edd57v3Jm3XqeddlqR65SZmRmZz5tvvukpKSnu7j59+nRPSUmJPLd3714/fPiwN2vWzN9///1C0z777LP+m9/8xt3dDx065J07d/atW7f6qlWrvHr16r5lyxbPzc31K664whcuXFhkPXv37nV394MHD3rbtm3922+/dXd3wJcuXeru7uPHj48s58Ybb/Tf//73ke23b98+/+yzz7x///7+ww8/uLv7nXfe6XPmzDlhnZ9//nlv3LhxZJn33Xefz5s3z93dv/vuO2/ZsqUfOHDAV61a5T/72c9KHKeo7XTs61ig4HHB/9CBAwd8//793qZNG1+/fr1v27bNExIS/MMPP3R39xtuuCGyvGO1a9fOMzIyInW4u2dlZXl2dra7u2/atMkLPpNWrVrltWvX9p07d/qhQ4f83HPP9UmTJrm7+5NPPuljx45197z/0X79+vmRI0d806ZN3rhxY8/Ozi7V+t99993+4osvurv74cOH/eDBgyfUfPz71t0dWOfFfKbGdYvj2E4Oo0EtgNj31ltvMWjQoMh+64Jv+wcOHOCdd97hhhtuiIx7+PDhyN8DBw6kSpUqtGnThl27dgHQpUsXbrvtNnJychg4cCAdO3YscdnJycmMHj2a3bt388orr3D99ddTtWrht9jbb7/NzTffTEJCAg0bNiQ5OZm1a9cWapUcLzMzk2HDhrF582bMjJycHABWrlzJqFGjIss4++yz+eSTT2jUqBFdunQB4MwzzwTgb3/7Gx9//DGLFi2KzHPz5s1Uq1aNrl270rx5cwBuvvlm3n77bQYPHnxCHVOnTmXJkiUA7Nixg82bN1O3bl2qVasW+fbbuXNn3nzzTQD+8Y9/MHfuXIDIt+958+bxwQcfROrLzs6mQYMGRa73lVdeydlnnx2pf+nSpUyZMgXIO+373//+d6HxixunqO1UkrfffptBgwZx2mmnAZCSksJbb73FddddR7NmzSL/B507d2b79u0nTN+zZ0+GDx/OjTfeSEpKCpB3Yezdd99Neno6CQkJbNq0KTJ+ly5daNQo79hdixYtuOqqqwC46KKLWLVqVWS8G2+8kSpVqtCyZUuaN2/OF198Uar17969O48++igZGRmkpKT8+NYGcb6ryt2XAcuSkpJGhl2LlJ+iTi88evQoderUIT09vchpqlevHvk778sW9OrVi9WrV7N8+XKGDh3K+PHjufXWW0tc9tChQ5k/fz4vv/xyZFfHsQrmHcQDDzxAnz59WLJkCdu3b6d3796ReR2/rkUNKxg+bdo0+vXrV2h4WlraCeMXNX1aWhorV65kzZo11KpVi969e0fO/U9MTIxMk5CQUOIxCXdn2LBhPPbYY4WGL1myhIceegiA2bNnA0Q+uAumW7x4Ma1atSo0XUHIlzROcdukpBqLc+z/SUJCAtnZ2SeMM3PmTN577z2WL19Ox44dSU9PZ9q0aTRs2JCPPvqIo0ePUqNGjSLnWaVKlcjjKlWqFNqWJ3udilv/Cy+8kG7durF8+XL69evH7Nmzufzyy0vaBCelg+MSV3r16sWSJUvIzs5m//79LFu2DMj75t2sWTMWLlwI5L3JPvrooxLn9dVXX9GgQQNGjhzJ7bffzvr16ws9f8YZZ7B///5Cw4YPH86TTz4JQNu2bYusb8GCBRw5coQ9e/awevVqunbtWmIdmZmZNG7cGMjb91/gqquuYubMmZEPl//85z+0bt2anTt3snbtWgD2799Pbm4u/fr1Y8aMGZHWyqZNmyL7299//322bdvG0aNHWbBgAZdeeimQFwgF42dmZnLWWWdRq1YtvvjiC959990Sawbo27cvM2bMAODIkSN8//339O3bl0WLFrF79+5IzV999RWDBg0iPT2d9PR0kpKSTphXv379mDZtWuRD/cMPPyz1OEVtJyj69YO81+jVV1/l4MGDZGVlsWTJEi677LKTrm+BLVu20K1bNx5++GHq1avHjh07yMzMpFGjRlSpUoV58+Zx5MiRUs+vwMKFCzl69Chbtmxh69atJwREceu/detWmjdvzj333MN1113Hxx9/HHjZx4vrFoeEr7xPI+7UqRNDhgyhY8eOXHDBBYXe8PPnz+fOO+/kkUceIScnh5tuuokOHToUO6+0tDQmT55MYmIip59+emS3S4G6devSs2dP2rVrx9VXX83kyZNp2LAhF154IQMHDixynoMGDWLNmjV06NABM+Pxxx/nnHPOKXGd7r33XoYNG8YTTzxR6JviiBEj2LRpE+3btycxMZGRI0dy9913s2DBAsaMGRM5mL5y5UpGjBjB9u3b6dSpE+5O/fr1efXVVwHo3r07EyZM4JNPPqFXr14MGjQIgNTUVNq3b0+nTp147rnnmDlzJu3bt6dVq1ZccsklJdYM8NRTT5Gamsof//hHEhISmDFjBt27d+eRRx7hqquu4ujRoyQmJvLMM89wwQUXlDivBx54gHHjxtG+fXvcnaZNm0ZOQS345l3cOMVtp9TUVK6++moaNWpUaJdQp06dGD58eCTQR4wYwcUXX1zkbqmijB8/ns2bN+Pu9O3blw4dOnDXXXdx/fXXs3DhQvr06VOoNVVarVq1Ijk5mV27djFz5sxCrZaS1n/BggW8+OKLJCYmcs455zBp0qTAyz6enUrTuaJJSkryU72Rk67jCObzzz/nwgsvDLuM0Bw8eJCLLrqI9evXUxF6LAjrOoBoWbx4MUuXLmXOnDlhl1Kmhg8fTv/+/Ys89hQNRb1vzewDdz+x+Yd2VYlEzcqVK2ndujVjxoypEKFR0S1dupT777+fO+64I+xSKh21OE5CLY5gKnuLQ6QiUotDRETKVFwHh7pVFxGJvrgODle36iIiURfXwSEiItGn4JC4c2yvtWlpaSd0BlfWXnjhBXbu3Bl5PGLECD777LPA8wmjdpHSUHBI3Il2d+dFKalbjeODY/bs2bRp06ZM6xEpTwoOiTvHdnc+fvx4Dhw4wODBg2ndujW33HJLpEuGDz74gOTkZDp37ky/fv34+uu8U6/T09O55JJLaN++PYMGDeK7774DoHfv3kycOJHk5GSeeuqpIqdftGgR69at45ZbbqFjx45kZ2fTu3dvCk4H/+tf/0qnTp3o0KEDffv2BfK6/OjRowcXX3wxPXr0YOPGjSFsNZHSU5cjEnd+97vf8emnn5Kenk5aWlrkngTnnnsuPXv25F//+hfdunVjzJgxvPbaa9SvX58FCxZw//3389xzz3Hrrbcybdo0kpOTmTRpEg899FCk/6l9+/bxz3/+k5ycHJKTk4uc/umnn2bKlCkn9Lm0Z88eRo4cyerVq2nWrFmkz6TWrVuzevVqqlatysqVK5k4cSKLF6snZoldCg6JqmHDhkW6NC8PaWlpJx2na9euNGnSBICOHTuyfft26tSpw6effsqVV14J5HXC16hRIzIzM9m3bx/JyclA3voc2xV7wQ2WNm7cWOT0JXn33Xfp1asXzZo1A/63e+/iuk0XiVUKDol7x3eFnZubi7vTtm1b1qxZU2jck13zU9A5XXHTl6S47r2L6zZdJFYpOCSq5syZE3qXI8V1l32sVq1asWfPHtasWUP37t3Jyclh06ZNtG3blrPOOou33nqLyy67jHnz5kVaH6Wdvrjld+/endGjR7Nt27bIrqqzzz672G7TRWJVXAdHtO8AKBXDsd2d16xZk4YNG54wTrVq1Vi0aBH33HMPmZmZ5ObmMm7cONq2bcucOXMYNWoUBw8epHnz5jz//POBph8+fDijRo2iZs2ahVok9evXZ9asWaSkpHD06FEaNGjAm2++WWy36SKxSp0cnoQ6OQxGnRyKVDzq5FBERMqUgkNERAJRcEjUVYbdnyLx4lTerwoOiaoaNWqwd+9ehYdIBeDu7N2794T7l59MXJ9VJeWvSZMmZGRksGfPnrBLEZFSqFGjRuQC2dJScEhUJSYmRq6MFpH4pF1VIiISiIJDREQCUXCIiEggCg4REQkkroPDzK41s1kn6/FURERKL66Dw92XuXtq7dq1wy5FRCRuxHVwiIhI9Ck4REQkEAWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigSg4REQkEAWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCSSug0N3ABQRib64Dg7dAVBEJPriOjhERCT6FBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigSg4REQkEAWHiIgEUjXsAqR0Xtn4ddglAJDSqlHYJYhIyNTiEBGRQNTiqIAmDb2+XJf38LzF5bo8EYltanGIiEgganFUQGoBiEiY1OIQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigSg4REQkEAWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCaTCBYeZNTezP5rZorBrERGpjMo1OMzsOTPbbWafHjf8p2a20cy+NLMJJc3D3be6++1lW6mIiBSnvO/H8QLwNDC3YICZJQDPAFcCGcBaM1sKJACPHTf9be6+u3xKFRGRopRrcLj7ajNretzgrsCX7r4VwMxeBga4+2NA/1NdlpmlAqkA559//qnORkREjhMLxzgaAzuOeZyRP6xIZlbXzGYCF5vZfcWN5+6z3D3J3ZPq168fvWpFRCq5WLh1rBUxzIsb2d33AqPKrhwRESlJLLQ4MoDzjnncBNgZUi0iInISsRAca4GWZtbMzKoBNwFLQ65JRESKETg4zOy0/DOhAjOzPwFrgFZmlmFmt7t7LnA38AbwOfBnd99wKvMvYnnXmtmszMzMaMxOREQoxTEOM6tCXivgFqALcBiobmZ7gBXALHffXJqFufvNxQxfkT+vqHL3ZcCypKSkkdGet4hIZVWaFscqoAVwH3COu5/n7g2Ay4B3gd+Z2c/LsEYREYkhpTmr6gp3zzl+oLv/B1gMLDazxKhXJiIiMemkLY6C0DCzR45/ruBYR1HBIiIi8SnIwfHGZhY5RmFmDYCV0S8penRwXEQk+oIExx1Aqpl1NbMuwD+AKWVTVnS4+zJ3T61du3bYpYiIxI3SnFU1F1gPfAiMBl4CcoGB7v5l2ZYnIiKxpjQtjjn5491GXmg0Bb4Dfm5mg8uuNBERiUUnbXG4+9+Bvxc8NrOqQBugA3AJoBsqVSKvbPw67BJIadUo7BJEKrXS7Koyd490Oph/pffH+T/zihpHRETiV6kuADSzMWZW6KYWZlbNzC43sznAsLIp78fRWVUiItFXmgsAf0re8Y0/mVkzYB9Qg7w79P0N+L27p5dVgT+GuhwpO5OGXl+uy3t43uJyXZ6IFK80xzgOAdOB6flXiNcDst19XxnXJiIiMSjQjZzcPcfMbnX3/ymrgqRiUAtApPI6lftxfGlmT+Yf4/i5mekTRESkEjmV4FhCXktlZ/5vXcshIlKJnEpw/An4C9ATuAo4K6oViYhITAscHO4+xN2Xu/tG4L+AZ6JfVnTodFwRkegrdXCY2dVm9p6ZbTSzP5tZd3ffBYwow/p+FHVyKCISfUFaHNOB/0NeNyOzgMlmdrO7Z5VJZSIiEpOCnI67y93/lf/3SjNbA7xH3jEPERGpJIK0OLab2SNmVi3/cQ6wvwxqEhGRGBYkOBxIAXaY2dvAl0CambUsk8pERCQmlXpXlbvfDGBmNYB25HWr3gGYbWbN3f28silRRERiSaAuRyDSd9W6/B8REalkTuUCwApD13GIiERfXAeHruMQEYm+uA4OERGJPgWHiIgEouAQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCSTwleMiYXtl49dhlwBASqtGYZcgEoq4bnHoynERkeiL6xaHuy8DliUlJY0MuxYpG5OGXl+uy3t43uJyXZ5ILIrrFoeIiERfXLc4JP6pBSBS/tTiEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBBLXwaG+qkREoi+ug8Pdl7l7au3atcMuRUQkbsR1cIiISPQpOEREJBAFh4iIBKLgEBGRQBQcIiISiO7HIXKKdO9zqazU4hARkUDU4hCJAt37XCoTtThERCQQtThEokAtAKlM1OIQEZFAFBwiIhKIgkNERAJRcIiISCAKDhERCUTBISIigcR1cOgOgCIi0RfXwaE7AIqIRF9cB4eIiESfgkNERAJRcIiISCAKDhERCUSdHIpUcLqhlJQ3BYdIHNF9QaQ8aFeViIgEohaHSBxRC0DKg1ocIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOEREJBAFh4iIBKLgEBGRQBQcIiISiK4cF5GoUGeLlYdaHCIiEohaHCISdeqlN76pxSEiIoGoxSEiUacWQHxTi0NERAJRcIiISCAKDhERCUTBISIigVS44DCzgWb2BzN7zcyuCrseEZHKplyDw8yeM7PdZvbpccN/amYbzexLM5tQ0jzc/VV3HwkMB4aUYbkiIlKE8j4d9wXgaWBuwQAzSwCeAa4EMoC1ZrYUSAAeO27629x9d/7fv86fTkQkIha6Pon3bk/KNTjcfbWZNT1ucFfgS3ffCmBmLwMD3P0xoP/x8zAzA34H/MXd1xe3LDNLBVIBzj///OisgIiIxMQFgI2BHcc8zgC6lTD+GOAKoLaZ/cTdZxY1krvPAmYBJCUleZRqFZEKRF2flI1YCA4rYlixH/TuPhWYWnbliIhISWIhODKA84553ATYGVItIhJHKksLoLzFwum4a4GWZtbMzKoBNwFLQ65JRESKUd6n4/4JWAO0MrMMM7vd3XOBu4E3gM+BP7v7higt71ozm5WZmRmN2YmICOV/VtXNxQxfAawog+UtA5YlJSWNjPa8RUQqq1jYVSUiIhWIgkNERAKJhbOqRETiSixcvQ5ldwV7XLc4dHBcRCT64jo43H2Zu6fWrl077FJEROKGdlWJiJSheOz2JK5bHCIiEn1qcYiIlKF47PZELQ4REQkkroNDZ1WJiERfXAeHzqoSEYm+uA4OERGJPgWHiIgEouAQEZFAFBwiIhKIgkNERAIxdw+7hjJnZnuAr8KuIwbUA74Nu4gYoW1RmLZHYdoecIG71y/qiUoRHJLHzNa5e1LYdcQCbYvCtD0K0/YomXZViYhIIAoOEREJRMFRucwKu4AYom1RmLZHYdoeJdAxDhERCUQtDhERCUTBISIigSg44pyZnWdmq8zsczPbYGZjw64pFphZgpl9aGavh11LmMysjpktMrMv8v9HuoddU5jM7L/y3yefmtmfzKxG2DXFIgVH/MsF/q+7XwhcAow2szYh1xQLxgKfh11EDHgK+Ku7twY6UIm3iZk1Bu4Bkty9HZAA3BRuVbFJwRHn3P1rd1+f//d+8j4YGodbVbjMrAnwM2B22LWEyczOBHoBfwRw9x/cfV+oRYWvKlDTzKoCtYCdIdcTkxQclYiZNQUuBt4LuZSwPQncCxwNuY6wNQf2AM/n77abbWanhV1UWNz9/wFTgH8DXwOZ7v63cKuKTQqOSsLMTgcWA+Pc/fuw6wmLmfUHdrv7B2HXEgOqAp2AGe5+MZAFTAi3pPCY2VnAAKAZcC5wmpn9PNyqYpOCoxIws0TyQmO+u78Sdj0h6wlcZ2bbgZeBy83sxXBLCk0GkOHuBS3QReQFSWV1BbDN3fe4ew7wCtAj5JpikoIjzpmZkbcP+3N3fyLsesLm7ve5exN3b0regc9/uHul/Fbp7t8AO8ysVf6gvsBnIZYUtn8Dl5hZrfz3TV8q8ckCJakadgFS5noCQ4FPzCw9f9hEd18RXkkSQ8YA882sGrAV+EXI9YTG3d8zs0XAevLORvwQdT1SJHU5IiIigWhXlYiIBKLgEBGRQBQcIiISiIJDREQCUXCIiEggCg4REQlEwSEiIoEoOERCYGZXmNm8sOsQORUKDpFwdCDvymSRCkfBIRKODsCHZlbdzF4ws9/m948kEvPUV5VIODoAu4E3gNnuXll76JUKSH1ViZSz/G7uvwW+Au5w9zUhlyQSiHZViZS/NsBa8npgPRJyLSKBKThEyl8H4B3y7gfyvJk1DLkekUAUHCLlrwPwqbtvAn4F/Dl/95VIhaBjHCIiEohaHCIiEoiCQ0REAlFwiIhIIAoOEREJRMEhIiKBKDhERCQQBYeIiATy/wEBQQ9uKI9JhwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -361,7 +361,7 @@ "text": [ "\n", "After trial and error, I found that 43402 trials were sufficient to arrive at an 1-sigma error of around 0.001,\n", - "with mean 0.486091 and standard deviation sigma 0.001000. Wolfram Alpha tells me it should be approximately I =\n", + "with mean 0.487692 and standard deviation sigma 0.000993. Wolfram Alpha tells me it should be approximately I =\n", "0.487595, so I am happy.\n", " \n" ] @@ -392,10 +392,10 @@ }, { "cell_type": "markdown", - "id": "605f8777", + "id": "e8ab8be3", "metadata": {}, "source": [ - "After hours of thinking, I tried $$f_Z(x) = \\sec^2(\\pi{}x) = \\frac{1}{\\cos^2(\\pi{}x)}$$ on $(0, 1)$.\n", + "After hours of thinking, I tried $$\\sqrt{1.2x}$$ on $(0, 1)$. It works for the first half, so let's consider $$f_Z(x) = \\sqrt{1.2\\left( \\frac{1}{2} - |\\frac{1}{2} - z| \\right)}$$ on $(0, 1)$.\n", "It has cumulative density $$F_Z(x) = \\int_0^x\\sec^2(\\pi{}t)dt = \\left[ \\frac{\\tan (\\pi{}t)}{\\pi} \\right]_{t=0}^x = \\frac{\\tan (\\pi{}x)}{\\pi}.$$ This density gives us $$X' := g(Z)f_U(Z)/f_Z(Z) = \\sin(\\pi{}Z(1-Z))\\sin^2(\\pi{}Z)$$ for $Z$ on $(0, 1)$. $Z$ can be sampled using the inversion method as $$F_Z^{-1}(p) = \\frac{\\arctan(\\pi{}p)}{\\pi}.$$\n", "\n", "Unluckily, I found that the function is way too non-constant." @@ -403,7 +403,45 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 138, + "id": "8b01dfc9", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABioklEQVR4nO3dd1gUxxvA8e/QkSKi2AtiQcTeu2KLLbHEWKIxRpNoTDTWGGOLLVFTrTEaSxI1ll+MGms0omLHjqJYEBU7IL3D/P44JJQDTil3wHye5x693dmddwe493Z3dkZIKVEURVEUQ2Ok7wAURVEURRuVoBRFURSDpBKUoiiKYpBUglIURVEMkkpQiqIoikFSCUpRFEUxSCpBKYqiKAZJJShFURTFIKkEpRgEIYSfEKLjK257VQjRLmcjyrJOZyHEBSFEmBBijA7l1wkh5ib9P0fjTbm/7LRjVvvOSS/bfkrhpBKUkqGkD7tYIUSJNMsvCiGkEMJRT6GlIqV0lVIezuNqPwMOSyltpJSLX2ZDXePVNdnk1PFrqy8X2/aV208bIcQMIUR4mld00u9p/xyIV9EDlaCUrNwBBr54I4SoDVjqLxyDUQm4qs8AhBAm+qw/m165/bQdt5RytpTS+sULKA1cBPYA27ITqKI/KkEpWfkdGJLi/bvAby/eCCEmCSH+TLmBEGKJEOJHbTsTQlQQQmwTQjwTQgQKIZamWF1PCHFZCBEihNgshLBIsZ2LEOKwECI46bLTGynWpfrmn1EdQoiyQog/k5bfyezSUhb1HQLcgKVJ39Sra9m+vhDifNIlrM1AymNJG+9kIcSDpLI+QogOQojfgYrA30l1fJZi28lCiMtAhBDCRMuZT2MhhLcQ4rkQYm2adpRCiKop3q8TQszNor6OWbVJirITM/oZZtZ+Ou471XFn8rOzBHYBEcCbUsq4jMoqBk5KqV7qpfUF+AEdAR/ABTAG7qP59isBR6AMmg8Cu6RtTICnQEMt+zMGLgE/AFZoPrRbpajrDFAWsAeuASOT1pkCt4AvADOgPRAGOKeMM7M60HwZOwfMSNqHE+ALvKYlzkzrSypzGHg/g3YzA+4C45L21ReIA+Zqidc5qU3LJr13BKqkLZfmZ3IRqABYatmfH3Alab09cPxFvUnrJVA1xft12uLS8jugS5tk+DPU0kbJ7fcS+0513Jm0/V7gBGCt778h9creS51BKbp4cRbVCbgOPHixQkr5CDgKvJW0qAsQIKU8p2U/TdB8eE2SUkZIKaOllMdSrF8spXwopQwC/gbqJS1vBlgD86WUsVLKQ2i+IQ8kvYzqaAw4SM2loFgppS+wChigZR8vU582zdB86P4opYyTUv4P8MygbAJgDtQUQphKKf2klLez2P9iKeV9KWVUBuuXJq0PAua9RNyZ0bVNMvoZ5tS+MzzupLOqTUApoKuUMlzHY1MMlEpQii5+B94GhpLi8l4KvwKDk/4/OKm8NhWAu1LK+AzWP07x/0g0H1qgSTj3pZSJKdbfBcq9RB2VgLJJl5CChRDBaL6xl9Kyj5epT5uywAMpZcq5bO5qKyilvAWMBb4EngohNgkhymax//svsf5uUjzZpWubZPQzzIl9Z3jcQggjNGeDzkBnKWWIDvUqBk4lKCVLUsq7aDpLdEP7DeftQB0hRC2gB7Ahg13dByq+ws39h0CFpA+hFyqS4kxOhzruA3eklHYpXjZSym7ZrE+bR0A5IYRIs71WUsqNUspW/HfpdMGLVRltkkX9FdLU+zDF+0igSIr3pXXcb3bbJDO67juz+H5CcybWUUoZkAMxKQZAJShFV8OB9lLKiLQrpJTRwP+AjcAZKeW9DPZxBs2H93whhJUQwkII0VKHuk+juc/1mRDCVGiey3kdzeUcXes4A4Qm3Wi3FEIYCyFqCSEaZ7M+bU4C8cCYpE4MfdBcekxHaJ4Hai+EMAeigSg0l/0AnqC5V/ayPhZClBdC2KM5S9ycYt1F4O2k4+8CtE2xLrP6stsmmcnWvoUQ3wNdgQ5Jl5yVAkIlKEUnUsrbUsqzmRT5FahNxpf3kFImoPngqQrcA/yBLJ9RkVLGAm+g+RAKAJYDQ6SU13WtI8XyemjOBgOAX4Ci2akvk3j7oLkk+jzpGDPq6mwOzE+q5zFQEk1SAfgamJZ0SXKiLnUn2Qj8g6YTiC8wN8W6T9G0QzAwCM3Z7wsZ1pfdNslMdvaddNY+Ds2Z4FWR+jmoJ2nOypR8RqS+TK4or0YIURFNB4rSUspQfcejKEr+p75dKNmW9C11PLBJJSdFUXJKfn4SXTEAQggrNPcu7qLpYq4oipIj1CU+RVEUxSCpS3yKoiiKQdLbJb4SJUpIR0dHfVWvKIqiGIhz584FSCkd0i7XW4JydHTk7NnMei0riqIohYEQQutIK+oSn6IoimKQVIJSFEVRDJJKUIqiKIpBUglKURRFMUgqQSmKoigGSSUoRVEUxSCpBKUoiqIYJDUWn6LksUP3DnHU/2jye2NhTJ9qfXAt4arHqBTF8KgEpSh5JCYhhm88v2Gzz2ZszWyxMLYAIDwunG23tjGh4QQGuQwi9US8ilJ4qQSlKHngbuhdJh6ZyPWg6wx1HcqYBmMwNTIFICQmhGnHprHAcwGejz2Z3XI2Rc3TzaOoKIWOugelKLlsn98++u/qz6OIRyxtv5QJjSYkJyeAouZFWdx+MZMaTeKo/1H67+rPlYAreoxYUQyDSlCKkktiEmKYe2ouk45MoopdFbb22ErbCm21lhVCMMR1CL92/ZVEmcg7e99hw7UNqOlwlMJMJShFyQX3Q+/zzp532OyzmXdrvsu6LusoY10my+3qONRh6+tbaVW2FfPPzGf84fGExqpJipXCSSUoRclhB+8epN+ufviH+7PYbTETG09MdUkvKy8u+U1oOAH3++70/7s/3oHeuRixohgmlaAUJYfEJcSx4MwCxh0eh6OtI1tf34pbRbdX2pcQgqG1hrK2y1piE2N5Z887bPHZoi75KYWKSlCKkgMehT9i6P6hrL+2nkEug/it62+Usy6X7f3WL1mfra9vpXHpxsw5NYfPPT4nMi4yByJWFMOnEpSiZNOxB8d4a9db3A6+zXdtv+PzJp9jaqz7Jb2s2FvYs7zjckbXH80+v30M3D2Q28G3c2z/imKoVIJSlFeUkJjAkgtLGHVwFKWKlGJzj810duycK3UZCSM+rPMhqzqtIjgmmIG7B7Lbd3eu1KUohkKnBCWE6CKE8BFC3BJCfK5lfVEhxN9CiEtCiKtCiPdyPlRFMRyBUYGMODiClZdX0qtqLzZ020Al20q5Xm+TMk3Y+vpWXOxd+Nzjc+acnENsQmyu16so+pBlghJCGAPLgK5ATWCgEKJmmmIfA95SyrpAO+A7IYRZDseqKAbhwtML9Pu7HxefXmR2i9nMbjkbCxOLPKu/ZJGSrH5tNe/Veo8tN7YwZO8QHoQ/yLP6FSWv6HIG1QS4JaX0lVLGApuAnmnKSMBGaAYRswaCgPgcjVRR9ExKye/evzNs3zDMTcxZ3209vav11kssJkYmjG84nkVui7gXeo9+f/dLNQCtohQEuiSocsD9FO/9k5altBRwAR4CXsCnUsrEtDsSQnwohDgrhDj77NmzVwxZUfJeRFwEE49MZKHnQtqUb8OmHpuoYV9D32HRvmJ7NvfYTBmrMnz878csubCEhMQEfYelKDlClwSlbWjltA9jvAZcBMoC9YClQgjbdBtJuVJK2UhK2cjBweElQ1UU/bgdfJuBuwdy8N5BxjUcx49uP2Jrlu7XW28q2FZgfbf19Krai5WXVzLq31E8j36u77AUJdt0SVD+QIUU78ujOVNK6T1gm9S4BdwB9P/1UlGyad8dTbfukJgQfun8C8NqDTPI6TAsTCyY3WI2M5vPxPOxJ/129cPrmZe+w1KUbNElQXkC1YQQlZM6PgwAdqYpcw/oACCEKAU4A745Gaii5KW4RM2oEJOOTsK5mDNbemyhcenG+g4rU0II+lbvy+9df8cII97d9y5bb2xVo08o+VaWCUpKGQ98AuwHrgFbpJRXhRAjhRAjk4rNAVoIIbyAf4HJUsqA3ApaUXJTQFQA7+9/P3lUiDWvraGUVSl9h6Uz1xKubO6xmSalmzD75GxmnJhBdHy0vsNSlJcm9PXtqlGjRvLs2bN6qVtRMnLh6QUmHJ5AeFw4M5vPpLtT99ypKCFFJ1chwMg456tITOCnSz/x8+WfcbF34Qe3H3Jk+CVFyWlCiHNSykbplqsEpSiaLuQbr2/kW89vKWtdlh/cfqB6seo5XQncPACHv4aH51OsEODcFdp9DmXq5mydwJH7R5jiMQUjIyMWtl5Ii3ItcrwORckOlaAUJQNR8VHMPjmbXb67aFe+HfNaz8vZXnpSwu1/wf1reHAW7CpCnQFgnPQse3QwXPgdokOgRg9oNwVK18q5+oF7ofcYe3gst57fYnT90bxf+32D7OyhFE4qQSmKFvfD7jPOfRw3nt9gVL1RfFjnQ4xEDg5R6XsE3L+C+6egaAVoMxHqvg0maQZaiQqG0yvg5DKICYWavTSJqmTOdYaNjIvky5NfsvfOXtpXaM+8VvOwNrPOsf0ryqtSCUpR0jjx4ASTjk5CIpnfej5tyrfJuZ3fPQnu88DPA2zKQpsJUH9I+sSUVtRzTZI69RPERkDtvppEVbxKjoQlpWT9tfV8d/Y7KtpW5Ee3H3Eq6pQj+1aUV6USlKIkkVKy+spqFp9fTNViVVnUbhEVbCtkvaEuHpyHQ3M1l/SsSkLrCdBwKJi+5Fh9EYFwYhGcWQXxMVBvILSdrLk8mAM8H3sy8chEYhJimNdqHh0qdsiR/SrKq1AJSlHQXOaadnwaB+4eoItjF2a1mEUR0yLZ3/ETb80Z0/VdYGkPrcZB4/fBLJv7Dn8Kx34Az9UgEzXJrs1EsCmd7ZAfRzxmrPtYrgZe5cM6H/JxvY9z9vKmouhIJSil0Lsfep8x7mPwDfFlbIOxDHUdmv2OAoG3Nb3yvP4H5rbQYjQ0GwnmNjkT9AshD8DjWzj/GxiZQpMPNEmwiH22dhuTEMPcU3PZfms7bcu35evWX2NjlsOxK0oWVIJSCrVjD47x2dHPMBJGLGyzkBZls9nVOuQBHFkAF9aDiTk0HalJTtlMGFkKugOH58PlzWBmDS0+geYfZyshSinZ7LOZBWcWUN6mPIvcFuFkp+5LKXlHJSilUJJSsvbqWn489yPVilVjkdsiytuUf/UdRgaBx3eae0MyERq9B60ngk0ejzTx9JrmXtf1XVCkuCaGRsNe/l5XCueenGP84fHEJMTwVauvaF+xfQ4GrCgZUwlKKXSi4qOYeXwme/32Zv9+U0w4nFoOxxdDXITmOaZ2n0Ox3J9FN1MPzsG/s8H3MNiWB7cpUHfgK49MkfK+1Kh6oxhRZ4S6L6XkOpWglELlQfgDxrqPxSfIh08bfPrqo5DHx8L5X+HIQoh4qnmQtv30HH0+KUf4HoGDX2pGqHCooYmxRnfNMEovKTo+mtknZ/O37990qNiBea3mYWVqlfMxK0oSlaCUQsPzsScTDk8gPjGehW0X0qpcq5ffSWIiXN0Gh+bAcz+o1Ao6fgkVDHhEcynh2k74dw4E3oTyTaDTLKj08vfbUj4vVbloZRa3X0wFmxzqiq8oaWSUoNS5u1KgbL6+mQ//+RA7Czs2dt/4asnptjusagd/Dtd0RBj0Pxi6y7CTE2jOlmr2hFGn4PXFEHIf1naFjQM096xealeCd2q+w08df+Jp5FMG7h7IqUencilwRdFOnUEpBUJcQhxfn/marTe20qZ8G+a3nv/y3aUfXYaDM+H2IShaEdpPg9pvgVE+/R4XG6kZPunYjxAbBvXeBrepYFv2pXbzonv+nZA7TGo8ibdrvK3G8VNylLrEpxRYQdFBjHMfx/mn53m/9vt8Uu8TjF+mk0DwPTg0T9N129IO2kzSPGRrYv7SscTGJxIUEUtAeAzBkXGERccRFh1PWEw8MfEJxMYnEhufSEJi6r87E2OBmbEx5qZGWJgYYW1hio2FCTYWJthbmVHcyhx7KzOMjV4hMST3PFwJwhiaj4KWn4JFUZ13EREXwRSPKbjfd6dPtT5MbToVM+Mshm1SFB2pBKUUSD5BPow5NIbA6EBmt5hNN6duum8cFQzHvodTKzTvm32kefjV0i7DTaLjErgXFMmdgAjuBETg/zySh8HRPAyO4lFINCFRcVlWa2osUiUaKSE+UaZLWmkJAcWtzChrZ0mZohaUtbOkon0RHEtY4VTCinJ2lpgYZ3K299xP0zXda6uma3rbzzXd5I1Ns4wZIFEmsvTCUlZ5raJ+yfr80O4HilsW12lbRcmMSlBKgfPv3X+ZcmwKNmY2LHZbjGsJV902jI+Fs2s0D9pGBWm6jLefBnb/dQKQUvIwJBov/2C8H4Xh8zgUn8dh3A2KJOWfjF0RU8oWtUxOGg425hS3NqOEtTl2lqbYpDgTsjA1xszYCKMMzoLiExKJTUgkKjaB8Jh4wqLjCY2O43lEHIERMQSEx/IsLDo5IT4IjiIyNiF5e1NjQRUHa5xL2+Bc2oaaZWypU94Oe6s0ZzoPL8A/0zUD2dpX0XSkqNFD5x5/++7sY/rx6RSzKMbi9oupYW9gPRqVfEclKKXAkFKyymsVSy4soU6JOvzo9iMORRx02RCu74YD0yHIFyq3gU5zoGw9ouMSuHQ/GE+/IM7fC+ayfzAB4bEAGAlwLG6Fc2kbqpeywcnBisolrHAsYYWthW5nH7lBSklAeCx+gZqzudvPwrnxOAyfx2E8DPlvivcK9pbUKW9Hw4rFaFLZHpcythgL4MZ+ODADAnygYgt4bR6Ua6BT3VcDr/LpoU8JjQ3lq1Zf0bFSx1w6SqUwUAlKKRCi46OZcWIGe+/spbtTd2a1mIW5sQ73ih5egP1T4e5xKOFMfMdZXDBvgsetQE7cCuCyfwixCYkAVC1pTd3ydtStUJQ65e2oUdoGC9Ocn5I9N4VExXH1YQhe/iFc9g/h4v1gHgRHAWBtbkLDSsVoVbUErarY4fxwO0aHv4KIZ1CnP3SYAUWzHm3jWeQzxrqP5XLAZUbXH80HtT9QnSeUV6ISlJLvPYt8xphDY7gaeJUxDcYwvNbwrD8QQx9qRlq49AeJlsU55zSKXyJbc+z2cyJiEzASULu8Hc2c7GlcyZ5GjsWwK1Iwb/4/DI7C0y8IT78gTvkGcetpOAAlrM3oVKUIw9lOlVu/atq0xWhoORbMM5/QMCYhhpknZrLbdzddK3dldovZWJi8+nBLSuGkEpSSr3kHejP60GjCYsOY33p+1uPExUbCiSUkHvsBmRDPdouezHzehXCKUM7OErcaDrSq6kBzp+IULaK/y3T69CgkiuO3AvG4+YzDPs8IiYqjknEAXxf9ixaR7iRalcKo4wzNDMCZdLV/Mb/WovOLqFOiDovaL6KEZYk8PBIlv1MJSsm3Dt49yBfHvqCoeVGWtl+Ks71zxoWl5PmZPzA99CXWMU/YndCE+fEDsS9Xnc6upengUhLnUjbqUlQa8QmJnL8XzL/Xn7D/ymPsgy4y3XQ99Y1uEWTrgtnr32BdrXWm+3jRaaWoeVGWtF+iOk8oOlMJSsl3Un0zd6jDIreMv5lHxMRz0uMfKpyZg3OsN1cSHdlY7CMqN+xM19qlKV8sByYlLCSklFx7FMaeyw+IubCZ96J/pawI4qy1G9HtZtKsft0Mu7NfC7zG6EOjCY0N5evWX6uZehWdqASl5CuxCbHMOjmLnbd30q1yN2a3nK21M8SVByHsPH4Bl6s/0FscJoiinKs6muqvjaCSg60eIi9YpJRcufOIwH++odnj9SAl6417EdF4NH2aVqOCffrE/yzyGZ+6f8qVgCuMazguZyaGVAo0laCUfON59HPGuo/l/NPzfFzvY0bUGZHqAy4mPoHdlx+x/vhNGj3ezBiT7ViIOJ7VGkbpHtMQLzFCgqK7uKC7BPz1OWXu7+GBLM7X8YOIqNKDoa2caF21RKrnu6Ljo5l2fBr7/fbTp1ofpjWdhqmODwQrhY9KUEq+4Bviyyf/fsKTiCfMbTWXrpW7Jq97FhbD76fusvH0XVwjPZlrsZ4KiQ+Iq9IJ067zoURVPUZeiPgdJ3bXJMwCrnIOV6bGvENcCReGtnDkzYblKWJmAmhGnlh2cRkrL6+kSekmfN/ue4qaqy8PSnoqQSkG7/Sj04w7PA5TI1MWt19MXYe6ANwNjGDlUV+2nvOndMIjFtltoX7USaR9FUSX+VC9s54jL4QSE+D8r8h/Z0NUCLstuvNF8OsYFynGuy0cebe5I8WSRrD4+/bfzDwxk3LW5VjWYRkVbSvqOXjF0KgEpRi0v27+xeyTs3Es6sjSDkspZ12OG0/CWPzvTfZ4PcLaKI4fyh3CLXATRkYm0HYSNBv1SgO6KjkoMgjc58HZNcSZ2bHe5j1m+9fHwtSUgU0qMrKtEyVtLTj35Bxj3ccCsMhtEQ1K6TZihVI4qASlGKREmciSC0v4xesXWpRtwbdtv+XRc1iUlJiszIz50vkevR8vxjj0PtTqC53nvPSUEUoue3QZ9kyC+6eIKtWAn4p8xDIfa0yMBIOaVmJkWyeiecrH/37Mg/AHzG45mx5OPfQdtWIgVIJSDE5MQgxTj01lv99++lbvy+CqY1l00Je/Lz+kiKkxYxua8W7Icsx8D4CDC3T7Bipn/iyOokdSwqVNmrEOIwIIqz2Eb+L6seFyKCZGgqEtHBnUogQzT33G2SdntXaAUQonlaAUg/I8+jljDo3h4rOLjKj1KU/vN2XjmfuYGAs+aF6WUaa7sTz1o2YqiHZToOkInaeFUPQsKhgOf62Zf6pIcZ41n8bX/nX569JDbMxN+LBtJe4Z/cpev130rNKTmc1nqh5+hZxKUIrB8AvxY9S/o3ga+ZQ2xcaw/3RJouMTGdC4ApOqPMDu8BTNaOOufTQjbKvLefnTo8uwezz4e0KlltxuMouvPOHf608pZWtO0wbncH+ynqalm/K92/fYmqnn1gorlaAUg3Dh6QVG/zua+ETgyTAePy3Fa66l+KJ1MSp5zoWr26B4Vej2LVRx03e4SnYlJsLF9ZppPWLCoMVozlb6gFn77uD1IIQqla8TaLkeR9tKLO+4nLLW6stIYaQSlKJ3+/32M8XjC0S8HUG+7+Jc3JEZ3ZxpEbxTM+J4fAy0maiZjlz1zitYIgI1SerierCrSGKXb/gz3JWF+30ISvDGttIGbC0s+anjcmoWr6nvaJU8phKUojdSSlZdXsuSiz+QEOmIScAwJnVqwNsVQzDePRYenAMnN+j+HRSvou9wldzkd1xz2e/ZdajZk4j281jiGcHq0yewqLAOU9Movm/3Le0qttV3pEoeUglK0YuExATG/PMlR59sJy60Nl1LjWNqxyqUOPsDnFwGRezhta+hdl+dpxxX8rn4WDixGI4s1Jwpd5jBzQpv8fnfJ7mWuAhji4eMcJ3EJ43f0XekSh7JVoISQnQBFgHGwC9SyvlayrQDfgRMgQApZaZfgVSCKvgeh4YyaMenPE08S5Go9izqPJ1miZdg11gIvgcNhkDHWZokpRQ+gbdh1zi4cwTKN0G+vojN90356uwXSMvr1LXuy+o3pmGez2YzVl5eRgkq41nI/tvQGFgGdAVqAgOFEDXTlLEDlgNvSCldgbdyImgl/9p64TqdNw3iScI5mti+x9GBX9Ds4hewvg8Ym8PQPfDGEpWcCrPiVWDIDuj9MwTeQvzchgGhWzjY/2fKGrXlUvj/aL12BOfuBug7UkVPskxQQBPglpTSV0oZC2wCeqYp8zawTUp5D0BK+TRnw1Tyi5DIOD784x++9PwIaebPhLpzWF2lPOYrmsGVbdDmMxh5DBxb6jtUxRAIAXUHwCeeUKsPHF1IyfWd2df6bbqWe5co89O8s/tDFuy/RFxCor6jVfKYLgmqHHA/xXv/pGUpVQeKCSEOCyHOCSGG5FSASv7hcfMZHZdt5ETkl1haRLOm5Ve8570etn0A9k4w0gPaTwVTC32HqhgaqxLQZyUM/hPiYxDrurIw5jFT6k3GxOo2v96ZzBs/7eXW0zB9R6rkIRMdymi7c532xpUJ0BDoAFgCJ4UQp6SUN1LtSIgPgQ8BKlZUIxoXFDHxCSzc58O6C/9gVWE9JSzsWFO+M05b3weZCF3mQ5MPwUjdS1CyULUjjDoJh+bA6Z9526cc5VuNYKzPb/jHLqT7T4FMe60tg5tWVEMkFQK6nEH5AxVSvC8PPNRSZp+UMkJKGQAcBeqm3ZGUcqWUspGUspGDg8OrxqwYEN9n4fRZfoJfL/+FVcV1VLUtw+ZoE5wOzoPyjTQfNs0+UslJ0Z25NXRdAMP/ATMr2uyZzq+WLhQtkkCRSj8xc99eRq4/R3BkrL4jVXKZLgnKE6gmhKgshDADBgA705TZAbQWQpgIIYoATYFrORuqYmj+POdPjyXH8E/Yh2W5TTSyLM2v189R6tFVeH0xvLMdijnqO0wlv6rQBEYchVbjqX11NxsfPaO0hRlFnX7h8L2jdFvkgadfkL6jVHJRlglKShkPfALsR5N0tkgprwohRgohRiaVuQbsAy4DZ9B0Rb+Se2Er+hQdl8Dnf15mwtYLlKx0gET7nXSSlqzwPolNpZYw6hQ0fFc916Rkn6kFdJwJH/xLRYvi/H7zMlWFCZYVfkNan2XAylOsOuqLvp7nVHKXelBXeSn3AiP5aMM5rj58Tt36+/CN9mBAeBSfh8Zi3HU+1B2oEpOSO+JjweNbwj2+Y2yZ0pw2M8JR9MfLuz5dXEuz8K062FqoUdHzo1d+DkpRXnC//pQeSzy49zyY5o034xvtwcfPg/nCti7GH5+Cem+r5KTkHhMzcPsC6/cPsjzOhtfCI/CTm2nfxIMD1x7zxpJj3HiievkVJCpBKVmSUrLM/RbDfvWkrH0itZyX4B12iRnPIxjZZh5i8P/UlBhK3ilbH7MPj7Cg+hD6h4bjGbabbo23Eh4bS+9lx9l/9bG+I1RyiLrEp2QqMjaez/53mV2XH9GnZhS+cj73ZDQLcKBTz3VQrJK+Q1QKMXnvNCv2fshyi0RamTjw9PlnnPNPYFzH6oxuXxUjI3VGnx+oS3zKS3sUEsVbK06y2+sRcxp4czVuJo8So/ipQk86DflXJSdF70TFpnz03gmmWtfkeNxTbCw/Y4RLID8cvMGoDeeJik3Qd4hKNqgEpWh15UEIvZYdJyAwgN9r/M4v4WuJMjZmTetvaNphHhipXx3FQJhZMeDNzSyoMZRLxomcj5nH6pruHPB+SP+VJ3kaGq3vCJVXpD5llHT+ufqYt1acpB43WGw/g0mJXpiZ2bDu9a24Vu2m7/AURauuzSayuPVC/MzN+TFmF7vKf0PEE196LTvO9ceh+g5PeQUqQSmprD1+h1HrzzDDegeDTObwcVEoYV2W33ptx6m4s77DU5RMta7SjZVd1hFoYc1oy0B+tp5Gx/gj9P3pJEdvPNN3eMpLUglKASAxUfL13mus2eXOP7ZfUcxoO5+WcsCpWHV+fX0zZazL6DtERdFJ/VINWNt9A3GW9gwvVYyBxstZYraUMeuO8NcFf32Hp7wElaAU4hISmbj1Eo88fueA5VQuFHnG5JIO1ClZn9Vd12FvoeZsUvIXZ3tnfu32O2ZWJRleviJFjT35x/ILftvyP1Yeva1GnsgnVIIq5CJj4xm11oMWV6az2GwZW8o78aVdEZqXbcGKTiuwMbPRd4iK8kocizryW9ffKG5TjhHlynHTzpz/mc8ieP8C5u26QmKiSlKGTiWoQiwkKo4ZKzbyxb0P6W18jJ/rdecb41A6VerE4vaLsTSx1HeIipItZazLsK7LOioVrcxoO1MOO7flM9PNuJ0ZwdxNh0hQScqgqQRVSAWERbNh8TS+ChxH6SKJLGozjKUhXrxR5Q0WtlmImbGZvkNUlBxR3LI4q19bjYu9CxPi/NjV5iOamN5mlM9Qlv+ykth4NVOvoVIJqhB6/OQx3ot6MSpqBSHlWvJ9qwGsvf8P/Z37M6flHEyMdJnHUlHyj6LmRVnZeSUNSjXgi/t72NF9OkbWJRn98DP+WfwRUVHqWSlDpBJUIfPk+kkSV7ShedwZ7jT8nMU167Dp9naGug5latOpGAn1K6EUTFamVizvsJyW5Voy68rP7O06mtsV3qRH6Cbu/tCeqMB7+g5RSUN9GhUWUvL88HKKbeqBkPHcfn0zK4oGs/32DkbVHcX4huPVFNpKgWdhYsFit8V0qNiB+ee+40jTVpxtuJAKMbeIXdaKqOsH9B2ikoJKUIVBTDgRf7xHscNTOE1tnr29n5/D9rDnzh7GNhjLR/U+UslJKTRMjU35pu03dHXsyg/nfsCzYiKnOvyPJ/E2mG96i9iD8yBRjeFnCFSCKuie+RC3oh0WN3awWLyN9bDNrLy/iIP3DjK58WSG1x6u7wgVJc+ZGpnydeuv6VmlJ8svLueq9UVu99zB9oRWmB1bSPzvfSFSTSevbypBFWRX/yJxpRvhz58wUkyjzbDZ/HxjJkf8jzC92XQG1xys7wgVRW+MjYyZ3XI2b1V/i1Veq7gmd2LR92emxr+PvONB4s9t4MF5fYdZqKkEVRAlxMP+qbB1KNcSyvMWCxj53mCWXfuCEw9OMKvFLPo599N3lIqid0bCiOnNpjPAeQBrr67FK2Y9jd8cR9/YGQSFxyDXvAbnftV3mIWW6k9c0EQEwNah4OfBdtNufBkziJ/ea8Dya1/g+diTOS3n0LNqT31HqSgGQwjBF02/wMTIhPXX1jPAOZ5Bvd+h058OrLdbhevfY+DhBei6UDPtvJJnVIIqSB5ehM2DkeFP+a7IWH4Ja8bKIbVZ6TOV80/P81Xrr+jh1EPfUSqKwRFC8FnjzzA1MmXt1bX0d5aM6zmI13dYs6LMbjqfWwtPvaHfb2BTWt/hFhrqEl9BcXkrrHmNRJnI50UX8nNIU5a87cqa29M4//Q8X7f6WiUnRcmEEIJxDccxrNYwNvts5p7YwOddXfjw0etsrDgL+dgLfm4L/mf1HWqhoRJUfpeYAAdmwLb3SSzbgE9tf2DLIwcW9K3BhrszuPD0AvNbz6ebk5poUFGyIoRgbIOxyUnqidkffNS2Ml/cqMZal1VgagFru8HFP/QdaqGgLvHlZ9Eh8Of7cPMfZKPhjA8dyN+XnzKrVzV2PpmdnJy6Vu6q70gVJd94kaQEgtVXVvNWdcnAJr2YfeY+Jp3WMsR/JmwfCU+uQMdZYKw+RnOLatn8KsgXNvaHIF9k9++Z/bgZ2y/7Ma6zI0dC5qvkpCjZIITg0wafArD6ymr6VTeiW+2uzDjwCJu+i+nt8BOcXArPrkPfNWBRVM8RF0zqEl9+5HccVnWAiGcwZAero91Ye9yPd1qU4XLc95x7co6vWn2lkpOiZMOLJDXUdShbbmymXJX9NHOyZ9K2a3hUmwSvLwLfw7C6MwTd0Xe4BZJKUPnNhfXwW08oUhze/5e/Q5yYu/sar9Wy57HFCjwfezK35Vy6O3XXd6SKku8JIRjfcDyDXQazyecPatU+QpWSVny0/jxXy/SGd/6CsMfwSwe4e1Lf4RY4KkHlF4mJcPBL2PExOLaE9w9wMtiOCVsu0dDRGkr+yulHp5jdcjavV3ld39EqSoHxogv62zXeZvONDbRqchprC2PeW+uJv10jeP9fsLCD396Ay1v0HW6BohJUfhAXDX8Oh2M/QMOhMOh/3AozZcTvZ6lQ3JySVbZy4tFxpjefTq+qvfQdraIUOEIIPm/yOf2q92PLzd/p1tqLqLgE3lvrSah1JXj/IFRoCts+gCPfgFQz9eYElaAMXUSg5pLe1W3QaTb0+JGgaMmwdZ6Ymkiq19rJsYdHmNJkCm9Vf0vf0SpKgSWEYGqzqfSq2outt9fQq9117gRE8PGG88Sb28HgbVBnALjPhZ2fQEKcvkPO91SCMmRBd2B1J80wK2+tg5afEpOQyIjfz/I4NJJGjQ5y7NG/TGw0kbdd3tZ3tIpS4BkJI75s/iXdKndj+91feL31TTxuBvDl31eRxqbQewW0/Vxzr3hDX4gO1XfI+ZpKUIbq4QVNcooKgnd3gmtvpJR8/qcXnn5BtG52lBNP9vNJvU941/VdfUerKIWGsZEx81rNo1OlThx8+gsdGvuy/tQ91h73AyHAbQr0XA53PGBddwh7ou+Q8y2VoAzRrX9hXQ8wsYBh/0DFZgAsc7/FXxf8adH4BGcC9zC81nBG1B2h52AVpfAxMTJhQesFtC7XGs/wVdR38WXubm/cfZ5qCtQfBG9vgcDbmi+aAbf0G3A+pRKUobm8BTb2g2KOMPwAOFQH4ID3E7795wa1a53GK/xvBtYYmPwgoaIoec/U2JTv231Po9KNuCNWU7GCL2P+uMDtZ+GaAtU6wtC/ITZCk6T8z+k34HxIJShDcmqFphdQxebw3h6wLQPAjSdhjN10gUqVPfFL2E6vqr34vMnnapp2RdEzCxMLlrRfgmsJV0Js1mJidZMPfj1LSFRSB4lyDWH4P2BhC7++Drfd9RtwPqMSlCGQEty/gn2ToUYPGPS/5KFTgiNj+eC3s5gVO0uQxZ90rtSZL5t/iZFQPzpFMQRWplYs77CcKkWdMCr9K/6R3oz54wIJiUldzYtXgWH7NVdFNvaDq9v1GW6+oj7l9C0xEfZMgiMLoP5geOtXzYjJQHxCIqP/uMCT+DMk2G+lZdmWzG89H2MjYz0HrShKSkXNi7Ki0wpKW5XEtvJveNy9zMJ91/8rYFMa3tsNZetrJhQ9t05foeYrOiUoIUQXIYSPEOKWEOLzTMo1FkIkCCH65lyIBVhCPOwYBZ6roMVoeGNpqpGRv/3nBicensC87CbqlazL9+2+x9TYVI8BK4qSkRKWJVjVeRVFLaywr7KOlSfPsMfr0X8FLItphkaq2gH+/hROLtNfsPlElglKCGEMLAO6AjWBgUKImhmUWwDsz+kgC6T4WM3oEJf+ALdp0GmOpotqkn1XHrHytDs2FTdQrVgVlnZYShHTInoMWFGUrJS1LsuqTquwNDPCzmktk7Yd5dbTsP8KmFnBgD/A5Q3Y/wUc/UZ/weYDupxBNQFuSSl9pZSxwCagp5Zyo4E/gac5GF/BFBcNW94B7+3QeR60nZQqOd16Gs7E7f9g4/grZaxLsqLTCmzNbPUXr6IoOnOyc2JFx58wNYvCuOwvvL/eg7DoFKNKmJhB37VQpz8cmgsHZ6mhkTKgS4IqB9xP8d4/aVkyIUQ5oDewIrMdCSE+FEKcFUKcffbs2cvGWjDERcGmgXBjH3T/Dlp8kmp1eEw8H2w8gCjzC3YWFqzs/DMlLEvoKVhFUV6FawlXFrn9iJF5AE8tf2LCVk9kyiRkbAK9VkCDd+HY97B/qkpSWuiSoLT1ZU7bkj8Ck6WUCZntSEq5UkrZSErZyMHBQccQC5C4KPhjoKarac9l0Pj9VKullEz433GeFlmChXkcKzv/TAWbCnoKVlGU7Ghetjnz23yNcZG7eIT8wM9Hb6YuYGSkmVOqyQg4tUxzyU8lqVR0mVHXH0j5KVkeeJimTCNgU9JzOSWAbkKIeCnl9pwIskCIjdScOfkegV7LoV76sfPWnLyBR9gCzCyD+anjSpztnfUQqKIoOaWLYxeCo4OZd3oeiy59TaNK39HI0f6/AkJA1wWaf08t1ySoLl+nuuRfmOlyBuUJVBNCVBZCmAEDgJ0pC0gpK0spHaWUjsD/gFEqOaXw4rKe7xHo9ZPW5HTpfiDfXZiOsYU/37RdSKPSjfQQqKIoOW1AjQEMrfkBJkXP8uGueQRFxKYuIAR0mQ9NP4LTP8G+KepMKkmWCUpKGQ98gqZ33jVgi5TyqhBipBBiZG4HmO/Fx8CmQZrk1HsF1BuYrkhIVCzDdk3B2Poa4xpMppNjRz0EqihKbhnfaDRuZXsQZ/MP72z5kcTENAlICM2ZU7OPNUnqwAyVpNDtEh9Syj3AnjTLtHaIkFIOzX5YBURCHGx9D27/q3nGqe6AdEWklLy9dTaxRU7yesV3GFZnkB4CVRQlNwkh+L7DHN7c9oTb4euZvK8c33QblLYQvDYPEmLgxGIwtQS3L/QTsIFQI0nkloR4zbh6Pruh27fQ4B2txT775xfuyR24WLdnXrtJeRykoih5xcTIhI1vLMHWqDJ7n3zL5sse6QsJAV2/0Ywqc2QBeHyf94EaEJWgckNiIuwcDVf/0jyA2+QDrcU2XTnA3kdLsZWu/N7zGzX4q6IUcFZmVmx6YxXGiXbMOzsJ72e+6QsZGcHri6FWX/h3lmYQ6UJKJaicJiX8Mw0ubdTMrNlyjNZil55489XZKRjFlWZjz2WYm5jlcaCKouhDRbuSfN1yEYlIhu75kKCooPSFjIw196xr9NAMIn1pc94HagBUgsppx77XPNPQZAS00z5s4eOIx7y//yMS4s2Z2/wHKhUrnsdBKoqiT91q1KF32WlEJgYy6O+PiEmISV/I2BTeXA2OrTVjdt4ofKPIqQSVk86uhX9nQ+23NN1GtVyyi4iLYMjuEUTFR9C91HTeqJVuWENFUQqBWZ17UClxOP5R3nx6cDKJMjF9IVMLGLARSrnCliFw92TeB6pHKkHllGt/w65xULWT5lkno/RNG58Yz5h/J/Ao0o/SMR8yt2tnPQSqKIohMDISrHnrfYye9+D4439ZfH6p9oIWtjDoTyhaHjb2hyfeeRuoHqkElRPun4E/39fMntnvN82puRbfeH7DmSfHSXjWk5V938bMRDW/ohRmpWwtWNBxDLHBjVh9ZRU7b+/UXtDaQTNVh6klbHgLQh9pL1fAqE/I7Aq8rflWY1MG3t4MZtqnxNh4bSMbr28kNrAVX7QejpODdR4HqiiKIepSqwxvlB1NfEQVZhyfydnHZ7UXtKsIg7ZA1HPY+BbEhGkvV4CoBJUdEYGwIWluxsF/gpX2UcePPzjOgjMLSAh3oYX9UAY1rZiHQSqKYuhmvlGH4hEfIOPs+dR9LPdD72svWKYu9PtVc5lv61DNYAAFmEpQryouGv4YAKEPNWdOxatoLeYb4svEIxMxTiiDxfMhLHyznnreSVGUVKzNTVjUrwXhd98lKjaBTw59QnhsuPbC1TpBj+/h1kHYM7FAD4mkEtSrkFIzZbP/Gej9M1RoorVYSEwIo/8dTXyCEc99B7OgT2McbMzzOFhFUfKDhpWK8XGrpoTcHYhfyF0+O/oZCYkZzGDUcCi0Gg/n1sHpn/MyzDylEtSrOLEYLm8Ct6ng2ktrkfjEeCYemciD8Ic8v/M2b9atTaeapfI2TkVR8pXRHapRw64+IqgXHg88WHR+UcaF208H5+6wfwrc+jfvgsxDKkG9LJ99cGAmuPaBNhmPnfeN5zecenQKy9B+lDStwYzX1fNOiqJkztTYiO/61SUyoCllRHvWXl3Ljls7tBc2MoI+K6FkTc2g1AE3tZfLx1SCehlPr8Gfw6FsPc2MuBncS/rr5l9svL6R6hbdeeRfhwV962Brob3ruaIoSko1StsytlM1bnh3wMm6HrNPzuZKwBXthc2tYeAfmkdbNvbX9PArQFSC0lV0iGZeJzMrzZPdGXQnv/zsMnNOzcHFriHnL7ZgUNOKtK5WCKe3VxTllX3Y2ol6FYpzx7s3xSyK86n7pwREBWgvbFcR+q+H4Lvw10jNYNUFhEpQupASto+C537w1jqwLau1WEBUAOPcx1HC0oFHN/tSvpg1X3RzydNQFUXJ/0ySLvXFxFjiEDmC0JhQJhyeQFxG3corNYfXvoIb++DYd3kbbC5SCUoXxxfB9V3QeQ5UaqG1SFxCHOMPjyc0NpTaJp9yP0Cw4M06WJnrNCekoihKKlUcrJnY2ZmT1yx4vdynnH96ngWeCzLeoMmHmik6Ds2D24fyLtBcpBJUVnyPaOZkce0NzUZlWGyh50IuPL3Ae86T+eu05O2mFWlRRfuDu4qiKLoY1qoydSvYsf1YKfpXe4fNPpsz7jQhBLyxGBxqwP+GQ3AGD/vmIypBZSbsMfxvGBSvCm8sybBTxC7fXWzy2cSgGu+wzcOBUrYWTOlaI4+DVRSloDE2EnzTtw5h0XE8uuNGk9JNmHNqDteDrmvfwMxKcz8qIQ62vpvvR5pQCSojiYnw1wiIjdAMAGtuo7XYzec3mX1yNg1KNsA4uDu3nobzVZ/a2Khee4qi5IDqpWwY3b4au72e0q3URIqaFWWc+zhCY0O1b1CiKvRcCg/Ogfu8vA02h6kbJBk5uQR8D0OPH6Ck9o4O4bHhjDs8DitTKz6sMZN3V12nT4NyuDmXzNtYk8TFxeHv7090dLRe6leUgs7CwoLy5ctjapq3X0A/aleFvVces2D3A75/ZwGjD3/AVI+pLGq/CCOh5TzDtRfcHgLHfgQnN3Bqm6fx5hQh9TSOU6NGjeTZsxmM2qtvD87D6k7g3BX6/a710p6UkvGHx+N+351VnX5hzp/RPAyO4uD4ttgV0c/07Xfu3MHGxobixYur8f4UJYdJKQkMDCQsLIzKlSvnef1XHoTwxtJj9G9ckVoul5l/Zj5j6o/hgzofaN8gNgJWtoPoUPjoBFgZ7szdQohzUspGaZerS3xpxYRpHsa1LgWvL87wvtPG6xs5eO8g4xqOw+t2cS77hzDjdVe9JSeA6OholZwUJZcIIShevLjerlDUKleU4a0q88eZe1Q170LXyl1ZenFpxtNzmFlB3zUQFQQ7Ps6Xg8qqBJXWvima5536rIIi9lqLXAm4wrdnv6Vd+Xa0L9OXb//xwc3ZgdfrlMnbWLVQyUlRco++/77GdapOOTtLvvjLi88bT6OCTQUmH51MUHSQ9g1K14ZOc+DGXji7Om+DzQEqQaV08yBc+B1afgqOLbUWCY0NZeKRiThYOjCn5Rxm7vRGSpjTq5bef3kVRSnYipiZMK93LW4/i+C344/5ru13BMcE84XHFyTKDEaQaDpCcx/qnxnw/G7eBpxNKkG9EB0Cf4/RPEPQborWIlJKZh6fyZOIJyxss5BjNyI5dP0pEzpXp3wx7UMfKYqi5KR2ziXpWa8sy91vYxxflslNJnP84XHWXFmjfQMhkh6TMYKdn+SrS30qQb2wfyqEPYKey8FE+5xNm3w2cfDeQT5t8CmVbVz5cqc3dcoX5b2WeX/DND+Kioqibdu2JCRkMMcNEBsbS5s2bYiPj8/Rulu00D4CSEq6xPcy2w8bNoySJUtSq1atV9pfZtK2U261m2KYpveoiaWZMV9su0Lfan3p4tiFpReWcv7Jee0b2FXQjIRz5yiczSCRGSCVoCD1pb3yDbUXeX6Tbz2/pVW5VgxxHcK3+30Iiojhq961MTZSl/Z0sWbNGvr06YOxsXGGZczMzOjQoQObN2/O0bpPnDiRZRld4nuZ7YcOHcq+ffteaV9ZSdtOudVuimEqYW3OlK41OOMXxJ/nHzCz+UzKWJVhiseUjJ+PajgUnNrBgfxzqU8lqOhQzaW9Es7Q9nPtReKj+ezoZ9iY2TC35Vy8/ENZf/ouQ5o7Uqtc0TwO2PBdu3aNNm3aUKdOHb755huqVq0KwIYNG+jZs2dyOTc3Nw4cOADAtGnTGDNmDAC9evViw4YNr1R3REQE3bt3p27dutSqVSv5A9va2hoAPz8/XFxc+OCDD3B1daVz585ERUWli+/WrVs4ODjg6OhIvXr1sLe3p0qVKoSGhmYYd9rja9OmDfb22jvavCxtdaZtp+y0m5L/9GtUgQYV7fh673Xi4sxY0GYBTyKfMPfkXLQ+PvTiUh/AztH54lKfelD38HwIfQjDD4CphdYiP5z7gVvBt/ip40/YmdszdPtxSlibM75z9TwOVnez/r6K98MMvkm9opplbZn5umumZeLj4xk0aBCrV6+mfv36fPTRR9SqVYvY2Fh8fX1xdHT8L8ZZs5gxYwZPnz7lwoUL7Ny5E4BatWrh6emZbt+tW7cmLCws3fJvv/2Wjh07ArBv3z7Kli3L7t27AQgJCUlX/ubNm/zxxx+sWrWKfv368eeff9KvX79U8VWtWpVWrVoxfvx4WrduTbt27ViyZAm2trZa49Z2fDkpo7ZK2U4ZtZtSMBkZCeb1rk2PJcdYuP86X/epw6h6o1hyYQmtyrfijSpvpN/IriJ0mgW7J8CVP6F237wP/CUU7gT15CqcXgEN34UKjbUWOep/lI3XNzLYZTCtyrXi1xN+eD0IYfHA+moSQi22bdtG3bp1qV+/PgA1a9akZMmSBAQEYGdnl6psmzZtkFLy/fffc/jw4eRLY8bGxpiZmREWFoaNzX9DTHl4eGRZf+3atZk4cSKTJ0+mR48etG7dOl2ZypUrU69ePQAaNmyIn5+f1viuXr2afP/o+vXrODs7Zxj3kydP0m2flY4dO/L48eN0y+fNm5fqTCyjOoFU7ZRRuykFl0sZW95r4cgvx+7Qt2EFhtcazomHJ5h3ah71HOpR0bZi+o0avgcX1mvuu1frDBa2eR+4jgpvgpIS9kzS/HA6zNRaJDAqkOnHp+NczJlxDcfxNDSab/f70KpqCYN45ikzWZ3p5JbLly8nf/gDXLlyhS5dumBpaZnuAUcvLy8ePXpEiRIl0n2gxsTEYGGR+oxWlzOo6tWrc+7cOfbs2cOUKVPo3LkzM2bMSFXe3Py/TjDGxsZERUWliy8qKoro6GiKFSvG/fv3KV68OGZmZhnGre34snLw4EGdy2bUVmnbSVu7KQXb2E7V2e31iGnbr/D3Jy2Z33o+fXb24XOPz/mt62+YGKX5mDcyhu7fwaoOcGQBvGa44/UV3ntQXlvh7nHo+KXWB3KllMw6OYvw2HDmt56PmbEZX++9Tkx8IrN7uqpnnjJQvHhxbty4AcDFixdZv349devWpVixYiQkJCR/iD969IhBgwaxY8cOrKys2L9/f/I+AgMDcXBwSDfemYeHBxcvXkz3epGcAB4+fEiRIkUYPHgwEydO5Pz5DHo1pZE2Pm9vb1xcNGMwXrt2Lfn/GcWddvuclFGdadspo3ZTCjZrcxNmvl6Ta49C+e3kXUpblWZG8xl4BXix2iuDh3PLNdRcOTr1EzzxztuAX0LhTFDRIfDPNM0Pqf4QrUV23t6J+313xjQYQ9ViVTlzJ4i/LjzgwzZOODlY53HA+cc777zD2bNnady4MWvWrMHR0REnJycAOnfuzLFjx4iMjKRPnz589913uLi4MH36dL788svkfbi7u9OtW7dXqt/Ly4smTZpQr1495s2bx7Rp03Te9kV8kPrynqWlJefPn8fb2zvTuFNuDzBw4ECaN2+Oj48P5cuXZ/Xql3+SP7O2SttO2Wk3JX97zbU0bao78MOBGzwLi6GLYxe6OnZlxaUVGU/N0WGm5grSnomG22FCSqmXV8OGDaXe7J0i5cyiUvqf07r6YdhD2WxDMzlkzxAZnxAv4+IT5Gs/HJHNvzooI2Li8jbWl+Dt7a3vEGRYWFjy/xcuXCinTp2a/P78+fNy8ODBWe6jd+/e8vr167kSX2Z0jS+3tn9ZadtJX+1W2BjC35k2t56Gyapf7JYTtlyUUkoZHB0s3Ta7yV7be8mY+BjtG51dK+VMWykvbcm7QLUAzkoteaLwnUE9vwtnVkL9wVCuQbrViTKR6SemkyATmNtqLsZGxqw/dZfrj8OY1qMmRcwK7207Xfzwww+4urpSr149/Pz8mD59evK6+vXr4+bmluWDur169UrukJCXdIkvN7d/GWnbSZ/tphiGKg7WDG/lxP/O+XPu7nOKmhflyxZfciv4FssuLtO+Uf0hUKYuHJoN8bF5G7AutGWttC+gC+AD3AI+17J+EHA56XUCqJvVPvV2BrVtpJRzSkoZ7K919cZrG2WtdbXkFh/NN4pnYdGy1sx9ctCqUzIxMTEvI31phvrNTlEKEkP+OwuPjpNN5x2U3RYdlfEJms+rmcdnyjq/1pEXnlzQvtHNg5qzqFM/512gafCqZ1BCCGNgGdAVqAkMFELUTFPsDtBWSlkHmAOszInkmeOeeMOlP6DJB1C0XLrVD8Mf8sO5H2hZtiV9q2meD1i47zpRsQl8+UZN1TFCURSDZmVuwtTuLlx9GMofZ+4BMKnxJEoXKc2MEzOISYhJv1GV9uDYGo4uhJjwPI44c7pc4msC3JJS+kopY4FNQKqHNKSUJ6SUz5PengLK52yYOeTQXM3U7a3Gp1slpWT2ydkAzGg+AyEEF+8Hs+WsP8NaVaZqSfVciaIohq9HnTI0dyrON/t9eB4Ri5WpFTObz+ROyB1WXtZy7iCEpsNExDNNrz4DokuCKgfcT/HeP2lZRoYDe7MTVK647wk+u6HlGK3dynf57uL4w+OMbTCWstZlNd3M/75KCWtzRrevqoeAFUVRXp4Qgplv1CQsOo4fD2oe+WhRrgVvVHmDNV5r8AnySb9RhcZQowecWAyRGcwtpQe6JCht17W09kkUQrihSVCTM1j/oRDirBDi7LNnz3SPMrukhINfgpUDNP0o3erAqEAWeC6gnkM9BtQYAMCOiw+5cC+Yz7o4Y6NGjFAUJR+pUdqWQU0rsf70PW480TzcPqnRJGzNbZlxYgbxiVpGvW8/TTOj+LHv8zjajOmSoPyBCinelwcepi0khKgD/AL0lFIGatuRlHKllLKRlLKRg4PDq8T7au4cgbvHoM1nYJ7+Gab5Z+YTGRfJrBazMBJGRMTE8/Xea9QuV5S+DQzzaqWiKEpmxneqjrW5CXN2eSOlxM7CjilNp+Ad6M167/XpNyjpAnUHwumVEPYk7wPWQpcE5QlUE0JUFkKYAQOAnSkLCCEqAtuAd6SUN3I+zGw6sRSsSkKD9A/lHvU/yj6/fYyoMwInO80DpSuO3OZJaAxfvlETIzWVhqIo+VAxKzPGdayGx80ADl57CsBrlV7DrYIbSy8u5UH4g/QbtZkICbHg+UseR6tdlglKShkPfALsB64BW6SUV4UQI4UQI5OKzQCKA8uFEBeFEGdzLeKX9fQ63DoATT5MN1p5TEIMX5/+mspFKzOs1jAA7gdF8vNRX3rWK0vDSjkzVYKiKIo+DGpWiWolrZm725uY+ASEEHzR9AuMhBELzixIv0HxKuDcVZOgYiPzPuA0dHpQV0q5R0pZXUpZRUo5L2nZCinliqT/vy+lLCalrJf0apSbQb+UU8vAxAIaDUu3ao3XGvzD/ZnadCqmxpr7TPP3XsdYCD7vWiOvI1UURclRpsZGzHi9JncDI1l73A+A0lalGVFnBO733TnqfzT9Rs0/gagguLwpb4PVomCPJBH+FC5t1lxXtSqeatX90Pv84vULXR270rRMUwDO+gWx2+sRI9o6UaaopT4iVvKAp6cnQgj10uNLyTutqznQoUZJlh26RWC45jmoITWH4FTUia9Of0V0fJoBjiu1gDL14OQySEzM+4BTKNgJynM1JMRA849TLZZS8vWZrzExMmFCowkAJCZK5uy+Rilbcz5s46SPaJU8cubMGb2NQaleyaPPKHloSjcXIuMS+PHgTQBMjU2Z2nQqD8IfsPpKmkGMhdCcRQXegpv/6CHa/xTcBBUXpbmOWr0LlKiWapX7fXc8Hngwqt4oSlmVAuDvyw+5dD+YiZ2d1Xh7BVh8fHzyvE6KUlhULWnNoKYV2XjmHjeTup03KdOErpW7ssZrDfdC76XewLUX2JaDk0vzPtgUCm6CurwZIgM03wRSiE2IZaHnQqraVeVtl7cBiI5LYOE+H1zL2vKm6lZeoB09elTrLLuKUtB92qEaRcyM+WrPteRlkxpNwtTYlG/Pfpu6sLEpNB0Bfh7w6FIeR/qfgpugTq+E0nXAsVWqxX9c/4MH4Q+Y1HgSpkaajhGrj93hQXAUU7u7qG7lOWTChAnUrFmT0aNHp1sXFRVF27ZtsxzVvE2bNsTHa3mgMBuuXbtGjRr/dYAZNmwYJUuWTJ77Ka379+/j5uaGi4sLrq6uLFq0KNX6tMeS2f5GjBjBkSNHMt3fy2rXrh1+fn56qTuzeHLr56e8uuJJo+K4+zzj6A3NQAkORRx4v/b7uN93x/OxZ+oNGrwLZtaaz1I9KZgJ6ok3PL2qee4pxQ3ZkJgQfr78My3LtaRF2RYAPAuLYbn7LTrVLEWLKiX0FXGB4uvry/Hjx/H29mbJkiXp1q9Zs4Y+ffpgbGyc4T7MzMzo0KEDmzdvztHY0t6gHzp0KPv27cuwvImJCd999x3Xrl3j1KlTLFu2DG/v/2YgTXssme3v9OnTVK1aNdP9ZYc+604rt35+Sva828KRCvaWzNt9jYREzb3AwS6DKVWkFN+d/Y5EmaJThKUduLwB13ZCvJZBZvNAwUxQV7eBMIKaqca0ZcWlFUTERTCh4YTkZYv/vUl0fKLqVp5DfHx8aNu2LXfv3qV+/fpERESkK7NhwwZ69vzvZ+Pm5saBAwcAmDZtGmPGjAGgV69ebNiwIcdiu3r1Kq6urqmWtWnTBnv7jJ93K1OmDA0aaOYNs7GxwcXFhQcP/nvAMe2xZLS/a9euUb16dcqVK5fp/rIjr+tu37499erVo169elhYWLB169ZU63P656dkn7mJMZO71MDnSRh/nvcHwMLEgk8bfMrVwKvsubMn9Qa13oSYULh1UA/RQsHrDSAlXPkTKrcB65LJi++F3mOTzyZ6V+1NtWKaThO3n4Wz8cw93m5SkSoFbRr3vZ/DY6+c3Wfp2tB1fqZFnJ2deffdd3F0dOT9999Ptz42NhZfX18cHR2Tl82aNYsZM2bw9OlTLly4wM6dmoFKatWqhaenZ7p9tG7dmrCwsHTLv/32Wzp27Jj8fuPGjZQsWTJ52fHjxxk+fLhOh6qNn58fFy5coGnTphkeS0b27t1Lly5dMt1fbsmtug8dOgTATz/9hLu7O3369GHZsv8mxsvo56foV/faZVhV4Q7f/3OD1+uUxdLMmO5O3fnd+3cWn19Mx4odsTBJGtTAqS1Y2ms+U2t0z/NYC94Z1KNLEOSryfwp/Hj+R0yNTPm43n9dzr/Z54OFiRFjOlRLuxclG7y8vKhbt67WdQEBAdjZ2aVa1qZNG6SUfP/992zatCn5cpmxsTFmZmbpkpGHhwcXL15M90qZnGJiYrC1tU31DT4hISHTy4qZCQ8P58033+THH3/E1tY2w2PJyP79+1MlCW37e6Fjx47UqlUr3WvHjh2vFPvL1P2y9f/222/s3buXDRs2pGvbjH5+in4JIfiiaw0eh0az5vgdAIyEERMbTeRRxCM2XEtx1mtsqrkS5bMXYtNfDcltBe8M6sqfYGSiGTo+yaVnlzhw9wCj6o3CoYhmkNpzd4PYd/Ux4ztVx8HGXF/R5p4sznRy04tLaQEBAUybNg2Au3fvEh8fz5YtW4iOTv1goJeXF48ePaJEiRLY2KSedysmJgYLi9RDVOlyBmVubk6PHj1YuHAhCQkJBAYGUqpUqVc6nri4ON58800GDRpEnz59kpdbWlqmOxZtIiMjCQ4OpmzZspnu74WDB3PucsrL1v0y9W/dupUNGzawY8cOTE21j/iv7een6F9Tp+J0dCnFT4dvM6BxBYpbm9OkTBPalW/HL16/8Ga1N7GzsNMUrtUHzq2FG/s1/89DBesMKjERrv4FVTqkmvPpp4s/Ucy8GO/WfBfQPKj71Z7rONiY837ryvqKtkAKCwvD1NSUIkWKUKJECVasWMGkSZOws7Pjzz//pFixYiQkJCR/sD969IhBgwaxY8cOrKys2L9/f/K+AgMDcXBwSPfhp8sZ1AvNmzfnxIkTHDhwgE6dOr308UgpGT58OC4uLowfn3qiy7THkhF3d3fc3Nyy3F9uyK26d+3axfLly9m2bVuGCSijn59iGD7v6kxkbDxLDt1KXvZpg0+JiIvgN+/f/itYqSVYl9J8+c9jBStB+XtCyP1Ul/cuP7vM8YfHGVprKEVMiwCw/+oTzt19zvhO1dVDuTnsypUrqbo537x5k+nTp7Nq1arky0mdO3fm2LFjREZG0qdPH7777jtcXFyYPn06X375ZfK27u7udOvWLVvx9OzZkx07dhAWFpbu7Axg4MCBNG/eHB8fH8qXL8/q1f89Vd+tWzc8PDz4/fffOXToUHKHgD17/ruR/OJYMttfyntAx48fz3R/2ZGXdb/77rv4+/vTsmVL6tWrl6rdXsiJn5+Se6qWtKF/44qsP3UXvwDN5buqxarSqVInNl7fSEhMiKagkTG49oabByA6JG+D1NdQJw0bNpQ5bvckKeeUlDIqJHnRyAMjZes/WsuI2AgppZRx8QnS7Vt32eG7wzIuPiHnY9Ajb29vfYeQyrVr1+TAgQNleHh4quXnz5+XgwcPznL73r17y+vXr2crhoSEBNm0aVO5bt26bO0nI7ocS/369WVsbGyu1N+2bVt5584dvdSdVTw58fMzRIb2d5YdT0KjpMv0vfLjDeeSl/kE+cha62rJJeeX/Ffw3mkpZ9pKeWFjrsQBnJVa8kTBOYNKTADv7VCtM1hovql7PfPi2INjDHEdknz2tPWcP77PIvjsNWdMjAvO4Rua0NBQ2rdvj4mJCRMmTGDkyJFERmqG769fvz5ubm5ZPqjbq1cvnJ2dsxWHkZERrq6uyZe5cpoux3L+/Hm9XebSV9059fNTcldJGwuGt6rMrsuPuPJAc3ZUvVh1OlbsyIZrG/47iyrfGIpW1DzCk5e0Za28eOX4GZTvUU2Gv7ItedGog6Nkyz9ayvBYzTf4qNh42WTeAdl72TGZmJiYs/UbgIL0zS4nBQUF6TuEXLN27Vr5/PlzfYeRzNDiyQ0F7e8sJCpW1pu1Xw7+5VTysuuB12WtdbXksgvL/iu4f5qUs+yljAjM8Rgo8GdQD85p/q3SHoCrAVc56n+Ud2u+i5WpFQDrTvjxJDSGyV1qqCH/C5FixYrpO4RcM3ToUJ27uucFQ4tHyZqthSkfu1XF42YAJ24FAOBs70z7Cu1Z772e0NhQTcGqHSAxPuefr8xEwUlQQbfBygEsigLw8+WfsTWzZWCNgQCERMax3P0Wbs4ONHUqntmeFEVRCpXBzSpRtqgFC/b7JE+HMrLuSMLiwth4baOmkH0Vzb9Bt/MsrgKUoO4kN6B/mD+H7x9mYI2BWJtpRohYcfQ2YTHxTHpNDWmkKIqSkoWpMWM7VefS/WD2X30MgEtxF1qVa8UWny3EJcZppt8wNtcMhJBHCk6CCrwN9pqJBrfc2IKRMOKt6m8B8CQ0mrXH79Czbllqlk3/5LyiKEph92aD8lQrac3C/T7EJ2gGjR1YYyDPop5x6N4hMDIC+8oQqBLUy4mNhLCHYO9ETEIMf938i/YV2ydPRrjk0E3iEyTjO6keRYqiKNoYGwkmvuaM77MItl3QDCLcsmxLylmXY9P1TZpC9k7qDOqlPffT/Fvcif1++wmOCaa/c38A7gdFsunMffo3rkDF4kX0F6OiKIqB61yzFHXLF2XRwZvExCdgbGRMP+d+nH1yllvPb2kS1PM7mlF78kDBSFAvbtrZO7Hp+iYqF61Mk9JNAPjx4E2MjQSj26sBYRVFUTIjhOYs6kFwFJvO3Aegd9XemBmZsclnkyZBxUdrrljlgQKSoDSnnFeJxSvAi/7O/RFCcOtpGH9d8GdI80qULqoGrFQURclKq6olaFrZnqXut4iKTaCYRTG6VO7C37f/JqJoOU2hPLrMV3ASVJESbLqzC0sTS96o8gYA3x+4gaWpMR+1q6rnABVFUfIHIQSTXnPmWVgMv570A6C/c38i4yP5OyIpMakE9RICbxNi78jeO3vp4dQDGzMbrjwIYY/XY4a3qoy9lZm+I1QUJRcJIXL0Vdg1crSnnbMDPx2+TWh0HLVL1KZm8ZpsvvsP0thM02s6DxSMBBV0h/3W1sQkxNDPuR8A3/3jQ1FLU95v46Tn4BRFyU23b9/m4sWLOToUmwITOzsTEhXHLx53EELQ37k/t0Ju4128kjqD0llcFIT642mcQEnLkjgXc+b8vee4+zxjRFsnbC3UXDSKUpBdvHgxwxmclVdXq1xRutYqzZpjdwiOjKV1udYAnLWxUwlKZ8/9kMC5uEAalm6IEIIfDtyguJUZ7zZ31Hd0ShZ8fX0ZPnw4ffv21XcoSjZs376dDz74gJ49e/LPP/8kL//111/x8PBg5MiR9O3bl59++inH61ZnPLlnbMfqRMTGs8rDF4ciDlSyrcRZU6EZuScPuprn/wQV5Ms9ExOexYXTqFQjztwJwuNmACPbVsHKXE1GqA8///wzZcqUSZ4U75133smwrJOTk9bJ7qKiomjbtm2WU3K0adOG+Pj4V4ozOjqaJk2aULduXVxdXZk5c2aG9Q8bNoySJUummowxpREjRnDkyBHc3NxwcXHB1dWVRYsWvVJcAO3atcPPzy/5fWb153TdusSVtu179erFqlWrWLduHZs3b04uf+7cOVq1asWKFSvYsmULZ8+eTRf78ePHU7X3/fv3dT6W58+fU7z4f2NrZvVzymrf2f2dKmicS9vQo05Z1h73IygiloalGnI+7jmJ8VEQ/jjX68//CSrwNucszAFoVKoRPxy4gYONOYObVdJzYIXX5cuXmTt3bvJU7L///jteXl706NEj1evp06cZ7mPNmjX06dMHY2PjDMuYmZnRoUOHVB+IL8Pc3JxDhw5x6dIlLl68yL59+zh16pTW+ocOHcq+ffsy3Nfp06epWrUq3333HdeuXePUqVMsW7YMb2/vV4otrczqz+26tcmo7efOncvHH38MQFxcHCYmJggh2LlzJ61ataJDhw7pYm/WrFmq9jYxMdH5WDw8PGjZsmXy+6x+TlntO7u/UwXRpx2qER2XwM9Hb9OoVCNCE2O4aWqaJx0l8n+CCvLlrLUt9hb2PAqw4aRvIKPaVcHSLOMPNiV3eXl5Ua9evVTLateuza5du1K9SpYsmeE+NmzYQM+ePZPfDxgwgP79+9O0aVMqVarE7t27Ac039w0bNrxSnEIIrK01gwnHxcURFxeX3IMrbf1t2rTB3t5e636uXbtG9erVKVeuHA0aNADAxsYGFxcXHjx48EqxpZVR/blVd0btnVLKtpdSMnnyZLp27Zocx9GjR2ndWnPf4o033uDEiROpflYvYjc2Nk7V3mXKlMnwWMLDw5kxY0byPuLi4jAz+6+XbmY/p6z2re24FKha0pqe9crx24m7OFrXBuCsZd4MGpv/r4EF+XLOwpwGJRvw48GblLI1Z2CTivqOSu8WnFnA9aDrObrPGvY1mNxkcpblrl69ynvvvYeRkRElSpTg4MGDGZYNDAxk6tSpXLhwga+//popU6YQGxuLr68vjo6OyeUuXbpEr1692Lx5M8eOHWP8+PF0796dWrVq4enpmW6/rVu3JiwsLN3yb7/9lo4dOya/T0hIoGHDhty6dYuPP/6Ypk2baq0/M3v37qVLly6plvn5+XHhwgWaNm2q0z5eVW7VnVF7p5Sy7ZcsWcLBgwcJCQnh1q1bjBw5kgMHDjBt2jQOHz7Mtm3biImJoVu3buliz6y90x6LpaUllpaW3Lhxg8qVK6dKTi8ro3bK6HeqMBvToRo7Lz1ku2ckZa3KcC4ymkEqQWXtYbAvD+0ELY1r8Jffc+b0dMXCVJ096cv9+/cpXbo0ly9f1ql88eLFWbFiRaplAQEBqSa9i4qKIiAgIPkeUc2aNXn+/DkAxsbGmJmZERYWho2NTfI2Hh4eOtVvbGzMxYsXCQ4Opnfv3ly5cgV7e/uXmnRv//79rF27Nvl9eHg4b775Jj/++CO2tqlHz+/YsSOPH6e/dj9v3rxUZ2z6rDuz9k4pZduPGTOGMWPGpFofHh6OtbU17dq1o127dhnGnvbnndmxGBsb06tXL3bs2EGjRo1SXd57GZm1U0a/U4VZ5RJW9K5fjvWn7tK9fQPOhT1CBt4it58Yy98JKi6ac3FBQHHOeNtRtqgF/RpX0HdUBkGXM53ccPnyZVxdXbO1D0tLS6Kjo5PfX7lyhWrVqmFhoRmu6vz586m6FcfExCSve0HXM6gX7OzsaNeuHfv27WP48OGp6s9MZGQkwcHBlC1bFtBccnrzzTcZNGgQffr0SVc+s7PJl5VbdWfV3ilpa/sXli5dqlPsz58/T9femR2Li4sL33zzDU5OTplezstIVu0EmR9XYTWmfTX+uvCAoKByBBnBneDb5PZTpvk7QQXf5ayFOVbCnCt+RZjTqyrmJursSZ+8vLzSJaiAgACmTZsGwN27d4mPj+fAgQMZ7qNYsWIkJCQQHR2NhYUFly5d4t69e0RHR5OQkMDMmTNZuHAhoLlE6ODggKlp6ufddDmDevbsGaamptjZ2REVFcXBgweZPHlyuvoz4+7ujpubG6C5DzN8+HBcXFwYP358lvVnV27VnVl7p5RR279s7GnbW5djcXBwICYm5qXr1WXf2Tmugqxi8SL0qV+OnVcDMXOEs9FPcJIScnHkjfzdSSKpB59VbDnKFC1Cv0bl9R1Roefl5UXNmjVTLStRogQrVqxg0qRJ2NnZ8eeff2a5n86dO3Ps2DFA84E5aNAg2rVrR+PGjfnoo4+SL+24u7unuq/xMh49eoSbmxt16tShcePGdOrUiR49eqSrH2DgwIE0b94cHx8fypcvn9w1PuU9oOPHj/P7779z6NCh5C72e/bseaXY0tJWf27VnVl7p5Sdtk977yxle+tyLG+99RbNmjVLt9+Mfk7dunXj4cOHOu07O8dV0H3SvirxMfbYSAvN81BhudzVPCeHB3mZV8OGDWV2PT0yX9ZaV0vWXjBW/nbiTrb3l995e3vrO4QM3bhxQw4cOFCGhYXpVP78+fNy8ODBUkopW7duLa9fv661XO/evTNclx0p689M/fr1ZWxsbI7X37ZtW3nnzh291J1Ze6eMKzttnzZ2Xds7L2R1XIb8d5YXJmy5KNss6S7b/1JDJvp65Mg+gbNSS57Q6QxKCNFFCOEjhLglhPhcy3ohhFictP6yEKJBjmdSLTyfXgTAXNRR954M2PXr15k5cyarVq1K7tadlfr16+Pm5kZCQgK3b9+mWrX083nFxsbSq1cvnJ1zfqbklPVn5vz583q7FJRbdWfU3illt+3Txq5re+e23PydKig+catKTERVnpqY4P/4XK7WlWWCEkIYA8uArkBNYKAQomaaYl2BakmvD4GcH89Ei5Mhd7BIhI9btlH3ngxUaGgo7du3x8TEhAkTJjBy5EgiIyN12nbYsGEYGxvz4MEDjIzS/6qamZkxZMiQnA45Xf36MHTo0JfqSZiTMmpv+C+u3Gh7fbb3C7n9O1UQOJawon7F9gAcf3g2i9LZo0sniSbALSmlL4AQYhPQE0j5aHdP4LekU7VTQgg7IUQZKeWjHI84iZSSC4lhVIorwoDGjrlVjZJNtra2PHyYN7NvFiRDhw7VdwhaGWpcSt76rGN7Lv8t8Qi6xYBcrEeXS3zlgPsp3vsnLXvZMgghPhRCnBVCnH327NnLxppKcGQk1tKUelbO6rknRVGUPFTZwZo68cVAmudqPbqcQWnrQ5h2+GBdyiClXAmsBGjUqFG2hiAuZmXFpg8vkpioRjJWFEXJa8veP4qRUe4+qqvLGZQ/kLIHQnkg7TUbXcrkitxuIEVRFCW9vPjs1SVBeQLVhBCVhRBmwABgZ5oyO4EhSb35mgEhuXn/ScmYVHPjKEquUX9feSvLS3xSynghxCfAfsAYWCOlvCqEGJm0fgWwB+gG3AIigfdyL2QlIxYWFgQGBlK8ePHkUbkVRckZUkoCAwPVEEh5SOjrG0GjRo1k2snLlOyJi4vD399f53HkFEV5ORYWFpQvX14Ng5TDhBDnpJSN0i7P32PxKamYmppSuXJlfYehKIqSI/L3WHyKoihKgaUSlKIoimKQVIJSFEVRDJLeOkkIIZ4Bd3NgVyWAgBzYT0Gj2iVjqm20U+2SMdU22uVUu1SSUjqkXai3BJVThBBntfX+KOxUu2RMtY12ql0yptpGu9xuF3WJT1EURTFIKkEpiqIoBqkgJKiV+g7AQKl2yZhqG+1Uu2RMtY12udou+f4elKIoilIwFYQzKEVRFKUAUglKURRFMUj5IkEJIboIIXyEELeEEJ9rWS+EEIuT1l8WQjTQR5z6oEPbDEpqk8tCiBNCiLr6iDOvZdUuKco1FkIkCCH65mV8+qRL2wgh2gkhLgohrgohjuR1jPqgw99SUSHE30KIS0ntUihmbRBCrBFCPBVCXMlgfe59/kopDfqFZoqP24ATYAZcAmqmKdMN2ItmZt9mwGl9x21AbdMCKJb0/66FoW10aZcU5Q6hmS6mr77jNpS2AewAb6Bi0vuS+o7bQNrlC2BB0v8dgCDATN+x50HbtAEaAFcyWJ9rn7/54QyqCXBLSukrpYwFNgE905TpCfwmNU4BdkKIMnkdqB5k2TZSyhNSyudJb0+hme24oNPldwZgNPAn8DQvg9MzXdrmbWCblPIegJSyMLSPLu0iARuhmWzNGk2Cis/bMPOelPIommPNSK59/uaHBFUOuJ/ivX/SspctUxC97HEPR/NNp6DLsl2EEOWA3sCKPIzLEOjyO1MdKCaEOCyEOCeEGJJn0emPLu2yFHABHgJewKdSysS8Cc+g5drnb36YD0rb1LBp+8brUqYg0vm4hRBuaBJUq1yNyDDo0i4/ApOllAmFbPZhXdrGBGgIdAAsgZNCiFNSyhu5HZwe6dIurwEXgfZAFeCAEMJDShmay7EZulz7/M0PCcofqJDifXk032BetkxBpNNxCyHqAL8AXaWUgXkUmz7p0i6NgE1JyakE0E0IES+l3J4nEeqPrn9PAVLKCCBCCHEUqAsU5ASlS7u8B8yXmhsvt4QQd4AawJm8CdFg5drnb364xOcJVBNCVBZCmAEDgJ1pyuwEhiT1JmkGhEgpH+V1oHqQZdsIISoC24B3Cvg34JSybBcpZWUppaOU0hH4HzCqECQn0O3vaQfQWghhIoQoAjQFruVxnHlNl3a5h+asEiFEKcAZ8M3TKA1Trn3+GvwZlJQyXgjxCbAfTU+bNVLKq0KIkUnrV6DphdUNuAVEovmmU+Dp2DYzgOLA8qSzhXhZwEdl1rFdCiVd2kZKeU0IsQ+4DCQCv0gptXYxLih0/J2ZA6wTQnihuaw1WUpZ4KfgEEL8AbQDSggh/IGZgCnk/uevGupIURRFMUj54RKfoiiKUgipBKUoiqIYJJWgFEVRFIOkEpSiKIpikFSCUhRFUQySSlCKoiiKQVIJSlEURTFIKkEpigFJehpf/V0qCipBKUqOEUK4CyE6Jf1/rhBisZYy44UQV5JeY5OWOQohrgkhlgPnST2uGUKITUKIzUKI00KIu0KI7nlwOIqidwY/1JGi5CMzgdlCiJJAfeCNlCuFEA3RDAPTFM1QOaeTZqt9jmZct/eklKO07LcusF1K2V8I0Qr4Htide4ehKIZBnUEpSg5JmthNAOOBAVLKhDRFWgF/SSkjpJThaAbxbZ207m7SZG+pCCEs0Yy2PitpkTdQLDfiVxRDoxKUouQQIURtoAwQI6UM01Ykk80jMlheC7gppYxOet8AzXTkilLgqQSlKDkgaYrrDWimv44QQrympdhRoJcQoogQwgrNjL4eWey6LlBRCGGRtM0s4IccDF1RDJZKUIqSTUlzJm0DJkgpr6GZluHLtOWklOeBdWgmuDuNZhqLC1nsvi6axHcYzZxFP0kpj+dU7IpiyNR0G4piwJJms/1ASumj71gUJa+pBKUoBkwI8QCoIKVM1HcsipLXVIJSFEVRDJK6B6UoiqIYJJWgFEVRFIOkEpSiKIpikFSCUhRFUQySSlCKoiiKQVIJSlEURTFIKkEpiqIoBun/dEtTsMPEvoIAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Test for the function I found.\n", + "\n", + "x = np.linspace(0, 1, 100)\n", + "g = lambda x: np.sin(np.pi*x*(1 - x))\n", + "f_Z = lambda z: np.sqrt(1.2*(1/2-np.abs(1/2 - z)))\n", + "F_Z_inv = lambda p: 3/(2*np.sqrt(1.2))*(1/2 - np.abs(p-.5))**(2/3)\n", + "\n", + "plt.figure()\n", + "plt.title(\"My choice of distribution for $Z$\")\n", + "plt.plot(x, g(x), label=\"$g(x) = \\sin(\\pi x(1-x))$\")\n", + "plt.plot(x, f_Z(x), label=\"$f_Z(x) = \\sqrt{1.2(1/2-|1/2 - z|)}$\")\n", + "plt.plot(x, F_Z_inv(x), label=\"$F_Z^{-1}(p) = 3(1/2 - |1/2 - p|)^{2/3}/(2\\sqrt{1.2})$\")\n", + "plt.xlabel('$x$ or $p$')\n", + "plt.legend()\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 134, "id": "5793420a", "metadata": { "deletable": false, @@ -425,8 +463,8 @@ "output_type": "stream", "text": [ "\n", - "After trial and error, I found that 43402 trials were sufficient to arrive at an 1-sigma error of around 0.001,\n", - "with mean 0.175938 and standard deviation sigma 0.000418. Wolfram Alpha tells me it should be approximately I =\n", + "After trial and error, I found that 1500 trials were sufficient to arrive at an 1-sigma error of around 0.001,\n", + "with mean 0.484070 and standard deviation sigma 0.001016. Wolfram Alpha tells me it should be approximately I =\n", "0.487595, so I am happy.\n", " \n" ] @@ -436,14 +474,20 @@ "def sample_nice_Z():\n", " '''Sample from the nice distribution Z'''\n", " p = rng.random()\n", - " return np.arctan(p*np.pi)/np.pi\n", + " #return 3/(2*np.sqrt(1.2))*p**(2/3)\n", + " #return 3/(2*np.sqrt(1.2))*(p%.5)**(2/3)\n", + " return 3/(2*np.sqrt(1.2))*np.abs(p - .5)**(2/3)\n", + " #return np.arctan(np.pi*p)/np.pi\n", + " #return 1/3*np.arccos(1-4*p)\n", " \n", "def sample_X_prime():\n", " '''Sample from X'.'''\n", " z = sample_nice_Z()\n", - " return np.sin(np.pi*z*(1 - z))*np.cos(np.pi*z)**2\n", + " return np.sin(np.pi*z*(1 - z))/np.sqrt(1.2*(1/2-np.abs(1/2 - z)))/2\n", + " #return np.sin(np.pi*z*(1 - z))*np.cos(np.pi*z)**2\n", + " #return np.sin(np.pi*z*(1 - z))/(3/2*np.sin(3/2*z)*np.cos(3/2*z))\n", " \n", - "n = 43402\n", + "n = 1500\n", "sample_mean, sample_stdev = estimate_expectation(sample_X_prime, n)\n", "print(\"\"\"\n", "After trial and error, I found that {} trials were sufficient to arrive at an 1-sigma error of around 0.001,\n", diff --git a/exercise_sheet_01.html b/exercise_sheet_01.html new file mode 100644 index 0000000..f2aa53d --- /dev/null +++ b/exercise_sheet_01.html @@ -0,0 +1,1361 @@ + + + + + +exercise_sheet_01 + + + + + + + + + + + + + + + + + + +
+
+
+

exercise_sheet_01 (Score: 97.0 / 100.0)

+
+
    + + + + + + + + + + + + + + + + + + + + + + + +
  1. Test cell (Score: 12.0 / 12.0)
    2. + + + + + + + + + + +
  2. Test cell (Score: 12.0 / 12.0)
    3. + + + + + + + +
  3. Coding free-response (Score: 12.0 / 12.0)
    4. + + +
  4. Comment
    5. + + + + + + + + + +
  5. Test cell (Score: 12.0 / 12.0)
    6. + + + + + + + +
  6. Coding free-response (Score: 11.0 / 12.0)
    7. + + +
  7. Comment
    8. + + + + + + +
  8. Coding free-response (Score: 14.0 / 14.0)
    9. + + +
  9. Comment
    10. + + + + + + +
  10. Coding free-response (Score: 14.0 / 14.0)
    11. + + +
  11. Comment
    12. + + + + + + +
  12. Coding free-response (Score: 10.0 / 12.0)
    13. + + +
  13. Comment
    14. + + + +
+
+
+
+
+
+ + +
+
+ +
+
+
+
+

Exercise sheet

Some general remarks about the exercises:

+
    +
  • For your convenience functions from the lecture are included below. +Feel free to reuse them without copying to the exercise solution box.
    • +
  • For each part of the exercise a solution box has been added, but you + may insert additional boxes. Do not hesitate to add Markdown boxes for +textual or LaTeX answers (via Cell > Cell Type > Markdown). But make sure to replace any part that says YOUR CODE HERE or YOUR ANSWER HERE and remove the raise NotImplementedError().
    • +
  • Please make your code readable by humans (and not just by the Python + interpreter): choose informative function and variable names and use +consistent formatting. Feel free to check the PEP 8 Style Guide for Python for the widely adopted coding conventions or this guide for explanation.
    • +
  • Make sure that the full notebook runs without errors before submitting your work. This you can do by selecting Kernel > Restart & Run All in the jupyter menu.
    • +
  • Each sheet has 100 points worth of exercises. Note that only the +grades of sheets number 2, 4, 6, 8 count towards the course examination. + Submitting sheets 1, 3, 5, 7 & 9 is voluntary and their grades are +just for feedback.
    • +
+

Please fill in your name here:

+ +
+
+ +
+
+
In [1]:
+
+ 
NAME = "Kees van Kempen"
+
+ +
+
+
+ +
+
+
+
+
+
+ +
+
+ +
+
+
+
+

Exercise sheet 1

+

Code from the lecture:

+ +
+
+ +
+
+
In [2]:
+
+ 
import numpy as np
+import matplotlib.pylab as plt
+
+rng = np.random.default_rng()
+%matplotlib inline
+
+def random_in_square():
+    """Returns a random position in the square [-1,1)x[-1,1)."""
+    return rng.uniform(-1,1,2)  
+
+def is_in_circle(x):
+    return np.dot(x,x) < 1
+
+def simulate_number_of_hits(N):
+    """Simulates number of hits in case of N trials in the pebble game."""
+    number_hits = 0
+    for i in range(N):
+        position = random_in_square()
+        if is_in_circle(position):
+            number_hits += 1
+    return number_hits
+
+def random_in_disk():
+    """Returns a uniform point in the unit disk via rejection."""
+    position = random_in_square()
+    while not is_in_circle(position):
+        position = random_in_square()
+    return position
+
+def is_in_square(x):
+    """Returns True if x is in the square (-1,1)^2."""
+    return np.abs(x[0]) < 1 and np.abs(x[1]) < 1
+    
+def sample_next_position_naively(position,delta):
+    """Keep trying a throw until it ends up in the square."""
+    while True:
+        next_position = position + delta*random_in_disk()
+        if is_in_square(next_position):
+            return next_position
+    
+def naive_markov_pebble(start,delta,N):
+    """Simulates the number of hits in the naive Markov-chain version 
+    of the pebble game."""
+    number_hits = 0
+    position = start
+    for i in range(N):
+        position = sample_next_position_naively(position,delta)
+        if is_in_circle(position):
+            number_hits += 1
+    return number_hits
+
+def naive_markov_pebble_generator(start,delta,N):
+    """Same as naive_markov_pebble but only yields the positions."""
+    position = start
+    for i in range(N):
+        position = sample_next_position_naively(position,delta)
+        yield position
+
+def sample_next_position(position,delta):
+    """Attempt a throw and reject when outside the square."""
+    next_position = position + delta*random_in_disk()
+    if is_in_square(next_position):
+        return next_position  # accept!
+    else:
+        return position  # reject!
+
+def markov_pebble(start,delta,N):
+    """Simulates the number of hits in the proper Markov-chain version of the pebble game."""
+    number_hits = 0
+    position = start
+    for i in range(N):
+        position = sample_next_position(position,delta)
+        if is_in_circle(position):
+            number_hits += 1
+    return number_hits
+
+def markov_pebble_generator(start,delta,N):
+    """Same as markov_pebble but only yields the positions."""
+    position = start
+    for i in range(N):
+        position = sample_next_position(position,delta)
+        yield position
+
+ +
+
+
+ +
+
+
+
+
+

Empirical convergence rate in the pebble game

(a) Write a function pi_stddev that estimates the standard deviation σ of the π-estimate using n trials by running the direct-sampling pebble game m times. Store this data for n=24,25,…,214 and m=200 in an array. Hint: you may use the NumPy function np.std. (12 pts)

+ +
+
+ +
+
+
In [3]:
+
Student's answer(Top)
+
+
+ 
def pi_stddev(n,m):
+    """Estimate the standard deviation in the pi estimate in the case of n trials,
+    based on m runs of the direct-sampling pebble game."""
+    
+    return np.std([simulate_number_of_hits(n)/n*4 for i in range(m)])
+    
+stddev_data = np.array([[2**k,pi_stddev(2**k,200)] for k in range(4,12)])
+stddev_data
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
Out[3]:
+ + + + +
+
array([[1.60000000e+01, 3.54522126e-01],
+       [3.20000000e+01, 2.87754643e-01],
+       [6.40000000e+01, 2.21671341e-01],
+       [1.28000000e+02, 1.50942352e-01],
+       [2.56000000e+02, 1.00401902e-01],
+       [5.12000000e+02, 6.94667450e-02],
+       [1.02400000e+03, 4.96484068e-02],
+       [2.04800000e+03, 3.59752212e-02]])
+
+ +
+ +
+
+ +
+
+
+
In [4]:
+
Grade cell: cell-90b5022604409d60 + + Score: 12.0 / 12.0 (Top) +
+
+
+ 
# If done correctly, your code should pass the following tests
+from nose.tools import assert_almost_equal
+assert_almost_equal(pi_stddev(2**3,1000),0.582,delta=0.03)
+assert_almost_equal(pi_stddev(2**4,1000),0.411,delta=0.03)
+### BEGIN HIDDEN TESTS
+assert_almost_equal(pi_stddev(2**5,1000),0.291,delta=0.03)
+assert_almost_equal(pi_stddev(2**6,1000),0.207,delta=0.03)
+### END HIDDEN TESTS
+
+ +
+
+ +
+
+ +
+
+
+
+
+

(b) Write a function fit_power_law that takes an array of (n,σ) pairs and determines best-fit parameters a,p for the curve σ=anp. This is best done by fitting a straight line to the data on a log-log scale (log⁡σ=log⁡a+plog⁡n). Hint: use curve_fit from SciPy. (12 pts)

+ +
+
+ +
+
+
In [5]:
+
Student's answer(Top)
+
+
+ 
import scipy.optimize
+
+def fit_power_law(stddev_data):
+    """Compute the best fit parameters a and p."""
+    
+    # Define the to be fitted function. We use the log, due to the calculation of the uncertainties.
+    log_f = lambda n, a, p: np.log(a) + p*np.log(n)
+    
+    n, 𝜎 = stddev_data.T
+    a_fit, p_fit = scipy.optimize.curve_fit(log_f, n, np.log(𝜎))[0]
+    
+    return a_fit, p_fit
+
+ +
+
+ +
+
+ +
+
+
+
In [6]:
+
Grade cell: cell-311260afe5628598 + + Score: 12.0 / 12.0 (Top) +
+
+
+ 
from nose.tools import assert_almost_equal
+assert_almost_equal(fit_power_law(np.array([[16.0,1.4],[32.0,1.1],[64.0,0.9]]))[0],3.36,delta=0.05)
+assert_almost_equal(fit_power_law(np.array([[16.0,1.4],[32.0,1.1],[64.0,0.9]]))[1],-0.31,delta=0.03)
+### BEGIN HIDDEN TESTS
+assert_almost_equal(fit_power_law(stddev_data)[0],1.82,delta=0.3)
+assert_almost_equal(fit_power_law(stddev_data)[1],-0.5,delta=0.08)
+### END HIDDEN TESTS
+
+ +
+
+ +
+
+ +
+
+
+
+
+

(c) Plot the data against the fitted curve σ=anp in a log-log plot. Don't forget to properly label your axes. Hint: use loglog from matplotlib. (12 pts)

+ +
+
+ +
+
+
In [7]:
+
Student's answer + Score: 12.0 / 12.0 (Top) +
+
+
+ 
a_fit, p_fit = fit_power_law(stddev_data)
+
+plt.figure()
+plt.title("Convergence of standard deviation of direct-sampled $\pi$ guesses over $n$ throws among $m = 200$ trials")
+plt.loglog(*stddev_data.T, label="data")
+n_range = np.logspace(1, 4, 10)
+plt.loglog(n_range, a_fit*n_range**p_fit,
+           label="fitted $\sigma = an^p$\n($a = {:.3f}, p = {:.3f})$".format(a_fit, p_fit))
+plt.xlabel("$n$")
+plt.ylabel("$\sigma$")
+plt.legend()
+plt.show()
+
+ +
+
+
+
+ Comments:

Good labeling! Use dots for the data.

+ +
+
+
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

Volume of a unit ball in other dimensions

In this exercise you will adapt the direct-sampling Monte Carlo method of the pebble game to higher dimensions to +estimate the d-dimensional volume of the d-dimensional ball of radius 1.

+

(a) Adapt the direct-sampling Monte Carlo method to the pebble game in a d-dimensional hypercube (−1,1)d with an inscribed d-dimensional unit ball. (12 pts)

+ +
+
+ +
+
+
In [8]:
+
Student's answer(Top)
+
+
+ 
def number_of_hits_in_d_dimensions(N,d):
+    """Simulates number of hits in case of N trials in the d-dimensional direct-sampling pebble game."""
+
+    def random_in_d_dimensional_box():
+        """Returns a random position in the box [-1,1)^d."""
+        return rng.uniform(-1,1,d)
+        
+    number_hits = 0
+    for i in range(N):
+        position = random_in_d_dimensional_box()
+        if is_in_circle(position):
+            number_hits += 1
+    return number_hits
+
+ +
+
+ +
+
+ +
+
+
+
In [9]:
+
Grade cell: cell-e8e97677ad30041e + + Score: 12.0 / 12.0 (Top) +
+
+
+ 
from nose.tools import assert_almost_equal
+assert_almost_equal(number_of_hits_in_d_dimensions(100,1),100,delta=1)
+assert_almost_equal(number_of_hits_in_d_dimensions(2000,3),1045,delta=50)
+### BEGIN HIDDEN TESTS
+assert_almost_equal(number_of_hits_in_d_dimensions(2000,2),1572,delta=50)
+assert_almost_equal(number_of_hits_in_d_dimensions(2000,4),619,delta=50)
+### END HIDDEN TESTS
+
+ +
+
+ +
+
+ +
+
+
+
+
+

(b) Compare your estimates for the volume Vd of the d-dimensional unit ball for N=10000 trials and d=1,2,…,7 to the exact formula Vd=πd/2Γ(d2+1) in a plot. Hint: the Gamma function is available in scipy.special.gamma. (12 pts)

+ +
+
+ +
+
+
In [10]:
+
Student's answer + Score: 11.0 / 12.0 (Top) +
+
+
+ 
import scipy.special
+
+def volume_d_ball(d):
+    """Returns the volume of a d-dimensional unit ball."""
+    return np.pi**(d/2)/scipy.special.gamma(d/2 + 1)
+
+N = 10000
+d_range = np.array(range(1, 8))
+V_d_mc = [number_of_hits_in_d_dimensions(N, d)/N*2**d for d in d_range]
+
+## Plotting
+plt.figure()
+plt.plot(d_range, V_d_mc, label="direct-sampling")
+plt.plot(d_range, volume_d_ball(d_range), label="exact")
+plt.xlabel("dimension $d$")
+plt.ylabel("volume $V_d$")
+plt.legend()
+plt.show()
+
+ +
+
+
+
+ Comments:

Use dots

+ +
+
+
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

Efficiency of the Metropolis algorithm and the 1/2-rule

In + the lecture we mentioned the "1/2 rule" for Metropolis algorithms: if +moves are rarely rejected then typically you do not explore the domain +efficiently, while if most moves are rejected you are wasting efforts +proposing moves. A good rule of thumb then is to aim for a rejection +rate of 1/2. + In this exercise you are asked to test whether this rule of thumb makes + any sense for the Markov-chain pebble game by varying the throwing +range δ.

+

(a) Estimate the mean square deviation E[(4hitstrials−π)2] from π for different values of δ ranging between 0 and 3, but fixed number of trials (say, 2000). For a decent estimate of the mean you will need at least 100 repetitions. (14 pts)

+ +
+
+ +
+
+
In [11]:
+
Student's answer + Score: 14.0 / 14.0 (Top) +
+
+
+ 
N = 2000 # number of trials
+m = 200  # number of repetitions
+delta_range = np.linspace(0, 3, 10) # throwing ranges
+
+p0 = np.array([0., 0.]) # starting position in origin
+hits = np.zeros((len(delta_range), m), dtype=int)
+
+for index_delta, delta in enumerate(delta_range):
+    for i in range(m):
+        hits[index_delta][i] = markov_pebble(p0, delta, N)
+
+## Print the results
+print("{:4}   {:6}   {:8}   {:<10}".format('N','delta','4*hits/trials','E[(4*hits/trials - pi)^2]'))
+for index_delta, delta in enumerate(delta_range):
+    pi_estimates = hits[index_delta]/N*4
+    print("{:4}   {:1.4f}   {:1.4f}          {:1.3f}"
+          .format(N, delta, np.mean(pi_estimates), np.mean((pi_estimates - np.pi)**2)))
+
+ +
+
+
+
+ Comments:

Nice formatting

+ +
+
+
+
+ +
+
+ + +
+ +
+ + +
+
---------------------------------------------------------------------------
+KeyboardInterrupt                         Traceback (most recent call last)
+Input In [11], in <cell line: 8>()
+      8 for index_delta, delta in enumerate(delta_range):
+      9     for i in range(m):
+---> 10         hits[index_delta][i] = markov_pebble(p0, delta, N)
+     12 ## Print the results
+     13 print("{:4}   {:6}   {:8}   {:<10}".format('N','delta','4*hits/trials','E[(4*hits/trials - pi)^2]'))
+
+Input In [2], in markov_pebble(start, delta, N)
+     70 position = start
+     71 for i in range(N):
+---> 72     position = sample_next_position(position,delta)
+     73     if is_in_circle(position):
+     74         number_hits += 1
+
+Input In [2], in sample_next_position(position, delta)
+     59 def sample_next_position(position,delta):
+     60     """Attempt a throw and reject when outside the square."""
+---> 61     next_position = position + delta*random_in_disk()
+     62     if is_in_square(next_position):
+     63         return next_position  # accept!
+
+Input In [2], in random_in_disk()
+     23 def random_in_disk():
+     24     """Returns a uniform point in the unit disk via rejection."""
+---> 25     position = random_in_square()
+     26     while not is_in_circle(position):
+     27         position = random_in_square()
+
+Input In [2], in random_in_square()
+      7 def random_in_square():
+      8     """Returns a random position in the square [-1,1)x[-1,1)."""
+----> 9     return rng.uniform(-1,1,2)
+
+KeyboardInterrupt:
+
+
+ +
+
+ +
+
+
+
+
+

(b) Measure the rejection rate for the same simulations as in (a). (14 pts)

+ +
+
+ +
+
+
In [12]:
+
Student's answer + Score: 14.0 / 14.0 (Top) +
+
+
+ 
# To calculate the rejection rate, we need to know the positions. I will copy
+# and adapt the above code and the code in the markov_pebble function as is
+# provided at the beginning of this document.
+
+N = 2000 # number of trials
+m = 200  # number of repetitions
+delta_range = np.linspace(0, 3, 10) # throwing ranges
+
+p0 = np.array([0., 0.]) # starting position in origin
+hits = np.zeros((len(delta_range), m), dtype=int)
+rejections = np.zeros((len(delta_range), m), dtype=int)
+
+for index_delta, delta in enumerate(delta_range):
+    for i in range(m):
+        prev_position = None
+        for position in markov_pebble_generator(p0,delta,N):
+            if np.all(position == prev_position):
+                rejections[index_delta][i] += 1
+            if is_in_circle(position):
+                hits[index_delta][i] += 1
+            prev_position = position
+
+## Print the results
+print("{:4}   {:6}   {:8}   {:<10}    {:<10}"
+      .format('N','delta','4*hits/trials','E[(4*hits/trials - pi)^2]', 'rejection rate'))
+for index_delta, delta in enumerate(delta_range):
+    pi_estimates = hits[index_delta]/N*4
+    print("{:4}   {:1.4f}   {:1.4f}          {:1.3f}                        {:1.3f}"
+          .format(N, delta, np.mean(pi_estimates), np.mean((pi_estimates - np.pi)**2), np.mean(rejections[index_delta]/N)))
+
+ +
+
+
+
+ Comments:

Nice formatting

+ +
+
+
+
+ +
+
+ + +
+ +
+ + +
+
---------------------------------------------------------------------------
+KeyboardInterrupt                         Traceback (most recent call last)
+Input In [12], in <cell line: 13>()
+     15 prev_position = None
+     16 for position in markov_pebble_generator(p0,delta,N):
+---> 17     if np.all(position == prev_position):
+     18         rejections[index_delta][i] += 1
+     19     if is_in_circle(position):
+
+File <__array_function__ internals>:5, in all(*args, **kwargs)
+
+File /opt/jupyter-conda/lib/python3.9/site-packages/numpy/core/fromnumeric.py:2450, in all(a, axis, out, keepdims, where)
+   2367 @array_function_dispatch(_all_dispatcher)
+   2368 def all(a, axis=None, out=None, keepdims=np._NoValue, *, where=np._NoValue):
+   2369     """
+   2370     Test whether all array elements along a given axis evaluate to True.
+   2371 
+   (...)
+   2448 
+   2449     """
+-> 2450     return _wrapreduction(a, np.logical_and, 'all', axis, None, out,
+   2451                           keepdims=keepdims, where=where)
+
+File /opt/jupyter-conda/lib/python3.9/site-packages/numpy/core/fromnumeric.py:86, in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs)
+     83         else:
+     84             return reduction(axis=axis, out=out, **passkwargs)
+---> 86 return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
+
+KeyboardInterrupt:
+
+
+ +
+
+ +
+
+
+
+
+

(c) Plot both the mean square deviation and the rejection rate as function of δ. How well does the 1/2 rule apply in this situation? (12 pts)

+ +
+
+ +
+
+
In [13]:
+
Student's answer + Score: 10.0 / 12.0 (Top) +
+
+
+ 
plt.figure()
+plt.title("Mean square deviation and rejection rate of $\\pi$ guesses by the Metropolis algorithm for {} trials of each {} repetitions".format(N, m))
+plt.plot(delta_range, np.mean((hits/N*4 - np.pi)**2, axis=1), label="$\mathbb{E}[(\\frac{4hits}{trials} - \\pi)^2]$")
+plt.plot(delta_range, np.mean(rejections/N, axis=1), label="rejection rate")
+plt.xlabel("throwing range $\\delta$")
+plt.legend()
+plt.rc('legend', fontsize=13)
+plt.show()
+
+print("The rejection rate is around 1/2 for delta = 1.3333, which also yields near lowest mean square deviation to pi.")
+print("The rule of thumb works quite nicely in this case, although it does not yields the best result.")
+
+ +
+
+
+
+ Comments:

Good labels. Due to the outlier at 0, we cannot see the minimum MSD in the plot.

+ +
+
+
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+ +
+ + +
+
The rejection rate is around 1/2 for delta = 1.3333, which also yields near lowest mean square deviation to pi.
+The rule of thumb works quite nicely in this case, although it does not yields the best result.
+
+
+
+ +
+
+ +
+
+
+ + +
+
+
+
+
+ + + \ No newline at end of file diff --git a/exercise_sheet_02.html b/exercise_sheet_02.html new file mode 100644 index 0000000..58fe1d3 --- /dev/null +++ b/exercise_sheet_02.html @@ -0,0 +1,1152 @@ + + + + + +exercise_sheet_02 + + + + + + + + + + + + + + + + + + +
+
+
+

exercise_sheet_02 (Score: 99.0 / 100.0)

+
+
    + + + + + + + + + + + + + + + + + + + + +
  1. Written response (Score: 5.0 / 5.0)
    2. + + + + + + + +
  2. Test cell (Score: 5.0 / 5.0)
    3. + + + + + + + +
  3. Written response (Score: 5.0 / 5.0)
    4. + + + + + + + +
  4. Test cell (Score: 5.0 / 5.0)
    5. + + + + + + + + + + +
  5. Test cell (Score: 15.0 / 15.0)
    6. + + + + + + + + + + + +
  6. Comment
    7. + + + +
  7. Test cell (Score: 9.0 / 10.0)
    8. + + + + + + + +
  8. Written response (Score: 15.0 / 15.0)
    9. + + + + + + + +
  9. Coding free-response (Score: 10.0 / 10.0)
    10. + + +
  10. Comment
    11. + + + + + + +
  11. Written response (Score: 15.0 / 15.0)
    12. + + + + + + + + + + +
  12. Coding free-response (Score: 15.0 / 15.0)
    13. + + + + +
+
+
+
+
+
+ + +
+
+ +
+
+
+
+

Exercise sheet

Some general remarks about the exercises:

+
    +
  • For your convenience functions from the lecture are included below. Feel free to reuse them without copying to the exercise solution box.
    • +
  • For each part of the exercise a solution box has been added, but you may insert additional boxes. Do not hesitate to add Markdown boxes for textual or LaTeX answers (via Cell > Cell Type > Markdown). But make sure to replace any part that says YOUR CODE HERE or YOUR ANSWER HERE and remove the raise NotImplementedError().
    • +
  • Please make your code readable by humans (and not just by the Python interpreter): choose informative function and variable names and use consistent formatting. Feel free to check the PEP 8 Style Guide for Python for the widely adopted coding conventions or this guide for explanation.
    • +
  • Make sure that the full notebook runs without errors before submitting your work. This you can do by selecting Kernel > Restart & Run All in the jupyter menu.
    • +
  • For some exercises test cases have been provided in a separate cell in the form of assert statements. When run, a successful test will give no output, whereas a failed test will display an error message.
    • +
  • Each sheet has 100 points worth of exercises. Note that only the grades of sheets number 2, 4, 6, 8 count towards the course examination. Submitting sheets 1, 3, 5, 7 & 9 is voluntary and their grades are just for feedback.
    • +
+

Please fill in your name here:

+ +
+
+ +
+
+
In [1]:
+
+ 
NAME = "Kees van Kempen"
+NAMES_OF_COLLABORATORS = ""
+
+ +
+
+
+ +
+
+
+
+
+
+ +
+
+ +
+
+
+
+

Exercise sheet 2

+

Code from the lecture:

+ +
+
+ +
+
+
In [2]:
+
+ 
import numpy as np
+import matplotlib.pylab as plt
+from scipy.integrate import quad
+
+rng = np.random.default_rng()
+%matplotlib inline
+
+def inversion_sample(f_inverse):
+    '''Obtain an inversion sample based on the inverse-CDF f_inverse.'''
+    return f_inverse(rng.random())
+
+def compare_plot(samples,pdf,xmin,xmax,bins):
+    '''Draw a plot comparing the histogram of the samples to the expectation coming from the pdf.'''
+    xval = np.linspace(xmin,xmax,bins+1)
+    binsize = (xmax-xmin)/bins
+    # Calculate the expected numbers by numerical integration of the pdf over the bins
+    expected = np.array([quad(pdf,xval[i],xval[i+1])[0] for i in range(bins)])/binsize
+    measured = np.histogram(samples,bins,(xmin,xmax))[0]/(len(samples)*binsize)
+    plt.plot(xval,np.append(expected,expected[-1]),"-k",drawstyle="steps-post")
+    plt.bar((xval[:-1]+xval[1:])/2,measured,width=binsize)
+    plt.xlim(xmin,xmax)
+    plt.legend(["expected","histogram"])
+    plt.show()
+    
+def gaussian(x):
+    return np.exp(-x*x/2)/np.sqrt(2*np.pi)
+
+ +
+
+
+ +
+
+
+
+
+

Sampling random variables via the inversion method

(35 Points)

+

Recall from the lecture that for any real random variable $X$ we can construct an explicit random variable via the inversion method that is identically distributed. This random variable is given by $F_X^{-1}(U)$ where $F_X$ is the CDF of $X$ and $U$ is a uniform random variable on $(0,1)$ and

+$$ +F_X^{-1}(p) := \inf\{ x\in\mathbb{R} : F_X(x) \geq p\}. +$$

This gives a very general way of sampling $X$ in a computer program, as you will find out in this exercise.

+

(a) Let $X$ be an exponential random variable with rate $\lambda$, i.e. a continuous random variable with probability density function $f_X(x) = \lambda e^{-\lambda x}$ for $x > 0$. Write a function f_inverse_exponential that computes $F_X^{-1}(p)$. Illustrate the corresponding sampler with the help of the function compare_plot above. (10 pts)

+ +
+
+ +
+
+
+
Student's answer + Score: 5.0 / 5.0 (Top) +
+
+
+

Reasoning from the PDF, we can find the CDF and invert that as follows.

+$$ +f_X(x) = \lambda{}e^{-\lambda{}x} +$$$$ +\implies F_X(x) + = \int_{-\infty}^x f_X(t)dt + = \int_0^x \lambda{}e^{-\lambda{}t}dt + = \left[ -e^{\lambda{}t} \right]_{t = 0}^x + = 1 - e^{\lambda{}x} + = \mathbb{P}(X \leq x) = p +$$

for $x \in [0, \infty)$, otherwise zero.

+

Now we seek $x$ as a function of $p$.

+$$ +1 - e^{\lambda{}x} = p +\iff -\lambda{}x = \ln{(1-p)} +\iff x = \frac{\ln{(1-p)}}{-\lambda} = F^{-1}_X(p) +$$

which works, as $1 - p \geq 0$ as $p \in [0, 1]$, allowing $\ln{0} = -\infty$.

+ +
+
+ +
+ +
+
+
In [3]:
+
Student's answer(Top)
+
+
+ 
def f_inv_exponential(lam,p):
+    return -np.log(1 - p)/lam
+
+f_X = lambda x, lam: lam*np.exp(-lam*x) if x >= 0 else 0
+
+# Input parameters as list for flexibility in testing.
+for lam in [1.5]:
+    pdf = lambda x: f_X(x, lam)
+    samples = [inversion_sample(lambda p: f_inv_exponential(lam, p)) for _ in range(100000)]
+    compare_plot(samples, pdf, -1, 4, 30)
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
In [4]:
+
Grade cell: cell-2022e00546cf1bb0 + + Score: 5.0 / 5.0 (Top) +
+
+
+ 
from nose.tools import assert_almost_equal
+assert_almost_equal(f_inv_exponential(1.0,0.6),0.916,delta=0.001)
+assert_almost_equal(f_inv_exponential(0.3,0.2),0.743,delta=0.001)
+### BEGIN HIDDEN TESTS
+assert_almost_equal(f_inv_exponential(4.3,0.02),0.0046983,delta=0.0000001)
+assert_almost_equal(f_inv_exponential(0.2,0.98),19.560,delta=0.001)
+### END HIDDEN TESTS
+
+ +
+
+ +
+
+ +
+
+
+
+
+

(b) Let now $X$ have the Pareto distribution of shape $\alpha > 0$ on $(b,\infty)$, which has probability density function $f_X(x) = \alpha b^{\alpha} x^{-\alpha-1}$ for $x > b$. Write a function f_inv_pareto that computes $F_X^{-1}(p)$. Compare a histogram with a plot of $f_X(x)$ to verify your function numerically. (10 pts)

+ +
+
+ +
+
+
+
Student's answer + Score: 5.0 / 5.0 (Top) +
+
+
+ $$ +f_X(x) = \alpha b^{\alpha} x^{-\alpha-1} +\\ +\implies F_X(x) = \int_{-\infty}^x f_X(t)dt + = \int_{b}^x \alpha{}b^\alpha{}t^{-\alpha-1}dt + = \alpha{}b^\alpha{} \int_{b}^x t^{-\alpha-1}dt + = \alpha{}b^\alpha{} \left[ \frac{-t^{-\alpha}}{\alpha} \right]_{t = b}^x + = b^\alpha (b^{-\alpha} - x^{-\alpha}) = 1 - b^\alpha x^{-\alpha} = p, +$$

for $x > b$, otherwise $F_X(x) = 0$.

+

To find $F_X^{-1}(p)$, we write $p$ as function of $x$.

+$$ +p = 1 - b^\alpha x^{-\alpha} +\iff b^\alpha x^{-\alpha} = 1 - p +\iff x^{-\alpha} - b^{-\alpha}(1-p) +\iff x = \frac{b}{(1-p)^{1/\alpha}} +$$

Thus, $F_X^{-1}(p) = \frac{b}{(1-p)^{1/\alpha}}$ for $p \in [0, 1]$.

+ +
+
+ +
+ +
+
+
In [5]:
+
Student's answer(Top)
+
+
+ 
### Solution
+def f_inv_pareto(alpha,b,p):
+    return b/(1-p)**(1/alpha)
+
+# plotting
+f_X = lambda alpha, b, x: alpha*b**alpha*x**(-alpha-1) if x >= b else 0
+
+# Input parameters as list for flexibility in testing.
+for params in [(3., 1.)]:
+    alpha, b = params
+    pdf = lambda x: f_X(alpha, b, x)
+    samples = [inversion_sample(lambda p: f_inv_pareto(alpha, b, p)) for _ in range(100000)]
+    compare_plot(samples, pdf, b-1, 4, 30)
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
In [6]:
+
Grade cell: cell-726b321246679d28 + + Score: 5.0 / 5.0 (Top) +
+
+
+ 
from nose.tools import assert_almost_equal
+assert_almost_equal(f_inv_pareto(1.0,1.5,0.6),3.75,delta=0.0001)
+assert_almost_equal(f_inv_pareto(2.0,2.25,0.3),2.689,delta=0.001)
+### BEGIN HIDDEN TESTS
+assert_almost_equal(f_inv_pareto(0.1,3.5,0.3),123.90,delta=0.1)
+assert_almost_equal(f_inv_pareto(3.0,0.5,0.3),0.56312,delta=0.0001)
+### END HIDDEN TESTS
+
+ +
+
+ +
+
+ +
+
+
+
+
+

(c) Let $X$ be a discrete random variable taking values in $\{1,2,\ldots,n\}$. Write a Python function f_inv_discrete that takes the probability mass function $p_X$ as a list prob_list given by $[p_X(1),\ldots,p_X(n)]$ and returns a random sample with the distribution of $X$ using the inversion method. Verify the working of your function numerically on an example. (15 pts)

+ +
+
+ +
+
+
In [7]:
+
Student's answer(Top)
+
+
+ 
def f_inv_discrete(prob_list,p):
+    assert np.isclose(np.sum(prob_list), 1), "The probabilities should sum to one."
+    
+    p_cum = 0
+    i = 0
+    while p_cum < p:
+        p_cum += prob_list[i]
+        i += 1
+    return i
+
+# plotting
+f_X = lambda prob_list, x: prob_list[np.rint(x).astype(int) - 1] if np.rint(x) in range(1, len(prob_list) + 1) else 0
+
+# Input parameters as list for flexibility in testing.
+for prob_list in [[.1, .3, .2, .4], [0.5, 0.5], [0.7, 0.1, 0.2]]:
+    alpha, b = params
+    pdf = lambda x: f_X(prob_list, x)
+    samples = [inversion_sample(lambda p: f_inv_discrete(prob_list, p)) for _ in range(100000)]
+    compare_plot(samples, pdf, .5, len(prob_list) + .5, len(prob_list))
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+ +
+ + + + +
+ +
+ +
+ +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
In [8]:
+
Grade cell: cell-140af6b31464fbef + + Score: 15.0 / 15.0 (Top) +
+
+
+ 
assert f_inv_discrete([0.5,0.5],0.4)==1
+assert f_inv_discrete([0.5,0.5],0.8)==2
+assert f_inv_discrete([0,0,1],0.1)==3
+### BEGIN HIDDEN TESTS
+assert f_inv_discrete([0,0,0.1,0.3,0.5,0.1,0],0.42)==5
+assert f_inv_discrete([0,0,0.1,0.3,0.5,0.1,0],0.02)==3
+### END HIDDEN TESTS
+
+ +
+
+ +
+
+ +
+
+
+
+
+

Central limit theorem?

(35 Points)

+

In this exercise we will have a closer look at central limits of the Pareto distribution, for which you implemented a random sampler in the previous exercise. By performing the appropriate integrals it is straightforward to show that

+$$ +\mathbb{E}[X] = \begin{cases} \infty & \text{for }\alpha \leq 1 \\ \frac{\alpha b}{\alpha - 1} & \text{for }\alpha > 1 \end{cases}, \qquad \operatorname{Var}(X) = \begin{cases} \infty & \text{for }\alpha \leq 2 \\ \frac{\alpha b^2}{(\alpha - 1)^2(\alpha-2)} & \text{for }\alpha > 2 \end{cases}. +$$

This shows in particular that the distribution is heavy tailed, in the sense that some moments $\mathbb{E}[X^k]$ diverge.

+ +
+
+ +
+
+
+
+

(a) Write a function sample_Zn that produces a random sample for $Z_n= \frac{\sqrt{n}}{\sigma_X}(\bar{X}_n - \mathbb{E}[X])$ given $\alpha>2$, $b>0$ and $n\geq 1$. Visually verify the central limit theorem for $\alpha = 4$, $b=1$ and $n=1000$ by comparing a histogram of $Z_n$ to the standard normal distribution (you may use compare_plot). (10 pts)

+ +
+
+ +
+
+
In [9]:
+
Student's answer(Top)
+
+
+ 
def sample_Zn(alpha,b,n):
+    assert alpha > 2
+    assert b > 0
+    assert n >= 1 and type(n) == int
+    
+    E_X = alpha*b/(alpha - 1)
+    Var_X = alpha*b**2/( (alpha - 1)**2*(alpha - 2) )
+    
+    inv_pareto_samples = [inversion_sample(lambda p: f_inv_pareto(alpha, b, p)) for _ in range(n)]
+    return np.sqrt(n/Var_X)*(np.mean(inv_pareto_samples) - E_X)
+
+# Plotting
+alpha = 4
+b = 1
+n = 1000
+pdf = gaussian
+samples = [inversion_sample(lambda p: sample_Zn(alpha, b, n)) for _ in range(1000)]
+compare_plot(samples, pdf, -5, 5, 100)
+
+ +
+
+
+
+ Comments:

Use less bins to prevent overfitting

+ +
+
+
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
In [10]:
+
Grade cell: cell-5d16b014bef9d86f + + Score: 9.0 / 10.0 (Top) +
+
+
+ 
assert_almost_equal(np.mean([sample_Zn(3.5,2.1,100) for _ in range(100)]),0,delta=0.3)
+assert_almost_equal(np.std([sample_Zn(3.5,2.1,100) for _ in range(100)]),1,delta=0.3)
+### BEGIN HIDDEN TESTS
+assert_almost_equal(np.mean([sample_Zn(7.5,0.7,100) for _ in range(100)]),0,delta=0.3)
+assert_almost_equal(np.std([sample_Zn(7.5,0.7,100) for _ in range(100)]),1,delta=0.25)
+### END HIDDEN TESTS
+
+ +
+
+ +
+
+ +
+
+
+
+
+

(b) Now take $\alpha = 3/2$ and $b=1$. +With some work (which you do not have to do) one can show that the characteristic function of $X$ admits the following expansion around $t=0$,

+$$ +\varphi_X(t) = 1 + 3 i t - (|t|+i t)\,\sqrt{2\pi|t|} + O(t^{2}). +$$

Based on this, prove the generalized CLT for this particular distribution $X$ which states that $Z_n = c\, n^{1/3} (\bar{X}_n - \mathbb{E}[X])$ in the limit $n\rightarrow\infty$ converges in distribution, with a to-be-determined choice of overall constant $c$, to a limiting random variable $\mathcal{S}$ with characteristic function

+$$ +\varphi_{\mathcal{S}}(t) = \exp\big(-(|t|+it)\sqrt{|t|}\big). +$$

(15 pts)

+ +
+
+ +
+
+
+
Student's answer + Score: 15.0 / 15.0 (Top) +
+
+
+ \begin{align} +\phi_{Z_n}(t) &= \mathbb{E}\left[ e^{itZ_n} \right] + \\ &= \mathbb{E}\left[ e^{itcn^{1/3}(\bar{X_n} - \mathbb{E}[X])} \right] + \\ &= \mathbb{E}\left[ e^{itcn^{1/3}(\frac{1}{n}\sum_{i=1}^n X_i - \mathbb{E}[X])} \right] + \\ &= \mathbb{E}\left[ \left( \prod_{i=1}^n e^{itcn^{-2/3}X_i} \right) e^{itcn^{1/3}\mathbb{E}[X])} \right] + \\ &= \left( \prod_{i=1}^n \mathbb{E}\left[ e^{itcn^{-2/3}X_i} \right] \right)\mathbb{E}\left[ e^{itcn^{1/3}\mathbb{E}[X])} \right] + \\ &= \left( \prod_{i=1}^n \phi_X(cn^{-2/3}t) \right)\mathbb{E}\left[ e^{itcn^{1/3}\mathbb{E}[X])} \right] + \\ &= \left( \phi_X(cn^{-2/3}t) \right)^n \mathbb{E}\left[ e^{itcn^{1/3}\mathbb{E}[X])} \right] +\end{align}

where we used the identity for products of indepedent expectation values https://hef.ru.nl/~tbudd/mct/lectures/probability_random_variables.html#equation-product-expectation, and the definition of $\phi_X(t) := \mathbb{E}\left[ e^{itX} \right]$.

+

Next, we will use the Taylor expansion around $t = 0$ as is given above, and, for the latter exponential, $\mathbb{E}(X) = 3$ for $\alpha = 3/2, b = 1$ as given.

+\begin{align} +\phi_{Z_n}(t) &= \left( \phi_X(cn^{-2/3}t) \right)^n \mathbb{E}\left[ e^{itcn^{1/3}\mathbb{E}[X])} \right] + \\ &= \left( 1 + 3 i cn^{-2/3}t - (|cn^{-2/3}t|+i cn^{-2/3}t)\,\sqrt{2\pi|cn^{-2/3}t|} + \mathcal{O}(t^2) \right)^n e^{3itcn^{1/3}} + \\ &= \left( 1 + \frac{1}{n} \left[ 3 i cn^{1/3}t - (|ct|+i ct)\,\sqrt{2\pi|ct|} \right] + \mathcal{O}(t^2) \right)^n e^{3itcn^{1/3}} +\end{align}

Taking the limit $n \to \infty$, the first set of parentheses can be rewritten in terms of an exponential using the identity $\lim_{n\to\infty} (1 + \frac{a}{n})^n = e^{a}$, we find a way to our desired expression.

+\begin{align} +\lim_{n\to\infty} \phi_{Z_n}(t) &= \lim_{n\to\infty} \left( 1 + \frac{1}{n} \left[ 3 i cn^{1/3}t - (|ct|+i ct)\,\sqrt{2\pi|ct|} \right] + \mathcal{O}(t^2) \right)^n e^{-3itcn^{1/3}} + \\ &= \lim_{n\to\infty} \exp{({3 i cn^{1/3}t - (|ct|+i ct)\,\sqrt{2\pi|ct|}})} \exp{(e^{-3itcn^{1/3}})} + \\ &= \exp{({-(|ct|+i ct)\,\sqrt{2\pi|ct|}})} +\end{align}

This matches $\phi_S(t) = \exp\big(-(|t|+it)\sqrt{|t|}\big)$ for $\sqrt{2\pi} c^{3/2} = 1 \implies c = (2\pi)^{\frac{-1}{3}}$.

+ +
+
+ +
+ +
+
+
+
+

(c) The random variable $\mathcal{S}$ has a stable Lévy distribution with index $\alpha = 3/2$ and skewness $\beta = 1$. Its probability density function $f_{\mathcal{S}}(x)$ does not admit a simple expression, but can be accessed numerically using SciPy's scipy.stats.levy_stable.pdf(x,1.5,1.0). Verify numerically that the generalized CLT of part (b) holds by comparing an appropriate histogram to this PDF. (10 pts)

+ +
+
+ +
+
+
In [11]:
+
Student's answer + Score: 10.0 / 10.0 (Top) +
+
+
+ 
from scipy.stats import levy_stable
+
+def sample_Zn(alpha, beta, c, n):
+    assert n >= 1 and type(n) == int
+    
+    E_X = alpha*b/(alpha - 1)
+    
+    samples = [inversion_sample(lambda p: f_inv_pareto(alpha, beta, p)) for _ in range(n)]
+    return c*n**(1/3)*(np.mean(samples) - E_X)
+
+alpha = 3/2
+beta  = 1
+c     = (2*np.pi)**(-1/3)
+n     = 1000
+pdf = lambda x: levy_stable.pdf(x, alpha, beta)
+samples = [sample_Zn(alpha, b, c, n) for _ in range(10000)]
+compare_plot(samples, pdf, -5, 5, 100)
+
+ +
+
+
+
+ Comments:

You are a bit overfitting.

+ +
+
+
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+
+
+

Joint probability density functions and sampling the normal distribution

(30 Points)

+

Let $\Phi$ be a uniform random variable on $(0,2\pi)$ and $R$ an independent continuous random variable with probability density function $f_R(r) = r\,e^{-r^2/2}$ for $r>0$. Set $X = R \cos \Phi$ and $Y = R \sin \Phi$. This is called the Box-Muller transform.

+

(a) Since $\Phi$ and $R$ are independent, the joint probability density of $\Phi$ and $R$ is $f_{\Phi,R}(\phi,r) = f_\Phi(\phi)f_R(r) = \frac{1}{2\pi}\, r\,e^{-r^2/2}$. Show by change of variables that $X$ and $Y$ are also independent and both distributed as a standard normal distribution $\mathcal{N}$. (15 pts)

+ +
+
+ +
+
+
+
Student's answer + Score: 15.0 / 15.0 (Top) +
+
+
+

The coordinate transformation $T$ is defined as $T(\Phi, R) = (R\cos{\Phi}, R\sin{\Phi}) = (X, Y)$. As $T$ is invertible differentiable, we can write the equality between the joint probability density in both coordinate pairs as follows, using the Jacobian. +$$ +f_{X,Y}(x,y) \Big|\frac{\mathrm{d}x}{\mathrm{d}\phi}\frac{\mathrm{d}y}{\mathrm{d}r}-\frac{\mathrm{d}y}{\mathrm{d}\phi}\frac{\mathrm{d}x}{\mathrm{d}r}\Big| += f_{X,Y}(T(\phi,r)) \Big|\frac{\mathrm{d}x}{\mathrm{d}\phi}\frac{\mathrm{d}y}{\mathrm{d}r}-\frac{\mathrm{d}y}{\mathrm{d}\phi}\frac{\mathrm{d}x}{\mathrm{d}r}\Big| += f_{X,Y}(T(\phi,r)) \Big|-r\sin{\phi}\sin{\phi}-r\cos{\phi}\cos{\phi}\Big| += f_{X,Y}(T(\phi,r)) r += f_{\Phi,R}(\phi,r) += \frac{1}{2\pi}\, r\,e^{-r^2/2} +\\ \implies f_{X,Y}(x,y) = \frac{1}{2\pi}\,e^{-r^2/2} = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}\frac{1}{\sqrt{2\pi}}e^{-y^2/2} = f_{X}(x)f_Y(y) +$$ +using that $r^2 = x^2 + y^2$ in the second to last step. We conclude that $X$ and $Y$ are independent, and the factorization shows that they are both distributed as a standard normal distribution $\mathcal{N}(0,1)$.

+ +
+
+ +
+ +
+
+
+
+

(b) Write a function to sample a pair of independent normal random variables using the Box-Muller transform. Hint: to sample $R$ you can use the inversion method of the first exercise. Produce a histogram to check the distribution of your normal variables. (15 pts)

+ +
+
+ +
+
+
+
+

For the sampling of $R$, we take its PDF, calculate its CDF, invert it, and use the function inversion_sample to pull values for $R$.

+$$ +f_R(r) = re^{-r^2/2} +\\ \implies F_R(r) = \int_{0}^r te^{-t^2/2}dt = 1-e^{-r^2/2} +$$

integrating over values $r > 0$ as it cannot be negative, and using a substitution with $z := t^2$.

+

Now the inversion.

+$$ +p := F_R(r) = 1-e^{-r^2/2} \implies r = \sqrt{-2\ln{(1-p)}} +$$

for $p \in [0, 1]$. Do note that we can also calculate $r = \sqrt{-2\ln{(p)}}$ as $p$ and $1 - p$ are identically distributed, also keeping the argument to $\ln$ positive, saving just one calculation, although I do not use this in the following.

+ +
+
+ +
+
+
In [12]:
+
Student's answer + Score: 15.0 / 15.0 (Top) +
+
+
+ 
def random_normal_pair():
+    '''Return two independent normal random variables.'''
+    phi = rng.random()*2*np.pi
+    r = inversion_sample(lambda p: np.sqrt(-2*np.log(1-p)))
+    x, y = r*np.cos(phi), r*np.sin(phi)
+    return x, y
+
+# Plotting
+pdf = gaussian
+samples = np.array([random_normal_pair() for _ in range(100000)])
+for index in [0, 1]:
+    compare_plot(samples[:,index], pdf, -5, 5, 100) 
+
+ +
+
+ +
+
+ +
+
+ + +
+ +
+ + + + +
+ +
+ +
+ +
+ +
+ + + + +
+ +
+ +
+ +
+
+ +
+
+
+ + +
+
+
+
+
+ +