Skip to content

Instantly share code, notes, and snippets.

@tomkreker
Created April 1, 2020 10:40
Show Gist options
  • Select an option

  • Save tomkreker/650c63a7fb9ba7adc9b551fb597dde5d to your computer and use it in GitHub Desktop.

Select an option

Save tomkreker/650c63a7fb9ba7adc9b551fb597dde5d to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.stats as sts\n",
"from scipy import signal\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def target(x):\n",
" '''Target distribution we want to sample from, copied from PCW'''\n",
" Z = 24.44321494051954\n",
" if abs(x) > 7:\n",
" return 0\n",
" elif abs(x) > 3:\n",
" return 3 * (1 - (x / 7) ** 2) ** 0.5 / Z\n",
" elif abs(x) > 1:\n",
" return (\n",
" (3 - abs(x)) / 2 -\n",
" 3/7 * 10**0.5 * ((3 - x**2 + 2*abs(x))**0.5 - 2)\n",
" ) / Z\n",
" elif abs(x) > 0.75:\n",
" return (9 - 8 * abs(x)) / Z\n",
" elif abs(x) > 0.5:\n",
" return (3 * abs(x) + 0.75) / Z\n",
" else:\n",
" return 2.25 / Z\n",
"\n",
"def proposal(x_null):\n",
" '''The proposal distribution, Gaussian centered on the current X with SD=2'''\n",
" return(sts.norm(loc=x_null,scale=2))\n",
"\n",
"def Metropolis():\n",
" '''The Metropolis algorithm '''\n",
" x_current = sts.norm(loc=0,scale=3).rvs() #start X0 somewhere\n",
" samples = []\n",
" num_steps = 100000\n",
" gap=100\n",
" rate = 0\n",
" for t in range(num_steps):\n",
" x_new = proposal(x_current).rvs() # sample from the proposal\n",
" r = target(x_new)/target(x_current) # calculate r \n",
" if (r > np.random.random()): # check if R is greater than a random number\n",
" new_sample = x_new # set new sample value\n",
" x_current = x_new # change current X status\n",
" rate+=1 # count acceptance rate\n",
" else:\n",
" new_sample = x_current # set the new sample to be the current x without changing it\n",
" \n",
" if t % gap == 0: # to take only the gap-th sample \n",
" samples.append(new_sample)\n",
" print('Acceptance rate: %.2f' %(rate/num_steps))\n",
" return(samples)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Acceptance rate: 0.71\n"
]
}
],
"source": [
"# run the algorithm to obtain samples\n",
"samples = Metropolis()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAswAAAD8CAYAAABjNPKeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xlc1NX+x/HXmWHYBGFQRFEUMbU0\n17ze8mp2q5uaZt2WW2YuaZamldrt51K2WJq2W9qilrdcSttui4ptVlpWahqKS5qiIi4o4ArDDHN+\nfwBeJZRBBg4z83k+Hj2Cme8Mb2jMN2fOorTWCCGEEEIIIUpnMR1ACCGEEEKI6kwKsxBCCCGEEOcg\nhVkIIYQQQohzkMIshBBCCCHEOUhhFkIIIYQQ4hykMAshhBBCCHEOUpiFEEIIIYQ4BynMQgghhBBC\nnIMUZiGEEEIIIc4hyHSAkmrXrq0TExNNxxBCCCGEEH5u7dq1h7TWsWVdV+0Kc2JiImvWrDEdQwgh\nhBBC+Dml1C5PrpMpGUIIIYQQQpyDFGYhhBBCCCHOwaPCrJTqrpTaqpTarpQaW8r9lyulflVKuZRS\nN592e1ul1CqlVKpSKkUpdas3wwshhBBCCFHZyizMSikrMAPoAbQA+iilWpS4bDcwEFhQ4vaTQH+t\ndUugO/CSUiq6oqGFEEIIIYSoKp4s+usIbNda7wBQSr0HXA9sKr5Aa51WdJ/79AdqrX8/7eMMpdRB\nIBbIqXByIYQQQgghqoAnUzLqA3tO+zy96LZyUUp1BIKBP8r7WCGEEEIIIUzxpDCrUm7T5fkiSql6\nwFzgTq21u5T771ZKrVFKrcnMzCzPUwshhBBCCFGpPCnM6UDCaZ83ADI8/QJKqZrAYuARrfVPpV2j\ntZ6pte6gte4QG1vm3tFCCFGtaV2uMQWfFSjfpxBCeFKYVwNNlVKNlVLBwG3Ap548edH1HwPvaK3f\nP/+YQgjhG1555RVatWqFw+EwHaVSvf322zRq1Mjvv08hhAAPCrPW2gWMAJYBm4FFWutUpdREpVRv\nAKXUX5RS6cAtwBtKqdSih/8LuBwYqJRaX/RP20r5ToQQohpYt24dqampvP3226ajVJr8/HwmTJjA\nnj17OHjwoOk4QghR6Tw6GltrvQRYUuK2R0/7eDWFUzVKPm4eMK+CGYUQwmdkZ2cDMHXqVAYNGkRQ\nkEf/m/Up8+fPZ8+ewrXg2dnZJCQklPEIUR0ljl1c4edIm9LTC0mEqP7kpD8hhPCi7OxswsLC2LFj\nB4sWLTIdx+sKCgqYMmUKYWFhwP9+QRBCCH8mhVkIIbwoOzubq6++mpYtWzJ58mTc7j9tDOTTPvzw\nQ37//XdGjhwJSGEWQgQGKcxCCOFF2dnZxMTEMG7cOFJTU/nss89MR/IarTWTJ0+mefPmDB48GJDC\nLIQIDFKYhRDCi7Kzs7Hb7dx6660kJSUxefJk05G8ZunSpfz222+MGzeOWrVqAVKYhRCBQQqzEEJ4\nidPp5Pjx49jtdoKCgrj77rv55Zdf8JcDmT7//HMiIyO5/fbbqVmzJkopKcxCiIAghVkIIbwkJycH\nALvdDkD79u0B2LBhg7FM3rRhwwZat26NzWbDYrEQHR0thVkIERCkMAshhJcUl8fiwtyqVSvAPwqz\n1poNGzac+p6g8PuUwiyECAT+t0GoEH7CG3ukguyTWpVKFua4uDhq167tF4U5PT2dI0eOSGEWQgQk\nGWEWQggvKVmYlVK0bt2alJQUk7G8ovh7aN269anbpDALIQKFFGYhhPCSkoUZCqdlpKam+vx+zMWj\n5BdffPGp26QwCyEChRRmIQKM1poVK1bQr18/4uPj2bx5s+lIfuNshfnkyZPs2LHDVCyv2LBhAwkJ\nCURHR5+6TQqz9w0ZMoTWrVvz8ssvy89WiGpE5jALv+GNOb/+PN/30KFDvPPOO8yaNYstW7YQEhKC\nw+EgJSWFiy66yHQ8v1BaYS6ewpCSksIFF1xgJJc3pKSknDEdA/5XmLXWKKUMJfMvK1asYOfOnTzw\nwAOMHTuWW2+9laFDh9KxY0f5GQthkIwwC+HHtNasXLmS22+/nfr16/Pggw9it9uZM2cOGzduBOTg\nCW/Kzs4mLCyMkJCQU7e1bNkSpZRPL/zLz89ny5YtZyz4g8LC7HQ6OXnypKFk/ic7O5uBAweydu1a\n+vXrx/vvv8+ll15Ku3bteP311zl27JjpiEIEJCnMQvght9PBsd++YN9/HqBLly4sWbKEe+65h5SU\nFH788UcGDhxI/fr1ASnM3lR8yt/pwsPDadKkiU8X5q1bt+JyuUotzCCvIW/RWp96DbVv35433niD\njIwMXnvtNQCGDRtGfHw899xzD+vXrzecVojAIoVZCD/iOnKA7OVvsffVgWQlvwzuAt544w327t3L\nyy+/fEbhKR4JLT5sQ1RcaYUZCqdl+HJhLs5e2pQMkMLsLbm5uTidzjNeQzVr1mTo0KGsW7eOn376\niZtvvpm5c+fSrl07rrjiCj7++GMKCgoMphYiMEhhFsIPOPZuJvPjyex9YwhHV/+X0IatiOvzNPUG\nTefuu++mRo0apT5OFm1519kKc6tWrdi2bZvPTl1ISUnBZrPRvHnzM26XwuxdxT/H0xdWFlNK8de/\n/pU5c+awd+9ennvuOdLS0rjxxhtp2rQpL774IkeOHKnqyEIEDCnMQvgo7S7g5NYf2T/33+yf9xB5\nu1Oo+dcbqT/0TWL/OZ7Qhq3KXCQkRxt717kKs9aaTZs2GUhVcRs2bODCCy/EZrOdcbsUZu8qbdFo\naex2Ow8++CDbt2/ngw8+oH79+owePZqEhARGjhzJrl27qiKuEAFFCrMQPsbtzOPYr4vJmDWUzP9O\npuBENvar76H+sDnYuw4kqGasx88lI8zeda7CDL57RHbJI7GLSWH2Lk8Lc7GgoCBuuukmVqxYwerV\nq+nduzczZsygSZMm9O/f/9TCXiFExcm2ckL4CLfjBMd+XczR1f/FnXuU4HrNqd11AOHNLkNZrOf1\nnHa7nQMHDng5aeAouZXhnv2ZfJyaw3clbtfuAlRQCA+8+glPbK1TlRHLreTWijk5OezZs+dP85dB\nCrO3Fa8n8LQwn65Dhw7MmzePp59+mhdeeIFZs2Yxd+5cevXqxZgxY+jcubO34woRUGSEWYhqruDk\nEbK/n0v6a4PI+f4dgus1Ja7vVOr2e44aF3Y+77IMMsLsTdpdgM7PxRIa8af7lMWKrXZDnJlpVR+s\ngopHxUsbYY6KikIpJa8hLynvCHNpEhISePHFF9m9ezcTJ07kp59+okuXLvztb38jOTkZrbW34goR\nUKQwC1FNuY4dJuvrWex9fRBHVy0irFEb6g54ibhbniC0QUuvHGIgc5i9x513HKDUwgxgi21Efqbv\nzS09V2G2WCxERUXJa8hLzrXor7xiYmKYMGECu3bt4pVXXiE9PZ0ePXpw2WWXsWTJEinOQpSTFGYh\nqpkDBw4watQo9r5xF8fWfkZ4s07ED36V2H+OJ6Sud0+Ks9vt5OTk4Ha7vfq8gaiswhwc2xj3yRwK\nTvhWudywYQPR0dE0aNCg1PvlXQrv8WZhLhYeHs6IESPYtm0bM2fO5MCBA/Ts2ZOOHTtycvvPUpyF\n8JDMYRaimjh06BDPPvss06dPx+FwUKPF34nqdBu26LqV9jXtdjtaa44ePerVv6QDkScjzAD5mbsI\nq3H+b7lXtZSUFFq1OvuOK1KYPVdyzntJWV/9igoOp8nDyee8ruQ8c08EBwczZMgQBg4cyNy5c5k0\naRKZa54kOK4JUV3uICypgxy9LcQ5yAizEIZlZ2fzyCOP0LhxY5599lluvPFGNm/eTO1rR1ZqWYb/\nzZWUw0sqruwR5kQAn5rHrLVm48aNpU7HKCaF2XvcjhNnff14i81mY9CgQWzZsoVa147E7ThB5gdP\ncGD+GPL2yK4aQpyNFGYhDMnNzeWZZ54hKSmJSZMmce2115KamsrcuXNp2rRplWQoHlWWwlNxZRVm\na41oLDWiyfehwrx7926OHj1a6g4ZxaQwe4877ziW0NIPGfI2m81GRKurib/rdWKuuRfXkf0cWDCW\nA4sew7F/e5VkEMKXSGEWoooVFBQwZ84cmjVrxpgxY+jUqRO//fYbCxcu5KKLLqrSLLItmPe4HSeA\nsxdmgODaiTh9aOHfuRb8FZPC7D2FhblyR5hLUtYgIttdS/zdM4m+4k7y9/3O/rdHkvnfKTiz9lZp\nFiGqM48Ks1Kqu1Jqq1Jqu1JqbCn3X66U+lUp5VJK3VzivgFKqW1F/wzwVnAhfI3WmsWLF9O2bVsG\nDRpEfHw8y5cvZ/HixeccwatMUpi9p3iE2XqOwmOLbYTz0C60u6CqYlVISkoKABdffPFZrykuzLJ4\nrOJMFOZiFlsoUX+9ifpDZxPV6TZyd6wh4817yfpqJgW5R41kEqI6KbMwK6WswAygB9AC6KOUalHi\nst3AQGBBicfGAI8BfwU6Ao8ppXxntYsQXrJ+/XquvPJKevXqRV5eHosWLeKnn37iiiuuMJpL5jB7\njzvvOCooGBUUfNZrgus0RrvyceXsr8Jk52/Dhg0kJiZSs2bNs15jt9vJz88nNze3CpP5J3feCSwh\nZgpzMUtIDaK73EH9u2cR0epqjv36ORlvDOHoLx+jXU6j2YQwyZMR5o7Adq31Dq11PvAecP3pF2it\n07TWKUDJvam6AV9qrbO01tnAl0B3L+QWwiccPHiQu+++m/bt27Nx40amT5/Opk2buOWWW6rFinQZ\nYfYeT0YHbUUL/3xlHvPZjsQ+nbyGvMftqLo5zGWxRtip1f0+6t35MsHxzcle/iYZbw7jxJaV8m6C\nCEieFOb6wJ7TPk8vus0TFXmsED4rPz+f5557jqZNmzJnzhxGjhzJtm3bGD58ODabzXS8UyIiIrBa\nrVJ2vMCdd7zM0UFbrQRQFpwH06omVAU4HA62bNkihbmK6AIn2ukwNiXjbIJjE4n710Tq3PIEKiiE\nQ59M4cCCseQf3Gk6mhBVypN9mEsbBvP010uPHquUuhu4G6Bhw4YePrUQ1Y/Wms8//5zRo0ezfft2\nevbsyfPPP0/z5s1NRyuVUkpO+/OSAg9GmC22EILs8eQfSquaUBWwZcsWCgoKypxfL4XZOzyZA1+s\nrP2cK0NY0iWEJrbleMqX5Hz/Dvv+8wCR7XuSnd2pQkd5C+ErPBlhTgcSTvu8AZDh4fN79Fit9Uyt\ndQetdYfY2FgPn1qI6mXnzp307t2b3r17Y7PZSE5O5vPPP6+2ZblY8Wl/omI83RIsOLaRT+zF7MkO\nGSCF2VvceWXvsmKasliJbNud+CFvENG2B8d+XUyzZs2YPXu2nBYq/J4nhXk10FQp1VgpFQzcBnzq\n4fMvA65RStmLFvtdU3SbEH7D4XAwefJkWrZsyfLly3nuuef47bff6Natm+loHpFtwbzD0x0ObLGJ\nuLL3487Pq4JU52/Dhg0EBweXuSe4FGbvOLWPt+FFf56whkVS65ph1BvwEs2bN2fIkCFceumlrF69\n2nQ0ISpNmYVZa+0CRlBYdDcDi7TWqUqpiUqp3gBKqb8opdKBW4A3lFKpRY/NAp6ksHSvBiYW3SaE\nX/jmm29o06YNDz/8MD179mTLli08+OCD1WqecllkSoZ3eFqYg2MbARrnoeq9H3NKSgoXXXRRma9l\nKcze8b+Db6rHoj9PBMclsWLFCubNm8eePXu49NJLeeCBBzh27JjpaEJ4nSdzmNFaLwGWlLjt0dM+\nXk3hdIvSHvsW8FYFMgpRZTydG1hwIoesb2ZxctN3BEXXpc7Nj7O6SQc6T/8N+I20KT0rN6gX2e12\ndu2q3uWtutPuAnT+SQ9HmBsDkJ+5i5D46jtdZ8OGDVx55ZVlXhcVFQVIYa4ot+PcJ0VWV0op+vbt\nS69evRg/fjyvvPIKH330EdOnT+f6668v+wmE8BFy0p8Q5aC15njqcjJmD+Pk1h+I6tSHeoNmENak\ng+lo503mMFecJ6f8FQuKjkPZQqv1POasrCz27t1b5vxlAKvVSlRUlBTmCvKFOcznEhUVxYwZM/jx\nxx+x2+3ccMMN3HjjjezdK6cFCv8ghVkID7mOHuTgB49z+PPnscXUJ37gK0R36YvFFmI6WoXISW0V\n97+308suO0pZsNVuWK33YvZ0wV8xmQdfcb40h/lcLr30UtauXcuUKVNYunQpF110Ea+++qosChQ+\nTwqzEGXQ2s2xXxeT8eZwHHtSsV99D3F9p2KrnVD2g31AdHQ0TqeTkydPmo7is8pTmKFwb1tnZlq1\n/SWluDB7emS7FOaKKzwpMgQV5DvrH87GZrMxZswYUlNTufTSSxk+fDhXX301O3fK3s3Cd0lhFuIc\nnNkZHFgwlqwvXyMk/kLiB8+g5iXXoSxW09G8RhZtVVx5C7OtTmPcuUer7cK/5cuXU6tWLerVq+fR\n9VKYK67Aw20JfUlSUhLLli1j1qxZrFmzhlatWvHaa6/JaLPwSVKYhSiF1ppj65awb859ODN3Ueva\nUdT510SCouJMR/O64sIs85jPX3nfTq9x0eVYQiPJ+uJVtK5e5eHLL7/ko48+Yvjw4R4f3y6FueK0\n44TPzl8+F6UUd911Fxs3bqRTp07ce++9/OMf/yAtLc10NCHKRQqzECW4jh3i4PuPkfXFq4TUb0G9\nwTOIaHWVx+XB18gIc8WVd4TZGh6F/e+DcKRv4vhvX1RmtHJxO/MYOnQozZo1Y9y4cR4/TgpzxXly\nUqQva9iwIcuWLWPmzJmsXr2aiy++mDfeeKPaTksSoiSPtpUTIlCc2PQdWV++hnY5ifnHMCLaXeu3\nRbmYFOaKK29hBqjR6mqOp35D9rdzCL/gr1gjzB8vfOSH9zi6YwfLly8nNDTU48dJYa44d95xgiJr\nmY5RbuU/pjuemndM4/CSlxk6dCijnnuLWj0eYM/Lt1dKPiG8RQqzEEBB7jGyvniVk1tWEBzfnNo9\nR2OLqW86VpWIjo4GpDBXhDvvOFht5doxRSlFrW4jyHhrBFlfvkbNS2+uxIRlKziRzdFfPmLQoEFc\nccUV5Xqs3W7H4XCQm5tLWFhY5QT0c+6841hiG5mOUSWCatahzq0TObb2M7K/nUPGW8NJvjaG7t27\nm44mxFlJYRYBLy89lUOfPkfBiSyiu/Sj5qU3+9WivrLICHPFufOOYz2Pt9NtMfWJ6nQrR1bM4+Tv\nP1ZCsvKxhEfx7LPPlvtxp7+GpDCfH09PivQXSlmo2eF6Qhu25tBnz9GjRw/uv/9+pk6dWq53N4So\nKlKYRcDS7gKO/PQ+R1YuICiqDnXveJaQes1Mx6pyxSe1yaK/81eRshN12a2ENmiJ25nn5VTlF1wn\niZiYmHI/7vTCHB8f7+1Yfq88J0X6m+A6janb/wX+6fqOl19+mW+++YYFCxZ4vAe4EFVFCrMISK5j\nhzn0+fM4dqcQflFXanUbjiUk3GvPX/55febISW0V53acf2FWShHa0LfLgbxLUTFuR+Ee6IFYmAEs\nthA+sV1DnZtj2LT0JVq3u4SYq+8hok23cq0hSZvSsxJTikAnu2SIgJP7x2r2zbmP/H1bqdXjAWpf\n92+vlmVfFB0dLWWnAtx5/rklmKekMFeMv5zyV1FhTToQf+d0QhMuJmvZdA5//jzu/FzTsYQAZIRZ\nBBDtLiDn+3c4+vOH2GITie09xm9O66so2eWgYtx5x7HVbmg6hjGyl3fFnM8uK/7KWiOaOrc8zpFV\niziycgGO/duJvWEswbGJpqOJACcjzCIgFJzI4cDCCRz9+UMi2nanXv8XpCyfxm63S9mpAHfecSwh\n/nVKW3nITisV43acAPC7k/7Ol7JYif5bH+rc9hRux3H2v/Mgxzd8bTqWCHAywiyMq+z5vo69W8j8\n79O4845R69pRRLS6qlK/ni+y2+1s3brVdAyfpLUbtyMwF2wVk8JcMTLCXLqwRm2IH/gKmZ89w+El\nL5K3ZwMx/xhWru0bhfAWGWEWfktrzbFfF7N/wViUNYi6dzwnZfksZA7z+StcsKUDuuwEBQURGRkp\nr6HzJHOYz84aYSfu1qeIuuxWTmz4igMLxuI6mmk6lghAUpiFX3I78zi8+AWyvnyNsMbtqDtwGsFx\nSaZjVVsyh/n8yehgIXkNnT95DZ2bsliJvrwfsTdOwJmVzr63R5GXnmo6lggwUpiF33EdzeTA/DGc\nSP2WqM59ib1pwnkdKhFI7HY7ubm5OBwO01F8jpSdQlKYz5/bcRwsQSiZanBO4U3/Sr1+L2AJCefA\nu+M5tm4JWmvTsUSAkDnMwq84MrZy8KMn0U4HsTc/SniTv5iO5BNO3+UgLi7OcBrfIoW5kBTm81d8\n8E159hwOVLbaCdTr/wKHPnuOrC9eJf/AH8RcPRQVZPPa1/DGuhrZE9r/yAiz8BsnNn3L/gVjsdhC\nqdvveSnL5SCLts6fFOZCUpjPX+E+3rJDhqcsoRHE3jSBmpf9i+O/LePAu+MoOCGvPVG5pDALn6e1\nm+zv3+HQZ88REt+cuv2eJziA98Q9H3LwxPmTBVuFpDCfv8JtCQP79VNeymLFfnl/al8/lvyDO9k3\n999s2rTJdCzhx2RKhvBp7vw8Di1+ntzfVxHR+hpirhmGsnrvrblAIQdPnD9/G2E+37ejszdkc+zg\nYRLHLpa3o8up8Gj1mqZj+KQaF3YmKKoOBz+cSKdOnfjoo4+48sorTccSfkhGmIXPKjiezYEFY8jd\n9jP2K4cQ0/0+KcvnqVatWgAcPHjQcBLf486TBVtQ+AuDdjnQLqfpKD6n4ORRrGGRpmP4rJB6zajX\n73nq169Pt27d+M9//mM6kvBDUpiFT3Jm7WX/vH/jzEovnMv2l+tlwUwFNGrUCKvVyvbt201H8Tmy\nYKtQ8Qh78Yi78Ix2OSk4cpAgez3TUXxaUFQcP/zwA1dccQV33nknEyZMkB00hFdJYRY+x7F3C/vn\nPYTbmUdcn6dlcZ8XBAcHk5iYyO+//246is8pLsyBTgrz+XHm7AM0QfZ401F8XnR0NEuWLGHw4ME8\n9dRT9O3bV7bKFF4jc5iFTzm5/WcOffIM1ogY6vxrIjYZlfGaZs2asW3bNtMxfE5hYZYdDooLc4EU\n5nJxZWcAYIupbziJf7DZbMyaNYsLLriAcePGcfDgQT7++GMiI2XKi6gYj0aYlVLdlVJblVLblVJj\nS7k/RCm1sOj+n5VSiUW325RSbyulNiilNiulxnk3vggkx9Ynk/nRJGy1G1L3jmelLHtZs2bN+P33\n3+VtzHIqXLAlI8ynRpgdUpjLw5m1FwCbjDB7jVKKsWPH8vbbb/Ptt99y5ZVXcujQIdOxhI8rszAr\npazADKAH0ALoo5RqUeKywUC21voC4EVgatHttwAhWutWwCXAPcVlWghPaa3JWTmfrGXTCW3cjrg+\nk7HWiDYdy+80a9aMEydOsG/fPtNRfIpMySgkUzLOjytrL5bwaHkNVYL+/fvz8ccfs3HjRrp06cKe\nPXtMRxI+zJMR5o7Adq31Dq11PvAecH2Ja64H3i76+APgKlW4AkYDNZRSQUAYkA8c9UpyERC01mR/\nM5sjP7xLjVZXU+fGCViCw0zH8kvNmjUDkHnM5eTOOy5HryOF+Xw5szOwxcjocmW57rrrWLZsGRkZ\nGfztb39jy5YtpiMJH+XJHOb6wOm/lqUDfz3bNVprl1LqCFCLwvJ8PbAPCAdGaa2zKhpaBAat3WQt\ne5XjvyUTeUlv7FcNCfidCM6Hp/vquo4Wbil305QPiUw+ccZ9sq9u6bR2F57SJodOYAkpnMcthbl8\nXFl7CU3qYDqGX7v88sv57rvv6NatG507dyY5OZkOHeRnLsrHkxHm0hpKyUmOZ7umI1AAxAONgQeV\nUkl/+gJK3a2UWqOUWpOZmelBJOHvtLuAw4tf5PhvydS87F9SlquANbI2KigYV9GcSlE27TgJaHk7\nHVDWIFRwmBTmcnA7TlJwIlsW/FWBtm3bsnLlSiIjI/n73//Ot99+azqS8DGeFOZ0IOG0zxsAGWe7\npmj6RRSQBdwOJGutnVrrg8APwJ9+rdNaz9Rad9Bad4iNjS3/dyH8inY5yfxkCidSlxN9eX/sl/eX\nslwFlLIQZI/HmS2F2VMFfnbKX0VZQiKkMJeD89QOGTIloyo0bdqUH374gYYNG3LttdfyzTffmI4k\nfIgnhXk10FQp1VgpFQzcBnxa4ppPgQFFH98MfKMLl9rvBq5UhWoAlwIygUicldvp4OBHT5H7+yrs\nVw4h6rJ/mY4UUGz2eJxZJX8fFmfjb8diV5QltIYU5nIofjdH9mCuOvHx8SxfvpykpCR69erF119/\nbTqS8BFlzmEumpM8AlgGWIG3tNapSqmJwBqt9afAm8BcpdR2CkeWbyt6+AxgDrCRwmkbc7TWKZXw\nfQg/4M7P4+CHT+DYvZGYbiOIbNvddKSAExQTz8ntP6PdBSiL1XScak8K85ksoTLCXB6FW8opgqJl\ni0xv8HS9BkDBleNxHXyYf3S/ltibHiUssW0lJhP+wKODS7TWS4AlJW579LSP8yjcQq7k446XdrsQ\nJbmdDjI/ehLHnlRq9RpNRMu/m44UkGz2+uAuwHXkoOxz7QEpzGeyhEbgypZtCT3lys7AWjMWiy3E\ndJSAYw2PIu62SRx472EyP5xI7I2PENa4velYohqTo7GFcdrlJPO/k8nblUKtHg9IWTYoqGjxkSz8\n84wU5jPJCHP5OLP3yoI/g4pLc1BMAw5++CS5O9aajiSqMSnMwiin00nmp1PJ27GWmG7DiWh1lelI\nAa34L29Z+OeZ4lPtpDAXsoRGyEl/HtJa48ySPZhNKyzNTxFcu2Hh+hkpzeIspDALY1wuF3379iV3\n20/Yr75H5ixXA5awmlhCasjCPw+5846DxYqyhZqOUi1YQiPQTgf5+fmmo1R77pNH0I4TBNllhNk0\na1hN6txaWJozP55MXnqq6UiiGpLCLIwoKChg4MCBvP/++9j/Poial1xnOpIAlFIExdSXKRkeKj4W\nW7Y9LFR84mF2drbhJNVf8bs4MiWjerCGRVLnXxOx1ozl4PtPkH/gD9ORRDUjhVlUObfbzT333MP8\n+fOZNGkSNTveaDqSOI0tpr7fdIGVAAAgAElEQVRMyfCQO++ETMc4jUUKs8dObSknhbnasIZHEXfr\nk1hCa3Bg0aM4D6ebjiSqESnMosqNGzeON998kwkTJjB+/HjTcUQJQfZ4Co5m4nY6TEep9tx5x+VY\n7NMU/yykMJfNmZUBliCCasphXdVJUM1Y4m59ClAcWDgB19GDpiOJakIKs6hS06ZN45lnnuHee+/l\niSeeMB1HlMJWdIiCK0e2BytL8ZQMUUhGmD3nzN5LUHRd2e+8GrLF1CfuXxNx55/kwMIJFJzIMR1J\nVANSmEWVWbhwIaNGjeLGG2/k5Zdflnmf1VTxW8ROmcdcJinMZ5LC7DlXVobMX67GguOSqHPzYxQc\nPcSBRY/KdolCCrOoGsuXL6d///507tyZ+fPnY7XKqEp1ZYtpAMqC8+BO01GqPSnMZ5LC7Bm304Ez\nKx1b7Yamo4hzCG3Qgth/jsd5aDcHP3oK7XKajiQMksIsKt1vv/3GDTfcQNOmTfnkk08IDZUtuKoz\nS3AotthGODK2mo5SrbndbtwOWfR3OinMnsk/8Ae4CwiJb246iihDWNIl1O45EseejRxOfhmttelI\nwhApzKJSpaWl0aNHD2rWrMnSpUux2+2mIwkPhMQ3x7Hvd7R2m45SbR07dgy0G2toDdNRqg1lDULZ\nQqUwlyG/6JfRkHpSmH1BjRZXENW5LydSl3Pkx/dMxxGGSGEWlSYrK4vu3buTm5tLcnIyCQkJpiMJ\nD4XEX4h2nMB1WOYxn01xKZQR5jNZQiOkMJfBkbEVa1Qc1ggZQPAVUZ1uo8bFV3Jk5XyOpy43HUcY\nIIVZVAqn08m//vUvdu7cyaeffkrLli1NRxLlUDzy5cjYYjhJ9SWFuXRSmMvmyNhKSL1mpmOIclBK\nUav7fYQ0bMXhpdPI27PRdCRRxaQwi0rx4IMP8vXXXzNz5ky6dOliOo4op6Ba9VEhNWQe8zlIYS6d\nFOZzcx07TMGxTELiLzQdRZSTstqI/efDBEXVJfOjSbKTUICRwiy8btasWbzyyiuMHj2aAQMGmI4j\nzoNSFkLqNcOxTwrz2UhhLp0U5nPLL/ozJQv+fJM1NII6tzwOFgsHP3icgtyjpiOJKiKFWXjV999/\nz7333kv37t155plnTMcRFRASfyHOzF2483NNR6mWpDCXzhIihflcHBlbwRpEcFwT01HEebJF16XO\njY/gOnqIzI8moQtku7lAIIVZeE1aWho33XQTTZo04d1335W9ln1cSHwz0G7y928zHaVaOlWY5Wjs\nM1hCa0hhPgdHxlaC6yShgmymo4gKCKl/EbWvHYkjPZWsr2ebjiOqQJDpAMK3JY5dDIA7P5f98x6i\n4Fguwf8cTdspPxhOJioqOL544Z9MyyhNdnY2KAsqOMx0lGrFEhrBsRMncDqd2GxSCk+n3QXk79tG\nRJtrTEcRXlCjRVfyD/zB0V8+Ijguicg23UxHEpVIRphFhWnt5tDnz+M8tJva14+R4179hDWsJkH2\neCnMZ5GdnY0lNEKOeC9BDi85O2dmGtrlkAV/fiS66wBCE9uR9eVrOPZuNh1HVCIpzKLCjqxaRO62\nn7BfOZiwxu1NxxFeFBLfnPyMrXK6VSmKC7M4kxTmsyv+5TNYFvz5DWWxUrv3/xEUWZvM/z5NwQl5\n3fsrKcyiQnLT1nNkxXxqtLiCyEt6m44jvCwkvjkFJ7LZvXu36SjVTmFhjjQdo9qxFv1MpDD/mSNj\nK5bwKIKi4kxHEV5kDYsk9p8P4847Tuanz6LdBaYjiUoghVmct/T0dA59+gy22gnEdBshb037oeCi\nt45XrVplOEn1k5WVhSVMRphLsoQVFuasrCzDSaofR8ZWQuKby/8r/VBwncbEXHMvjt0p5KycbzqO\nqARSmMV5KT7JTxc4ib1hHJbgUNORRCUIjk1EBYfx7bffmo5S7WRlZZ0aTRX/UzwlQwrzmTIyMnBl\npRNSv4XpKKKSRLS6mojW13B01SIWL15sOo7wMinM4ryMHz+eVatWUav7fdhqJZiOIyqJsgYR2qgN\nS5culXnMJRSOMEthLklGmEuXnJwMQFiSrPPwZ/ar78FWJ4l+/fqxZ88e03GEF0lhFuWWnJzMc889\nx9ChQ6lx0eWm44hKFpZ0Cbt372bzZlkBXqygoICcnBxZ9FcKS0gNQApzSUuXLsUaEYMttrHpKKIS\nWWwhxF4/BqfTSd++fXG5XKYjCS+RwizKZf/+/QwYMICLL76YF154wXQcUQXCki4BCv/CF4VycnIA\nZNFfKZTFSnR0tBTm07hcLr788ktCG18i85cDgC2mPq+++iorVqzgqaeeMh1HeIlHhVkp1V0ptVUp\ntV0pNbaU+0OUUguL7v9ZKZV42n2tlVKrlFKpSqkNSimZ7Oqj3G43/fr149ixY7z33nuEhcmBDYEg\nqGYdWrRoIYX5NMVlUKZklC4mJkYK82lWrVrFkSNHTv3yKfxfv3796NevH08++STfffed6TjCC8os\nzEopKzAD6AG0APoopUquWhgMZGutLwBeBKYWPTYImAcM1Vq3BK4A5NB1H/XCCy/w1Vdf8dJLL9Gy\nZUvTcUQV6tGjBytWrOD48eOmo1QLxWXQKlMySiWF+UzJyclYrVbCEtuajiKq0IwZM0hKSuKOO+6Q\nbRb9gCdHY3cEtmutdwAopd4Drgc2nXbN9cDjRR9/AExXhe87XQOkaK1/A9BaH/ZSblFBxUdaeyr/\n4E72vTOOsGaXMemPeCaX8/HCt/Xo0YPnn3+e5cuXc91115mOY9zhw4X/K5MpGaWLiYk59TMShdOZ\nOnXqxG75BSugREZG8u6773LZZZcxYsQI5s+X7eZ8mSdTMuoDpy/1TC+6rdRrtNYu4AhQC2gGaKXU\nMqXUr0qp/6t4ZFHVtMvJoc+fxxIaQS3Zbzkgde7cmRo1asi0jCIyJePcatWqJSPMRfbv38+6devo\n0aOH6SjCgA4dOjBhwgQWLFjAokWLTMcRFeBJYS6tHZXcX+ps1wQBnYG+Rf/+p1Lqqj99AaXuVkqt\nUUqtyczM9CCSqEo5K+fjzEyjVvf7sYZHmY4jDAgJCeGqq66S7eWKSGE+N5mS8T/F28lJYQ5c48eP\np2PHjgwdOpSMjAzTccR58qQwpwOnb7TbACj5X/zUNUXzlqOArKLbv9NaH9JanwSWAH/ahFJrPVNr\n3UFr3SE2Nrb834WoNHnpmzj684dEtL6G8As6mo4jDOrRowdpaWls3brVdBTjThXmoi3UxJliYmLI\nzs7G7XabjmLc0qVLqVevHm3atDEdRRgSFBTE3LlzycvLY/DgwTLo4KM8KcyrgaZKqcZKqWDgNuDT\nEtd8Cgwo+vhm4Btd+IpYBrRWSoUXFemunDn3WVRjbqeDw0tewhpVB/uVd5mOIwwrHiFbsmSJ4STm\nZWVlER0djbJYTUeplmJiYtBac+TIEdNRjHI6nXzxxRd0795dprIFuGbNmvHss8+SnJzMf/7zH9Nx\nxHkoszAXzUkeQWH53Qws0lqnKqUmKqV6F132JlBLKbUdGA2MLXpsNvAChaV7PfCr1lpWi/mIIz8s\nwJWdQa3u92EJCTcdRxjWqFEj2rZty8KFC01HMS4rK4uYmBjTMaqt4p9NoE/L+PLLL8nJyeH66683\nHUVUA8OGDePyyy9n9OjR7Nu3z3QcUU4e7cOstV6itW6mtW6itZ5UdNujWutPiz7O01rforW+QGvd\nsXhHjaL75mmtW2qtL9Zay6I/H+HY9ztHf/mYiDbdZCskccrtt9/OL7/8wrZt20xHMUoK87lJYS40\nf/587Ha7zF8WAFgsFmbPnk1eXh7Dhg2TqRk+Rk76E3+iC5wcXvoy1hp27H8fZDqOqEb69OmDUop3\n333XdBSjpDCfmxRmOHHiBP/973+55ZZbCA4ONh1HVBNNmzblySef5JNPPpFdM3yMFGbxJ0d++gBn\nZhox3YbLoiZxhgYNGtC1a1fmz58f0KMjUpjPTQozfPLJJ5w8eZK+ffuajiKqmZEjR/KXv/yFESNG\nyH7lPkQKsziDMzuDI6sWEX5hF9kVQ5Sqb9++/P7776xdu9Z0FGOkMJ+bFObC6RgJCQl07tzZdBRR\nzQQFBfHmm2+SnZ3N2LFjTccRHpLCLE7RWpP1xWsoqw37VUNMxxHV1M0330xwcHDAnlrldrvJzs6W\nwnwOdrsdCNzCnJmZybJly7j99tuxWOSvWfFnrVq1YtSoUcyePZsff/zRdBzhAfmTLE45uWUleWnr\niL68H0ERUgZE6aKjo+nZsyfvvfceBQUFpuNUuaNHj+J2u6Uwn4PNZiMyMjJgC/OiRYsoKCiQ6Rji\nnB577DESEhIYNmwYLpfLdBxRBinMAgC34yTZ38wiOK4Jke2uNR1HVHN9+/Zl//79fPPNN6ajVLni\nEiiF+dwC+bS/+fPn06pVK1q1amU6iqjGIiIimDZtGikpKbz88sum44gySGEWQOHx1wXHs4npNlwO\nYxBl6tmzJ1FRUQG5Ab8UZs8EamHeunUrq1at4vbbbzcdRfiAG264gZ49e/Loo4+yd+9e03HEOUhh\nFjgP7+HYr58T0aYbIfWamY4jfEBoaCgDBw7k/fffD7gN+KUweyZQC/Mrr7xCcHAwd955p+kowgco\npXj55ZdxOp2MHz/edBxxDlKYBdnfvIkKCiG6yx2mowgfct999+FyuXjttddMR6lSUpg9E4iFOScn\nh//85z/06dOHuLg403GEj0hKSmL06NG88847/PLLL6bjiLOQwhzgcnesJXfHGqL/dhvWGtGm4wgf\n0qRJE3r16sXrr79OXl6e6ThVRgqzZwKxML/11lucOHGC+++/33QU4WPGjx9P3bp1GTlyZEDvcV+d\nSWEOYLrARdbXswiyxxN5yXWm4wgfdP/995OZmcnChQtNR6kyxSWweOs0Ubriwhwof/kXFBQwffp0\nOnfuTPv27U3HET4mMjKSyZMns2rVqoA/SbW6ksIcwI6tX4orKx37lYNRVpvpOMIHXXXVVbRs2ZJp\n06YFTDE6fPgwERERctxxGWJiYnC5XBw7dsx0lCrx+eefs3PnTh544AHTUYSPGjBgAO3bt2fMmDGc\nPHnSdBxRghTmAOV2nOTID+8S2qgNYU3kRD9xfpRS3H///axbt46VK1eajlMl5JQ/zwTaaX/Tpk0j\nISGBG264wXQU4aMsFgsvvfQS6enpvPLKK6bjiBKkMAeoo798jDv3KNFdB6KUMh1H+LA77rgDu93O\ns88+azpKlZDC7JlAKsxr165l+fLlDB8+nKCgINNxhA/r0qULvXr1YsqUKWRnZ5uOI04jhTkAHTx4\nkKOrPya8eWdC6jU1HUf4uPDwcEaPHs1nn33Gzz//bDpOpZPC7JlAKsyPPPIIMTExDBs2zHQU4Qcm\nTZrEkSNHmDp1quko4jRSmAPQpEmT0K582UZOeM0DDzxA7dq1eeSRR0xHqXRSmD0TKIV55cqVJCcn\nM2bMGGrWrGk6jvADrVu3pm/fvkybNk0OM6lGpDAHmLS0NF577TUiWv8DW60GpuMIPxEZGcm4ceP4\n6quv+Pbbb03HqVRSmD0TCIVZa83DDz9M3bp1GTFihOk4wo888cQTFBQU8OSTT5qOIopIYQ4wjz32\nGFarlahOfUxHEX5m2LBh1K9fn4cffthvd8zQWkth9lAgFOYvv/yS77//nkceeYTw8HDTcYQfSUpK\n4p577mH27Nls27bNdByBFOaA8scffzBv3jyGDRtGUM3apuMIPxMWFsaECRP48ccfWbp0qek4leL4\n8eO4XC5q1aplOkq1FxoaSnh4uN8W5uLR5UaNGjFkyBDTcYQfevjhh7HZbDz99NOmowikMAeUp59+\nGpvNxkMPPWQ6ivBTgwYNIikpiYceeoj8/HzTcbxOTvkrH38+7W/BggWsWbOGxx57TPbkFpWibt26\n3H333cydO5e0tDTTcQKeFOYAsWvXLt5++23uuusu6tWrZzqO8FM2m41p06axadMmnn/+edNxvE4K\nc/n4a2HOyspi9OjRdOzYkf79+5uOI/zYQw89hMViYcqUKaajBDwpzAFi6tSpKKUYM2aM6SjCz/Xq\n1YubbrqJiRMn8scff5iO41VSmMvHXwvz2LFjOXz4MG+88QZWq9V0HOHHGjRowJ133smcOXNIT083\nHSegSWEOAHv37uXNN99k4MCBJCQkmI4jAsC0adOw2WwMHz7crxYASmEuH38szD/88AOzZs1i5MiR\ntG3b1nQcEQDGjh2L2+3mmWeeMR0loElhDgDPPfccBQUFjB071nQUESDq16/PpEmTWLZsGQsXLjQd\nx2ukMJePvxVmp9PJPffcQ8OGDXn88cdNxxEBIjExkX79+jFr1iz2799vOk7AksLs57Kzs5k5cya3\n3347SUlJpuOIAHLvvffyl7/8hREjRvjNW4nF5c9utxtO4huKC7O/vMvwyCOPkJqayvTp04mIiDAd\nRwSQcePG4XA4mD59uukoAUsKs5+bNWsWJ0+eZPTo0aajiABjtVqZN28eDoeDPn364HK5TEeqsKys\nLMLCwggLCzMdxSfExMTgcDjIzc01HaXClixZwjPPPMM999zDddddZzqOCDBNmzbluuuu4/XXX/eL\nP0++yKPCrJTqrpTaqpTarpT60/v6SqkQpdTCovt/Vkollri/oVLquFLq396JLTzhdDp55ZVXuOKK\nK2SunTCiWbNmzJw5k5UrV/Loo4+ajlNhcmhJ+fjL4SXp6en079+fNm3a8OKLL5qOIwLUqFGjOHz4\nMHPnzjUdJSCVWZiVUlZgBtADaAH0UUq1KHHZYCBba30B8CIwtcT9LwL+eZJBNfbhhx+Snp7OqFGj\nTEcRAaxPnz4MGTKEp59+muTkZNNxKuTw4cNSmMuh+Gd1+PBhw0nOn8vlok+fPjgcDhYtWiTvLghj\nunbtStu2bXnppZf8ZpqTLwny4JqOwHat9Q4ApdR7wPXAptOuuR54vOjjD4DpSimltdZKqRuAHcAJ\nr6UWZdJa8+KLL3LBBRfQq1cv03FEgJs2bRo//fQTffr0YeXKlbRs2dJ0pHKbO3cuixcv5p///Kfp\nKNVe4tjFAORn7gcUf+15O7E3PYI1rGa5nidtSs9KSOc5rTXDhw9n5cqVzJ8/n2bNmhnNI3xH8Z+B\nijr9z4BSilGjRjFgwACWLVtG9+7dvfI1hGc8mZJRH9hz2ufpRbeVeo3W2gUcAWoppWoAY4AnKh5V\nlMeqVav45ZdfeOCBB7BYZKq6MCssLIzPPvuMsLAwunfvzp49e8p+UDWhtebpp5+mf//+dOnShVmz\nZpmO5DOCYxOpff0YHPu3sX/e/+E6csB0pHJ5/PHHmTlzJg8//DC333676ThCcNttt1G3bl2ZGmSA\nJ01KlXJbyfcCznbNE8CLWuvj5/wCSt2tlFqjlFqTmZnpQSRRlpdeeono6GgGDhxoOooQADRq1Ijk\n5GSOHTtGt27dfOJteofDwdChQxk/fjx9+/YlOTmZqKgo07F8So0LOxN365O4T+awf+6/cezdYjqS\nR1599VUmTpzI4MGDefLJJ03HEQKA4OBghg8fzhdffEFqaqrpOAHFk8KcDpx+2kUDIONs1yilgoAo\nIAv4K/CMUioNGAmMV0qNKPkFtNYztdYdtNYdYmNjy/1NiDPt37+fjz/+mMGDB8vWR6Jaad26NZ98\n8gk7duygR48e1bo0p6en07VrV2bOnMm4ceN45513CA4ONh3LJ4UmXEzdvs+igoLZv2Asx9YtqdZz\nMN9++21GjBhxalcCpUobExLCjKFDhxIcHMwbb7xhOkpA8aQwrwaaKqUaK6WCgduAT0tc8ykwoOjj\nm4FvdKEuWutErXUi8BIwWWstmwhWsrfffhuXy8WQIUNMRxHiT7p27cqiRYtISUmhc+fO7N6923Sk\nP/n6669p3749qampfPDBB0yePFmmNlWQrXYCdQdOIyyxLVlfvMrhxS/gzs8zHesMWmumTp3KwIED\nufLKK3nvvfcICvJkqY8QVad27drceOONzJ07V7aYq0Jl/g1QNCd5BLAM2Aws0lqnKqUmKqV6F132\nJoVzlrcDowE5Us4QrTWzZ8+mS5cuNG/e3HQcIUrVu3dvli1bxr59+7jsssvYuHGj6UgAnDx5kpEj\nR3L11VdTu3ZtVq9ezU033WQ6lt+whkYQe/OjRHXuy4nUb9n3n/vIS99U9gOrgNvtZtSoUYwdO5Zb\nb72VxYsXEx4ebjqWEKW66667yMnJ4eOPPzYdJWB4NGSitV6itW6mtW6itZ5UdNujWutPiz7O01rf\norW+QGvdsXhHjRLP8bjW+jnvxhclff/992zfvp277rrLdBQhzqlr166sWLECrTWdOnVi3rx5RvOs\nWrWKdu3aMW3aNIYPH87q1au58MILjWbyR0pZiP5bH+L6TEK73RyYP4bs5W/hdjqMZTpw4AA9e/Zk\n2rRpPPDAAyxYsICQkBBjeYQoy9///ncaN27M7NmzTUcJGPJek5+ZPXs2UVFR3HzzzaajCFGmVq1a\n8fPPP9OnTx/69evHsmXLmDFjBjVrnn37MW9v15SRkXFqjnLDhg356quvuOqqq7zyNcTZhTZsTfyd\nr5D97Vsc/eUjTmxZgf2KOwm/sEuVzhlOTk5mwIABHD16lFdffZWhQ4eW+vW99boTwhssFgt33XUX\nDz/8MNu3b+eCCy4wHcnvyaQ8P5Kdnc0HH3xA37595a1E4TMSEhL49ttveeKJJ1iwYAGtW7fmvffe\nq/RFYTk5OUycOJGmTZvy3nvvMWbMGDZs2CBluQpZQsKp1W0EcX0mYwmN4NCnz3Bg/hjy9mys9P/+\ne/fuZdCgQfTo0YO4uDjWrFnDsGHDZIGf8BkDBw7EYrHw1ltvmY4SEKQw+5H58+eTl5cn0zGEzwkK\nCuLRRx9lxYoVREVF0adPHy699FK+++47rxcn17FDZC9/i4SEBB577DG6devGpk2bmDJlyjlHtkXl\nCW3YmnoDXiKm+304szM4sGAs++f9m08++QS32+3Vr5WTk8MjjzxC06ZNmT9/Pv/3f//HL7/84pOH\n6YjAFh8fT8+ePZkzZw5Op9N0HL8nhdmPvPnmm7Rv35527dqZjiLEeenUqRO//vorc+bMYe/evVxx\nxRW0bduWGTNmkJOTc97P687P5Xjqcg4seoy9rw3i6Or/ct1117Fu3To++ugjmjRp4sXvQpwPZbES\n2aYb9YfOJuYfQyk4kcMNN9xAYmIiY8aMISUl5byfW2vNqlWruPPOO4mPj2fSpEnccMMNbNmyhalT\npxIaGurF70SIqnPXXXexf/9+lixZYjqK31PVbS/MDh066DVr1piO4XM2b95MixYtmDZtGvfff3+Z\n18t8PFEe3jqiuDyvO7czjxOpyzm+Ppn8A3+ANYiQ+AsJbdSG0ISWBEXHY42woyzWMx6nC1y4jh3C\nlbUXx97N5O3ZiCNjKxQ4sdasQ40WXYlo0429rw+u0u9HlI92F/DcXxy88847JCcn43K5aNSoEZdf\nfjldu3albdu2NG7cGLvd/qdpFMePH2fXrl2sXbuWe5+dS97u3yg4mokKDqPGRV2JbHctwXFJhr4z\nIbxHuwtIn96P0MS2xPb+v1KvMX3EfHWnlFqrte5Q1nWy6M9PvP/++yilZLGf8BsWWyiRbXsQ2bYH\njv3bObn5e/J2p3Bk5QKOFB82agkqLM3KUnhLgZOCEzmgi97GVxaC45oQ2b4n4c0uI6T+RSglb6z5\nAmWxcsstt3DLLbdw6NAhPvjgA7766iuSk5OZO3fuqetq1qx5qjRbLBZycnLIyso6db8lrCahCRcT\n2uk2alzYBUuIrO8Q/kNZrIQ378SJ1G9xO/Ow2OTdksoihdlPLFy4kC5duhAfH286ihBeF1L3AkLq\nFq4CL8g7Tn7GVlxHD+I6cpCC44dBa1AWUBaCImsRFBVHUFQcwXUvkILkB2rXrs3QoUMZOnQoWmu2\nbt3K5s2b2blzJzt37uTYsWNorXG73URGRtKoUSMaNWrERRddRO8Fu+WXJOHXwi/swvH1yeT+sYYa\nF3Y2HcdvSWH2A6mpqWzatInp0+UQReH/rKERhCVdYjqGMEQpxYUXXujxHtnq3fRKTiSEWaEJF2MJ\nj+bklpVSmCuRFGYfU9qcyZwV80FZmLI5imdlTqWoBDJXV5ggrzshynZqWsaGr3Hn52EJlmkZlUHe\np/JxWmtObFlBSMLFWCPspuMIIYQQoorVuLAL2uUg94/VpqP4LSnMPs55aBeurHR5G0YIIYQIUCEN\nWmCtYefklhWmo/gtKcw+7uTmFaAshDfrZDqKEEIIIQwonJbxN3J3rMGdn2s6jl+SwuzDtNac2LqS\n0IatsNaINh1HCCGEEIaEX9gZ7cond/svpqP4JSnMPsyVsw9X1l7Cm11mOooQQgghDApp0AJLeBS5\nO+Twt8oghdmH5e1cB0BoohyFLYQQQgQypSyENmpLbto6dPHhTcJrZFs5H5abtg5rVBxBdjmsRIjy\nki3LhBD+JqxxO05u/g5nZhrBdeT4d2+SEWYfpQtc5O1KISyxHUop03GEEEIIYVhoYlsAcneuN5zE\n/0hh9lGOfb+j808S2limYwghhBACgiJrY6vdkLy0daaj+B0pzD4qb+c6UBZCG7UxHUUIIYQQ1URo\nYjvy9mzE7XSYjuJXpDD7qNy0Xwmu2xRraITpKEIIIYSoJsIS20GBE0d6qukofkUKsw9y5x0nf982\nwmQ6hhBCCCFOE5JwMViDyEuTeczeJIXZB+XtSgHtlvnLQgghhDiDJTiU0AYtyN35q+kofkUKsw/K\nTfsVFRxGSL3mpqMIIYQQopoJTWyPMzONguPZpqP4DSnMPigvbT2hDVujrLKNthBCCCHOdGp7uV0y\nLcNbpDD7GNexw7hy9hPaqLXpKEIIIYSohoLjkrCERuDYs9F0FL8hhdnH5B/cAUBw3QsMJxFCCCFE\ndaSUheC4pFOdQVScRz1HaiUAAAvVSURBVIVZKdVdKbVVKbVdKTW2lPtDlFILi+7/WSmVWHT7P5RS\na5VSG4r+faV34wce54Giwhzb2HASIYQQQlRXtjpJODN34XK5TEfxC2UWZqWUFZgB9ABaAH2UUi1K\nXDYYyNZaXwC8CEwtuv0QcJ3WuhUwAJjrreCBKv/AHwTZ62EJCTcdRQghhBDVVHBcE7Qrn61bt5qO\n4hc8GWHuCGzXWu/QWucD7wHXl7jmeuDtoo8/AK5SSimt9TqtdUbR7alAqFIqxBvBA1X+wZ0E10ky\nHUMIIYQQ1VhwncJ3otevl4V/3uBJYa4P7Dnt8/Si20q9RmvtAo4AtUpccxOwTmstZzWepyNHjuDK\n2UdwXBPTUYQQQghRjdlqJYDVxrp160xH8Que7EumSrlNl+capVRLCqdpXFPqF1DqbuBugIYNG3oQ\nKTClpKQA//utUQghhBCiNMpiJTg2UUaYvcSTEeZ0IOG0zxsAGWe7RikVBEQBWUWfNwA+Bvprrf8o\n7QtorWdqrTtorTvExsaW7zsIIMW/JdpkSoYQQgghyhBcpzHr1q1D65LjnKK8PCnMq4GmSqnGSqlg\n4Dbg0xLXfErhoj6Am4FvtNZaKRUNLAbGaa1/8FboQLV+/Xos4dFYI2JMRxFCCCFENRcc14SsrCzS\n09NNR/F5ZRbmojnJI4BlwGZgkdY6VSk1USnVu+iyN4FaSqntwGigeOu5EcAFwASl1Pqif+p4/bsI\nEOvWrSO4TmOUKm0GjBBCCCHE/xS/Iy3zmCvOo7OVtdZLgCUlbnv0tI/zgFtKedxTwFMVzCiA/Px8\nUlNTCWtfcoMSIYQQQog/C66TiFKK9evX07t377IfIM5KTvrzEZs2bcLpdMqCPyGEEEJ4xBIcRtOm\nTWWE2QukMPuI4lWusqWcEEIIITzVrl072SnDC6Qw+4h169YRHh5OkL2e6ShCCCGE8BFt27YlLS2N\n7Oxs01F8mhRmH7F+/XratGmDslhNRxFCCCGEj2jXrh0Av/32m+Ekvk0Ksw9wu92sX7+etm3bmo4i\nhBBCCB9S3B1kHnPFSGH2Aenp6Rw9epTWrVubjiKEEEIIHxIXF0dcXBwbNmwwHcWnSWH2ATt37gSg\nSRNZ8CeEEEKI8klKSiItLc10DJ8mhdkHFL/IExMTjeYQQgghhO9JTEyUwlxBUph9QPGLvGHDhmaD\nCCGEEMLnJCYmsmfPHlwul+koPksKsw9IS0sjPj6ekJAQ01GEEEII4WMSExNxuVxkZGSYjuKzpDD7\ngLS0NBo3/v/27i9GrrKM4/j3l9YtoikQG2O63bglFrVgRQIVJf4talFSbtS0G0yjGxsJVCAabSGS\naEoi0IhNxIsGalQaKpSqjUEBg/GuhVpAqLW6aTfdXTBAFLQx0MzyeHFOzXY7OztnZrrvO+zvczVz\n5kz765PO5Jkz7zyvd/gzMzOz6k70EF6W0To3zF1geHjY65fNzMysJSd6CDfMrXPDnLlarcbIyIgb\nZjMzM2vJid9AuWFunRvmzI2NjTE+Pu6G2czMzFoyb948Fi5c6Ia5DW6YM+eRcmZmZtYuj5Zrjxvm\nzLlhNjMzs3a5YW6PG+bMDQ8PI4m+vr7UUczMzKxLeRZze9wwZ84zmM3MzKxdnsXcHjfMmfNIOTMz\nM2uXR8u1xw1z5rxpiZmZmbXLm5e0xw1zxjyD2czMzDqhr68PSW6YW+SGOWOewWxmZmad4FnM7XHD\nnLEjR44AHilnZmZm7evv7/9/b2HVuGHOmGcwm5mZWad4FnPr3DBnzDOYzczMrFM8i7l1bpgzNjw8\nTG9vLz09PamjmJmZWZfr7+9nfHycsbGx1FG6TlMNs6SVkg5JGpK0oc7j8yT9onx8r6T+CY9tLI8f\nkvSZzkV/4/MMZjMzM+sUz2Ju3bQNs6Q5wF3AFcBSYI2kpZNOGwT+FRHvAu4EbiufuxRYDZwPrAR+\nXP551gTPYDYzM7NO8Szm1jVzhXk5MBQRhyPiOLADuGrSOVcBPy1v7wRWSFJ5fEdEvBYRR4Ch8s+z\nadRqNUZHR32F2czMzDrCs5hbN7eJc3qBkQn3R4EPTnVORNQkvQK8rTy+Z9Jze1tOexotW7Ysq1Er\nEeEZzGZmZtYxPT099Pb2smnTJjZv3pw6zkm2bt3KmjVrUseYUjMNs+ociybPaea5SFoHrCvvHpN0\nqIlcp8MC4KVEf3ddg4ODDA4Opo4xlezqlTnXqxrXqxrXqxrXqxrXq7osaqbb6h8/duzYzAaZxsDA\nwIKBgYEU9XpnMyc10zCPAhPnmi0CnpvinFFJc4GzgH82+VwiYiuwtZnAp5OkfRFxceoc3cL1qsb1\nqsb1qsb1qsb1qsb1qs41qyb3ejWzhvkJYImkxZJ6KH7Et3vSObuBteXtzwOPRUSUx1eXUzQWA0uA\nxzsT3czMzMzs9Jv2CnO5Jvk64GFgDrAtIg5I+h6wLyJ2A/cAP5c0RHFleXX53AOS7gf+AtSAayNi\n/DT9W8zMzMzMOq6ZJRlExEPAQ5OO3TLh9qvAF6Z47q3ArW1knEnJl4V0GderGterGterGterGter\nGterOtesmqzrpWLlhJmZmZmZ1eOtsc3MzMzMGnDDPImkCyXtkfSUpH2SvNHKNCStL7c+PyDp9tR5\nuoGkb0oKSQtSZ8mZpDsk/VXSnyX9UtLZqTPlSNLK8jU4JGlD6jw5k9Qn6Q+SDpbvWdenztQNJM2R\n9KSk36TOkjtJZ0vaWb53HZT0odSZcibpxvK1+Kyk+ySdkTpTPW6YT3U78N2IuBC4pbxvU5D0CYod\nHZdFxPlAXpPQMySpD/gUcDR1li7wKHBBRCwD/gZsTJwnO5LmAHcBVwBLgTWSlqZNlbUa8I2IeC9w\nKXCt69WU64GDqUN0iS3A7yLiPcD7cd2mJKkX+DpwcURcQDFcYnXaVPW5YT5VAPPL22dRZ260neQa\n4PsR8RpARLyQOE83uBP4FnU28bGTRcQjEVEr7+6hmOVuJ1sODEXE4Yg4Duyg+BBrdUTE8xGxv7z9\nH4pmJssdaHMhaRHwOeDu1FlyJ2k+8FGK6WFExPGIeDltquzNBd5c7uNxJpn2XW6YT3UDcIekEYqr\npb6i1dh5wEck7ZX0R0mXpA6UM0mrgLGIeDp1li70FeC3qUNkqBcYmXB/FDeATZHUD3wA2Js2SfZ+\nSPEh//XUQbrAucCLwE/KJSx3S3pL6lC5iogxil7rKPA88EpEPJI2VX1NjZV7o5H0e+AddR66GVgB\n3BgRD0r6IsWnxMtnMl9upqnXXOAciq82LwHul3RuzOLxK9PU6ybg0zObKG+N6hURvy7PuZniq/Tt\nM5mtS6jOsVn7+muWpLcCDwI3RMS/U+fJlaQrgRci4k+SPp46TxeYC1wErI+IvZK2ABuA76SNlSdJ\n51B8I7YYeBl4QNLVEXFv2mSnmpUNc0RM2QBL+hnFWi2AB/BXUNPV6xpgV9kgPy7pdWABxSfsWWmq\nekl6H8WbwtOSoFhesF/S8oj4xwxGzEqj/18AktYCVwIrZvMHsQZGgb4J9xeR6VeauZD0JopmeXtE\n7EqdJ3OXAaskfRY4A5gv6d6IuDpxrlyNAqMRceJbi50UDbPVdzlwJCJeBJC0C/gwkF3D7CUZp3oO\n+Fh5+5PA3xNm6Qa/oqgTks4DeoCXkibKVEQ8ExFvj4j+iOineGO9aDY3y9ORtBL4NrAqIv6bOk+m\nngCWSFosqYfiBzO7E2fKlopPq/cAByPiB6nz5C4iNkbEovI9azXwmJvlqZXv5yOS3l0eWkGx27HV\ndxS4VNKZ5WtzBZn+SHJWXmGexleBLeXi81eBdYnz5G4bsE3Ss8BxYK2vAloH/QiYBzxaXpXfExFf\nSxspLxFRk3Qd8DDFL8y3RcSBxLFydhnwJeAZSU+Vx24qd7Q164T1wPbyA+xh4MuJ82SrXLayE9hP\nsezuSTLd8c87/ZmZmZmZNeAlGWZmZmZmDbhhNjMzMzNrwA2zmZmZmVkDbpjNzMzMzBpww2xmZmZm\n1oAbZjMzMzOzBtwwm5mZmZk14IbZzMzMzKyB/wGI9/mAV+b38QAAAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x280945fed30>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# plot samples and target distribution\n",
"plot_x = np.linspace(-8,8,200)\n",
"plt.figure(figsize=(12,4))\n",
"plt.plot(plot_x, [target(_) for _ in plot_x], color='black') #plot the distribution\n",
"plt.hist(samples,density=True, bins=30) #plot the samples\n",
"plt.show() #very similar to the PCW plot, hurray"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment