Skip to content

Instantly share code, notes, and snippets.

@adrn
Created October 5, 2019 13:42
Show Gist options
  • Save adrn/565f876cfc052d39ea48eee486e398c8 to your computer and use it in GitHub Desktop.
Save adrn/565f876cfc052d39ea48eee486e398c8 to your computer and use it in GitHub Desktop.
hq/notebooks/Simple-Linear-Nonlinear-Inference.ipynb
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "import matplotlib as mpl\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport numpy as np\nfrom scipy.stats import multivariate_normal",
"execution_count": 1,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "ndata = 16\nnpars = 3\n\ny = np.random.uniform(size=ndata)\nM = np.vstack((np.ones((1, ndata)),\n np.random.uniform(-1, 1, (npars-1, ndata)))).T\nC = np.diag(np.exp(np.random.uniform(-1, 0, size=ndata)))\nCinv = np.linalg.inv(C)\n\nth = np.array([1., 2., 3.])\nmu = np.random.uniform(size=npars)\nV = np.diag(np.exp(np.random.uniform(-1, 0, size=npars)))\nVinv = np.linalg.inv(V)\n\nZ = C.dot(M).dot(Vinv).dot(mu)\nX = C.dot(Cinv + M.dot(Vinv).dot(M.T)).dot(C)\n\nprint(y.shape, M.shape, C.shape)\nprint(th.shape, mu.shape, V.shape)\nprint(Z.shape, X.shape)",
"execution_count": 41,
"outputs": [
{
"output_type": "stream",
"text": "(16,) (16, 3) (16, 16)\n(3,) (3,) (3, 3)\n(16,) (16, 16)\n",
"name": "stdout"
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "eval_y = np.random.normal(y, size=(128, ndata))\neval_th = np.random.normal(th, size=(128, npars))",
"execution_count": 42,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "norm_y = multivariate_normal(mean=M.dot(th), cov=C)\nnorm_th = multivariate_normal(mean=mu, cov=V)\ntruth = norm_y.logpdf(eval_y) + norm_th.logpdf(eval_th)",
"execution_count": 43,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "Sig = np.linalg.inv(M.T.dot(Cinv).dot(M) + Vinv)\nphi = Sig.dot(M.T.dot(Cinv).dot(y) + Vinv.dot(mu))",
"execution_count": 44,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "norm_th2 = multivariate_normal(mean=phi, cov=Sig)\nnorm_y2 = multivariate_normal(mean=Z, cov=X)",
"execution_count": 45,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "test = norm_th2.logpdf(eval_th) + norm_y2.logpdf(eval_y)",
"execution_count": 46,
"outputs": []
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "plt.figure(figsize=(5, 5))\nplt.scatter(truth, test)\nplt.xlim(-220, 0)\nplt.ylim(-220, 0)",
"execution_count": 47,
"outputs": [
{
"output_type": "execute_result",
"execution_count": 47,
"data": {
"text/plain": "(-220, 0)"
},
"metadata": {}
},
{
"output_type": "display_data",
"data": {
"text/plain": "<Figure size 360x360 with 1 Axes>",
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAr0AAAJyCAYAAADJmni5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdf3ScV33v+48UK7GUyFZMLSm2MQmb/LCkcEIMNA6/E/k6bmnVFgOBmoTDKZflQrilh9NySld85V5YDS1NuIGrEziF+uDDj3uTA8YsjBslQBJsVsGBJhrbcdhgTKRINjiyVDRKlEj3j5mpR6P58Twzz+/n/VoryzPSnmd2JEw+s5/v/u6mhYUFAQAAAEnWHPYEAAAAAL8RegEAAJB4hF4AAAAkHqEXAAAAiUfoBQAAQOIRegEAAJB4hF4AAAAkHqEXAAAAiUfoBQAAQOIRegEAAJB4hF4AAAAkHqEXAAAAiUfoBQAAQOIRegEAAJB4hF4AAAAk3rKwJ4AcY8wbJL1f0mskLeS//IikT1trHwltYgAAAAnASm/IjDEXGGM+L+l+ST+UdIW1dp2kKyU9KulBY8x/M8a0hDlPAACAOGtaWFioPQq+MMYsk/Q1SW+W9AFr7WfKjPmgpE9J+l+S3matfSHYWQIAAMQfK73h2qVc4D0paajCmM9IekrSH0m6PaB5AQAAJAqhNyTGmKsl/UX+6ZettfPlxuVXdr+cf/pRY0xPEPMDAABIEkJveHZJOi//+OEaYwsb2c7Lvw4AAAAuEHpDYIx5iaTfL/rSozVecrjo8R8YY9Z5PysAAIDkIvSG460697OfttY+XW2wtXZU0m/yT8+TtM3HuQEAACQOoTccNxU9PunwNb8oerzVw7kAAAAkHqE3HK8uenza4Wt+VfT4tz2cCwAAQOJxIlvAjDGXSGov+tIzDl96pujxSmNMl7V2otzAV77ylT+fn5+/dH6+bEMIR5YtW1arzhgAAECS9Pzzz19b72ubm5vV3Nx84kc/+tFlXs6pFKE3eJeUPH/W4etKx62RVDb0SlrR1NTkalKlrrrqqrr/x+ulkZERSVJfX1/IM0HQ+N2nF7/7dOL3Hm+F3189mpqaND8/f6l3symP0Bu8C0ueP+fwdaXjSq9T7MSGDRtW7dmzx/msIsoYI0lKwr8L3OF3n1787tOJ33t6XXPNNYG8DzW9wZsreX5e2VFLlY4rvQ4AAAAqIPQGb7Lk+fkOX3dBjesAAACgAkJv8H4haaHoebUyhWJtRY8XtLiFGQAAAKog9AbMWpuV9GTRl1Y7fGnxuGPW2lnvZgUAAJBshN5wPFz0uMvha7qLHj/k4VwAAAASj9Abjq8XPV5njFlebbAxplXS2qIvfc2XWQEAACQUoTccBySN5h83S7qyxvirJBUa756UNOzTvAAAABKJ0BsCa+2cpH8o+tL1NV5S/P1PWmtf8H5WAAAAyUXoDc+nJR3LP76pxtjC949IGvJtRgAAAAlF6A2JtfY5Se+QNCPpd40xl5YbZ4y5TNJWSb+RdHN+lRgAAAAuEHpDZK39iaQ/kDQr6SvGmPbi7xtjVkj6qnLB+M3W2seDnyUAAED8LQt7Amlnrb3fGPNqSf+PpOPGmN3KHTzxEkm3SspIeqW19niI0wQAAIg1Qm8EWGuPSHqjMaZX0mslvUjSzyTdYK09GurkAAAAEoDQGyHW2oxyK7sAAADwEKEXkWatDXsKCAm/+/Tid59O/N7Tq6+vT8eOHXvU7/dhIxsAAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASb1nYEwAAJM/oZFbDRyY0lZ3TitYWbe7p0pqO1rCnBSDFCL0AAM+cmp7Vzr0ZHciMa37h3NcH92W0pbdbgwO96mxfHt4EAaQWoRcA4IlT07PaNnRIJ8/MLPne/IK0f2RcmbEp3bfjeq1uvyCEGQJIM2p6AQCe2Lk3UzbwFjt5Zka37x0JaEYAcA6hFwDQsNHJrA5kxh2NPZAZ19hk1ucZAcBihF4AQMOGj0wsquGtZn5BGj464e+EAKAEoRcA0LCp7Jyr8Wdn3I0HgEYRegEADVvR2uJq/Mo2d+MBoFGEXgBAw/p7utTc5Gxsc5PUv6HL3wkBQAlCLwCgYWs7WrWlt9vR2C293RxUASBwhF4AgCcGB3q1flVb1THrV7Vp10BfQDMCgHMIvQAAT3S2L9e9OzZpa1/3klKH5iZpa183B1MACA0nsgEAPNPZvlxD2zdqbDKr4aMTOjszp5VtLerf0EVJA4BQEXoBAJ5b09GqWzZdGvY0AODfUd4AAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPEIvQAAAEg8Qi8AAAASj9ALAACAxCP0AgAAIPE4htgBY0yPpJXW2kM1xjVJ+kNJ37HWPhPI5AAAAFATodeZT0t6gzHm25K+KumQpF9KWpDUIekKSW+Q9A5JWUn3O72wMeY8Sdsl3SKpV9ILkmYkfVvSJ621Jzz7twAAAEgpyhucac7/8zuSdks6rly4nZU0LukhSX+T/9pWa+20k4saY9ZJ+r6kjysXrNdZa9dKeqOkFZKOGWPe6+m/CQAAQAqx0uuNOUl3Svpra+2ckxcYYy6R9D1J6yRtstY+WvietXZU0q3GmIskfdYYc5G19k4f5g3AhdHJrIaPTGgqO6cVrS3a3NOlNR2tYU8LAOAAode5e5QraeiX9FvKlTaMSnpE0hettSedXihf0vBVSS+V9D+LA2+Jv5L0R5L+3hjzuLV2uIH5A6jTqelZ7dyb0YHMuOYXzn19cF9GW3q7NTjQq8725eFNEABQE6HXuXFr7cckfcyDa71P0uvyj/dUGmStfcIY80NJr1JuxXeDtfZZD94fgEOnpme1beiQTp6ZWfK9+QVp/8i4MmNTum/H9VrdfkEIM0wvVt4BuEHoDZgx5gJJt+efzks6WOMljygXei+T9F7lan8BBGTn3kzZwFvs5JkZ3b53REPbNwY0q3Rj5R1APdjIFry3SurKP/6ptXaqxvjDRY8/4M+UAJQzOpnVgcy4o7EHMuMam8z6PCMUVt73jywOvNK5lfdtQ4d0epqbYgAWI/QG7+aix8cdjC8ec6Ux5hqP5wOgguEjE0uCVSXzC9Lw0Ql/JwRXK+8AUIzyBheMMecr10/3DyW9OP/Pc5KelPTPku6x1lb8r54xpkXSDUVf+oWDtz1R8nyrpJ84nzWAek1lHTVj+XdnZ9yNhzv1rLxT4wuggJVe53okPSxpuaS3WWtfbq29WNKAcj/HQUk/M8b871Wu8XJJxf8PfNrB+55Rrva34DpXswZQtxWtLa7Gr2xzNx7usPIOoBGEXufeKul2a+2nrbW/KXzRWvsD5VZvfyKpTdI9xpgPVbjGFSXPax5VbK19QVJx3W/pNQD4pL+nS81NzsY2N0n9G7pqD0TdWHkH0AhCr3P3W2sPlPuGtXZW0keKvvQJY8wrygy9pOS5050WxePWOHwNgAat7WjVlt5uR2O39HZzK91nrLwDaAQ1vc48q9zpadXcL+lXyh1csUzSX0t6S8mYC0ueP+fw/YvHlV6jrJGRERljHF5+KWtt3a8FkmRwoFeZsamqm6fWr2rTroG+AGe1WJT71Xo5t/6eLg3uyzgqcWDlHQhWI5mjvb1dy5b5H0kJvc78hXKnsVVkrZ3PHySxNf+lNxtjWq21xT2MSu+1nefw/YvHcb8OCFBn+3Ldu2NT2b6wzU25Fd5dA32hHEwR5X61fsytsPK+f6T2ZjZW3gGUimToNcasUG7jWBh+XHrqmbX2Xx2+9mdFj8+XdK2k7xd9bbJk/PkOr1v8X9PSa5TV19enPXsqHvYGwIXO9uUa2r5RY5NZDR+d0NmZOa1sa1H/hvBWVKN8Upyfc4vDyjuQRo3cId6+fbuOHTvm4WzKi2ToVS4sfiek975MS9uEOTVd8rz03lrpdR2VKii3Qa7g524mBMA7azpadcumS8OehqRonxTn59yivPIOINqiGnoflbQppPd+uvQLxph2SW+SNGOtHa7y2vmS56X37kpXjFfXmowx5kItbnNGj14g5aLcrzaIuUVx5R1A9EUy9OaP5v1B2POQJGPMyyQ9IGl9/vnd1toPVhheunL7q+In1tqnjTFWUqHa28kui9Kt4w85eA2ABKunX21QK9RBzi1KK+8Aoo+WZbV9QvnAm3drlbGl7cQeLzPm60WPL3fw/sVjnpP0LQevAZBgUe5XG+W5AUg3Qm9tryl5Plpl7FVFj39srV1SKiHp88XjjTG1Wt9vKHp8b34VHECKRblfbZTnBiDdCL21lW4iK7vSa4xZJam36EufLjfOWntE0r7804u1OCiXc33+z3lJf1drsgCSz4+T4kYns9p98ITufuBJ7T54QmOT2ZqvCWpuAOCFSNb0RswRSa9Wrsb4RmttpS3Jb9O5DxEPSdpd5ZoflrRZuY1uN0k6Wm6QMeYC5TbQSdI/WmvZxAbA0361XvfTpZcugKhipbe2e/J/PlUp8ObD6Z/nn56W9E5r7QuVLmitPS6psBnufcaYSh8+3i7pRcoF7w+5nTiA5Boc6NX6VW1Vx9TqV1vop7t/ZHzJ5rNCP91tQ4d0etrpienezQ0AvEbore0Lkr4qaYsxZkkpQj6wfk65DWePSbrOWlut7leSZK39nKSPSrpS0l3GmEW/C2PMBkl3STomabO19jeN/osASI5Cv9qtfd1Lygmam6Stfd01D39w00836LkBgNcob6jBWrtgjPljSX8m6VvGmGHlujLMSlon6R2SOiX9raS/qVL+UO7aHzfGPC7pTkk/MsZ8TdIzkq6W9E5JX5L0YWtt6aEXANBQv1q/++nSSxdA1BB6HciXKnzSGPMpSa9XblV3taQzkv5U0iPW2tk6r73PGPMt5bpEXCOpXdLDknZaa539FwlAqtXTrzaofrr00gUQFYReF6y1z0t6MP+Pl9d9QbnNbxw8ASAQ9NMFkDbU9AJACtFPF0DasNILABE1OpnV8JEJTWXntKK1RZt7vKuH7e/p0uC+jKMSB/rpAkgCQi+AVPMzWNbL69655dBPF0DaEHoBpFIQwbLeeW0bOlS2lVihd25mbMqTll+DA73KjE1VbVtGP10ASUFNL4DU8etQBi/41Tu3HPrpAkgTVnoBpI6bYDm0fWNAs/K/d2459NMFkBaEXgCpEkawdCqo3rnl0E8XQNJR3gAgVeoJlkGhdy4A+IeVXgCpEnSwdNMdgt65AOAfQi+AVAkqWNbTHYLeuQDgH8obAKRKf0/Xkk4FldQbLOvtDlHonesEvXMBwB1CL4BUCSJYNtJ2bHCgV+tXtVV9Lb1zAcA9Qi+A1PEzWNbTHaIYvXMBwB/U9AJInUKwLFdz29yUW+HdNdBXV7D0ou0YvXMBwHuEXgCp5Few9LI7BL1zAcA7hF4AqeZ1sKTtGABEEzW9AOChILpDAADcI/QCgIdoOwYA0UToBQCP0XYMAKKH0AsAHqPtGABEDxvZAMAHtB0DgGgh9AKAj2g7BgDRQHkDAAAAEo/QCwAAgMSjvAFAIoxOZjV8ZEJT2TmtaG3R5h5qZwEA5xB6AcTaqelZ7dyb0YHMuOYXzn19cF9GW3q7NTjQq8725eFNEAAQCYReALF1anpW24YO6eSZmSXfm1+Q9o+MKzM2RXswj7GqDiCOCL0AYmvn3kzZwFvs5JkZ3b53REPbNwY0q+RiVR1AnLGRDUAsjU5mdSAz7mjsgcy4xiazPs8o2Qqr6vtHFgde6dyq+rahQzo9/Ww4EwSAGgi9AGJp+MjEkvBVyfyCNHx0wt8JJZybVXUAiCJCL4BYmsrOuRp/dsbdeJzDqjqAJCD0AoilFa0trsavbHM3Huewqg4gCQi9AGKpv6dLzU3OxjY3Sf0buvydUIKxqg4gCQi9AGJpbUertvR2Oxq7pbeblloNYFUdQBIQegHE1uBAr9avaqs6Zv2qNu0a6AtoRslUz6r66GRWuw+e0N0PPKndB09Q5wsgdPTpBRBbne3Lde+OTWV7xzY35VZ4dw30cTBFgwqr6vtHam9me+OVq/U33zxCL18AkUPoBRBrne3LNbR9o8Ymsxo+OqGzM3Na2dai/g2cEualwYFeZcamqrYtW9vRqifG/02jZVZ1OSEPQNgobwCQCGs6WnXLpkt1242X65ZNlxJ4PVZYVd/a172k1KG5Sdra160ruy8qG3iL0csXQFhY6QUAOFJtVX1B0uvueNDRdQq9fPlgAiBIhF4gRUYnsxo+MqGp7JxWtLZocw8lAHCvsKpebPfBE657+ZZeI0j8XQDSh9ALpMCp6dmym73YXASvxKWXL38XgPQi9AIJd2p6VtuGDpXdgMTmInglDr18+bsApBsb2YCE27k3U3XHvcTmIjQuDifk8XcBSDdWeoEEG53M6kCmdm9Vic1FaIybXr5hnJDn198FaoOB+CD0Agk2fGQiVpuLEG9OevmGdUKe138XqA0G4ofyBiDB4rK5CMngpJdvWPWyXv5dKNQG7x8ZXxKkC7XB24YO6fT0s/VMFYBPWOkFEiwOm4uQLFE9Ic/LvwtuaoOHtm909b4A/MNKL5BgcdhchGRa09GqGzd0aUVri87OzOn+IxMaq3Fam5+8+rtQT20wgGhgpRdIsKhvLkIyRbHe1au/C9TJA/HFSi+QcIMDvVq/qq3qmLA2FyF5olzv6sXfBerkgfgi9AIJF+XNRUieKPfC9eLvAnXyQHxR3gCkQFQ3FyFZ4tAXutG/C/09XRrcl3FU4kCdPBAthF4gRdZ0tFJfCN9Eqd611qER9f5doE4eiC9CLwDAE1Godw1iE12UD+EAUBk1vQAAT4Rd7xrUJjrq5IF4YqUXAOCJsOtdgzw0gjp5IH4IvQAAT4RZ7xrWJjrq5IH4oLwBAOCZsPpC17OJDkC6EHoBAJ4Jq941CpvoAEQb5Q0AAE+FUe8a9iY6ANFH6AUA+CLIetewN9EBiD5CLwAErNbBCXCPQyMA1ELoBYCABHFwQppxaASAatjIBgABCOrghDTj0AgA1bDSCwABCPLghDTj0AgAlRB6AcBnYR2ckGYcGgGgFOUNAOAzDk4AgPARegHAZxycAADhI/QCgM84OAEAwkfoBQCf9fd0LekmUAkHJwCAPwi9AOCzwsEJTnBwAgD4g9ALAAEYHOjV+lVtVce0tpynD974soBmBADpQugFgAB0ti/Xf3vXtWptOa/imOzcC3rfFx/lgAoA8AGhFwACcvcDP1V27oWqYwoHVAAAvEXoBYAA1HNABQDAO4ReAAgAB1QAQLgIvQAQAA6oAIBwEXoBIAAcUAEA4SL0AkAAOKACAMJF6AWAAHBABQCEi9ALAAFxckDF+lVt2jXQF9CMACA9CL0AEJDO9uW6d8cmbe3rXlLq0Nwkbe3r1n07rtfq9gvCmSAAJNiysCcAAFE1OpnV8JEJTWXntKK1RZt7uhouO+hsX66h7Rs1NpnV8NEJnZ2Z08q2FvVvaPzaAIDKCL0AUOLU9Kx27s3oQGZ8UW/dwX0Zbent1uBArzrblzf0Hms6WnXLpksbmygAwDFCLwAUOTU9q21Dh3TyzMyS780vSPtHxpUZm6IMAQBihppeACiyc2+mbOAtdvLMjG7fOxLQjAAAXiD0AkDe6GRWBzLjjsYeyIxrbDLr84wAAF4h9AJA3vCRiUU1vNXML0jDRyf8nRAAwDPU9AJA3lR2ztX4szPOx/vRCQIA4ByhFwDyVrS2uBq/sq32+CA6QQAAaqO8AQDy+nu6lhwaUUlzk9S/oavqmEIniP0j40vKJgqdILYNHdLp6WfrnDEAwClCLwDkre1o1Zbebkdjt/R21yxPoBMEAEQHoRcAigwO9Gr9qraqY9avatOugb6qY+gEAQDRQk0vABTpbF+ue3dsKluH29yUW+HdNdBX82CKejpB1HtCWxw3ycVxzgDijdALACU625draPtGjU1mNXx0Qmdn5rSyrUX9G5wHMz87QRTEcZNcHOcMIBkIvQBQwZqO1rpXX/3oBFEsjsclezVnv1eJWYUGkonQCwA+6O/p0uC+jKMSh2qdICoFMDeb5Ia2b6znX8Fzjc7Z71ViVqGBZCP0AoAPCp0g9o/U3sxWrhNEtQD2+itW63tPnHY0j8ImubBXKuvZ2Fc8Z79XtuO4cg7AndR0bzDGvM0Y43BbSdXr/IEx5pvGmKeNMU8ZY04YYz5vjLk6StcEEL56O0HU6u/73SdOy+n/mUXluORGj3j2u/0b7eWA5EtN6JX0rkZebIzpMMZ8U9JuSV+T9BJr7TpJ10r6taSfGGN2hn1NAPUbncxq98ETuvuBJ7X74ImG24gVOkFs7etecuhFc5O0ta+77MqhkwDmRj2b5LzWyMY+v9u/0V4OSIdUlDcYY94j6c0NvP4iSf8s6VWSfs9a+83C96y1ZyT9F2PMeZL+T2PMxdbaPwvjmgDq42ctp9tOEG4CmFNuN8n5oZGNfX63fwuyvRyA8CR6pdcYs84Y8zFJn23wUp9VLpx+vzicltglKSvp/zDG/MeQrgnApaCOCi50grjtxst1y6ZLK9bYuglgTjg5LjkIjRzx7Hf7tyDaywEIX+JCrzHmK/ma2DFJI5I+JOm8Bq73O5LekX+6p9I4a+2kpEJ4vdMYszrIawKoT9RqOd0GsFqcHJcchEaOePa7/Zvf1wcQDYkLvdbam621l1pr11hrOyT9aYOX/FjR44drjH0k/+dKSX8Z8DUBuBTFWk63AawaJ8clB+XU9KyenXuh5rhyc25kldgJv68PIBoSF3q9ZIx5raRr8k9nJR2p8ZLDRY//kzFmyfKKH9cEUJ9GOwr4wW0Ae+MVq11tkgtDoYTkwRpt1m64qrPsnBtZJXbC7+sDiIZUbGRrwM1Fj39qra31n8fjRY87JG2R9PUArgmgDlGs5XTb37eR45KDOnnMaTeKC5Y1VwzpgwO9yoxNVb1OIyvbfl8fQPgIvdXdVPT4F7UGW2tPG2NmJBUac27V0oDqxzUB1CGqtZxuA5jb45KDPHms0UMpCgrt38rNu7kp9wFg10Bf3Svbfl8fQPgIvRUYY14kyRR9ydnxR9KvJK3PP77O72sCqJ9XRwV7zc8AFvTJY25LSO47/JRuu/Hyst932/7NLb+vDyBchN7Krih5/ozD153RuYB6uTGmqaiEwY9rljUyMiJjTLUhVVlr634tEBeNHhXsJ78CmJtuFUPbN9b9PgVuS0j+4f7jOvL0VNXV5uKV7dHJrO73uETD7co5ADWUOdrb27Vsmf+RlNBb2SUlz5026Swe16pcHW4h3PpxTQANiHotp5cBzKtSAzfclpAsyNlqc5AlGgCSgdBb2YUlz59z+LrScRfqXED145pl9fX1ac+eii2AAeSlqZYzjJPH3JSQFKu22hx0iQaA2hq5Q7x9+3YdO3bMw9mUR+itrPSenNMDLkrHzVV47NU1ATQoLbWcYXSrcFNCUqrSanPQJRoAkqHh0GuMWSGpx4O51OPH1trGzgatbLLk+fkOX1e6rHDW52sC8EjSaznD6lbhpISknHKrzWGUaABIBi9Weq+V9B0PrlOPyySd8OnapdctLU2opK3o8dPW2lmfrwkAjoTVraK4hOTbI+NyU+lQutocRokGgGTwIvQ+KmmTB9epx9M+Xvunkor74652+LricT8J4JoA4EiY3SoKJST/9wNP6h/uP177BXmlq81RPFAEQDw0HHqttVOSfuDBXCLFWvu8MeaQpBvzX6q55GGMaZG0quhLD/l9TQBwI+xuFW/ZuE53DR+ve7U5qgeKAIi+5rAnEHHFJ5+V75a+mNHin+nXAromgAgbncxq98ETuvuBJ7X74AmNTWZDm0uh1GBrX7eamxZ/r7lJ2trX7WvXg8JqsxPlVpv7e7qWzLuSIA8UARB9dG+o7suSPqFcb9xLjDEd1trSzWjFNhQ9PmitfSKgawKIoKj2kg27W0Ujq81RPlAEQLSx0luFtfbXkv570Zdq1S5fX/T4jqCuCSB6Cr1k94+ML7mVX+glu23okE5P+9WAprZCt4rbbrxct2y6NLCA2Ohq8+BAr9avaiv7vYIwDxQBEE2s9NY2KOkdkn5L0k2S9lcZe1P+z/uttd8I+JoAIoRestU1stqcpgNFAHiH0FuDtfbXxph3SfqmpHcZY/7aWjtdOs4Y83pJfZLGJd0a9DUBRAe9ZJ2rtzdy2CUaAOInDeUNncVPjDGdlQZWYq39tqT3SGqX9D+MMYsOlTDGXCLpn5RrodZvra3ZSs2PawLwn5NNafX0kkV9wirRABA/iVvpNcZ8RdJ1+adNktaVDLHGmF8XPb/LWntXretaa/+HMeankj4jKWOM+bJyK7BXSHq3pPsl3WatdXzWph/XBOAPN5vS4t5LdnQyq+EjE5rKzmlFa4s297B6CiD+Ehd6rbU3+3jtg8aYayW9UtKrJXVI+ldJ11hrT0TlmgC8VdiUVq5Gt7ApLTM29e+br+LaSzaq3SYAwAuJC71+s9YuSPph/owW4AMAACAASURBVJ/IXhOAd9xuSgvruN9GuA32ABA3aajpBYC61bMprdEDGMLgJtgDQBwRegGgino3pcWpl2w9wR4A4obQCyBxvDz2t95NaWEf9+sG3SYApAE1vQASw4+NWI1sSotLL9m4d5sAACcIvQASoZGNWNVadHmxKa3eAxiC4me3CdqfAYgKQi+ARKjn2F8nK8OFTWn7R2rXvEZlU5pbfnSboP0ZgKihphdA7NWzEauwMrx/ZHxJ2CusDG8bOqTT08/GalNaPbzuNuHmZwsAQSH0Aoi9ejZiuVkZjtOmtHp5GexpfwYgiihvABB7bjdi/fLMjOuV4TUdrbHYlFavQrAvV5LQ3JRb4d010Fcz2Nez6p6Enx+A6CP0Aog9txuxRiezrleGCxvRor4prRFedJuoZ9U9qT9PANFC6AUQe243Yq11ubKYthZdjQR72p8BiCpqegHEntuNWOsurl67WspNi66087P9GQA0gtALIBHcbMTq7+lasiGtEqctupDDzxZAVBF6gQTy8hjeuHDTYcHrFl04h58tgKiiphdIkLQfCOBmI9bgQK8yY1NVW2vFufdumPjZAogiVnqBhOBAgHMKG7Fuu/Fy3bLp0rKriWnovRsWfrYAooiVXiAh6jmGN+28aNGF8vjZAogaQi+QABwI0Jgk994NGz9bAFFBeQOQAPUcCAAAQJoQeoEE4EAAAACqo7wBSAAOBMgZncxq+MiEprJzWtHaos091I8CAHIIvUACuD2GN2kHAjz21KQ+ct9jOvr0tIp/BGlp1QYAqI3yBiAB0nogwKnpWb3nn/5Fv//p7+tISeCV0teqDQBQGaEXSAg3x/AmQaEv8YPHTtccW2jVBgBIL0IvkBBpOxDgI/c9VrMvcbFCqzYAQDpR0wskSFoOBHjsqUlHK7zFCq3a6BkLAOlE6AUSKOkHAnzkvsfqeh2t2gAgvQi9AGJldDKro09P1/XapLZqi6KotI+LyjwAhI/QCyBWho9MLOnS4EQSW7VF0anpWe3cm9GBzPiiFnpBt4+LyjwARAehF0CsuD19riBJrdqiqtBRo9wGw0L7uMzYlO8bKqMyDwDRQvcGALHi9vQ5Kf6t2kYns9p98ITufuBJ7T54IrJdKHbuzdTsqBFE+7iozANAtLDSCyBW3Jw+J0nXmxfpUze/IpYrenG6RT86mdWBzLijsYX2cX6svEdlHgCih5VeALHi5vS561/6In3pvdfFNvBuGzqk/SPjSwJ+FE+aGz4y4fiDSKF9XJLnASB6CL0AYsfp6XOfescrApqR9+J2i95trbVf7eOiMg8A0UPoBRA7ST99rp5b9GEanczq6PiUq9c4bR/ntp7Zbc03beyA9KCmF0As+XH6XFR6utZziz6Mw0gq1RzX4qR9XL31zG5qvmljB6QLoRdArHlx+lzUNozF4RZ9tbZgtdRqH9dIy7FCzff+kdor5bSxA9KF8gYAsddIS68obhiLwy16JzXH5ThpH9doPbPTmu84t7ED4B4rvQBiy4sVWjcBa2j7Ri+mXVPUb9G7qTkuaG7KrazuGuhbtDpbWlJy9bqVDbccK9R8l/vfRqV5AEg+Qi+AWPLi1K2o9nSN+i16NzXHkvS7V3fro7/bs2ielT6wNEmOj5muVs/sR803gHgj9AKIJS9WaKO8YWxwoFeZsamq/45h3aJ3W3N8VfeKJYG30gcWF1laUu16Zi9qvgEkAzW9AGLHq5ZeUd4wFuW2bI3WHNdbD+zk2gBQCSu9AGLHqxXaqG8Yi+ot+kZqjuupB3Z6bQCohtALIFZGJ7P63vHTrl5TaYU26hvGCqJ2i76RmmO39cBurg0A1RB6AcRCvQchSJVXaKO+YSzK6q05dltS4ubaAFANNb0AIq9aL91aaq3Q0tO1PvXWHLstKSm5dOj1zADii5VeAJHXyManWiu09HStXz01x25LSu7bcb0eHz0bmXpmAPFF6AUQaY1sfHK6QhvVDWO1lB7ssLknnPm6qTl2W1LyivUX6xXrL25whgBA6AUQcfVsfKp3hTZqG8Yq8eIkujBFuQcxgOSiphdApLnd+HTDVZ165C9v0ND2jYksSahW31w4iW7b0CGdnn42nAk6EOUexACSi5VeAJHmduPTG69cHemShEZ5cRJdFMS1pARAfBF6AURaXHrpBqGek+iiHiDjUlICIP4obwAQaYWNT04kvZduPSfRAQByCL0AIo9eujlu65srnUQHAGlE6AUQeWx8ynFb39xc+sMCgBSjphdALLDxyV19syR9/8lf6f1vepm/kwKAmCD0AoiVNG98WtvRqtdfsVrffeK0o/EHf/ZrV5vZonLYBQD4gdALADFyRedFjkOvJN336FO67YbLq46J+2EXAOAEoRcAQlLPyupPT/2bq/f48cnJqt8vHHZRrvdv4bCLzNhUKmqmASQboRcAAtbQymqTt5vTknLYBQDUQugFkEph1a82urL6H17coQePnXL8fq94cUfF7yXxsAsAqITQCyBVwq5fbXRlddvGdbrz/uOO3+8tG9dV/F49h12kdRMhgPijTy+A1Cissu4fGV8S9gqrrNuGDun09LO+vH89K6ul1na06oarVju6xg1XdVZdmeWwCwBpQugFkBpuVln94NUxwn/7lpdrbY0yg7UdrbrjLS+vOsbtYRcr29yNB4AoIfQCSAUvVlmrXXv3wRO6+4EntfvgiYqv9WpltbN9ub72/uurnlD39fe/pma3hf6eriWvr6S5Serf0OVsMABEEDW9AFLBj/pVt/XBXq6senFC3dqOVm3p7db+kdofBrb0drOJDUCsEXoBpILX9av1dGFwc4xwk6Sr166sOa7RE+oGB3qVGZuqWvaxflWbdg301f0e9eKEOABeorwBQCp4Xb9aT31wYWXViQVJbxk6qB17DuvU9Kyj19Sjs3257t2xqWqpRNAHU5yantWOPYf1ujse1M5vZPTJ+49r5zcyeu0dD/r+8wCQXKz0AvBF1Fbp3Kyy1qpfbaS/rZOV1YKgTkTzolTCK5wQB8AvhF4Angq7D24l9dSvVgrujdQHF1ZWy/2MKgnqRLRGSyW8wAlxAPxC6AXgmaiv0jmtX/3gjZdrx57DFYP7S1a1uXrf0vrgwsrqo794Rm8ZOign+TkNJ6JxQhwAP1HTC8AzYffBrcVJ/eo979qo933xcNUDLL76o1+6et9K9cGPj551FHgL712pb29SeNXHGADKYaUXgCfCWqVzWztcq351x57DNYP7My5OJqtWH8yJaIvx8wDgJ0IvAE/40Qe3mkZrh8vVr7oJ7k5V62/LiWiL8fMA4CfKGwB4IshVukLtcLUShG1Dh3R6+llX13UT3CXp4hqhq1Z/W05EW4yfBwA/EXoBeCLIVTq/aofdBve3verFDfW3ddO3Nw0novHzAOAnyhuAFPKjh66XfXCr8bN22G1wX9vRqv+6dUND/W2jfCJaGPh5APALoRdIET976NbTB7ceftYO1xvcG+lvW61vb3NT7me1a6AvNQcx8PMA4BdCL5ASQfTQDWKVzs/a4aCCe6konYgWBfw8APiB0AukRBAnXQWxSud37XCYt9ejcCJalPDzAOAlQi+QAkH20PV7lc7v2mFurwNAMhF6gRQIuoeu5N8qXRAlCNxeB4DkIfQCKZCkk65OTc/q2edfqDnOixIEbq9Hgx/dRgCkD6EXSIGknHRVbTNesRuuWq073vIfIluCUBzi5hcW1NTUpCYploHOz0DqZ7cRAOlD6AVSIKgeun5zshlPki5Ydl4kA2+lEFcsLoHO70AaRLcRAOlC6AVSwKs62EZX9Rp5fZCb8fzgdJU6DoEuiEAaRLcRAOlC6AVSopFWXI2u6nmxKhjGZjwvOV2lLohyoPM7kMb9Aw6AaGoOewIAglFoxbW1r1vNTYu/19wkbe3rLrsyV1jV2z+y9JZ8YVVv29AhnZ5+tuz7Nvr6gjhvxnMT4ooVAl2U1BNI3arnAw4A1MJKL5Ai9bTianRVz6tVwThvxnMT4orNL0h/svuH+qf3vDoy9b1BrLjH+QMOgOhipRdIoUIrrttuvFy3bLq0ag1vI6t6Xq4K9vd0LVmhriRqm/HchrhiR56edrQSHpQgAmmcP+AAiC5CL4CKGr3N7OVt6sJmPCfqPZTCL25DXKnCSngUBBFI4/wBB0B0EXoBVNToqp7Xq4KDA71av6qt6hgvDqXwmpsQV0lU6nuDCKRx/oADILoIvQAqanRVz+tVwXo344XNTYirJCobtoIKpHH9gAMgutjIBqCiRg+18ONQjHo240WBk5ZxtTSyYcvLk9MaaX/nVOEDTrlWd81NuUC9a6Avch9wAEQXoRdARY0eauHVoRjlFDbjxUW1EOdUPfWxXp6cVhyc3/bKdfrRL57RQ8dP+xZI4/oBB0A0EXodMMb0SFpprT1UY1yTpD+U9B1r7TOBTA7wWaOrekGsCsZFaYj75ZkZfe7hnzt6bT31sV6dnFYpODc3Sa+/fLVeddnFmp+Xb4E0bh9wAERTakKvMeZtkr5qra1nO8mnJb3BGPNtSV+VdEjSLyUtSOqQdIWkN0h6h6SspPtdzOs8Sdsl3SKpV9ILkmYkfVvSJ621J+qYL+CZRm8zc5t6qeIQ99QzWV9WwiVveiTXCs7fPX5aP/vVbyJZSw0AxVITeiW9q4HXNuf/+Z38P5UclrTVWjvt5KLGmHWS7pX0YkkfkLTPWvu8MWatpI9LOmaMuc1a+7kG5g40rNHbzNymrsyvlXCvjvL1+8hhAAhKKkKvMeY9kt7s41vMSbpT0l9bax3tNDHGXCLpe5LWSdpkrX208D1r7aikW40xF0n6rDHmImvtnT7MG3Cl0dvM3KZeyq+VcC9OTvMqOANAFCQ69OZXUndI+ksPLnePciUN/ZJ+S7nShlFJj0j6orX2pIt5nadcmcRLJf3P4sBb4q8k/ZGkvzfGPG6tHW5g/gAiyo+VcC96JAdx5DAABCVxodcY8xVJ10k6X1Jb/s/zPLj0uLX2Y5I+5sG13ifpdfnHeyoNstY+YYz5oaRXKbfiu8FaG42zSAF4zsuVcC96JAdx5DAABCVxoddae3Pxc2PMuyV9IZTJlGGMuUDS7fmn85IO1njJI8qF3sskvVe5TXUAUshNr10veiQHceQwAAQlcaE3Bt4qqfBfl59aa6dqjD9c9PgDIvQCqVNPr10veiR7dbiIlwdjAEC9CL3BK16JPu5gfPGYK40x11hrf+LxnABEVCO9dhvtDNFocPbyYAwAaBSh1wVjzPnK9dP9Q+XajL1Y0nOSnpT0z5LusdZOVHl9i6Qbir70Cwdve6Lk+VZJhF4gJRppGeZFZ4h6g7NXB2MAgFcIvc71SHpY0hclvc1a+xtJMsZcJ+kfJA1K+ktjzIestZ+tcI2XSypeCjnt4H3PKFf725x/fl0dcwcQQ160DPOix3I9wZn+vgCihtDr3Fsl3WStPVD8RWvtD4wxNyh3Sts1ku4xxlxYoa/uFSXPax5VbK19wRgzpdzJb+WuASChvGwZ1khnCLfBmf6+AKKI0Ovc/aWBt8BaO2uM+YhyRwdL0ieMMd+11v64ZOglJc+dth8rHrfGyQtGRkZkjHF4+aWstXW/FoA3otYyzGlwpr8vkD6NZI729nYtW+Z/JG2uPQTKhc7v1hhzv6Rf5R8vk/TXZcZcWPL8OYfvXzyu9BoAEiquLcOiFtYBQGKl16m/UO40toqstfP5gyS25r/0ZmNMq7U2WzSs9P/ZnR6aUTzO0X8d+vr6tGdPxXMvAMSAVy3DghbXsA6gfo3cId6+fbuOHTvm4WzKazj0GmNWKLfJKww/DuKEMmvtvzoc+rOix+dLulbS94u+Nlky/nyH1y3eIVJ6DQA+CrPHrBe9dsMQ17AOINm8WOm9VtJ3PLhOPS7T0pZeYZoueV76/+QnSp47LVVoK3r8czcTAlCfqPSYbbTXbhjiGtYBJJsXofdRSZs8uE49ng7iTYwx7ZLeJGnGWjtcZeh8yfPS/yKWrhivdvDeF2pxmzN69AI+i1KPWS967YYhjmEdQLI1HHrzx+j+wIO5RJIx5mWSHpC0Pv/8bmvtBysML125/VXxE2vt08YYK6mwxdHJPb3ukucPOXgNgAZErcdso712wxDXsA4gudjIVtsnlA+8ebdKqhR6S9uJPV5mzNcl/ef848sdvH/xmOckfcvBawDUKco9ZhvptRuGOIZ1AMlF6K3tNSXPR6uMvaro8Y+tteXKLz6vc6H3KmNMk7W22naPDUWP782vrAPwCT1mvRe3sA4gmejTW1vpJrJbyw0yxqyS1Fv0pU+XG2etPSJpX/7pxVoclMu5Pv/nvKS/qzVZAI2hxywAJBMrvbUdkfRq5eqWb7TWVir0e5vOfYh4SNLuKtf8sKTNym10u0nS0XKDjDEXKLeBTpL+0VrLJjbAZ/SYRZht6gD4h5Xe2u7J//lUpcCbD6d/nn96WtI7rbUvVLqgtfa4ztUFv88YU+nDx9slvUi54P0htxMH4F5/T5eam5yNpcdsspyantWOPYf1ujse1M5vZPTJ+49r5zcyeu0dD2rHnsM6NT0b9hQBNCANobez+IkxprPSwAq+IOmrkrYYY5aUIuQD6+eU23D2mKTrrLXV6n4lSdbaz0n6qKQrJd1ljFn0uzDGbJB0l6RjkjZba3/jct4A6lDoMesEPWaTo9Cmbv/I+JKa7kKbum1Dh3R62vfzkAD4JHHlDcaYr0i6Lv+0SdK6kiHWGPProud3WWvvqnQ9a+2CMeaPJf2ZpG8ZY4aV68owm7/2O5QL1n8r6W+qlD+Uu/bHjTGPS7pT0o+MMV+T9IykqyW9U9KXJH3YWlt66AWQSkHddqbHbPpErU0dAO8lLvRaa2/24ZovSPqkMeZTkl6v3KruaklnJP2ppEestXXd97LW7jPGfEu5LhHXSGqX9LCkndZaZ32TgIQL+nQ0L3rMUhcaH1FuUwfAO4kLvX6y1j4v6cH8P15e9wXlNr9x8ARSxUkwDOt0tHp7zEbl+GI4R5s6IB0IvQAC5yYYhn3b2U2P2SgdXwznaFMHpEMaNrIBiBA3G4bque0cJjcBHdFBmzogHQi9AALlJhjWc9s5LHEL6DiHNnVAOhB6AQTGbTB86hnHzVAkhXvbOU4BPe5GJ7PaffCE7n7gSe0+eKLhDxC0qQPSgZpeAIFxGwxHXYaZMG87UxfqPz83CdKmDkg+VnoBBMZtMFzb0Rqb287UhfrL78MjCm3qtvZ1L/nfXHOTtLWvmw2IQMyx0gsgMG6D4YtXtWlLb7f2j9QuiQj7tnN/T5cG92UcrWR7FdDT1As4iC4e9bapAxAPhF4AgaknGN7U1x2L286FutAgAnraegEHfXiEmzZ1AOKD8gYAgalnw1CcbjsPDvRq/aq2qmMaDeh+3+aPIjYJAvACK70AAlXPhqEwbjvXUzrgxfHFtbi9zZ+EEgg2CQLwAqEXQKAaCYZB3HZutHTAz4Du5jb/t0fG9e4v/IseOn469iUQbBIE4AVCL4DARXXDkJfHCPsR0N3c5l+Q9N0nTi/5ehyPQw5jkyCA5CH0AvCc01vqUdswFESHgEa4vc1fTZj/Hm4FuUkQQHIRegF4Js5dBYLuEFAPt7f5awnr36MeHB4BoFF0bwDgibh3FYhDh4D+ni7Hh3U4EadOB3Hq4gEgmljpBeCJqJcG1BKHDgFubvM7FadOB1GtBQcQD4ReAA2LQ2lALXHpEODkNr8bcex0ELVacADxQHkDgIbFoTSgmlPTs/rOsVOOx4fZIaDWbf43XrHacQkEnQ4ApAkrvQAaFofSgEqqtSmrJOwOAbVu8+/Yc5hOBwBQgtALoGFxKQ0ox0ktcrEodQiodJufTgcAsBTlDQAa5qarQJRuqbupRZakN165OhYdAuh0AABLsdILoGFxPTzATS2yJN1wVWdsgiKdDgBgMUIvAE/E8ZZ6nGuRnaLTAQDkUN4AwBNxvKXe5PKghyjVIgMA3GGlF4Bn4nRL/dT0rL70Lycdjw+jFnl0MqvhIxOays5pRWuLNvdE7+cIAHFB6AXguTjcUt+5N6OxyVnH44OsRT41PaudezM6kFl8pPPgvoy29HZrcKBXne3LA5kLACQFoRdA6rjt2rBm5fLAapGr9Q2eX5D2j4wrMzYVuVKRSlitBhAVhF4gYQgZtbnt2vDH170ksIDppG/wyTMzun3viIa2bwxkTvVgtRpA1BB6gYSIe8gIMqy77dow7yYhN+CxpyYdtX2TpAOZcY1NZiP5gSZpq9UAkoHQCyRAnENGGGE9iifInZqe1a2f/xfH4+cXpOGjE5GsnU7KajWAZKFlGZAAbkJGlBTC+v6R8SXlBoWwvm3okE5PP+vp+0bxBLmdezN6xmUf4Cj2DXZTL11YrQaAIBB6gZiLc8gIK6wXTpBzIoiuDW431hVEsW+wm3rpwmo1AASB0AvEXFxDRthhfXCgV+tXtVUdE9QJcm431knh9A12Ig2n3AGIJ0IvEHNxDRlhh/UonSDn9ncoBds32I0o1ksDgMRGNiD24hoyohDWo3KCnNvf4cVtLYH1DXarv6dLg/syjj7QRHW1GkAysdILxFwUN2U5EaWwXjhB7rYbL9ctmy4NfAXVze9Qkr74n347cl04CqJWLw0ABYReIObiGjLiGtb94OZ3uLWvW31rV/o8o8ZEqV4aAAoIvUACxDFkxDWs+yWOv8NKolQvDQAF1PQCCVAIGeUOeWhuyoXGXQN9kQsZgwO9yoxNVW1bFpeg16i4/g4riUq9NAAUEHqBhIhjyEha0GtUHH+HtRTqpQEgbIReIGHiFjKSGPQaFbffIQDEAaEXQCQQ9NwZncxq+MiEprJzWtHaos096f2QAABOEHoBIEZOTc+WLQcZ3JfRlt5uDQ70qrN9eXgTBICIIvQCQEycmp7VtqFDZTf+zS9I+0fGlRmbojMCAJRB6AWQGEm/5b9zb6ZqpwtJOnlmRrfvHdHQ9o0BzQoA4oHQCyD20nDLf3QyqwOZcUdjD2TGNTaZTVTgB4BGEXoBLBK31dK03PIfPjKxKNBXM78gDR+dYGMgABQh9AKQFN/V0ijf8vfyA8RUds7V+LMz7sYDQNIRegHEdrU0qrf8/fgAsaK1xdX4lW3uxgNA0jWHPQEA4XOzWhol9dzy91vhA8T+kfElcyt8gNg2dEinp591dd3+ni41Nzkb29wk9W/ocnV9AEg6VnqBlIvqaqkTYd7yr1S64Fe5xdqOVm3p7db+kdq/qy293ZH5HXklbrXmAKKH0AukXJw3SIVxy79a6cLrL1+th5487eg69XyAGBzoVWZsqmqoXr+qTbsG+hxfM+riWmsOIHoobwBSLs4bpIK+5V+rdOG7x0/7Wm7R2b5c9+7YpK193Uv+vZubpK193ZGru26EX6UiANKJlV4g5eK8QSroW/5OShfcqOcDRGf7cg1t36ixyayGj07o7MycVra1qH9D8m73R7kzB4D4IfQCKdff06XBfRlHK5RR3CAV1C1/N7XPTjXyAWJNR2tkykz8EOdacwDRRHkDkHKF1VInorhByq9b/qOTWe0+eEJ3P/Ckdh88ofsOP+W4dMGJKH6AiJIoduYAEG+s9AKI/QYpL2/5V9o45bUofoDwUqPdFuJcaw4gmgi9AP59tbRc2GtuygW0XQN9kd8g1egt/2qHdHgpyh8gGuVVt4U415oDiCZCLwBJ6dogVYlXG9WaJL3hytV6qKSbQ5w+QNTDy5P94l5rDiB6CL0AFkn6BqlKvNyodlNfdyo/QHjZbSHth3EA8B6hFwDkbuNULX/6ppdJStcHCD+6LcS91hxAtNC9AQDkfuNUNT8++Yxn14oLP7otpO0wDgD+YqUXAOR+41Q1aewk4Fe3BWrNAXiF0AsAcrdxqpY0dhLwu9tCmkpFAPiD8gYAkLtDOqpJayeB/p6uJSUIlaT1ZwQgXIReAMgbHOjV+lVtDV0jrZ0E4n6yH4DkI/QCQF7xximHi5aLpL2TgJMPDWn/GQEID6EXAIoUNk59/yM36L9suUKXd15UMwDTSSCHbgsAooyNbABQxpqOVr3/TZfr/W+6fEnngKvXrtTjo2fpJFAG3RYARBWhFwBqKNc54BXrLw5nMjFBtwUAUUN5AwAAABKPlV4AiTU6mdXwkQlNZee0orVFm3u4xQ4AaUXoBZA4p6ZntXNvRgcy44sOmxjcl9GW3m4NDvSqs315eBMEAASO0AsgUU5Nz2rb0CGdPDOz5HvzC9L+kXFlxqboIgAAKUNNL4BE2bk3UzbwFjt5Zka37x0JaEYAgCgg9AJIjNHJrA5kxh2NPZAZ19hk1ucZAQCigtALIDGGj0wsquGtZn5BGj464e+EAACRQegFkBhT2TlX48/OuBsPAIgvQi+AxFjR2uJq/Mo2d+MBAPFF6AWQGP09XWpucja2uUnq39Dl74QAAJFB6AWQGGs7WrWlt9vR2C293RxUAQApQp9eIMWSeGLZ4ECvMmNTVduWrV/Vpl0DfQHOCgAQNkIvkEJJPrGss3257t2xqey/X3NTboV310AfB1MAQMoQeoGUScOJZZ3tyzW0faPGJrMaPjqhszNzWtnWov4N8V/JBgDUh9ALpIybE8uGtm8MaFb+WNPRqls2XRr2NAAAEcBGNiBFOLEMAJBWhF4gRTixDACQVoReIEU4sQwAkFaEXiBFOLEMAJBWhF4gRTixDACQVoReIEU4sQwAkFaEXiBlBgd6tX5VW9UxnFgGAEgaQi+QMoUTy7b2dS8pdWhukrb2dcf6YAoAAMrhcAoghTixDACQNoReIMU4sQwAkBaUNwAAACDxCL0AAABIvESWNxhjLpD0dknbJHVKepmkeUljkh6S9P9aax+p47rnSdou6RZJvZJekDQj6duSPmmtPRGFawIAAGCxxK30GmNeKel+Sa2S3mmtvc5a+1uS3iTpuMdLCQAAGrRJREFUJ5Juk/SwMWafMWadi+uuk/R9SR+X9GlJ66y1ayW9UdIKSceMMe91OVfPrwkAAIClErXSa4zpkfR5SW+21p4s/p61NiPp3caYJyX9X5LeLOlqY8wma+3TNa57iaTvSVonaZO19tGi645KutUYc5GkzxpjLrLW3ulgrp5fEwAAAOUlbaX385J2lgbeEh9XbsVXkl4i6T5jTMWDWfPlB1+V9FJJ/19xOC3xV/k//94Y019tkn5cE4ii0cmsdh88obsfeFK7D57Q2GQ27CkBAFIqMaHXGPNqSb+tXGlDRdbaBUlfKPrSJkn/W5WXvE/S6/KP91S57hOSfqjcz/Sz+briIK8JRMap6Vnt2HNYr7vjQe38RkafvP+4dn4jo9fe8aB27DmsU9OzYU8RAJAyiQm9kq7J/7newdgflTz/3XKD8iHz9vzTeUkHa1y3sDnuMklla3H9uCYQJaemZ7Vt6JD2j4xrfmHx9+YXpP0j49o2dEinp58NZ4IAgFRKUuhdkf/z28aY240xG6qM/VXJ85dVGPdWSV35xz+11k7VmMPhoscfCPCaQGTs3JvRyTMzVcecPDOj2/eOBDQjAACSFXpt/s8XSxqU9Jgx5vcqjG0veb5QdpR0c9Hj4w7mUDzmSmPMNWXG+HFNIBJGJ7M6kBl3NPZAZpwaXwBAYJIUeg8o14e3YJmkP6kwdm3J86OlA4wxLZJuKPrSLxzM4UTJ861+XxOIkuEjE0tKGiqZX5CGj074OyEAAPISE3qttTOS/qOk4nKBwxWGF3dCWFD5zWQvV67Xb8FpB9M4o1ydbsF1AVwTiIyp7Jyr8Wdn3I0HAKBeierTa639Z2PMZcp1Rpi01n6vdEx+I1lxicFnrLU/KR0n6YqS5884eP8XjDFTkjoqXMOPa5Y1MjIiY4yToZXet+7XIr1WtLa4Gr+yzd14AEA0NZI52tvbtWyZ/5E0UaFXkqy1ZyTtrTJkUNLq/OP7JP15hXGXlDx3utW8eNyaAK4JREZ/T5cG92UclTg0N0n9G7pqDwQAwAOJC73VGGO2SPpw/umXJd1irX2+wvALS54/5/BtiseVXsOPa5bV19enPXsqtgAGfLG2o1Vberu1f6T2ZrYtvd1a09FacxwAIPoauUO8fft2HTt2zMPZlJeYmt5ajDGbJX09//S/SvrjKoFXkkqLDc9z+FbF40qv4cc1gUgZHOjV+lVtVcesX9WmXQN9Ac0IAAAPVnqNMSsk9Xgwl3r82Fpbs0TAGLNd0j9KmpR0s7X2Ow6uPVny/HyHcyo+Na30Gn5cE4iUzvblunfHJu3cm9GBzOIDKpqbciu8uwb6tLqdAwYBAMHxorzhWklOQqQfLtPSll6LGGNuV66O9yHlAu/TDq9del1HZQWSipe4fh7ANYHI6WxfrqHtGzU2mdXw0QmdnZnTyrYW9W/ooqQBABAKL0Lvo5I2eXCdelQMsPmeuJ+TdIukT0j6K2vtCy6u/a8lz1eXHbX4PS/U4pZkpV0h/LgmEFlrOlp1y6ZLw55GRaOTWQ0fmdBUdk4rWlu0uYdQDgBJ1XDozR+j+wMP5uIZY0ybcp0ZXi3p962136wy9nWSnrfWHir+urX2aWOMlVToweFkm3l3yfOH/L4mAPdOTc+WLb8Y3JfRlt5uDQ70qrN9eXgTBAB4LnEb2fIrvP9L0kslvaJa4M17u6RLK3zv60WPL3fw9sVjnpP0rYCuCcChU9Oz2jZ0SPtHxpe0VptfkPaPjGvb0CGdnnbaURBRNzqZ1e6DJ3T3A09q98ETHH8NpFQSW5Z9XLlTy/qstU85GH+tpH0Vvvd5Sf85//gqY0yTtbZaB9INRY/vza+CB3FNAA7t3JvRyTMzVcecPDOj2/eOaGj7xoBm5R6lGbWxog+gWKJCrzFmg3KHTdzpJPAaY7qVK4EoezKatfaIMWafpN+TdLGkqyQdrXLJ6/N/zkv6u6CuCSSdVwFvdDKrA5naPYQl6UBmXGOT2cgFSYKcM4UV/XIfcAor+pmxKd2343o6iQApkajQK+lPlCvZOK/aCqoxZplyZQOfUa4H7pkq1/ywpM2Slku6SRUCav544zfln/5jhaON/bwmkDheB7zhIxOOTouTcsFo+OhEpDbiEeScS8qKPgDvJK2m93fyf/6ZpOeNMWX/Ua429ojOBcqyK72SZK09LumD+afvywfmct4u6UX5636o2iT9uCaQNH7U3k5l3Z3tcnYmWmfBuAlyaVbPij6A5EtM6DXGLFeuVKCgWblV3HL/NBWNW1CNAx+stZ+T9FFJV0q6yxiz6OeWL6u4S9IxSZuttb+pNV8/rgkkiR8Bb0Vri6s5rGxzN74Wpxuqyo0jyDlXz4o+gORLUnlDaWsvp6ac9O+11n7cGPO4pDsl/cgY8zXlVoivlvROSV+S9GFr7bTTN/bjmkAS+FV729/TpcF9GUeBqLlJ6t/gpKtgbU7LNKqNu6q7PdalGUGK+4o+AH8kJvRaa09o8QquH++xzxjzLUmvkXSNpHZJD0vaaa119l/oAK4JxJ1ftbdrO1q1pbdb+0dq/9Xa0tvtySY2p3W497zrWr3vi49WHHfkaXeffdMc5MJe0QcQTYkJvUHJrwo/JA8PifDjmkCc+blSNzjQq8zYVNXSifWr2rRroM/VHCpxWqbx7i/8UBNT3vUGTnOQC2tFH0C0JaamF0By+LlS19m+XPfu2KStfd1qLrk31Nwkbe3r9qz7gZsyDS8Db9qDXGFF3wmvVvQBRB8rvQAix++Vus725RravlFjk1kNH53Q2Zk5rWxrUf8Gbw94cFOm4SWCXPAr+gCij9ALIHKCqr1d09Hq62Yvt2UaXiDI5RRW9MttDGxuyv3vZtdAX+r7GQNpQugFEElJWKlzW6bhVM8lK3RsfMpRkEvzccVBregDiAdCL4BISsJKnZsyDaeam6T/fusrJeW6VvzyzIxGJ7Na29GqdRe3ae6FeUkcV1zM7xV9APFA6AUQWXFfqXNTptG14gJHm9kK5Rynpmd1yP66bKh94xWr9cTEv2m0zAEVHFcMIK0IvQAiL84rdU7LNKr16S0et2ugr2bv3wefOF1zXoXT7Ia2b3T2LwIAMUfLMgDwkdMWaRsuWem4lZqT3r9OpP24YgDpwkovAPjMaZmGk3Fuev/WkvbjigGkC6EXAALitEyj2jive/+m+bhiAOlCeQMAxIjXvX/TfFwxgHQh9AJAjHjZ+zftxxUDSBdCLwDESH9P15KNbvXiuGIAaULoBYAYKfT+bVTUT7MDAK8RegEgZgYHerV+VVtdry1tfwYAaUH3BgCImULv34/c95gePFb7IIqL21r09le+WGsubo3NaXYA4DVWegEgpn566jeOxj0zM6dfnJnRLZsuJfACSC1CLwDEkNtT2Th9DUDaEXoRacYYGWPCngZCwO++snpOZSucvhYH/O7Tid97eo2M/P/t3XmwXGWZx/HvzU1IbjDJJUAWkwDxQbIBA5hBBJRFmAhGQWWTAplxnLFEoMShcJQxGmekCqaUWBMGEdGiXIApNWKYAYQBhsVQA0RmyEKEBwIYCCAhC2S9uT1/vKfnntvp5fRyuvue/n2qUn3O6bffPqm+z9vPefs977uCvr6+o9J+HyW9IiJDTK2rsmn1NRHpZEp6RUSGmFpXZdPqayLSyZT0iogMMbWsyqbV10Sk02nKMhGRNrZu4zbuW/Uam7ftYmzPCE6dPZFTZk9k4dKVVQ1x0OprItLplPSKiLSh17ds5xt3rOSelesHJbcLl65k3pxJnHjI/ty/pvIcvaDV10REQEmviEjbeX3Lds66YVnRKcn6c3DXivVM6e1hSm8P6ypMQ3byzAlc86nDtfqaiHQ8Jb0iIm0myRy86zZu4+SZ+3P41HF79AZ3AbMmj+Gas/6Mw6aMS/dkGyQ+jENEJA1KekVE2kg1c/A+uOYNHvnKyXx9/mzuW/0am7buYtzoEUNqqeFiwzi6o+e+8NMnWXjGHCaMGdWy8xOR7NDsDSIibaSaOXjzC07kgFz0mlwN8/e2Sn4Yx10r1hf9P9+1Yj1n3bCMN7bsaP7JiUjmqKdXRKSNVPvz/k+Wvcg3f7Oy6M1u7d5LmmQYx0sbtrLgjhXccMH7mnRWIpJV6ukVEWkj1c7B++zrb+/RS5q/2a2de0mrGcZxz8r1vFLhhj0RkUqU9IqItJFTZk9kWFdj6sr3krajWoZxiIjUQ0mviEgbmdLbw7w5kxpWX7v2klY7jGPTVs3qICL1UdIrItJmFp4xhwPGj25IXe3aS1rtMI5xo6tfellEJK4rN5Ru9ZVE5s6d++bIkSPHT58+vdWnUrcVK8JPs4ceqtWkOk2nf/Z9u/t5ZdP2oj2iI4d3s6Nvd+K6JowdxYQ2W5xi1+5+1qzfUvzJt14Oj/tM+/9DMyaNYUT3nv00u3b3s3lbH7tzObq7uhjbM7xoOWl/nR7znWz16tXkcjmWL1/eoMFdxSnpzaC5c+e+AIwF1rb4VESkTrmu7hG5ET29ua6u7q5cbnfXrm0bcyN6evtHjplW+dXBsB1bXh628+1kaxY30e6e8e/JDR/ZW6lcV9+Ojd3bNjwfP5brGja8f1TvAcVe39W3Y+Ow7Rtf6sr19zXyfEUkNQcBm5944olUe+uU9IqIiIhI5uk3IBERERHJPCW9IiIiIpJ5SnpFREREJPOU9IqIiIhI5g1v9QmIiNTLzGYD49x9WYVyXcAngAfc/a2mnJyIJKZYljQp6RWRtmBm5wC3u3st8zQuBk4ws7uB24FlwMtADugFDgFOAD4NbAPureK8uoELgM8Ac4DdwFbgbuA77r62hvMVkeIUy5IaJb3SMmY2EjgXOAuYABwM9AOvAA8B/+buj9RQb8MbNjWWTXFhHa8dFv07PfpXypPAae5eYlWEwcxsKvALYBpwCbDU3fvMbApwNfCMmV3q7jfVce4dr84Lnng9ZwKfA95HiNE+4H7gOnd/ul3qlLIUyx3GzE4AvggcR7i4AXgEWFxLDlCOxvRKS5jZXMIVeg9wvrsf4+77AScBTwGXAg+b2dKosUpa71TgUUIjthiY6u5TgBMJC3Y8Y2Z/U+W5NrxOGczMPgvMT/EtdgHXAh9w90SLNJjZZOC/gCOBj7n7EnfvA3D3de5+EfDvwA/M7PKUzrtT1HPBg5n1mtmdwC3AEuBAd58KHAW8CTxlZt9odZ3SEIrljDCzkWb2I0Iu8DhwSBRjM4DlwP1m9n0za9ga5FqcQpouGrN1GzDf3V8qUeYq4J+i3RcJDdyrFeqdTLg6nBqVX16kzC+BTwJfdvfrEpxrw+uUAdEFxReArwDdALX09pnZg8AzhJ9BTwH2I/QYrCN8fj8p9bdWor5u4AHgg8DP3P2CEuVmRO/bD8xz9/uqPfdOF13w3Aw1f/bvIvS8/jkhobmzSJnvApcD33P3L7WiTklGsdwZzGw44WJyPnCJu19fpMxlwPeAXwHnuHvytddLUNIrTWdmjwHXuPuSMmW6CFd6R0SHlgHHuXvRP9g0GjY1lukws9uAY4C9gNHRY0/++TqS3gfd/ZsNOseLgXwjfJq7312m7H8TkqMXgFnuvqMR55B1Dbzg+TlhfOej7n58iTK9hGFTPcBn3f3Hza5TklEsdwYzuxr4KvASMN3d+4uU6QbWEjqdvuXudf+youEN0lRmdjTwfircfBAlt/EvkQ8Af1HmJZ8nJKcAPy1T7xrCzyjDCD9ljWxynR3P3c9z94Pc/d3u3gtc3Opzios+vwXRbj/wuwovyY85mw5omEsZZnabma01s1eAFYSe0u466judkJxC+RjdCOR7a68zs/2bWae0hmK5PZnZYcCV0e6txRJegKhn99Zo96roV+K6KOmVZsv33B6QoOwTBfsfLVYojYZNjWVHOxuYGG0/5+6bK5R/MrZ9STqnlA0pXPB8O7b9cIWy+RgdR+hdbmad0hqK5fb0LQYudpPGWHf0uroo6ZVmGxs93m1mC8xsVpmyfyrYP7hEuTQaNjWWneu82PYfEpSPl5lhZkeULCkNY2bHM3ARvR1YVeEl8Rj9azPrKSyQRp3SUorlNmNmBwIfjx3a4z6ZAvEYO7OaG9uL0ZRl0mwePU4DFgJfN7NPuvvSImXHFOyXGoBed8Pm7k81oU5JmZntRZhW7hOEv7FpwE7gWeC3wI3u/lqZ148ATo4dejHB264t2D+NMAOJpCseo8+VGu8fE4/RXmAe8Osm1Ck1UCxn1tkMdLhuqXSDuruvM7N3gL0Jvb1nAYtqfXP19Eqz3UO4+SNvOGEOzGKmFOyvLizQwIYt1TqlKWYTfiobRbjT93B33wc4g9DWLQSeN7O/LVPH4cRuqgOSTIm0gTAEJu+Yqs5aavWR2HbFGI2mt9oaO1QsRtOoU6qnWM6ueIwlnYUjHot1xZiSXmkqd98K/BUQHy7wZInip8S2cxS/qSSNhk2N5dB0NrDA3Re7+zv5g+7+GOEi5inCbBE3lpmL85CC/YrLm0Y3W8T/ngvrkAYzs30Bix1KNF8rg4dMDYrRNOqUmimWs+vo2HYtMfb+et5cSa80nbv/lnDD15nAie6+x+D06Eay+E+N15cYLpBGw6bGcmi6193vKfaEu28H/j526FozO7JI0ckF+0mnLIqXe3fC10jtqo7RyIbY9nujqRHTrFNqo1jOoGje+/iwxVpibJyZTSxZsgKN6ZWWcPcNwB1liiwE8lMA/RL4colyaTRsaiyHnh2EFZfKuZfQY7Afoe37B+BTBWX2LtjfmfD94+UK65DGa0SM9hDG4ea/eNOoU6qnWM6uRsQYhO/XkuO5y1FPr7QdM5sHXBHt3gqc5+67ShRPo2FTYzn0XAl8v1yBaC7Ix2OH5he5277w7yzpHLLxcqX+VqVxGhGjhfWkUadUT7GcXS2PMfX0Slsxs1MZuPv5q4SV28rdQZ1Gw5apxtLMxhJuDGmF3zdjVSN3/5+ERZ+Pbe8FHAU8Gju2saD8XgnrjS9IUliHNF4jYrSwnjTqbHvt1j4oljOt5TGmpLdDtFvDVoyZXQDcTGhoznP3BxLUnUbDlrXG8ijCcsqtMJ09Z7ZopS0F+4Vjw9YW7CftURgd236hmhNK01CI+xo1IkYBNqVc51AwVNuHjorljGhUjNX8/aqkt3O0dcNmZgsI43gfIiS8ZefuiymstxENWxp1ttJywjLOrZD0c6yLmY0BTgK2uvt9ZYoWLnc5qmC/sJep4tKyZrY3g2f7aKd5Pds67utQWG8tMfpqdFNUmnUOBW3VPiiWM+1FwkxM+Zs9a4mxHMmmES1KSW/naKuGLS+aE/cmwiTk1wJfi2ZCSCqNhi1TjWW0mtxjrT6PtJjZwcB/Ei1tbWb/4u6XlShe2MgOWvXP3V81M2dg6qokdwlPKth/KMFrmqUt474BniPMj5v/MqwYo0XKFcZoGnW2vXZqHxTL2ebu28zsWQZmSqklxp6p58JSSW+HaKeGLc/MRhNmZjga+Li731mm7AeBPndfFj+eRsOmxnLIuZboSzJyEVDqi7JwVo2ni5T5NfB30fZ7E7x/vMxO4D8SvKYp2jHuG8Hd+8xsGfDh6FDFGI0usMfHDhXGfcPrlKoplrPvYQaS3qRTj8W/X+uKMc3eIC0RfVn8CngPcGS5hDdyLnBQiefiy342qmFLo05Jx3EF++vKlJ0Z2/59iWE0P4qXTzDv6qzY9i+iRFPSV22MGoO/85Y0qU5JTrGcffEYm2pmhcNSBolm5YivzlpXjCnplVa5mrB60YfdPclShEcxeILquDQaNjWWQ0fhWOqLihUys/HAnNihxcXKufsqYGm0uw+Dv1yLOTZ67Af+udLJSsPcCmyLtiebWW+F8vEY/Z27r2lSnZKcYjn77mHgYmYYMKNC+ZkMjAF+CSg3zrsiJb3SdGY2i7DYxA/d/Y8Jyk8iDIEoOuF7Gg2bGsshZVX0+BhwqLs/XqLcOQy0eQ8Bt5Sp8wogP27sI6UKRSsHnhTt3lxi1UBJgbu/CfwwdqjS2OVjY9vXNKtOqYpiOeOiOfe/Gzt0bKmyRZ7/TpX3/OxBSa+0wucIf3vd5XpQzWx4lCD/nDBPX6meXkinYVNjOTTcGD3+0d23FisQfUb5Vf3eAM4v13i6+x8YGEv4eTMrdf/DucC+hC/ry6s9canbQgZuYCoZowXP3+vuv2lynZKMYrkzLAaeibaTxtgq4IZ631hJr7TC6dHjl4A+Myv6jzA2dhUDCWXJpT3TaNjUWDbNhPiOmU0oVbCEHwO3A/PMbI8e+ehzu4kwRvN/gWPcvdxYQQDc/SbgKsLPb4vMbFB7GV2QLSI03qe6+ztVnrfUKeqZvRDYDVwYTXe1BzP7EHAosJ4SP5mnWackpljuAO6+E/g0YbaUj5rZQcXKmdl04DTgHcqvzJpYVy5XbrErkcaKBq1vq1hwTzlgRKWfNszsa8C3geuBy6LlKvPPzSKs2PMaYSzxKwnPueF1djIzu40wnhvCWK2pDL4Afxt4M7a/yN0XVaizm3AR9UXCmK+nCb30UwmN6wRCL8E/lupBKlP3x4DrgM2EmyjeAg4Dzif8CnGFuxdOlC8JmNmVDB4WMNHdX6+hns8QFra5Ezg3+lLNPzeZEKOjCAnNylbVKZUpljtHtALrEmAFIY62xJ4bS/j8ZxJmd3qwEe+ppFeaKrqiq2Xhhk3uXummkvx7NLxhU2M5NEQ9QR8i9ATtS+ghWAk8Us/cjtEX8XHAEYT5QV8G7nP39XWfdAdJ44InVvexhAvTdxFuSFtPmBrpL4F7gUur/bzSqFOSUSx3BjObDfwroRf+FsLCEwcSfj1ZCVwc/eraEEp6JZPSaNjUWIq0t+gegbmEG1/HERbIeMDd17ZTnSIymJnNAY4nXOD8CXjY3Vc3+n2U9IqIiIhI5ulGNhERERHJPCW9IiIiIpJ5SnpFREREJPOU9IqIiIhI5inpFREREZHMU9IrIiIiIpmnpFdEREREMk9Jr4iIiIhknpJeEREREck8Jb0iIiIiknlKekVEREQk85T0ioiIiEjmKekVERERkcxT0isiIiIimaekV0REREQyT0mviIiIiGSekl4RERERyTwlvSIiIiKSeUp6RURERCTzlPSKiIiISOb9HzxuzXGSz4q2AAAAAElFTkSuQmCC\n"
},
"metadata": {
"image/png": {
"width": 350,
"height": 313
},
"needs_background": "light"
}
}
]
},
{
"metadata": {
"trusted": true
},
"cell_type": "code",
"source": "",
"execution_count": null,
"outputs": []
}
],
"metadata": {
"kernelspec": {
"name": "conda-root-py",
"display_name": "Python [conda env:root] *",
"language": "python"
},
"language_info": {
"name": "python",
"version": "3.7.3",
"mimetype": "text/x-python",
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"pygments_lexer": "ipython3",
"nbconvert_exporter": "python",
"file_extension": ".py"
},
"gist": {
"id": "",
"data": {
"description": "hq/notebooks/Simple-Linear-Nonlinear-Inference.ipynb",
"public": true
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@davidwhogg
Copy link

This is sweet! Does it verify our claim?

Oh and you can make things more structured by making the M quadratic and generating the data from a line plus noise.

@davidwhogg
Copy link

I guess it doesn't verify...

@adrn
Copy link
Author

adrn commented Oct 5, 2019

@davidwhogg Hm, shouldn't the plot at the bottom be a straight line if the claim is right?? And sure, I'll make the data more realistic!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment