二项分布和Beta分布

时间:2021-12-11 10:59:42

原文为: 二项分布和Beta分布

二项分布和Beta分布

In [15]:

%pylab inline
import pylab as pl
import numpy as np
from scipy import stats
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.

二项分布

在概率论和统计学中,二项分布是n个独立的[是/非]试验中成功的次数的离散概率分布,其中每次试验的成功概率为$p$。举两个例子就很容易理解二项分布的含义了:

  • 抛一次硬币出现正面的概率是0.5($p$),抛10(n)次硬币,出现k次正面的概率。

  • 掷一次骰子出现六点的概率是1/6,投掷6次骰子出现k次六点的概率。

在上面的两个例子中,每次抛硬币或者掷骰子都和上次的结果无关,所以每次实验都是独立的。二项分布是一个离散分布,k的取值范围为从0到n,只有n+1种可能的结果。

scipy.stats.binom为二项分布,下面用它计算抛十次硬币,出现k次正面的概率分布。

In [16]:

n = 10
k = np.arange(n+1)
pcoin = stats.binom.pmf(k, n, 0.5)
pcoin
Out[16]: array([ 0.00097656, 0.00976563, 0.04394531, 0.1171875 , 0.20507813,
0.24609375, 0.20507813, 0.1171875 , 0.04394531, 0.00976563,
0.00097656])In [17]: pl.stem(k, pcoin, basefmt="k-")
pl.margins(0.1)

aaarticlea/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXYAAAD5CAYAAAAzzx7cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFpNJREFUeJzt3WFsU9fdx/HfBVuJBBqFiqZgZ0obp8ShJDAlRKwCeSuJC1GjkfEiG1NRm0URE4RNffWUBy3pIzqiqdoSshfRVKpWRTTvCCKQptHq0kETa6MV02jVzCKrCQtDG1kpGQmx7vOCzW1KcJzi+JLj70ey5Ht9js//xNc/XV9f31i2bdsCABhjgdMFAABSi2AHAMMQ7ABgGIIdAAxDsAOAYVxODm5ZlpPDA8C8leiERsf32G3bTtvt5z//eVrHS/eN+c3fm8lzY36pv83E8WAHAKQWwQ4AhsmoYA8EAk6XMKeY3/xl8twk5pdulp3MAZu5GtyykjpeBAD4wkzZmVF77ACQCQh2ADAMwQ4AhiHYAcAwBDsAGIZgBwDDzBjsPT09KiwsVEFBgVpaWu54/MiRIyopKVFxcbGeeOIJnT9/Pv5YXl6eiouLtW7dOq1fvz61lQMAppXwPPZYLKZVq1apr69PHo9HZWVlOnr0qPx+f7zN+++/r6KiIi1ZskQ9PT1qampSf3+/JOmRRx7RH//4Ry1btmz6wTmPHQBm7Z7OYw+Hw/L5fMrLy5Pb7VZtba26urqmtNmwYYOWLFkiSSovL9elS5emPE5wA0B6Jbxs7/DwsHJzc+PLXq9XAwMDd23/yiuvaOvWrfFly7K0efNmLVy4UA0NDaqvr7+jT1NTU/x+IBC4736ai8zT3X1abW29Gh93KStrUo2Nlaqq2uR0WchgoVBIoVAo6fYJg30210t/5513dPjwYZ05cya+7syZM1qxYoWuXr2qiooKFRYWauPGjVP6fTnYAad1d5/W3r1vKRI5EF8XieyTJMIdjvnqTm9zc3PC9gkPxXg8HkWj0fhyNBqV1+u9o9358+dVX1+v48ePa+nSpfH1K1askCQtX75c27ZtUzgcTmoSgFPa2nqnhLokRSIHdOjQ2w5VBMxewmAvLS3V4OCghoaGNDExoc7OTlVXV09p8+mnn6qmpkZvvPGGfD5ffP3Y2JiuX78uSbpx44Z6e3u1Zs2aOZgCkDrj49N/iL15c2GaKwG+voSHYlwul9rb2xUMBhWLxVRXVye/36+Ojg5JUkNDg1588UVdu3ZNu3btkiS53W6Fw2GNjIyopqZGkjQ5OakdO3aosrJyjqcD3JusrMlp12dnx9JcCfD1cdle4EumO8aen/+CWluf4hg77htctheYhaqqTWptDSoY3C9JCgb3E+qYd9hjB+7CsiQ2T9yP2GMHgAxDsAOAYQh2ADAMwQ4AhiHYAcAwBDsAGIZgBwDDEOwAYBiCHQAMQ7ADgGEIdgAwDMEOAIYh2AHAMAQ7ABiGYAcAwxDsAGAYgh0ADEOwA4BhCHYAMAzBDgCGIdgBwDAEOwAYhmAHAMO4nC4AmEl392m1tfVqfNylrKxJNTZWqqpqk9NlpYTJc4NzCHbc17q7T2vv3rcUiRyIr4tE9knSvA9Ak+cGZ3EoBve1trbeKcEnSZHIAR069LZDFaWOyXODswh23NfGx6f/UHnz5sI0V5J6Js8Nzpox2Ht6elRYWKiCggK1tLTc8fiRI0dUUlKi4uJiPfHEEzp//nzSfYGZZGVNTrs+OzuW5kpSz+S5wWF2ApOTk3Z+fr598eJFe2Jiwi4pKbEvXLgwpc3Zs2ft0dFR27Zt+9SpU3Z5eXnSfWcYHrBPnHjXzs9/wZbs+C0//3/sEyfenfOx53rzdHJumN9mys6Ee+zhcFg+n095eXlyu92qra1VV1fXlDYbNmzQkiVLJEnl5eW6dOlS0n2BmVRVbVJra1DB4H5JUjC4X62tTxnx5aLJc4OzEp4VMzw8rNzc3Piy1+vVwMDAXdu/8sor2rp166z6NjU1xe8HAgEFAoFka0eGqKrapKqqTbIsqafn/5wuJ6VMnhtSJxQKKRQKJd0+YbBblpX0E73zzjs6fPiwzpw5M6u+Xw52AMCdvrrT29zcnLB9wmD3eDyKRqPx5Wg0Kq/Xe0e78+fPq76+Xj09PVq6dOms+gIAUivhMfbS0lINDg5qaGhIExMT6uzsVHV19ZQ2n376qWpqavTGG2/I5/PNqi8AIPUS7rG7XC61t7crGAwqFouprq5Ofr9fHR0dkqSGhga9+OKLunbtmnbt2iVJcrvdCofDd+0LAJhb1n9OnXFmcMuSg8NjnrGs2ycFmjheuueG+W2m7OSXpwBgGIIdAAxDsAOAYQh2ADAMwQ4AhiHYAcAwBDsAGIZgBwDDEOwAYBiCHQAMQ7ADgGEIdgAwDMEOAIYh2AHAMAQ7ABiGYAcAwxDsAGAYgh0ADEOwA4BhCHYAMAzBDgCGIdgBwDAEOwAYhmAHAMMQ7ABgGIIdAAxDsAOAYQh2ADAMwQ4Ahpkx2Ht6elRYWKiCggK1tLTc8fjHH3+sDRs2KDs7Wy+//PKUx/Ly8lRcXKx169Zp/fr1qasaAHBXrkQPxmIx7d69W319ffJ4PCorK1N1dbX8fn+8zYMPPqhDhw7p2LFjd/S3LEuhUEjLli1LfeUAgGkl3GMPh8Py+XzKy8uT2+1WbW2turq6prRZvny5SktL5Xa7p30O27ZTVy0AYEYJ99iHh4eVm5sbX/Z6vRoYGEj6yS3L0ubNm7Vw4UI1NDSovr7+jjZNTU3x+4FAQIFAIOnnB4BMEAqFFAqFkm6fMNgty7qnYs6cOaMVK1bo6tWrqqioUGFhoTZu3DilzZeDHQBwp6/u9DY3Nydsn/BQjMfjUTQajS9Ho1F5vd6ki1mxYoWk24drtm3bpnA4nHRfAMDXkzDYS0tLNTg4qKGhIU1MTKizs1PV1dXTtv3qsfSxsTFdv35dknTjxg319vZqzZo1KSobAHA3CQ/FuFwutbe3KxgMKhaLqa6uTn6/Xx0dHZKkhoYGjYyMqKysTJ999pkWLFig1tZWXbhwQX//+99VU1MjSZqcnNSOHTtUWVk59zMCgAxn2Q6etmJZFmfNIGmWJaVzc0nneOmeG+a3mbKTX54CgGEIdgAwDMEOAIYh2AHAMAQ7ABiGYAcAwxDsAGAYgh0ADEOwA4BhCHYAMAzBDgCGIdgBwDAEOwAYhmAHAMMQ7ABgGIIdAAyT8D8oAdPp7j6ttrZejY+7lJU1qcbGSlVVbXK6LCSB1y4zEOyYle7u09q79y1FIgfi6yKRfZJEQNzneO0yB4diMCttbb1TgkGSIpEDOnTobYcqQrJ47TIHwY5ZGR+f/kPezZsL01wJZovXLnMQ7JiVrKzJaddnZ8fSXAlmi9cucxDsmJXGxkrl5++bsi4//wXt2VPhUEVIFq9d5iDYMStVVZvU2hpUMLhfkhQM7ldr61N8+TYP8NplDsu2bduxwS1LDg6Pe2RZUjpfPpPHM3luSL2ZspM9dgAwDMEOAIYh2AHAMAQ7ABiGYAcAw8wY7D09PSosLFRBQYFaWlruePzjjz/Whg0blJ2drZdffnlWfQEAqZfwdMdYLKZVq1apr69PHo9HZWVlOnr0qPx+f7zN1atX9de//lXHjh3T0qVL9fzzzyfdl9Md5zfTT9HjdEfcr+7pdMdwOCyfz6e8vDy53W7V1taqq6trSpvly5ertLRUbrd71n0BAKmX8LK9w8PDys3NjS97vV4NDAwk9cTJ9m1qaorfDwQCCgQCST0/AGSKUCikUCiUdPuEwW5Z1tcuJNm+Xw52AMCdvrrT29zcnLB9wkMxHo9H0Wg0vhyNRuX1epMq5F76AgC+voTBXlpaqsHBQQ0NDWliYkKdnZ2qrq6etu1XD+TPpi8AIHUSHopxuVxqb29XMBhULBZTXV2d/H6/Ojo6JEkNDQ0aGRlRWVmZPvvsMy1YsECtra26cOGCFi9ePG1fAMDc4uqO+NpMP0WP0x1xv+LqjgCQYQh2ADAMwQ4AhiHYAcAwBDsAGIZgBwDDEOwAYBiCHQAMQ7ADgGEIdgAwDMEOAIYh2AHAMAQ7ABiGYAcAwxDsAGAYgh0ADEOwA4BhCHYAMAzBDgCGIdgBwDAEOwAYhmAHAMMQ7ABgGIIdAAxDsAOAYQh2ADAMwQ4AhiHYAcAwBDsAGGbGYO/p6VFhYaEKCgrU0tIybZvGxkYVFBSopKREH3zwQXx9Xl6eiouLtW7dOq1fvz51VQMA7sqV6MFYLKbdu3err69PHo9HZWVlqq6ult/vj7c5efKk/vKXv2hwcFADAwPatWuX+vv7JUmWZSkUCmnZsmVzOwsAQFzCYA+Hw/L5fMrLy5Mk1dbWqqura0qwHz9+XDt37pQklZeXa3R0VFeuXFFOTo4kybbthAU0NTXF7wcCAQUCga8xDQAwVygUUigUSrp9wmAfHh5Wbm5ufNnr9WpgYGDGNsPDw8rJyZFlWdq8ebMWLlyohoYG1dfX3zHGl4MdAHCnr+70Njc3J2yfMNgty0pq0Lvtlf/+97/XypUrdfXqVVVUVKiwsFAbN25M6jmRvO7u02pr69X4uEtZWZNqbKxUVdUmp8sC2DYdkjDYPR6PotFofDkajcrr9SZsc+nSJXk8HknSypUrJUnLly/Xtm3bFA6HCfYU6+4+rb1731IkciC+LhLZJ0m8geAotk3nJDwrprS0VIODgxoaGtLExIQ6OztVXV09pU11dbVef/11SVJ/f78eeOAB5eTkaGxsTNevX5ck3bhxQ729vVqzZs0cTSNztbX1TnnjSFIkckCHDr3tUEXAbWybzkm4x+5yudTe3q5gMKhYLKa6ujr5/X51dHRIkhoaGrR161adPHlSPp9PixYt0quvvipJGhkZUU1NjSRpcnJSO3bsUGVl5RxPJ/OMj0//Et68uTDNlQBTsW06J2GwS9KWLVu0ZcuWKesaGhqmLLe3t9/R79FHH9WHH354j+VhJllZk9Ouz86OpbkSYCq2Tefwy9N5rrGxUvn5+6asy89/QXv2VDhUEXAb26ZzCPZ5rqpqk1pbgwoG90uSgsH9am19ii+n4Di2TedY9ky/IJrLwS1rxh8wIXmWJaXzz8l483OsTBjPdDNlJ3vsAGAYgh0ADEOwA4BhCHYAMAzBDgCGIdgBwDAEOwAYhmAHAMMQ7ABgGIIdAAxDsAOAYQh2ADAMwQ4AhiHYAcAwBDsAGIZgBwDDEOwAYBiCHQAM43K6ABN1d59WW1uvxsddysqaVGNjJf/nEUgD3nu3Eewp1t19Wnv3vqVI5EB8XSRy+z+1Z+IGBqQL770vcCgmxdraeqdsWJIUiRzQoUNvO1QRkBl4732BYE+x8fHpPwTdvLkwzZUAmYX33hcI9hTLypqcdn12dizNlQCZhffeFwj2FGtsrFR+/r4p6/LzX9CePRUOVQRkBt57X7Bs27YdG9yylI7h0/1NeXf3aR069Lbeeuv/FAzu1549FWn58saypHS+mow3P8cyebx0v/ecOgtnxuy0Z3Dq1Cl71apVts/nsw8ePDhtmz179tg+n88uLi62z507l3TfJIa/ZydOvGvn579g396sbt/y81+wT5x4d87HTsP0GM+Q8Uyem6njOZstiSeY8NHJyUk7Pz/fvnjxoj0xMWGXlJTYFy5cmNKmu7vb3rJli23btt3f32+Xl5cn3TcdwV5ZuW/KH/6/t2Dwf+d8bBM3Zsab/2MxXmo4my2JJ5jwGHs4HJbP51NeXp7cbrdqa2vV1dU1pc3x48e1c+dOSVJ5eblGR0c1MjKSVN904JtyAHPhfs6WhD9QGh4eVm5ubnzZ6/VqYGBgxjbDw8O6fPnyjH0lqampKX4/EAgoEAjMdg4JTf2m3Irfe/ddybKaUzrWdCxr5jaMx3jpHovxUuXLGXL7mPdcnIUTCoUUCoWSbp8w2K0k/zL2PXwr8uVgnwuNjZWKRPb954cLt+vMz39Bra1PZdyv0QCkznS/dL19Fs5TKR/rqzu9zc2Jd0oTBrvH41E0Go0vR6NReb3ehG0uXbokr9erW7duzdg3Hf4b3ocO7dfNmwuVnR3Tnj2EOoB7c19nS6ID8Ldu3bIfffRR++LFi/b4+PiMX56+//778S9Pk+k7w/AAgGnMlJ0J99hdLpfa29sVDAYVi8VUV1cnv9+vjo4OSVJDQ4O2bt2qkydPyufzadGiRXr11VcT9gUAzK2M+IESAJhkpuzkkgIAYBiCHQAMQ7ADgGEIdgAwDMEOAIYh2AHAMAQ7ABiGYAcAwxDsAGCYjAr22Vz2cj5ifvOXyXOTmF+6EewGYX7zl8lzk5hfumVUsANAJiDYAcAwjl/dEQAwe4miO+H12Ocal+wFgNTjUAwAGIZgBwDDEOwAYJiMCfaenh4VFhaqoKBALS0tTpeTUtFoVN/5zne0evVqPf7442pra3O6pDkRi8W0bt06Pf30006XklKjo6Pavn27/H6/ioqK1N/f73RJKfWLX/xCq1ev1po1a/TDH/5Q4+PjTpd0T5577jnl5ORozZo18XX//Oc/VVFRoccee0yVlZUaHR11sMIMCfZYLKbdu3erp6dHFy5c0NGjR/XRRx85XVbKuN1u/epXv9Kf//xn9ff36ze/+Y1R8/uv1tZWFRUVGXc21d69e7V161Z99NFHOn/+vFH/9H1oaEi//e1vde7cOf3pT39SLBbTm2++6XRZ9+TZZ59VT0/PlHUHDx5URUWFPvnkEz355JM6ePCgQ9XdlhHBHg6H5fP5lJeXJ7fbrdraWnV1dTldVso8/PDDWrt2rSRp8eLF8vv9unz5ssNVpdalS5d08uRJ/fjHPzbqbKp//etfeu+99/Tcc89Jklwul5YsWeJwVanzjW98Q263W2NjY5qcnNTY2Jg8Ho/TZd2TjRs3aunSpVPWHT9+XDt37pQk7dy5U8eOHXOitLiMCPbh4WHl5ubGl71er4aHhx2saO4MDQ3pgw8+UHl5udOlpNTPfvYz/fKXv9SCBWZtshcvXtTy5cv17LPP6lvf+pbq6+s1NjbmdFkps2zZMj3//PP65je/qZUrV+qBBx7Q5s2bnS4r5a5cuaKcnBxJUk5Ojq5cueJoPWa9S+7CtI/ud/P5559r+/btam1t1eLFi50uJ2VOnDihhx56SOvWrTNqb12SJicnde7cOf3kJz/RuXPntGjRIsc/xqdSJBLRr3/9aw0NDeny5cv6/PPPdeTIEafLmlOWZTmeORkR7B6PR9FoNL4cjUbl9XodrCj1bt26pe9///v60Y9+pO9973tOl5NSZ8+e1fHjx/XII4/oBz/4gX73u9/pmWeecbqslPB6vfJ6vSorK5Mkbd++XefOnXO4qtT5wx/+oG9/+9t68MEH5XK5VFNTo7NnzzpdVsrl5ORoZGREkvS3v/1NDz30kKP1ZESwl5aWanBwUENDQ5qYmFBnZ6eqq6udLitlbNtWXV2dioqK9NOf/tTpclLupZdeUjQa1cWLF/Xmm2/qu9/9rl5//XWny0qJhx9+WLm5ufrkk08kSX19fVq9erXDVaVOYWGh+vv79e9//1u2bauvr09FRUVOl5Vy1dXVeu211yRJr732mvM7V3aGOHnypP3YY4/Z+fn59ksvveR0OSn13nvv2ZZl2SUlJfbatWvttWvX2qdOnXK6rDkRCoXsp59+2ukyUurDDz+0S0tL7eLiYnvbtm326Oio0yWlVEtLi11UVGQ//vjj9jPPPGNPTEw4XdI9qa2ttVesWGG73W7b6/Xahw8ftv/xj3/YTz75pF1QUGBXVFTY165dc7RGRy8CBgBIvYw4FAMAmYRgBwDDEOwAYBiCHQAMQ7ADgGEIdgAwzP8DmPscwfQ+c5AAAAAASUVORK5CYII=">

下面是投掷6次骰子,出现6点的概率分布。

In [18]:

n = 6
k = np.arange(n+1)
pdice = stats.binom.pmf(k, n, 1.0/6)
pdice
Out[18]: array([ 3.34897977e-01, 4.01877572e-01, 2.00938786e-01,
5.35836763e-02, 8.03755144e-03, 6.43004115e-04,
2.14334705e-05])In [19]: pl.stem(k, pdice, basefmt="k-")
pl.margins(0.1)

aaarticlea/png;base64,iVBORw0KGgoAAAANSUhEUgAAAW8AAAD5CAYAAADodLT+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAExFJREFUeJzt3V9M3Xf9x/HX6TiB2C51ztrZc0gwByynm6NMWtzFmrOpnCopsbMXzOmWji5EUyheLbRrpHPVodkFf4why1yydGG9mBlNGQfW1VOaLoUpNZrgn3oC8YB1yRrRrbOwc/r9XfRXVgY9HJDDl/fh+UhIzvd7Pj2fV/+9+PI5fz4ex3EcAQBMWeN2AADAwlHeAGAQ5Q0ABlHeAGAQ5Q0ABuUsxyQej2c5pgGArHOrFwQu25W34zjL9vWjH/1oWedb7bktZ7ea23J2q7ndyJ4KyyYAYBDlDQAGZWV5h0IhtyMsitXckt3sVnNLdrNbzS2trOweZ76FlaWYxOOZd/0GADBTqu7MyitvAMh2lDcAGER5A4BBlDcAGER5A4BBlDcAGER5A4BBlDcAGER5A4BBy/KRsEitu7tfra19mpzMUW5uQvX1Faqs3OF2LAAr2LxX3pFIRMXFxSoqKlJzc/Mtx73zzjvKycnRr3/96yUNmO26u/t14ECv+vqe1ZkzTerre1YHDvSqu7vf7WgAVrCU5Z1MJrV//35FIhENDw+rs7NTf/rTn+Yc99RTT2nnzp18hskCtbb2KRY7OuNcLHZUbW1vupQIgAUpy3twcFCFhYUqKCiQ1+tVdXW1urq6Zo1ra2vTnj17tGHDhowFzVaTk3OvXF29etsyJwFgSco17/HxceXn508f+/1+DQwMzBrT1dWl06dP65133rnllmdNTU3Tt0Oh0Ir6aEU35eYm5jyfl5dc5iQA3BaNRhWNRtMam7K809l7sqGhQc8999z0Rxfeatnk5vLGx+rrKxSLHZqxdBIIHFRd3U4XUwFwwycvbI8cOXLLsSnL2+fzKR6PTx/H43H5/f4ZY373u9+purpakvTee++pp6dHXq9XVVVVi8m+6tx4VUlb22H19v5Y4fBh1dXt5NUmAFJKuRlDIpHQ5s2b9dZbb2nTpk3avn27Ojs7FQwG5xy/d+9e7dq1Sw8//PDMSdiMIS0ej8QfE4AbUnVnyivvnJwctbe3KxwOK5lMqqamRsFgUB0dHZKk2trapU8LAJgX26CtIFx5A7gZ26ABQJahvAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAzKcTvAUuru7ldra58mJ3OUm5tQfX2FKit3uB0LAJZc1pR3d3e/DhzoVSx2dPpcLHZIkihwAFkna5ZNWlv7ZhS3JMViR9XW9qZLiQAgc7KmvCcn5/4h4urV25Y5CQBkXtaUd25uYs7zeXnJZU4CAJmXNeVdX1+hQODQjHOBwEHV1X3dpUQAkDlZU96VlTvU0hJWOHxYkhQOH1ZLy06erASQlTyO4zgZn8Tj0TJMc9N80jJOt2Ss5gaQGam6M2uuvAFgNaG8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADKK8AcCgecs7EomouLhYRUVFam5unnV/V1eXSkpKVFpaqi9/+cs6ffp0RoICAD6W8h2WyWRSmzdv1qlTp+Tz+bRt2zZ1dnYqGAxOj7ly5YrWrl0rSfrjH/+o3bt3629/+9vMSXiHZVqs5gaQGYt+h+Xg4KAKCwtVUFAgr9er6upqdXV1zRhzo7gl6YMPPtBnP/vZJYgMAEgl5U464+Pjys/Pnz72+/0aGBiYNe71119XY2OjLl26pL6+vjkfq6mpafp2KBRSKBRaXGIAyFLRaFTRaDStsSmXTV577TVFIhG98MILkqRjx45pYGBAbW1tc44/e/as9u3bp7/85S8zJ2HZJC1WcwPIjEUvm/h8PsXj8enjeDwuv99/y/EPPPCAEomELl++vMioAIB0pCzvsrIyXbx4UaOjo5qamtLx48dVVVU1Y0wsFpv+zjA0NCRJuvPOOzMUFwAgzbPmnZOTo/b2doXDYSWTSdXU1CgYDKqjo0OSVFtbq9dee00vv/yyvF6v1q1bp1dffXVZggPAasZmDCuI1dwAMoPNGAAgy1DeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABlHeAGAQ5Q0ABuW4HQC2dXf3q7W1T5OTOcrNTai+vkKVlTvcjgVkvXnLOxKJqKGhQclkUvv27dNTTz014/5XXnlFP/vZz+Q4jm6//Xb98pe/1L333puxwFg5urv7deBAr2Kxo9PnYrFDkkSBA5nmpJBIJJxAIOCMjIw4U1NTTklJiTM8PDxjzNtvv+1MTEw4juM4PT09Tnl5+azHmWeaJbfM0y0Za7krKg45kjPrKxx+2u1oQFZI1Z0p17wHBwdVWFiogoICeb1eVVdXq6ura8aY+++/X+vXr5cklZeXa2xsLFPfZ7DCTE7O/YPb1au3LXMSYPVJuWwyPj6u/Pz86WO/36+BgYFbjn/xxRf1zW9+c877mpqapm+HQiGFQqGFJcWKk5ubmPN8Xl5ymZMA2SEajSoajaY1NmV5ezyetCf9zW9+o1/96lc6d+7cnPffXN7IDvX1FYrFDs1Y8w4EDqqubqeLqQC7Pnlhe+TIkVuOTbls4vP5FI/Hp4/j8bj8fv+scX/4wx/05JNP6sSJE7rjjjsWERkWVVbuUEtLWOHwYUlSOHxYLS07ebISWAae/18Un1MikdDmzZv11ltvadOmTdq+fbs6OzsVDAanx/z973/XQw89pGPHjukrX/nK3JN4PEoxzZLzeK4/dWaN1dyS7ezASpWqO1Mum+Tk5Ki9vV3hcFjJZFI1NTUKBoPq6OiQJNXW1uqZZ57Rv/71L33/+9+XJHm9Xg0ODi7xbwEAcLOUV95LNglX3mmxmluynR1YqVJ1J2+PBwCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCDKG8AMIjyBgCD5i3vSCSi4uJiFRUVqbm5edb9f/7zn3X//fcrLy9Pzz//fEZCAgBmykl1ZzKZ1P79+3Xq1Cn5fD5t27ZNVVVVCgaD02PuvPNOtbW16fXXX894WADAdSmvvAcHB1VYWKiCggJ5vV5VV1erq6trxpgNGzaorKxMXq83o0EBAB9LeeU9Pj6u/Pz86WO/36+BgYFFTdTU1DR9OxQKKRQKLepxACBbRaNRRaPRtMamLG+Px7MUeSTNLG8AwGyfvLA9cuTILcemXDbx+XyKx+PTx/F4XH6//39PCAD4n6Qs77KyMl28eFGjo6OamprS8ePHVVVVNedYx3EyEhAAMJvHmad1e3p61NDQoGQyqZqaGjU2Nqqjo0OSVFtbq3/+85/atm2b/vOf/2jNmjW6/fbbNTw8rHXr1n08icezrOXu8UgWv5dYzS3Zzg6sVKm6c97yznSAzMxns0is5pZsZwdWqlTdyTssAcAgyhsADKK8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADEq5GQOQrbq7+9Xa2qfJyRzl5iZUX1+hysodbscC0kZ5Y9Xp7u7XgQO9isWOTp+LxQ5JEgUOM1g2warT2to3o7glKRY7qra2N11KBCwc5Y1VZ3Jy7h84r169bZmTAItHeWPVyc1NzHk+Ly+5zEmAxaO8serU11coEDg041wgcFB1dV93KRGwcJQ3Vp3Kyh1qaQkrHD4sSQqHD6ulZSdPVsIU9rBcQazmluxmt5obqwN7WAJAlqG8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADKK8AcAgyhsADGIPS8AYNk+GRHkDprB5Mm5g2QQwhM2TcQPlDRjC5sm4Yd7yjkQiKi4uVlFRkZqbm+ccU19fr6KiIpWUlOjChQtLHhLAdZY3T+7u7lc4/LRCoSaFw0+ru7vf7UhpW5HZnRQSiYQTCASckZERZ2pqyikpKXGGh4dnjOnu7na+8Y1vOI7jOOfPn3fKy8tnPc480yy5ZZ5uyVjN7Th2s1vLffLkGScQOOhc37zt+lcg0OicPHnG7WgpzZ374IrP7TjuZk/VnSmvvAcHB1VYWKiCggJ5vV5VV1erq6trxpgTJ07o8ccflySVl5drYmJC7777bqa+1wCrmtXNky2v1a/U7ClfbTI+Pq78/PzpY7/fr4GBgXnHjI2NaePGjTPGNTU1Td8OhUIKhUL/Q+z5eOTxZPDhM8hqbsludqu5pWfV2yv19j7rdpA0HZ11prdX8ngs5L85+/UNgTPxPEM0GlU0Gk1rbMry9qT5r9r5xO7Gc/26m8s70z6ZB4C7wuGn1dc3u6TD4cOKRH7sQqL03Sp7Jp5n+OSF7ZEjR245NuWyic/nUzwenz6Ox+Py+/0px4yNjcnn8y00M4AsVl9foUDg0IxzgcBB1dV93aVE6Vup2VNeeZeVlenixYsaHR3Vpk2bdPz4cXV2ds4YU1VVpfb2dlVXV+v8+fP69Kc/PWvJBMDqdmNNvq3tsK5evU15eUnV1a38tXpp5Wb3OPOsMfT09KihoUHJZFI1NTVqbGxUR0eHJKm2tlaStH//fkUiEa1du1YvvfSS7rvvvpmTeDwsZQDAAqXqznnLO9MBAABzS9WdvMMSAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAyivAHAIMobAAzKyvJO9yMVVxqruSW72a3mluxmt5pbWlnZKe8VxGpuyW52q7klu9mt5pZWVvasLG8AyHaUNwAYtGyfKggAWLhbVXTKzRgyPTkAYHFYNgEAgyhvADCI8gYAg7KuvCORiIqLi1VUVKTm5ma346TtiSee0MaNG/WlL33J7SgLEo/H9eCDD+ruu+/WPffco9bWVrcjpeXq1asqLy/X1q1btWXLFjU2NrodacGSyaRKS0u1a9cut6OkraCgQPfee69KS0u1fft2t+MsyMTEhPbs2aNgMKgtW7bo/Pnz7gZyskgikXACgYAzMjLiTE1NOSUlJc7w8LDbsdLS39/vDA0NOffcc4/bURbk0qVLzoULFxzHcZz333/f+eIXv2jmz/zKlSuO4zjORx995JSXlztnz551OdHCPP/88853vvMdZ9euXW5HSVtBQYFz+fJlt2MsymOPPea8+OKLjuNc/zczMTHhap6suvIeHBxUYWGhCgoK5PV6VV1dra6uLrdjpeWBBx7QHXfc4XaMBbvrrru0detWSdK6desUDAb1j3/8w+VU6fnUpz4lSZqamlIymdRnPvMZlxOlb2xsTG+88Yb27dtn7tVc1vJK0r///W+dPXtWTzzxhCQpJydH69evdzVTVpX3+Pi48vPzp4/9fr/Gx8ddTLS6jI6O6sKFCyovL3c7SlquXbumrVu3auPGjXrwwQe1ZcsWtyOl7Yc//KF+/vOfa80aW/+FPR6Pvva1r6msrEwvvPCC23HSNjIyog0bNmjv3r2677779OSTT+rDDz90NZOtv/l58GYg93zwwQfas2ePWlpatG7dOrfjpGXNmjX6/e9/r7GxMfX396+oz61I5eTJk/rc5z6n0tJSc1ex586d04ULF9TT06Nf/OIXOnv2rNuR0pJIJDQ0NKQf/OAHGhoa0tq1a/Xcc8+5mimrytvn8ykej08fx+Nx+f1+FxOtDh999JG+/e1v67vf/a6+9a1vuR1nwdavX6/Kykr99re/dTtKWt5++22dOHFCX/jCF/TII4/o9OnTeuyxx9yOlZbPf/7zkqQNGzZo9+7dGhwcdDlRevx+v/x+v7Zt2yZJ2rNnj4aGhlzNlFXlXVZWposXL2p0dFRTU1M6fvy4qqqq3I6V1RzHUU1NjbZs2aKGhga346Ttvffe08TEhCTpv//9r958802Vlpa6nCo9P/nJTxSPxzUyMqJXX31VDz30kF5++WW3Y83rww8/1Pvvvy9JunLlivr6+sy8uuquu+5Sfn6+/vrXv0qSTp06pbvvvtvVTMvy9vjlkpOTo/b2doXDYSWTSdXU1CgYDLodKy2PPPKIzpw5o8uXLys/P1/PPPOM9u7d63aseZ07d07Hjh2bfvmXJP30pz/Vzp07XU6W2qVLl/T444/r2rVrunbtmr73ve/pq1/9qtuxFsXKcuG7776r3bt3S7q+DPHoo4+qoqLC5VTpa2tr06OPPqqpqSkFAgG99NJLruZZlg+mAgAsraxaNgGA1YLyBgCDKG8AMIjyBgCDKG8AMIjyBgCD/g87Gs73YTxyPAAAAABJRU5ErkJggg==">

Beta分布

对于硬币或者骰子这样的简单实验,我们事先能很准确地掌握系统成功的概率。然而通常情况下,系统成功的概率是未知的。为了测试系统的成功概率$p$,我们做n次试验,统计成功的次数k,于是很直观地就可以计算出$p = k/n$。然而由于系统成功的概率是未知的,这个公式计算出的$p$只是系统成功概率的最佳估计。也就是说实际上$p$也可能为其它的值,只是为其它的值的概率较小。

例如有某种特殊的硬币,我们事先完全无法确定它出现正面的概率。然后抛10次硬币,出现5次正面,于是我们认为硬币出现正面的概率最可能是0.5。但是即使硬币出现正面的概率为0.4,也会出现抛10次出现5次正面的情况。因此我们并不能完全确定硬币出现正面的概率就是0.5,所以$p$也是一个随机变量,它符合Beta分布。

Beta分布是一个连续分布,由于它描述概率$p$的分布,因此其取值范围为0到1。 Beta分布有$\alpha$和$\beta$两个参数,其中$\alpha$为成功次数加1,$\beta$为失败次数加1。

连续分布用概率密度函数描述,下面绘制实验10次,成功4次和5次时,系统成功概率$p$的分布情况。可以看出$k=5$时,曲线的峰值在$p=0.5$处,而$k=4$时,曲线的峰值在$p=0.4$处。

In [20]:

n = 10
k = 5
p = np.linspace(0, 1, 100)
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="k=5", lw=2) k = 4
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="k=4", lw=2)
xlabel("$p$")
legend(loc="best");

aaarticlea/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXYAAAENCAYAAAAbu05nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4jFf7wPHvZLOTBA0i1iiJJYmdWmKrnSq/oqWqKEpVq95S7Ut3VbW0lsba2nlpUVupSu1iibWqobYEUUskQtY5vz9Ok1BEJDPPJJP7c11zmcgzz308knvOnHOe+5iUUgohhBB2w8HWDRBCCGFZktiFEMLOSGIXQgg7I4ldCCHsjCR2IYSwM5LYhRDCzqSb2OPi4qhXrx7+/v74+voyevTohx43bNgwKlWqhJ+fH6GhoVZpqBBCiIxxSu+befPmZdu2beTPn5+kpCQaNWrEzp07adSoUeoxGzZs4PTp04SFhbFv3z4GDx7M3r17rd5wIYQQD/fYoZj8+fMDkJCQQHJyMu7u7vd9f+3atfTp0weAevXqERUVRWRkpBWaKoQQIiPS7bEDmM1matasyZkzZxg8eDC+vr73fT8iIgIvL6/Ur0uXLk14eDgeHh6pf2cymSzYZCGEyD0yUxzgsT12BwcHDh8+THh4ONu3byc4OPixgR+WyJVS8lCKsWPH2rwN2eUh10KuhVyL9B+ZleFVMUWKFKF9+/YcOHDgvr/39PTk4sWLqV+Hh4fj6emZ6QYJIYTImnQT+7Vr14iKigLg7t27bNmyhYCAgPuO6dSpEwsWLABg7969uLq63jcMI4QQwljpjrFfvnyZPn36YDabMZvN9O7dmxYtWhAUFATAwIEDadeuHRs2bMDb25sCBQowf/58QxqeUwUGBtq6CdmGXIs0ci3SyLXIOpPKykBORoOYTFkaLxJCiNwos7nzsatihBDCWtzd3bl586atm2Fzbm5u3Lhxw2Lnkx67EMJmJDdoj7oOmb0+UitGCCHsjCR2IYSwM5LYhRDCzkhiF0KIhyhXrhxbt261dTMyRRK7EEI8hMlkylKdq8DAQPLly0ehQoUoVKgQPj4+Fmxd+iSxCyGEFZhMJqZPn05MTAwxMTGcPHnSsNiS2IUQ4jFOnjxJhQoVWL58+RO9zlZLOSWxCyFEOg4dOkSbNm2YNm0a3bt3p0OHDri5uT300alTp/teO3r0aIoXL06jRo347bffDGuz3KAkhLCZx+UGS23lkJn0U758eV5++WXmzZvH4sWLadKkyRO9PiQkhKpVq+Li4sLSpUsZOnQohw8fpkKFCg8ca+kblCSxC0PFJsSy++Jubty9QXxyPAnJCRTJU4T6pevjVcTr8ScQdiU7J/Zy5coRFxdHYGAgy5Yty3Ib2rZtS/v27Rk6dOgD37N0YpdaMcLqrsZeZcmxJWwI28D289uJT45/6HFehb1oXLYxr9V8jSZlm8jOWyJTCdlSTCYTQUFBjB8/nrfffptJkyYBOkHv3Lnzoa9p0qQJ69evN7KZDyU9dmE1sQmxTNoziQm7J3A74TYAJkzULlWbcq7lyOOUhzyOeYiIiWDPxT3cir+V+tpnvJ7hvcbv0da7rSR4O5adc0P58uWZO3cutWrVokWLFrRq1YrPP/88Q6+9desWe/fupWnTpjg5ObF8+XIGDhzI4cOH8fb2fuB46bGLbE8pxYIjCxi9dTSXb18GoK13W16q/hLPVnyW4gWKP/AaszJz4uoJfjj5A1+HfM2ui7tov6Q9jcs0ZmGXhZR1LWv0P0MIQO8et2XLFpo1a4aLiwsffvjhY1+TmJjIBx98wB9//IGjoyM+Pj6sWbPmoUndGqTHLiwqPime1ze8zrzQeQDULlWbL1t9SWC5wAyfIyY+hqCDQXy5+0uuxl6lSJ4iBHUIonu17hZpY3Q07N4NBw9CYmLa35cuDU2aQKVKlhvbFemT3KDJ5KnIti7FXKLriq7sDd9LPqd8TG83nT7+fXAwZW5V7bU71+i/tj9rTq0B4BX/V5jZfiZ5nfI+8blu3YIZM2DVKggNBbP50cd6eEDz5vDmm1CvXqaaLjJIcoMmiV1kS0cjj9JmURsu375MmSJlWN19NQElAx7/wsdQShF0MIi3fn6LuKQ4WlZoyeruqyngUiBDr792DaZMgWnTdHIHcHKC2rWhQQMoXFj/ndkMJ0/C9u1w9Wra61u2hPff1z156cVbnuQGTRK7yHbCrofRaH4jrsZepUnZJqz8v5UPHUfPiiNXjtB6UWsiYyN5xusZ1r+4niJ5izzyeKVgwQJ44w2IidF/16wZjBgBgYFQ4BHvC0rBqVPw/fcwfXraazt3hjlzoFgxi/6zcj3JDZokdpGtXLh1gcbzG3Ph1gVaVWjFTz1/Io9THqvE+vP6n7Rc0JKL0RepWbImP/f6mWL5H8y0UVEweDCkLD1+9ln473/hmWeeLN7Nm/DNN/DVV3pcvmRJWLgQWrSwwD9GAJIbUkhiF9lG5O1ImnzXhD+v/0lDr4Zs7rU5w0MkmXU+6jwtF7bk9I3T1PWsy7Y+28jvnD/1+wcOQLducP687pV/8w288krWhlHOn4devWDnTn2ed9+FTz8FBynIkWWSGzTZGk9kC/FJ8XRc2pE/r/+Jn4cf619cb/WkDlDWtSzbX9lOOddyhESE0PvH3piVngkNDtbDLefP6zH00FDo2zfrY+Nly8K2bfDRRzqZjx+vz5uUlPV/jxDWID12kSlDNwxl+v7plHMtx95+e/Eo6GFo/N///p2GcxtyK/4WIxuOJDBxAl27QlwcvPgizJ8PLi6Wj/vLL/DccxAbC88/D0uWQB7rjDzlCpIbNBmKETa3/PhyeqzqgYujC7te3UXtUrVt0o6tf22lzeI2JJmTcNwQRHLIawwcqJc1WnOYZM8eaNtWr7Jp3Rp++AHy53/868SDJDdoMhQjbOrUtVP0/6k/AJNbT7ZZUgdoUaEFwyp8C0By6yG89J8QZs60/th3gwZ62Kd4cfj5Zz3+npxs3ZjCeLI1nsgV7ibepdv/unE74TY9qvVgcO3BNm3P77/D3CH9YO8wcExij2dPYhKiDYnt76+Tu6sr/PijnlAV9iWrW+OlCAsLI2/evPTu3dsCrcoYSewiwz787UOOXz1O5aKVmdVhlk2Lc0VGQvv2ejjkuYIT8Pfw56+bfzFo3SDDPtr7+uo7WZ2c9JLIb781JKzIYYYMGULdunUN/X2RxC4y5NDlQ0zcPRETJr5/7nsK5Slks7bcvatvGDp3DurWhcXf52FZt2UUcC7A0uNL+e7wd4a1pXlzmDVLPx86VA/NCPuT2a3xli1bhpubGy1atDB0LiHdxH7x4kWaNWtG1apVqVatGl9//fUDxwQHB1OkSBECAgIICAjgk08+sVpjhW0kmZPov7Y/ySqZYfWGUa+0bQuovPUW7NunlyGuXasnLisXq8y0dtMAGLpxKH9e/9Ow9vTtC++9p8fZe/aEixcNCy0MkNmt8aKjoxk7diyTJ082fII43bK9zs7OTJ48GX9/f27fvk2tWrVo1aoVPj4+9x3XtGlT1q5da9WGCtv5avdXhF4JpWyRsnzS3LZv3KtWQVCQXsq4erUu2JWij18ftvy1hSXHltBvbT9+e+W3TBcge1IffwxHjsD69fDSS/Drr3qIRmSN6UPLDF+osZlLrL/99tsDW+OtW7cuQ6/94IMP6N+/P6VKlTJ82DLdn/oSJUrg7+8PQMGCBfHx8eHSpUsPHCfLlexX2PUwxv02DoBZHWdR0KWgzdpy/jz01wtymDhRT2Dey2Qy8U3bb/Ao4MHOCzuZdXCWYW1zcNBr50uWhB079J2pImdTShEUFMQzzzzzxPudHj58mK1btzJ8+PDUcxkpw32Kc+fOERoaSr1/1TE1mUzs3r0bPz8/PD09mThxIr6+vg+8fty4canPAwMDCQwMzHSjhXHe3PQmcUlxvOz3Ms9WfNZm7UhK0j3hqCjo2FGPZz+Mez53vmn7DS+sfIH/bPkPHZ7uQOnCpQ1pY/HiupZMq1b6LtVmzXRVSJF5me1pW0JWtsYLDg7m3LlzlClTBoDbt2+TnJzMyZMnOXDgwCNjBgcHExwcnPXGqwyIiYlRtWrVUj/++OMD34uOjlaxsbFKKaU2bNigKlWq9MAxGQwjspmfT/+sGIcq/HlhFXk70qZtGTtWKVCqVCml/v47/WPNZrPqvLSzYhyq45KOymw2G9LGFKNH67aWLq3UjRuGhs5xsnNuKFeunNq6dauKiopStWrVUqNGjcrwa+/cuaMiIyNVZGSkunLlinrnnXdUt27d1LVr1x56/KOuQ2avz2MHIBMTE+natSu9evXiueeee+D7hQoVIv8/t921bduWxMREbty4kfV3HGFTyeZkRmweAcCYxmN4qsBTNmvLsWN6aMNkgkWLHl8612QyMb3ddArnKcxPf/7EihMrjGnoPz78UG/QER4u69vtQcrWeBs3bmTs2LEZek2+fPl46qmneOqpp/Dw8KBgwYLky5ePokWLWrm1/0gv65vNZtW7d281fPjwRx5z5cqV1B7Rvn37VNmyZS32riNsZ9aBWYpxqHJTyqm7iXdt1o6kJKXq1tU94Ndff7LXBh0IUoxDlZxYUkXHRVungY9w4oRSzs663cHBhobOUSQ3aI+6Dpm9PunWitm5cydNmjShRo0aqbO6n332GRcuXABg4MCBTJ8+nZkzZ+Lk5ET+/PmZNGkS9evXv+88Ug8iZ4mJj8H7G2+uxl5lebflvFD1BZu15euv9RZ1np76TtOUHY8ywqzMNJjbgJCIEEY1GsXnLTK2w7yljBune+9PP61XzOR98h397J7kBk2KgAmrG/PrGD7b8RkNSjdg16u7bHaH6YUL+u7O2Fi9tLFz5yc/x77wfdSfWx8XRxdOvH4Cb3djdokHiI+HgAC95d6YMSC3eDxIcoMmRcCEVV25fYXJeyYDMKn1JJsldaX0LkixsdC1a+aSOkC90vV42e9lEpITeGfzO5Zt5GPkyQOzZ+vnX3yh5wqEMIIkdnGfCbsmcDfpLs9VeY76pes//gVWsnYtbNgARYroXZCyYnyL8RR0KciaU2vYcmaLZRqYQc88o9+gkpL0Ek3pnN7Pzc0ttdhWbn64ublZ9LrKUIxIdTnmMhW+rkBcUhyhA0PxL+H/+BdZQXw8VK0KZ87A1KkwbFjWzzl+53hGbx2NTzEfjgw6grOjc9ZPmkFRUeDtDdevw8qV+hOIEBkhQzEiy77Y9QVxSXF0qdLFZkkdYNo0ndSrVNG9XUt4q/5bVHSryMlrJ5l/eL5lTppBrq665ADAyJF6lychrEkSuwB0bz3oYBAAY5tmbK2uNVy9qu/aBJg0CZwt1LHO45SHT5vr+/w//O1D7iTescyJM2jAAKhWDc6ehSlTDA0tciFJ7AKA8bvGE5cUx/M+z+NXws9m7fjvfyE6Wm8917atZc/9f1X/j5ola3Ip5hLf7MviwP0TcnKCyXpOmk8/hStXDA0vchkZYxdcirlEhakViE+O58igI9TwqGGTdhw9qpcHmkx6Bcm/iohaxOYzm2m9qDWueV35a9hfuOWz7KTV43TurCeGX30V5s41NLTIgWSMXWTalL1TiE+O53mf522W1AFGjQKzGV5/3TpJHaBVhVY0L9+cqLgovtj1hXWCpGPiRD28NH8+nDhheHiRS0hiz+Wi46NTx9ZHPTPKZu3YsQM2boSCBfVwjLWYTKbUO1Cn7ptKRHSE9YI9RKVK8NpretnjBx8YGlrkIpLYc7lZB2cRHR9N07JNqeNZxyZtUErvQATw9tuPL/KVVXU969LVpytxSXF8ssP420HHjIF8+fQm2Pv3Gx5e5AKS2HOxhOQEpuzVSzRGNhxps3b8/DPs3Anu7jqxG+HjZh9jwsTcQ3O5eMvYvexKlkxbmz9mjKGhRS4hiT0XW3psKRExEfgW96VtJQsvQckgszmttz5qlL7T1Ag+xX3oXq07ieZEPt9pbHEwgP/8Rxc027IFtm0zPLywc5LYcymlFBP3TATgnQbvGLY36L/98AOEhupe7JAhxsb+oMkHutceanyv3d1d36wEutcui8aEJUliz6U2nd7E8avHKVmwJC9Wf9EmbUhOTpso/eAD+Ge/FsP4Fvele7XuJCQnMH7XeGODo8sRFy8Oe/boujhCWIok9lxq0l69f+Ob9d4kj1Mem7Thhx90SdsyZaBfP5s0IbXXPufQHMN77YUKpe2w9PHH0msXliOJPRc6+fdJfvnrF/I75+e1Wq/ZpA1mc1p98lGjwMXFJs3At7gvL1R9wWa99kGDoGhR2LcPtm41PLywU5LYc6Hp+6cD0LtGb8PvvEyxbp2+07RUKejb1yZNSHVvr/1yzGVDYxcokLYSSDbiEJYiiT2XiY6P5vsj3wMwpI7Bs5X/UCotiY0cafst46o+VZUuPl1ISE5g8t7JhscfOlRXgPztN32jlhBZJYk9l/n+8PfcTrhN07JNqe5R3SZt2LxZ35hTvLi+CzM7GN1oNAAzD8zk5t2bhsYuXDhtXbv02oUlSGLPRczKzLT90wB4o+4bNmmDUmm1yUeMMH4lzKPULlWblhVacjvhdupQlZHefFOXU9i8GUJCDA8v7Iwk9lzkl79+4c/rf1K6cGk6V8nkJqJZtGMH7NoFbm662Fd2ktJrn7pvquH12t3d09bxf/qpoaGFHZLEnotMC9G99UG1BuHk4GSTNkyYoP984w293C87aVauGXU963LtzjXmHJpjePy33tIbYK9dC3/8YXh4YUcksecS56LOse7Pdbg4ujCg1gCbtOHECVi/Xk+WDh1qkyaky2QypfbaJ+6eSEJygqHxPTygTx/9/KuvDA0t7Iwk9lxibuhcFIquPl15qsBTNmnDRF3BgFdf1ROn2VGnyp3wLe7LxeiLLD221PD4I0bojUYWLIDLxq68FHZEEnsukGROYu4hvV3PwFoDbdKGiAhYvBgcHIyr4JgZDiaH1EqXX+35yvCdv55+Gp57DhIS4Btjd+8TdkQSey6w/s/1XL59maeLPk2Tsk1s0oavv4bEROjaFSpWtEkTMqxntZ6ULFiSY1ePseWvLYbH/89/9J8zZ0JMjOHhhR2QxJ4LzDo0C4DXar6GyWQyPH50NHz7rX4+0nZl3zMsj1MehtXTC8sn7p5oePz69aFRI4iKgjnGz+EKOyCJ3c5duHWBTac34eLoQh//PjZpw6xZOrk3bQp1bLNJ0xMbWGsgBZwLsOWvLRy5csTw+ClvgJMn6086QjyJdBP7xYsXadasGVWrVqVatWp8/fXXDz1u2LBhVKpUCT8/P0JDQ63SUJE580LnYVZmnvd5nmL5rbzn3EMkJelhGMgZvfUUbvnc6FdTl5z8ao/xS1Q6dIAqVeDiRV0FU4gnkW5id3Z2ZvLkyZw4cYK9e/cyffp0Tp48ed8xGzZs4PTp04SFhTFr1iwGDx5s1QaLjEsyJzE3VE+aDqhpmyWOP/6ok9PTT0Nb22zSlGnD6w3HweTA0uNLCY8ONzS2g4O+GxVgyhRDQws7kG5iL1GiBP7+/gAULFgQHx8fLl26dN8xa9eupc8/i2/r1atHVFQUkZGRVmqueBKbTm8iPDocb3dvAssF2qQNKUnpzTd1sspJyruVp5tvN5LMSXy97+GfVq2pd299h+7evfohREZl+PbDc+fOERoaSr169e77+4iICLy8vFK/Ll26NOHh4Xh4eNx33Lhx41KfBwYGEhgYmLkWiwxL6a33D+hvk63v9u+H3bt15cKXXzY8vEWMaDCCFSdWMPvQbP7b9L8UdCloWOwCBXSRtC++gKlT9aSqsG/BwcEEBwdn/UQqA2JiYlStWrXUjz/++MD3OnTooHbu3Jn6dYsWLdTBgwfvOyaDYYQFRd6OVE4fOSnHDx3VpehLNmnDSy8pBUqNHGmT8BbTcG5DxTjUtH3TDI994YJSjo76cfGi4eGFjWU2dz62G5eYmEjXrl3p1asXzz333APf9/T05OLFtC3FwsPD8fT0zPo7jsiSxUcXk2ROoo13G0oWKml4/EuXYPlycHTMnuUDnsTwesMBXRzMrMyGxvbygm7d9P6w040vOilyqHQTu1KKfv364evry/Dhwx96TKdOnViwYAEAe/fuxdXV9YFhGGEspRTzD88HoK+/bbYnmjFDr4jp0kXvaZqTdfHpQpkiZQi7EcbGsI2Gx0+ZRA0KgjvGFp0UOVS6iX3Xrl0sWrSIbdu2ERAQQEBAABs3biQoKIigoCAA2rVrR4UKFfD29mbgwIHMmDHDkIaLRzt0+RDHrh6jaL6idKzc0fD4d++m3ZD0iP5AjuLk4JRav94WOyzVrw9168LNm7BwoeHhRQ5k+mccx7pBTCbDa27kZkM3DGX6/ukMqzeMqW2mGh7/u+/0PqY1a8KBA7qoVU4XFRdF6UmliU2M5eigo4bvPrVkCbz0ElSrpveKtYdrKh4vs7kzhy1AE48TlxTHkmNLANsMwygF03TZd954w34SkGteV/oG6Os5dZ/xb5bduumyvsePw/bthocXOYwkdjvz06mfuBl3E/8S/viX8Dc8fkgIHDwIRYtC9+6Gh7eqYXWHYcLEoqOL+Dv2b0Nju7ik7Q8rk6jicSSx2xlbT5qm9Nb794d8+WzSBKupVLQS7Z9uT3xyPLMPzTY8/sCBepXRDz/oMshCPIokdjty5fYVfj7zM84OzrxY/UXD40dGwooV+g7TQYMMD2+IlEnUGftnkJhsbHUuT094/nm99PGftQtCPJQkdjuy5NgSzMpMu0rtbFLwa84cvUFEx45Qrpzh4Q3RqkIrqhSrQkRMBKv/WG14/JR7AoKCID7e8PAih5DEbkcWHtVr4XrX6G147KQkvTEE5PwbktJjMplSe+1fhxhfP6ZxY6heHa5ehVWrDA8vcghJ7HbiWOQxDl85jGteVzo83cHw+GvW6HHfypWhRQvDwxvqZb+XKZynMDsv7OTQ5UOGxjaZ0t44U+YzhPg3Sex2IqW33r1qd/I45TE8fsp9aa+/bj9LHB+loEtB+gXoWu3fhBi/MelLL0GRIrBnDxw+bHh4kQNIYrcDyeZkFh9bDNhmGObUKfj1V8ifH/rYZpMmww2pMwQTJpYeW2r40scCBdKqZaYMfwlxL0nsdmDbuW1cirlEBbcKNPRqaHj8lPIBL76oe5K5QUX3ijZd+piy6mjxYr3toBD3ksRuB+6dNDV6s+o7d3QJAYDctnlWyiTqzAMzSTInGRrb11fvIRsbC4sWGRpa5ACS2HO42IRYVv2ul0f0qtHL8PjLl0NUlC5SVbOm4eFtqmWFljxd9GnCo8NZe2qt4fFT3khnztSlHIRIIYk9h1v9x2piE2NpULoB3u7ehsdPGePNbb11AAeTA0Pr6CUqtphE7dIlrX7Mrl2GhxfZmCT2HC5l0tQWvfWDB/X2d25u9lcXJqP6+PehoEtBgs8Fc/zqcUNju7hAP704RyZRxX0ksedgV2OvsvnMZpwcnHih6guGx09JJq+8Yn91YTKqcJ7CvOynl6hM3298da7XXtPLS1eu1DctCQGS2HO0FSdWkKySaV2xteElBG7dgqVL9XN7rQuTUUPqDAFgwZEFRMVFGRq7bFlo316XckiZxBZCEnsOljIM81L1lwyPvWiRXhHTvDk8/bTh4bMV3+K+NC/fnDuJd/ju8HeGx095Y501C8zGbskqsilJ7DnUmRtn2Bu+lwLOBehUuZOhsZVKW7s+cKChobOtlEnU6funG77hdZs2el/ZM2dg61ZDQ4tsShJ7DrX0uB4Hea7KcxRwKWBo7D179EqMp56C554zNHS21bFyR7wKe3H6xmm2nNliaGxHR13/HqScr9AksedASimbDsOkJI++ffXKDKE3vB5cW6/5nLbf+Opc/frpBL9mDVy5Ynh4kc1IYs+BQq+E8se1PyievzitKrYyNPaNG3ozDUjbqk1o/Wv2x8XRhfV/rufszbOGxi5VCjp10uWT580zNLTIhiSx50ApvfXu1brj5OBkaOwFCyAuDp59FipUMDR0tle8QHF6VOuBQjHzgPELy1PmO2bN0rssidxLEnsOk2xOZtnxZQC8WM3Y7e+UShuGkUnTh0tZ+jjn0BzuJN4xNHarVlC+PJw/D5s3GxpaZDOS2HOYHRd2cCnmEuVdy1O/dH1jY++AP/6AEiX09nfiQXU961KnVB1uxt1MfQM2ioMDDBign8skau4miT2HSVkN06NaD8MrOc6apf/s1w+cnQ0NnaMMrauXPk4LmYYyuDpX377g5ATr1ukdrUTuJIk9B0lITuB/J/4HwIvVjR2GuX5d37ZuMqUtrRMP90LVFyiWvxihV0LZE77H0NglSuglqMnJMomam0liz0E2n9nMzbibVHuqGtWeqmZo7AULID4eWreGcuUMDZ3j5HXKS/+a+t3PVvVjAObMkUnU3EoSew6y5NgSwDaTpinDMDJpmjGDag3CweTA/078j8jbkYbGbtFCr1i6cEEmUXOrdBP7q6++ioeHB9WrV3/o94ODgylSpAgBAQEEBATwySefWKWRQm+osebUGkCPrxspZdK0ZEldcEo8XlnXsnSq3IlEc6LhW+fJJKpIN7H37duXTZs2pXuCpk2bEhoaSmhoKO+//75FGyfSrD21ljuJd6hfuj7l3cobGlsmTTMnpX7Mtwe+JTE50dDYMomau6Wb2Bs3boybm1u6JzB61j+3SlkNY/QwzL2TpimbOoiMaV6+OVWKVSEiJiL105ZRPDxkEjU3y9JtiyaTid27d+Pn54enpycTJ07E19f3oceOGzcu9XlgYCCBgYFZCZ2r3Lh7g02nN+FgcjB8Q42USdM2bWTS9EmZTCaG1BnCGxvfYFrINLr5djM0/muv6TflOXPgvfd0LRmRvQUHBxMcHJz1E6nHOHv2rKpWrdpDvxcdHa1iY2OVUkpt2LBBVapU6aHHZSCMSMfsg7MV41AtF7Q0NK7ZrFSVKkqBUj/8YGhou3Er7pYq+FlBxTjUschjhsZOTlaqQgX9/7dhg6GhhYVkNndmaVVMoUKFyJ8/PwBt27YlMTGRGzduZP3dRtwnZRimZ7Wehsa9d9K0QwdDQ9tKqYK9AAAgAElEQVSNwnkK08evD6BvWDKSTKLmXllK7JGRkalj7CEhISilcHd3t0jDhHY55jLbzm7DxdGF532eNzS2TJpaRkr9mIVHF3Lz7k1DY8skau6UbmLv2bMnDRs25NSpU3h5eTFv3jyCgoII+uftf+XKlVSvXh1/f3+GDx/OsmXG1sbIDVacWIFC0da7La55XQ2LK5OmluNT3IeWFVraZOs8mUTNnUxKWX9Zi8lkktUzmVR/Tn32RexjWddldK/W3bC4kyfD22/rSdONGw0La7fW/LGG55Y/R0W3ivz5xp84mIy7N/CXX3TlxzJl4K+/ZBI1J8ls7pQ7T7OxszfPsi9iHwWcC9CxsnHlFO8tzyubaVhGh6c7UM61HGdunmFjmLHvlM2bp92J+vPPhoYWNiKJPRtLKfvauUpn8jvnNyzujh1w6pRMmlqSo4Mjr9d+HYBvQr4xNPa9k6gp8ybCvkliz8ZstRom5Zf/1Vdl0tSS+tXsR16nvPx85mdOXTtlaGyZRM1dJLFnUyeunuDY1WO45XXj2YrPGhb32jWZNLUW93zuqZuPG1318d5J1LlzDQ0tbEASezaV0lvv6tsVF0cXw+J+/31aed7yxpakyRXeqPsGAN8d/o7o+GhDY6dU5pRyvvZPEns2pJSySW0YKc9rfX4l/GhStgkxCTF8f/h7Q2M3bw4VK8LFi7LSyd5JYs+GQiJC+OvmX5QsWJImZZsYFjc4GP78E0qVkklTaxpWdxigJ1HNymxYXAeHtFVOcieqfZPEng3du6+po4Nxi45Tftn799cTbcI6OlfpTJkiZQi7Ecam0+mXxba0V17RE+IbNuieu7BPktizmWRzMstPLAeMXQ1z9Sr88IPu1cmeptbl5OCUWqt96r6phsZ+6ino2hXMZj3WLuyTJPZsJvhcMFduX6GiW0Vql6ptWNz58yExUe+Q5OVlWNhcq1/NfuRzysfmM5s5+fdJQ2PfO4malGRoaGEQSezZTOra9eo9MZlMhsQ0m2XS1Gju+dx52e9lwPgblpo2hcqV4dIlva5d2B9J7NlIfFI8K39fCRi7GmbrVl1DpEwZXRtGGGNYPT2J+v2R74mKizIsrskkk6j2ThJ7NrLp9CZuxd/Cz8MPn+I+hsWdOVP/OWCAFIgykm9x39Sqj3MOGTvg3acP5Mmja8f89ZehoYUBJLFnI7YoIRARAWvX6lUwMmlqvOH1hgN6OCbJbNyAd9Gi0L37/QXfhP2QxJ5NxMTHsPbUWkAvczTK7Nn6LsQuXaBECcPCin+0rdSWykUrc+HWBX44+YOhsQcP1n/Om6fvNhb2QxJ7NrH6j9XcTbpL4zKNKeta1pCYiYk6sQMMGmRISPEvDiYH3qz3JgCT9kwyNHa9euDnl1YfSNgPSezZxOJjiwF4sbpxk6br1umVEZUrQ7NmhoUV//Ky38u45XVjX8Q+9lzcY1hckymt1/7tt4aFFQaQxJ4NRN6O5Je/fsHJwYn/8/0/w+KmTJoOGqR/yYVtFHApwMDaep3p5L2TDY390ktQqBDs3AnHjhkaWliRJPZsYMWJFSSrZNp4t6Fo/qKGxAwLgy1bIF8+vUJC2NbQOkNxcnBi1clVnI86b1jcggWhd2/9XHrt9kMSezaQMgyTUqvbCCkrIXr0ADc3w8KKR/As7MkLVV/ArMyG37CUMr+ycCHExBgaWliJJHYbO3PjTOq+pp0qdzIk5p07aTvWy6Rp9vFW/bcAmHVwFrfibhkWt3p1aNRIJ/VFiwwLK6xIEruNLTm2BIAuPl0M29d02TK4eRPq1IG6dQ0JKTKgdqnaNC3blJiEGGYfmm1o7CFD9J/Tp+u17SJnk8RuQ0opw4dhlIJp0/TzoUMNCSmewMiGIwGYsncKCckJhsV9/nl9H8OJE/Dbb4aFFVYiid2GDl0+xKnrpyievzgtyrcwJObevRAaCsWKwQsvGBJSPIG2ldriU8yHiJgIVpxYYVhcF5e0+jHTjd2OVViBJHYbWnh0IaDvNHV2dDYkZkpvvX9/yJvXkJDiCTiYHBjRYAQAX+7+EmXguMjAgbq0xI8/Qni4YWGFFUhit5Ekc1JqbZjeNXobEjMyEv73P72ZhkyaZl+9avTCo4AHRyOP8stfvxgWt1QpXVoiOTmtjLPImSSx28iWM1u4GnuVykUrG7ahxuzZuoxAx45Q1piqBSIT8jjlSS3pO3HPRENjp8y7zJoFCcYN8QsLk8RuIynDML1r9DZkQ42kpLS16ykrIET2Nbj2YAo4F2Dzmc0cvnLYsLiNG0O1avrTndSPybnSTeyvvvoqHh4eVK9e/ZHHDBs2jEqVKuHn50doaKjFG2iPYuJjWP3HagBeqmHMapiUcdPKlaGFMfO0Igvc8rkxoNYAAMbvHG9YXJMprdf+jbH3SQkLSjex9+3bl02bHr2L+oYNGzh9+jRhYWHMmjWLwSkVhUS6Vp1cxd2kuzQp24RyruUMiTn1nz2T33hDj7GL7G9EgxE4Ozjzv9//x+kbpw2L26sXuLrqFVQhIYaFFRaU7q9448aNcUvnfvO1a9fS559CI/Xq1SMqKorIyEjLttAO3TsMY4SDB2HXLihSROrC5CSlC5emV41emJWZL3d/aVjcAgX0blqQ1iEQOYtTVl4cERGB1z1b2pcuXZrw8HA8PDweOHbcuHGpzwMDAwkMDMxK6BwrPDqcbWe3kccxD918uxkSM+WXs18/XfRJ5BzvPvMu3x3+ju8Of8fYpmMpVaiUIXGHDIGvvoIVK+DLL/WKGWF9wcHBBAcHZ/k8WUrswAPrbB81EXhvYs/NFh9djELRqXInXPO6Wj3elSu6hICDg9xpmhNVLlaZ532eZ9XJVUzZO4UJrSYYErdsWb30cdUqXd75448NCZvr/bvT++GHH2bqPFkabfX09OTixYupX4eHh+Pp6ZmVU9o1pRTfHfkOMG4Y5ttv9RLHTp2gfHlDQgoLG9VoFAAzD8zk5t2bhsV9U2/sRFAQxMUZFlZYQJYSe6dOnViwYAEAe/fuxdXV9aHDMEILiQjhj2t/4FHAg7aV2lo9Xnx82mYaKb+kIuepXao2LSu05HbCbaaFTDMsbqNGEBAAf/8NS5caFlZYQLqJvWfPnjRs2JBTp07h5eXFvHnzCAoKIuifBdHt2rWjQoUKeHt7M3DgQGbMmGFIo3Oq+YfnA9DbrzdODlkeBXus5cvh6lWoUQOaNrV6OGFFYxqPAfQOSzHxxhRNN5nSOgRTp0rVx5zEpAwoRmEymQyteZEd3U28S8mvSnIr/hbHBh+j2lPVrBpPKd3bOnIE5s6FV1+1ajhhZUopmn7XlB0XdvB5i89Th2esLT5ej7dHRsLWrdC8uSFhxT8ymztlRbNB1pxaw634W9QuVdvqSR30L+GRI+Dhofe1FDmbyWTigyYfAPDVnq+ITYg1JG6ePGmT7hONrW4gskASu0FShmFe8XvFkHhffaX/fOMN/cspcr6WFVpSv3R9rt25xrcHjNugdPBgvTfuxo26XrvI/iSxGyA8OpwtZ7bg4uhCz+o9rR7v+HHYtAny55cqjvbEZDLx3yb/BXRJ3zuJdwyJW7Qo9O2rn0+aZEhIkUWS2A2w8MhCFIrOlTvjns/d6vFSfvn69tW/lMJ+tPFuQ+1StYmMjWT2QeO2z3vrLT2ZumiRvjdCZG+S2K1MKZU2DOP/itXjXb6sf/lMJhg+3OrhhMHu7bWP3zXesF67tzd07qxL+U4zbsWlyCRJ7Fa2/fx2wm6EUapQKZ6t+KzV402bpm9I6tJF/zIK+9Ph6Q7UKlmLK7evMHP/TMPivvOO/nPmTIg1Zu5WZJIkditL2W2+r39fq69dj4mBlFsJRoywaihhQyaTiY+b6Xv8x+8ab9i69oYNoV49uHFDL6EV2Zckdiu6cfcGK3/XuxX0C+hn9XjffgtRUfqOwYYNrR5O2FAb7zY09GrItTvX+Hrf14bENJng3Xf184kTZYel7EwSuxUtOrqI+OR4WlVoRXk36xZqiYtLmzQdPdqqoUQ2YDKZ+KTZJ4BeIWNUDZnOncHHBy5ehCVLDAkpMkESu5UopVKHYQbUHGD1eN9/r1cr+PlBW+uXoRHZQLPyzWhevjm34m8xaa8x6xAdHNJ67ePH642vRfYjid1KQiJCOH71OMXzF6dzlc5WjZWUBBP+qeY6apT+yCxyh5Sx9il7p/B37N+GxHzxRShTBk6dgtWrDQkpnpAkditJ6a338e+Di6OLVWOtWAF//aVXwfzf/1k1lMhmGno1pF2ldtxOuM0nOz4xJKazM4wcqZ9//rkUB8uOJLFbQUx8DMuOLwOsP2mqlP5IDPCf/4Cjo1XDiWzo8xafY8LEzP0zOXPjjCExX30VihfX2y7+8oshIcUTkMRuBYuOLiI2MZbGZRpTpVgVq8b66Sc4dkxvXfbyy1YNJbKpGh416OPfh0RzImN+HWNIzPz5026A+8SYDwriCUhitzClFNP3TwdgcO3BVo4FKTsOjhwpxb5ys48CPyKvU16Wn1jO/oj9hsQcMgRcXWH7drDANp3CgiSxW9j289s58fcJPAp40NW3q1VjrV0LoaFQogQMHGjVUCKb8yrixbB6wwAYuWWkIfsfFCkCb7+tn48dK2Pt2YkkdgubcUDf+jmg1gCrTpre21sfNUqXVRW52+hGo3HP585v539jfdh6Q2IOGya99uxIErsFXY65zA8nf8DR5MjAWtbtQq9dC4cPQ8mS8NprVg0lcgjXvK6pW+iN3DKSxOREq8e8t9c+bpz02rMLSewWNOvgLJLMSXSu0pnShUtbLY701sWjDKkzBG93b/649kfqXI+1DRsGbm7Sa89OJLFbSGJyIkEH9SbfQ+oMsWqsNWvSeusDrH9Tq8hB8jjlYXLryQCMCx5nyE1LMtae/Uhit5DVf6zm8u3L+BTzoVm5ZlaLk5wMH+itL6W3Lh6qfaX2tK7Ymlvxt3h/2/uGxEzpte/YAZs3GxJSpEMSu4V8HaIr7L1e53VMVrynf/FivfVdmTIyti4ezmQyMbn1ZJwcnJh9cDaHrxy2eszChXVHA/SfZrPVQ4p0SGK3gJCIEHZe2EmRPEXo49fHanHi4+G/evMcPvoI8ua1WiiRw/kU92Fo3aEoFMM2DjNk+eMbb4Cnpx4mXLHC6uFEOiSxW8BXe74CYGDtgRTKU8hqcb79Fs6fh6pVoVcvq4URdmJs07EUz1+cHRd2sPDoQqvHy5cvbVL//felXrstSWLPonNR51j5+0qcHJx4o+4bVosTHZ126/Znn0lNGPF4rnldmfjsRABGbB7B9TvXrR7zlVegcmU4cwbmzLF6OPEIktizaOq+qZiVmR7Velh1ieNXX8G1a/DMM9Cxo9XCCDvTu0ZvmpVrxrU713j3l3etHs/JCT79VD//6CO4fdvqIcVDSGLPgqi4KOYc0t2St+u/bbU4ly/rxA66kqPUWxcZZTKZmNl+Ji6OLswNncvOCzutHvP556FuXYiMTPu5FcaSxJ4Fsw/O5nbCbZqVa0ZAyQCrxXnvPb0rfOfOej9TIZ5E5WKVefcZ3VsftG4QCcnWHfw2meDLL/XzL77Q2+gJY0liz6SE5ITUJY4jGoywWpz9++G778DFRW8gLERmvNf4PbzdvTnx9wm+2PmF1eM1aaI3fbl7N20ZpDDOYxP7pk2bqFKlCpUqVeKLLx78gQgODqZIkSIEBAQQEBDAJ7mkOPPCIwsJjw7Hp5gPbStZZ5NRpdJqXg8frndIEiIz8jrlZVaHWQB8vP1jjkYetXrMCRN0KeklS2DPHquHE/dS6UhKSlIVK1ZUZ8+eVQkJCcrPz0/9/vvv9x2zbds21bFjx/ROox4TJsdJTE5UFaZWUIxDLTqyyGpxlixRCpTy8FDq1i2rhRG5yOvrX1eMQwV8G6ASkhKsHm/MGP0zXKeOUsnJVg9ndzKbO9PtsYeEhODt7U25cuVwdnamR48erFmz5mFvDlZ628melh5byl83/8Lb3Zvu1bpbJUZsrN7qDvTyxsKFrRJG5DJftPyC8q7lCb0Syuc7P7d6vFGj9O5e+/fDQusvpRf/cErvmxEREXh5eaV+Xbp0afbt23ffMSaTid27d+Pn54enpycTJ07E19f3gXONS7lzAQgMDCQwMDBrLbeRZHMyn+7Q67nea/QeTg7pXsJM++wzCA+HmjX12mAhLKGgS0HmdZ5Hs++b8fH2j+lUuRP+JfytF6+gXsn18svw7rvQqZOuKSMeLjg4mGBLlMhMrzu/cuVK1b9//9SvFy5cqIYOHXrfMdHR0So2NlYppdSGDRtUpUqVLPZxIjtadmyZYhyq3JRyVvsoe/y4Uk5OSplMSu3ebZUQIpcbsn6IYhyq+ozq6m7iXavGSk5WqlEjPSQzcKBVQ9mdzObOdIdiPD09uXjPWqWLFy9SuvT9N+EUKlSI/PnzA9C2bVsSExO5ceNG1t9xsiGzMvPJDj05PLrRaJwdnS0fw6yLeyUlwaBB0KCBxUMIwfiW46nkXoljV48xcstIq8ZycICgIHB21n/u3m3VcILHrIqpXbs2YWFhnDt3joSEBJYvX06nTp3uOyYyMjJ1jD0kJASlFO7u7tZrsQ2t/mM1x68ep3Th0lYr9jVnjv7BL1FCD8cIYQ0FXQqytOtSnB2cmRYyjbWn1lo1nq9v2pzRwIGQaP3NnXK1dBO7k5MT06ZNo3Xr1vj6+tK9e3d8fHwICgoiKEhvKrFy5UqqV6+Ov78/w4cPZ9myZYY03GhJ5iTG/Kq3HRvdaDR5nPJYPMaVK3ocEmDqVL2XpBDWUqtULca3HA9A3zV9iYiOsGq8MWP0kt3jx+WeDGszKWX9JS0mkynHr5yZfXA2r617jYpuFfl9yO9W2ai6e3dd7rRtW1i/XkoHCOszKzPtl7Rn0+lNBJYL5Jfev+DoYL0Kc7/8Aq1a6ZLTR47A009bLZRdyGzulDtPM+BO4h3GBo8F4NPmn1olqS9bppN6gQIwfbokdWEMB5MD33X+Do8CHgSfC079VGotLVvqFTJxcdCnj55LEpYniT0DpuydwuXbl6lVshb/V/X/LH7+S5fg9df180mToHx5i4cQ4pE8CnqwrNsyHE2OfLHrC1b+vtKq8aZO1Rty7N2r704VlieJ/TGu3bnGF7t0KYUvWn6Bg8myl0wp6NcPbt7UQzCyObWwhcBygXzZSlfuemX1K5y4esJqsVxdYf58/XzsWAgNtVqoXEsS+2N8uuNTouOjaV2xNS0qtLD4+WfNgk2bwN0d5s6VIRhhO8PrD+fF6i8SmxhLl+VdiIqLslqsVq1gyBA9FPPyy3rbR2E5MnmajhNXT+Af5E+yOZmDrx20eGneP/6AWrXgzh1YvhxeeMGipxfiid1JvEODuQ04GnmUZys+y7qe66xyvwboshkBARAWBm+9pYchxf1k8tTClFK8vuF1ksxJvFbrNYsn9dhY6NZNJ/WXXpKkLrKH/M75Wd19NcXzF2fzmc0MWj/Iap2yAgV0/RgnJ5g8GX780SphciVJ7I+w+Nhitp/fTrH8xfishWXvFFJKT5aeOAFVquhNqoXILsq7lWfdi+vI55SPeaHzUmsjWUO9emkTqK+8ovdKFVknif0houKieGfzOwBMaDkB93yWvZN27lxYsADy54eVK3WhJCGyk7qedVnadSkmTHyw7QMWHFlgtVjDh+vt9KKj9eYccXFWC5VrSGJ/iP9u+y+RsZE09GpIH3/Llg4IDYWhQ/Xzb7+FqlUtenohLKZzlc583VbvEvbqmlf58aR1xkpMJpg3DypW1L8fw4bpT7Ui8ySx/0tIRAjT90/HweTAjHYzLLq88dIlXbY0Pl4X+urd22KnFsIqhtYdynuN3yNZJdN9ZXfW/bnOKnGKFNGfXvPmhdmz4euvrRIm15DEfo/YhFh6/dALszLzdoO38SvhZ7lzx0LHjrrGeqNG8oMrco5Pmn3C2w3eJtGcSNcVXdl8ZrNV4vj765476FUyP/1klTC5giT2e/znl/8QdiOMak9V4+NmH1vsvMnJ8OKLcOiQ/rj54496L0ghcgKTycTEVhMZUmcICckJdF7WmU2nN1klVs+e8NFHeiimRw/9OyOenCT2f2wM28iM/TNwdnBmYZeF5HXKa5HzKgVvvw1r1+qdYzZsgGLFLHJqIQxjMpn4uu3XDKg5gLikODou7cjSY0utEuv99/VNS3fu6E+5Fy5YJYxdk8QOXL9znVfXvgrAx80+tuhWYR98oIddnJ11T12q2YmcysHkwLcdvmVEgxEkmZN46YeXmLF/hsXjmEz6juymTfW8VPPmEGHdisJ2J9cn9mRzMn1W9+HK7Ss0KtOIdxq+Y7Fzf/wxfPopODrC0qX6B1WInMzB5MCXrb5kfIvxKBRDNgzh/V/fx6zMFo2TJw+sXq3vzD5zBlq00PsViIzJ9SUF3tv6Hp/v/Bz3fO4cfO0g5VzLWeS8EyboTTMcHGDxYj1eKIQ9mXNoDgPXDcSszHSu3JmFXRZSKE8hi8a4cQOaNYOjR/XS4G3boHhxi4bI1jKbO3N1Yl9+fDk9VvXA0eTIz71+tkiRL6Xgww/1w2SC777T44VC2KOfT/9Mj1U9iIqLomrxqqzpsYaK7hUtGuPvvyEwEH7/XW+xt3EjlClj0RDZltSKeUKHrxym75q+AExqPckiST0pSa9P//BD3VOfM0eSurBvrb1bE9I/BJ9iPpz4+wS1Z9e2eD334sVh61bdY//9d2jYEI4ds2gIu5Mre+znos7ReH5jwqPD6evfl7md5mLKYr3c2Fg93LJuHeTLp3dE+te+30LYrej4aHr/2Dt1U+x+Af2Y0mYKBV0sVy/j5k147jnYvh0KF4Y1a3RP3p7JUEwGXbx1kabfNeVs1Fme8XqGrS9vzfLG1GfO6EqNhw/ruurr1kGDBhZqsBA5hFKKGftnMGLzCOKT46nkXon5nefzTJlnLBYjLg569YJVq/RKs0mTdF13e93HQBJ7BlyKuUTgd4GE3QijTqk6bOm9hSJ5i2TpnD/+qKvSRUfrHdh/+klXbBQitzp+9Tg9V/Xk+NXjAAyoOYDxLcdbrJhecjKMHKlL/YL+pDxrFhSy7LxttiCJ/TEioiNoubAlf1z7g4ASAWx9eStu+dwyfb67d2HMmLQfrq5dddXGIll7nxDCLsQlxfHZjs8Yv3M8ieZEiucvzoRWE+hdozeODo4WibFihd5W8vZtqFxZLykOsOy2CTYniT0dBy4doNPSTly+fZnqT1Xn1z6/Uix/5m//DA7We5OePq03CfjyS3jzTfv9OChEZp38+ySD1g9i+/ntAFR7qhqfNf+MDk93yPK8FsCpU3oY9Phxfb/IO+/ofVTz5cvyqbMFSeyPsOLECvqs7kNcUhxNyjZh1QurMp3Ur1+H0aN19TmAatV0L71uXQs2WAg7o5Ri0dFFfLDtA87fOg9AQ6+GjG40mnaV2mW5guqdO/rT89SpermxtzfMnAktW1qi9bYlif1f4pLiGBs8lgm79PYs/QL6MaP9DFwcXZ74XLdvw5QpumceHa0nbd5/H0aNApcnP50QuVJ8UjwzD8zk0x2fcu3ONQB8ivnwTsN3eLH6i1muz7Rvnx6aOXFCf/3ss/DZZ/ru1ZxKEvs9dl3YRb+1/Th1/RQOJgcmtprI8PrDn/ijX1SULiM6YQJERuq/a91az8T7+lqh4ULkAtHx0cw+OJsp+6YQHh0OgGteV16s/iKv+r9KzZI1Mz1Mk5AAX30F48frThjoXZlGjoQ6dSz1LzCOJHbgcsxlPtnxCTP3z0Sh8Cnmw9xOc2ng9WRrD0+ehGnT4Pvv9fp00Hszfv65vr1ZCJF1CckJLD++nKn7pnLw8sHUv/ct7svzPs/zfJXn8S/hn6kkf/26Tu7ffKM3tgG9BHnYML0NX075pJ2rE3t4dDhf7PqC2QdnE58cj6PJkVGNRvF+k/cz/PHu3Dk9y758+f01oFu00BOjHTpYZnI0ODiYQHu/qyKD5Fqkye3X4siVI8w/PJ+FRxdy4+QNKKf/vmyRsrSq2IoW5VvQrFwzPAp6PNF5w8P12PucOfoTOICrq07uPXrojpqTk2X/LZaU6dypHmPjxo2qcuXKytvbW40fP/6hx7zxxhvK29tb1ahRQx06dOiB72cgzBOLTYhVy48vV12WdVHOHzkrxqEYh+qyrIs6cuXIY19/86ZSa9cq9eabSlWrppSedtGPwoWVeu01pY4ds3iz1dixYy1/0hxKrkUauRZafFK86jWslxr400BVYmKJ1N/rlEelryupl1a9pKbsmaJ2nN+hbty5kaHz3r6t1LffKlWjxv2/666uSnXpotS0aUr9/rtSyclW/gc+oczmznTfq5KTkxk6dCi//PILnp6e1KlTh06dOuHj45N6zIYNGzh9+jRhYWHs27ePwYMHs3fv3id/h3mM+KR4Dl0+xK6Lu9h1cRdbzmwhNlGPk5gw0b1qd8Y0HkN1j+r3ve72bX1n6OnTemnU4cO6R37mzP3nL1BAF/Xv3h3atNF7LwohjOXi6EJFt4qM6zCOGe1ncOjyIbb+tZWtZ7ey88JOwm6EEXYjjMXHFqe+pmTBkvgW96Wie0XKu5anvGt5vIp4UbJgSUoWKklep7wUKAADB+rHyZP6k/myZTon/PijfoAuVRAQADVr6nm0SpX0KpuSJXX9p5wi3cQeEhKCt7c35cqVA6BHjx6sWbPmvsS+du1a+vTpA0C9evWIiooiMjISD4/7PzIt3Hog9blCoZQZM2bMZjMJ5ngSzfHEJ8cTmxRDdMINbiXe4EZ8JBF3/iLizhmuxJ3HTPJ95yzvUpeazj2pkvwC5kOlmPmr/rh19SpcvqwfN28+/N+WJ4+eLW/RQi+Lql8/54y7CZEbOJgcqF2qNrVL1ebdRu+SkJzA8avHCYkIISQihCORRzj590ku377M5duX2Xp26+0IBI0AAAYySURBVEPPUzhPYYrmK0rR/EUpmq8oRfIWoVDNQrRtUIhWdwoScS4f58Ly8def+bh1PQ+/3XDht59dYKMzmJ0g2RlnJ0fc3Rxwd3WkaFEHihRyoFAhE4UKOlAgv4m8eSFvXhN58oCzkwlnZz3E4+howsGBtIfJBCY9rHvv0K7Fb4FJrzv/v//9T/Xv3z/164ULF6qhQ4fed0yHDh3Url27Ur9u0aKFOnDgwAMfJ+QhD3nIQx5P/rD4UExGZ6PVvwb3//26f39fCCGE9aQ7auTp6cnFixdTv7548SKlS5dO95jw8HA8PT0t3EwhhBAZlW5ir127NmFhYZw7d46EhASWL19Op38VGe/UqRMLFiwAYO/evbi6uj4wvi6EEMI46Q7FODk5MW3aNFq3bk1ycjL9+vXDx8eHoKAgAAYOHEi7du3YsGED3t7eFChQgPnz5xvScCGEEI+QqZH5R7DEmnd78bhrsWjRIlWjRg1VvXp11bBhQ3XkyOPX3udUGfm5UEqpkJAQ5ejoqFatWmVg64yVkWuxbds25e/vr6pWraqaNm1qbAMN9Lhr8ffff6vWrVsrPz8/VbVqVTV//nzjG2mAvn37qqeeekpVq1btkcc8ad60WGJPSkpSFStWVGfPnlUJCQnKz89P/f777/cds379etW2bVullFJ79+5V9erVs1T4bCUj12L37t0qKipKKaV/wHPztUg5rlmzZqp9+/Zq5cqVNmip9WXkWty8eVP5+vqqixcvKqV0crNHGbkWY8eOVaNGjVJK6evg7u6uEhMTbdFcq9q+fbs6dOjQIxN7ZvKmxZbc37vm3dnZOXXN+70etebd3mTkWjRo0IAi/+zKUa9ePcLDw23RVKvLyLUA+Oabb+jWrRvFixe3QSuNkZFrsWTJErp27Zq6SKFYsczvG5CdZeRalCxZkuh/KnlFR0dTtGhRnLLz/f+Z1LhxY9zcHr3pT2bypsUSe0REBF5eXqlfly5dmoiIiMceY48JLSPX4l5z586lXbt2RjTNcBn9uVizZg2DBw8GMr7MNqfJyLUICwvjxo0bNGvWjNq1a7Nw4UKjm2mIjFyLAQMGcOLECUqVKoWfnx9Tp041upnZQmbypsXe/iy15t0ePMm/adu2bcybN49du3ZZsUW2k5FrMXz4cMaPH59a8OjfPyP2IiPXIjExkUOHDrF161bu3LlDgwYNqF+/PpUqVTKghcbJyLX47LPP8Pf3Jzg4mDNnztCqVSuOHDlCIXvc3PQxnjRvWiyxy5r3NBm5FgBHjx5lwIABbNq0Kd2PYjlZRq7FwYMH6dGjBwDXrl1j48aNODs7P7C0NqfLyLXw8vKiWLFi5MuXj3z58tGkSROOHDlid4k9I9di9+7djBkzBoCKFStSvnx5Tp06Re3atQ1tq61lKm9aagIgMTFRVahQQZ09e1bFx8c/dvJ0z549djthmJFrcf78eVWxYkW1Z88eG7XSGBm5Fvd65ZVX7HZVTEauxcmTJ1WLFi1UUlKSio2NVdWqVVMnTpywUYutJyPX4q233lLjxo1TSil15coV5enpqa5fv26L5lrd2bNnMzR5mtG8abEeu6x5T5ORa/HRRx9x8+bN1HFlZ2dnQkJCbNlsq8jItcgtMnItqlSpQps2bahRowYODg4MGDAAXzvcrisj1+K9996jb9+++Pn5YTabmTBhAu7u7jZuueX17NmT3377jWvXruHl5cWHH35IYmIikPm8achGG0IIIYyTgyoMCyGEyAhJ7EIIYWcksQshhJ2RxC6EEHZGErsQQtgZSexCCGFn7K+ijhAZcOjQIdasWYOXlxclSpTg1KlTjBgxwtbNEsIipMcucqW7d+9SqFAhSpUqRYcOHdiwYYOtmySExUhiF7nSM888w759+2jSpAlKKa5cuWLrJglhMZLYRa51/fp1ChYsyK+//krnzp1t3RwhLEbG2EWudObMGZKSkvjpp5/Yv38/H374oa2bJITFSK0YkSstXLgQk8lEr169bN0UISxOhmJErnP58mXmzJmT7q5WQuRk0mMXQgg7Iz12IYSwM5LYhRDCzkhiF0IIOyOJXQgh7IwkdiGEsDOS2IUQws5IYhdCCDsjiV0IIezM/wPbtBCjjLlbtgAAAABJRU5ErkJggg==">

下面绘制$n=10, k=4$和$n=20, k=8$的概率分布。可以看出峰值都在$p=0.4$处,但是$n=20$的山峰更陡峭。也就是说随着实验次数的增加,$p$取其它值的可能就越小,对$p$的估计就更有信心,因此山峰也就更陡峭了。

In [30]:

n = 10
k = 4
p = np.linspace(0, 1, 100)
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="n=10", lw=2) n = 20
k = 8
pbeta = stats.beta.pdf(p, k+1, n-k+1)
plot(p, pbeta, label="n=20", lw=2)
xlabel("$p$")
legend(loc="best");

aaarticlea/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXcAAAENCAYAAAD0eSVZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVFX/wPHPsIksAu6KuIImLoCaqKWiZiYqVlZq+bj0aGbl0mJl9Sssn7IySyt7sHJrca3cQs1UskfFlbRcUVEBFTVFxI0Bzu+PEyiKMMAsMHzfr9e8nJl77r3fbsx37px77vcYlFIKIYQQdsXB1gEIIYQwP0nuQghhhyS5CyGEHZLkLoQQdkiSuxBC2CFJ7kIIYYdMSu5ZWVmEhITQp0+ffJePGTOGgIAAgoKCiIuLM2uAQgghis6k5D5t2jQCAwMxGAy3LYuOjubw4cPEx8czc+ZMRo0aZfYghRBCFE2hyT0pKYno6GiGDx9Ofvc7LV++nCFDhgAQGhpKamoqKSkp5o9UCCGEyZwKa/D888/z4YcfkpaWlu/y5ORk/Pz8cl/XqVOHpKQkatSokftefmf8QgghClfcIgIFnrmvXLmS6tWrExISUuAObl2WXzJXSslDKd566y2bx1BaHnIs5FjIsSj4URIFJvfNmzezfPlyGjRowMCBA1m/fj2DBw/O08bX15fExMTc10lJSfj6+pYoKCGEECVTYHJ/9913SUxMJCEhgQULFtC1a1fmzZuXp01ERETue7GxsXh7e+fpkhFCCGF9hfa53yynuyUqKgqAkSNHEh4eTnR0NP7+/ri7uzN79mzzR2lHwsLCbB1CqSHH4gY5FjfIsTAPgyppx44pOzEYStx/JIQQ5U1JcmeRztyFEKIoKleuzIULF2wdRqnn4+PD+fPnzbpNOXMXQliMfPZNc6fjVJLjJ7VlhBDCDklyF0IIOyTJXQgh7JAkdyGEsEOS3IVVZatsluxbwtu/vc34teMZuXIkY1aNYU/KHluHJkSRGI1GHnnkERo0aICDgwO//fbbbW1eeeUVqlatStWqVXn11VetGp8MhRRWE5sUy+hVo9lxcsdtyz7d9im9Anox4d4J3FP3HhtEJ0TRderUieeff55HH330tppaUVFRLFu2jD179IlL9+7dadCgASNHjrRKbDIUUljchasXeH7N88zdPReA2p61GRw0GB9XHypVqMTes3v5etfXXM28CsCTIU8ys/dMHB0cbRm2MIPS/NmvX78+o0ePZt68eRw/fpwHHniAuXPnUqFChSJvy8/Pj++++45OnTrlvtehQweefPJJhg8fDsDs2bOZOXMmW7ZsuW19SwyFlDN3YVFXjFfo9X0vtiRtwcXRhZc6vMSEeyfg4eKRp92bnd5k+rbpTN0ylVlxs7iWeY25D87FyUH+RIVlGAwGFi9ezJo1a6hQoQL33HMPc+bMITw8nBYtWtyxVPkXX3zBgAEDCt3+vn37CAoKyn3dsmVL9u7da7b4CyOfHGExmdmZDFgygC1JW/Cr5Me6wesIqBKQb9tq7tV4p8s73N/wfsK/D+f7P7/HmGXku4e/w9nR2cqRC2sx11QPxf1xMGbMGGrWrAlAnz59+OOPPxg5ciSpqakljik9PR0vL6/c15UqVSI9Pb3E2zWVXFAVFqGUYtTPo1hxaAU+rj6sGbTmjon9Zh3rdeSXQb9QqUIlFu9bzIAfBpCtsq0QsSiPchI7QMWKFc2afD08PPJMcnTx4kU8PDwKWMO8JLkLi5i0cRJf7foKVydXVj6+kqbVmpq8bnu/9qwbvA5vV29+3P8jn237zIKRCltSyjwPc0pMTMTDwwNPT898H/PnzzdpO82aNeOPP/7Ifb17926aN29u3mALIMldmF3cqTgif4vEweDAwkcW0sGvQ5G30aZ2G+b0nQPAK7++wv6z+80cpRD58/PzIz09nUuXLuX7GDhwYG7b69evc+3atdueAwwePJipU6dy8uRJkpOTmTp1KkOHDrXaf4ckd2FWWdlZPP3z02SrbEa3HU1Ek4hib6vvXX0ZGjyUa5nX+NdP/8KYZTRjpELkZTAYijzfc5MmTXBzc+PkyZP06NEDd3d3Tpw4Aej5Lvr06UOLFi1o2bIlffr04amnnrJE6PmSoZDCrP6747+M+nkUtT1rs//Z/VSqUKlE20u7nkbLL1py/OJx3ur8FpFhkeYJVFiFfPZNI1UhRamWkp7ChHUTAJj2wLQSJ3aAShUqMefBORgwMGnjJHae3FnibQpRHkhyF2bz0tqXSL2WSk//nvRr2s9s2w2rH8bYdmPJUlmMXztezgSFMIF0ywiz2JK4hQ6zOuDq5MreZ/bS0KehWbefei2VhtMacuHaBX4Z9AvdG3U36/aFZchn3zRW75a5du0aoaGhBAcHExgYyIQJE25rExMTg5eXFyEhIYSEhDBp0qRiBSLKtnf/9y4AL7R/weyJHcDb1ZtX7nkFgAnrJsjYdyEKUeiZ+5UrV3BzcyMzM5N7772XKVOmcO+99+Yuj4mJYerUqSxfvvzOO5Fvb7u2+/RugqOCqehUkePjjlPNvZpF9nPFeAX/6f6cSj/FokcW8WizRy2yH2E+8tk3jU0uqLq5uQGQkZFBVlYWlStXvq2N/M8r3yZvmgzAU62fslhiB3BzduOtzm8B8Pr612VopBAFKLS2THZ2Nq1ateLIkSOMGjWKwMDAPMsNBgObN28mKCgIX19fpkyZclsbgMjIyNznYWFhhIWFlTh4YXuHzx9m0d5FODs482L7Fy2+vydDnmTKlinEn49nzh9zGNF6hMX3KYS1xMTEEBMTY5ZtmXxB9eLFi/To0YPJkyfnScyXLl3C0dERNzc3Vq1axdixYzl06FDenchPM7v11Iqn+HLXlzwZ8iRfR3xtlX0u+GsBA38YSF2vuhwZc0QqR5Zi8tk3jU3HuXt5edGrVy927Mg70YKnp2du103Pnj0xGo2cP3++WMGIsiU5LZk5f+gx6DkXO63hsWaP4V/ZnxMXT7D0wFKr7VeIsqTA5H7u3Lnc0pdXr15l7dq1hISE5GmTkpKS+82ybds2lFL59ssL+/PJ1k8wZht5tNmjNK7S2Gr7dTA4MDZ0LAAfx35stf0KcbPY2Fi6d+9OlSpVqF69Oo899hinT5/O08aW0+wVmNxPnTpF165dCQ4OJjQ0lD59+tCtWzeioqKIiooCYMmSJbRo0YLg4GDGjRvHggULrBK4sK3rmdeZHTcbwCp97bcaGjwUrwpebE7czLbkbVbfvxCpqak8/fTTHD9+nOPHj+Pp6cmwYcNyl988zd6ePXtYsWJFbt60CmUFVtqNsKIFfy5QRKKCvghS2dnZNonhpV9eUkSiBi4ZaJP9i8KV5s9+vXr11JQpU1TLli2Vl5eX6t+/v7p27Vqxt7dz507l6emZ+7p9+/bqyy+/zH09a9Ys1a5du3zXvdNxKsnxk/IDoli+3PUlACNajShyJT1zGd12NI4GRxbvW0xSWpJNYhBl183T7CUkJLBnzx7mzJlDYmIi3t7e+Pj45Pu4U+/Exo0b89Rrl2n2RJlz9MJR1iWsw9XJlcdbPG6zOOp61aVfYD8W7V3E59s/571u79ksFlE8honmOTFQbxVvRIm5ptnbs2cP77zzTp6bOWWaPVHmfB2nhzw+GvgoPhV9bBrLuNBxAETtiOKK8YpNYxFljzmm2Tt8+DDh4eFMnz6de+65J/d9W0+zJ2fuokgyszNzL6SOaGX7G4ja+7WnrW9btiVvY8m+JQwOGmzrkEQRFPeM25ISExNp2rTpHbsbZ86cmTsb0/Hjx+nevTtvvvkmTzzxRJ52OdPstWnTBrD+NHuS3EWRRMdHcyr9FE2qNOHeuvcWvoIVjGg1gm3J25gVN0uSuyixnGn2CpOcnEzXrl157rnn8p1hKWeavfDwcJRSTJ06lbFjx1oi5HxJt4wokpwLqcNbDbfZhdRbPdbsMdyc3fjt+G8cPn/Y1uGIMqqo0+x99dVXJCQkEBkZmTt5dqVKNyaokWn2RJlx5vIZan1UCweDA8kvJFPdvbqtQ8o1dOlQ5u6ey+sdX2dSVyk7XVrIZ980Ms2esKkl+5aQrbLp0ahHqUrsoAuKAcz5Yw5Z2Vk2jkYI25PkLky2cO9CAPo362/jSG7XsW5H/Cv7k3wpmbVH19o6HCFsTpK7MElyWjK/H/+dCo4V6HtXX1uHcxuDwcCwYH3r96y4WTaORgjbk+QuTLJ432IUivCAcCpVqFT4CjYwOGgwDgYHlh5Yyrkr52wdjgB8fHxyL1TK484PHx/z3y8iyV2YpDR3yeSoU6kOPRr1wJht5Ls939k6HAGcP38epZQ8CnlYoky6JHdRqGOpx4hNisXN2Y3ejXvbOpwCDQ0eCsD3f31v20CEsDFJ7qJQi/YuAqBP4z64u7jbOJqC9W7cGw8XD7Ylb+PI+SO2DkcIm5HkLgpVFrpkcrg5u9G3ib7gu+AvmVtAlF+S3EWB4v+OZ9epXXi6eNIzoKetwzHJwOa67seCvZLcRfklyV0U6If9PwDQ966+uDq52jga03Rv1B0fVx/+OvMXf535y9bhCGETktxFgZYf1PWpH7rrIRtHYjoXRxceCXwEkK4ZUX5Jchd3lJKeQmxSLBUcK3B/o/ttHU6RDGg+AID5f82X2iaiXCowuV+7do3Q0FCCg4MJDAxkwoQJ+bYbM2YMAQEBBAUFERcXZ5FAhfWtPLQShaJbw254uFhvkgFz6FyvM7U8anH0wlG2n9xu63CEsLoCk7urqysbNmzgjz/+YM+ePWzYsIH//e9/edpER0dz+PBh4uPjmTlzJqNGjbJowMJ6lh1cBpA7+qQscXRw5LFmjwHSNSPKp0In63BzcwMgIyODrKwsKleunGf58uXLGTJkCAChoaGkpqaSkpJCjRo18rSLjIzMfR4WFkZYWFgJQxeWdMV4JbcAV2m/celOBjQfwLSt01i4dyFT7p+Cg0F6IUXpFhMTQ0xMjFm2VWhyz87OplWrVhw5coRRo0YRGBiYZ3lycjJ+fn65r+vUqUNSUlKByV2Ufr8e/ZVrmddo69uW2p61bR1OsYT6hlLXqy4nLp5ga9JW2vu1t3VIQhTo1hPfiRMnFntbhZ7KODg48Mcff5CUlMTGjRvz/Va59YJVaZmhRxRfTpdMROMIG0dSfAaDgYebPgzcGNIpRHlh8u9ULy8vevXqxY4dO/K87+vrS2JiYu7rpKQkfH19zRehsLqs7CxWHFwBQESTspvcAfo17QfAj/t/lFEzolwpMLmfO3eO1NRUAK5evcratWsJCQnJ0yYiIoJ58+YBEBsbi7e3921dMqJs2Zq8lbNXztLAuwHNq1tvtnZLaF+nPTXca5CQmsDulN22DkcIqykwuZ86dYquXbsSHBxMaGgoffr0oVu3bkRFRREVFQVAeHg4DRs2xN/fn5EjRzJjxgyrBC4sJ+fGpYgmEWW+i83RwZEH73oQkK4ZUb7IBNniNs1mNGPf2X2sG7yOrg262jqcElt7ZC33f3s/gdUC2fvMXluHI4TJZIJsYTYnLp5g39l9eLp4cm/de20djlmE1Q/Dx9WHfWf3ceDcAVuHI4RVSHIXeaw+vBqAbg274eLoYuNozMPZ0Tn3wvCP+3+0cTRCWIckd5HHqsOrAOjpXzbK+5oqZ0ikJHdRXkhyF7kysjJYd3QdAA/4P2DjaMzr/kb34+7szs5TOzmWeszW4QhhcZLcRa7NiZu5lHGJwGqB1PWqa+twzMrVyZVejXsBsPTAUhtHI4TlSXIXuey1SyZHTgG0nKGeQtizQmvLiPIj52KqpZN7VhZcuQJGo344O8Mt9egsoqd/T5wcnNh4fCMXrl7Ap6KP5XcqhI3ImbsAIDktmT0pe3B3djf7EMjLl2HJEnjpJejYETw9oVIlqFIFatbU/9aoAd26wdixsGoVZGaaNQQAfCr60KleJ7JUFtHx0ebfgRCliCR3Adw4a+/aoCsVnCqYZZv79+tk7esLjz4KH30E//sfXL0K7u7g7Q3VqoGHB5w5A+vXw/TpEB6u1xk7FnabuWJAbtfMIemaEfZNkrsAYPUR83XJHDkCEREQGKiT9cWLEBoKb74JP/8MZ89CejpcuKCTeloaHD+ul731FjRpot+fPh2Cg6F/fzh8uMRhATcKoa2KX8X1zOvm2agQpZCUHxBkZmdS9YOqXLx+kaNjjtLAp0GxtnPtGrz/Prz3Hly/rs/On3gCnn4abqk3VyClYNcumDMHvvxSb8vJCUaOhHfeAZ8SdpUH/TeIPSl7WDNoTZmbG1aUL1J+QJTI1qStXLx+kcZVGhc7scfFQYsWEBmpk/G//qXP4KOiipbYAQwGaN0aPv0U4uPhySchOxs+/xyCguC334oVYq6cs/ecmvVC2CNJ7iJ3Or3insV+9x106KC7TgIDISYG5s3TF0lLys8Pvv4a9uzRXTuJidClC7z+uh5pUxw5E5AsP7hcflEKuyXJXeQm9+4NuxdpvcxMeP55GDRId8k8+STs3AmdO5s/xmbN4Pff4f/+T5/Zv/uuTvLnzhV9W61rt6a2Z22S0pKIOx1n/mCFKAUkuZdzF69dZGvSVhwNjoTVDzN5vevX4eGH4ZNP9Dj1L76Ar74CV1fLxersDG+/rX8Z1KkDmzbpXwxHjhRtOw4GB/o07gPIDU3CfklyL+dijsWQpbJoV6cdlSpUMmmda9d0Yl+xQt98tGGDvmhqrXk9OnaE2Fjd/x4fD+3bw9atRduG3K0q7J0k93KuqF0yV69C374QHQ1Vq+qx6ffcY8kI8+frq7tpevTQQyu7dNFfMqbq0qALbs5uxJ2OIzkt2XKBCmEjktzLudzk3qjw5J6RoRP7L7/om482bNBnz7bi6al/PQwdqr90evfWN0mZwtXJlfsa3gfAz/E/Wy5IIWykwOSemJhIly5daNasGc2bN2f69Om3tYmJicHLy4uQkBBCQkKYNGmSxYIV5nXi4gkO/X2IShUq0da3bYFtldLjzNeuherVdb9381Iwd7azsx5NM2SIrlcTHm56F01Ov/vKQystGKEQtlFg4TBnZ2c+/vhjgoODSU9Pp3Xr1nTv3p2mTZvmade5c2eWL5e+y7Jm7RF91t6lfhecHAquIffuu/qmIjc33SUTGGiFAE3k4KATfEYGzJ+vu2rWr4dWrQpeLzwgHIBfj/7KVeNVKjpXtEK0QlhHgWfuNWvWJDg4GAAPDw+aNm3KyZMnb2snY4XLJlP72+fPhzfe0BdMv/9e32BU2jg66rH1/frpcge9esGJEwWvU9uzNq1rteZq5lU2HCtCh70QZYDJJX+PHTtGXFwcoaGhed43GAxs3ryZoKAgfH19mTJlCoH5nNZFRkbmPg8LCyMsLKzYQYuSy1bZrEvQsy4V1N++davu0wZd+KtvXysEV0xOTvrL54EH9PWAnD74SgUMAurduDc7T+1k5aGVuWfyQthKTEwMMTEx5tmYMsGlS5dU69at1U8//XTbsrS0NHX58mWllFLR0dEqICDgtjYm7kZY0c6TOxWRqLof11XZ2dn5tvn7b6Xq1VMKlHr6aaXu0KzUOX9eqSZNdNw9eihlNN657fbk7YpIlN9UvzseByFspSS5s9DRMkajkX79+jFo0CAefPDB25Z7enri5uYGQM+ePTEajZw/f9483zzCYnL627s37I4hnwHqSsGwYbpa4913w7Rp1hvHXlI+Pvq6QLVqsGYNjBlz57atarWipkdNEtMS+fPMn9YLUggLKzC5K6X497//TWBgIOPGjcu3TUpKSm6f+7Zt21BKUdka0+qIEsnpkskZDnirTz6B5ct1zfWFC8HFxZrRlVzDhrBsGVSooO+enT07/3YOBgd6Bei5VVccXGHFCIWwrAKT+6ZNm/j222/ZsGFD7lDHVatWERUVRVRUFABLliyhRYsWBAcHM27cOBYsWGCVwEXxXc+8zv9O6AHhXRt0vW351q3w8sv6+ezZ0KB4hSJtrn17ndgBnnkG/vgj/3a9G/cGYGW8DIkU9kPquZdDMcdi6DK3Cy2qt2DPqD15ll2+rG9MOnIExo2Djz+2UZBmNGKErnvTsCHs2HF7Pfj0jHSqflCVjKwMTr90muru1W0TqBC3kHruokhyumS6Nex227IJE3Rib9lST7xhDz79VI95P3pU3+yUnZ13uYeLB2H1w1Co3OkGhSjrJLmXQ+uO/pPcG+RN7jExOhE6OcHcuWWvn/1OXF31BN3e3rpcQX6/RnL63eVuVWEvJLmXM2nX09iWvA1HgyOd6nXKfT89XddjBz0Rxj/3rtmNBg30FxbAa6/d3v/eq7FO7muOrMGYVcxZQIQoRSS5lzMbj28kS2Vxt+/deUr8vvoqJCTo/vbXXrNhgBYUEaFLE2dkwOOP61o0ORr6NKRp1aakXU9jU+Im2wUphJlIci9ncvvbb+qS2bRJz09qb90x+fnoI7jrLti/H8aPz7ss5+xdumaEPZDkXs7c2t+emQmjRullr75q2xK+1uDmpksUODvDjBm6Dz5H7wA9JFJKAAt7IMm9HDlz+Qx/nvkTVydX2vu1B/QF1D//1H3S9todc6uQEF3lEvQwyb//1s87+HXAq4IXB84d4Mj5Is7dJ0QpI8m9HNmQoCsf3lv3XlydXElOhjff1Ms++wwqlqOKty+8oKfrS0mBsWP1e86OzvTw7wHI2bso+yS5lyO39re/8IIeJfPQQ3qSi/LEwQFmzdJfaN99p0sVwI0hkZLcRVknyb0cuTm5r10LixbpPuhPPrFxYDbi7w/vvaefP/00nD8PPf17YsBAzLEY0jPSbRugECUgyb2cOJZ6jKMXjuJVwYuW1VqRUwfuzTehbl3bxmZLo0fr7pnTp3X3TDX3aoTWCSUjK4Nfj/5q6/CEKDZJ7uVETn97WP0wZn3tyL590KgR3KHYZ7lxc/fMt9/C6tUyakbYB0nu5UROl0z7ml1zL6J+8IEuiVve+fvD22/r56NGQVe/f/rdD/0sBe9EmSXJvRxQSrE+YT0AB1Z15dw56NRJX0gV2rhxuuTCsWPww4wgfD19OZV+irjTcbYOTYhikeReDhz8+yCn0k9RxbU6333SDICpU8vOzErW4OQEM2fqbppPPjYQWvnG2bsQZZEk93Ig56zd7XRXjBkGBg+G1q1tHFQpdPfd+gJrVhbsWSwTeIiyTZJ7OZDT3564sSsVK964O1Pc7p13oE4dOLy2K05UYHvydlLSU2wdlhBFJsndzmWr7NyRMiR0ZexY8PW1bUylmaenvlsXozskdEWhWHV4la3DEqLIJLnbud2nd3Ph2gVIrYuXapg7N6q4s4gI6NULMvfJ3aqi7CowuScmJtKlSxeaNWtG8+bNmT59er7txowZQ0BAAEFBQcTFyeiC0mTdP/3tJHRjwquG2+YPFbczGGD6dHA5rpN79ME1ZGRl2DgqIYqmwOTu7OzMxx9/zN69e4mNjeXzzz9n//79edpER0dz+PBh4uPjmTlzJqNy6seKUuH7Lbq/3ftCV0aPtnEwZUjDhvDaM/UhpTlXsi4Rc/R/tg5JiCIpMLnXrFmT4H/mW/Pw8KBp06acPHkyT5vly5czZMgQAEJDQ0lNTSUlRS5AlQZXrxv548JGAF55rAtubjYOqIx55RXwPqvP3v+zWEbNiLLFydSGx44dIy4ujtDQ0DzvJycn4+fnl/u6Tp06JCUlUaNGjTztIiMjc5+HhYURFhZWvIiFySbN3o5yvozzxSa8+JpcRS0qV1d4/dHejN//Pr+nrOTUqanUqmXrqIQ9i4mJISYmxizbMim5p6en88gjjzBt2jQ8PDxuW37rLdqGfO6OuTm5C8szGmHG6nUQAp38uuLsbOuIyqZxj7TjtbcqY6wczzNvHuSnL5vYOiRhx2498Z04cWKxt1XoaBmj0Ui/fv0YNGgQDz744G3LfX19SUxMzH2dlJSEr4y1s7m5cyHVR/e3j+zerZDW4k6cHJwID9DF7pfuW0FsrI0DEsJEBSZ3pRT//ve/CQwMZNwdygdGREQwb948AGJjY/H29r6tS0ZYV0YGvP3eFaizBQMGujXqYuuQyrTHW/fRTxqvYMwYyM62bTxCmKLAbplNmzbx7bff0rJlS0JCQgB49913OXHiBAAjR44kPDyc6Oho/P39cXd3Z/bs2ZaPWhRozhxIZBM4ZRBSsxWVK1a2dUhlWo9GPXBycCKz7ia2LzzP3LmVGTbM1lEJUTCDskJNU4PBIKVTreT6dQgIgMQmr8K97zO+w3g+6P6BrcMq8+6bd58u4/DDt1RPeYJDh8DLy9ZRCXtXktwpd6jamVmzIDERXAPzzpcqSqZPY901U/WeFZw5A5Mm2TggIQohyd2OZGT8Myeo6wWuV96Js4Mz99a919Zh2YU+TXRyv1ZnNTgamTYNDh+2cVBCFECSux355ht91l7n3t9QKNrVaYe7i7utw7ILDX0aElgtkPTMi/QY8TtGI4wfb+uohLgzSe52IjPzRinfwF7SJWMJOV0zde9bgbs7LF0K69fbOCgh7kCSu51YsACOHtXzgSY6/5PcG0pyN6ec5L4+eQUTJuiLXOPG6ck9hChtJLnbgexs+M9/9PNRL59k/7n9uDu709a3rW0DszPt6rSjSsUqHLlwhJ6D91OvHvz5J3z1la0jE+J2ktztwA8/wIEDUK8eVG2jJ+boWK8jLo4uNo7Mvjg6ONKrsS4k9svx5Xz4oX7///4P0tJsGJgQ+ZDkXsYpdeOs/dVXIeaE9Ldb0oNNdAmOpQeW8sgjcM89cPasTF0oSh9J7mXczz/D7t1QuzYMGaL49eivgCR3S7m/0f24OrmyNXkrp9NPMXWqfv/jjyEhwbaxCXEzSe5lmFL/jGsHXnwRTlw+RGJaItXcqhFUM8i2wdkpdxd3ujfsDsDyg8tp2xaeeELfYzBhgo2DE+ImktzLsN9/h82boXJleOopWHt0LaBHyTgY5H+tpfRt0heAZQeXAfoL1tUVFi7U/z+EKA0kA5RhOWfto0eDhwe5XTI5Z5bCMvo06YMBA+sS1nHp+iX8/OCll/Sy55+XqpGidJDkXkbt2gWrV4O7u07umdmZbDimR8rc1/CvkJ14AAAgAElEQVQ+G0dn36q7V6eDXwcysjJYfXg1oKfkq1kTtm3TZ/BC2Jok9zJq8mT971NPQZUqsD15O2nX02hcpTF1veraNrhy4MG7/hk1c3ApoH855RQTe/VVuHrVVpEJoUlyL4MOHYIlS8DZGV54Qb+X098uXTLWkdPvHh0fjTHLCMDQodCyJZw4AdOm2TA4IZDkXiZ9+KEeKTN4MNSpo9+T5G5dAVUCCKwWSOq1VDYe3wiAoyN89JFe/u67cOaMDQMU5Z4k9zLm1CmYNw8MBnj5Zf3epeuXiE2KxdHgSFj9MJvGV57knL3/dOCn3Pfuuw969YJLl+Ctt2wVmRCS3MucTz7RY6ofeggaN9bv/Xb8NzKzM2nr2xYvV5keyFoeuushQCf3bHVjiMyHH+qz+JkzYe9eW0UnyjtJ7mVIaip88YV+/sorN97P7ZJpJF0y1tSmdhvqetXl5KWTbE3amvt+06YwcqQeEik134WtFJrcn3zySWrUqEGLFi3yXR4TE4OXlxchISGEhIQwSeYfs5j//lf/3O/SBdreVPAxZ3z7fQ1kCKQ1GQwGHm76MABL9i/JsywyEipVglWr4JdfbBCcKPcKTe7Dhg1j9erVBbbp3LkzcXFxxMXF8cYbb5gtOHHDtWu6SwbynrUnpyWz7+w+PFw8aFennW2CK8ceafoIAD/s+yHPRMbVqkHOR+HFF6Xmu7C+QpN7x44d8fHxKbBNcWfnFqabNw9SUiA4GO6//8b7a46sAaBrg644OzrbKLryq71fe2p51OL4xePsOrUrz7LRo6F+ffjrLz1xuRDW5FTSDRgMBjZv3kxQUBC+vr5MmTKFwMDA29pFRkbmPg8LCyMsLKykuy43srLIrR3+yit6pEyOnDskH2j0gA0iEw4GBx5q+hAzts/gh/0/0Lp269xlrq7w/vvQv78+ix8wADw9bRisKPViYmKIiYkxz8aUCRISElTz5s3zXZaWlqYuX76slFIqOjpaBQQE3NbGxN2IO1i8WClQqkEDpYzGG+8bs4zKe7K3IhJ15PwR2wVYzq07uk4RiQqYHqCys7PzLMvOVqp9e/3/b8IEGwUoyqyS5M4Sj5bx9PTEzc0NgJ49e2I0Gjl//nxJNyv+oRR88IF+/tJL4HTTb61tydtIvZZKQOUAGvo0tE2Agk71OlHVrSrx5+P568xfeZYZDLrWO8DUqXD8uA0CFOVSiZN7SkpKbp/7tm3bUEpRuXLlEgcmtJgY2L4dqlbVt7ffLLdLxl+6ZGzJycEpt9bMD/t/uG15aCgMHAjXr+u6M0JYQ6HJfeDAgXTo0IGDBw/i5+fHrFmziIqKIioqCoAlS5bQokULgoODGTduHAsWLLB40OXJ++/rf8eMgX9+IOXKuZgqyd32+jXtB8CSfUvyXT55su6DX7AAtmyxZmSivDIoZfmhLgaDQUbUFMPu3Xp0jJubLkZVpcqNZeeunKP6h9VxdnTm/MvncXdxt12ggoysDGpMqUHqtVT2PrOXwGq3Dyp44w09321oqJ7Uw0FuIRSFKEnulD+vUiynr3348LyJHWDtkbUoFJ3qdZLEXgq4OLrk3tC0cG/+Bd1ffVXXfN+6VZ/BC2FJktxLqWPH9KQPjo43yvrebPURGQJZ2gxsPhCA+X/Oz/dsy8NDn7mDTvRXrlgzOlHeSHIvpT76SI9vHzAA6tXLuyxbZbPmsPS3lzZh9cOo7l6d+PPxxJ2Oy7fNkCG6qy0xEaZMsXKAolyR5F4KnT0LX3+tn99caiDHnpQ9pFxOwdfTN9++XWEbTg5OPBr4KAAL/sq/38XR8UYZifffh6Qka0UnyhtJ7qXQZ5/padrCwyG/em2r4lcB0MO/B4abb1cVNjeg+QBA97vfXAb4Zp07wyOP6G4ZGRopLEWSeylz+bJO7pD/WTvAikMrAOgd0NtKUQlTdfDrgF8lP05cPEFsUuwd2334IVSoAN99J0MjhWVIci9lvvoKzp+Hdu2gY8fbl5+9fJbYpFhcHF2kfnsp5GBwoH/z/sCdu2ZAFxR76SX9fOxYXftdCHOS5F6KGI36FnW4vUBYjuj4aBSKLvW74OHiYd0AhUkGNNNdM4v2LiIr+861fl99FWrX1ncgf/ONtaIT5YUk91JkwQJ9s9Jdd0FERP5tVsavBKB3Y+mSKa1a1WqFf2V/Ui6nEHMs5o7tPDz0naugv8zT0qwTnygfJLmXEtnZN0oNvPxy/ncvZmRl5A6BlOReehkMhtwx79/sKfiU/IknoH17Xav/7betEZ0oLyS5lxIrV+rJlOvU0R/4/Gw8vpFLGZdoXr059b3rWzU+UTT/avkvQNeaSc9Iv2M7Bwd9Ad1ggGnTYP9+a0Uo7J0k91JAKXjvPf38xRfBxSX/disP6S6ZPo37WCkyUVwBVQLo4NeBy8bL/LT/pwLbtmoFI0ZAZiaMG6f/HoQoKUnupcDGjRAbq+vHjBiRfxul1I0hkNIlUyYMCRoCwNzdcwttO2kSeHvrybSXLbN0ZKI8kOReCuSctY8ZA+53qAF24NwBjl44SlW3qoT6hlovOFFsjzV7jAqOFVifsJ7Ei4kFtq1WDd55Rz9//nmpOyNKTpK7je3aBWvW6KT+3HN3bpdz1h4eEI6jg6OVohMl4e3qTd+7+qJQfLvn20LbP/00BAXponE5X/hCFJckdxvLGQo3ciQUNIGV3JVaNg1uORiAeXvmFVqX28kJZszQzz/4AA4dsnR0wp5JcrehgwdhyRJwds6/rG+O0+mn2XRiEy6OLvTw72G9AEWJ9fDvQXX36hw4d4DtJ7cX2r5DBxg2DDIy9C85ubgqikuSuw29957+8A4bBr6+d2639MBSFIruDbtTqUIl6wUoSszJwYknWuixraZcWAV9v4OPD6xdq7/8hSgOSe42kpAA336rS8DeqUBYjpxJl3Pm6RRly9DgoQB8/+f3XDEWfqW0WrUbfe7jxsGlSxYMTtitApP7k08+SY0aNWiRX93Zf4wZM4aAgACCgoKIi8t/ggJxu/ff15NxPPEENGx453bnr55nQ8IGHA2ORDS5Q00CUaq1rNGStr5tSb2WyuK9i01aZ/hwuPtuOHkS3nzTwgEKu1Rgch82bBirV6++4/Lo6GgOHz5MfHw8M2fOZNSoUWYP0B4lJcHs2fquxAkTCm67/OByslQWXRp0oYpblYIbi1LrqVZPATBz10yT2js6QlSUvoN1+nTYudOS0Ql7VGBy79ixIz4+Pndcvnz5coYM0TdqhIaGkpqaSkpKinkjtENTpugLZo88oouEFUS6ZOzDgOYD8HTxZHPiZv4685dJ64SE6G6Z7Owbd7AKYSqnkqycnJyMn59f7us6deqQlJREjRo1bmsbGRmZ+zwsLIywsLCS7LrMOnMGZv5z8vb66wW3Tbuexi9HfsGAgQfvetDywQmLcXdxZ1DLQXyx4wtm7pzJ9J7TTVpv4kR9UTUuDj79VN/gJOxXTEwMMTExZtlWiZI7cNvY3TtN+3Zzci/PpkzRU+j17q1vWCnIz4d+JiMrg451O1LTo6Z1AhQWM7L1SL7Y8QXzds9j8n2TcXN2K3QdDw899r13b/i//4N+/aBuXSsEK2zi1hPfiRMnFntbJRot4+vrS2Lijduqk5KS8C1oTF85d+YMfP65fv7WW4W3ly4Z+xJUM4i2vm25eP2iyRdWAXr1gkcf1VMwjholY9+FaUqU3CMiIpg3bx4AsbGxeHt759slI7QPP9Q1Q/r0gTZtCm57xXiFVYf1RNgPN33YCtEJaxjZeiQAUTujirTetGm6sFh0tJ53VYjCGFQB90QPHDiQ3377jXPnzlGjRg0mTpyI0WgEYORI/Uf63HPPsXr1atzd3Zk9ezatWrW6fScGQ6G3Xtu7M2f0vJlXr8KOHdC6dcHtF+1dRP8l/bm79t1sG7HNKjEKy7uccZlaH9XiUsYl4kbGEVwz2OR158zRN7z5+MC+fVBTeursXolyp7ICK+2mVHvxRaVAqT59TGsfMT9CEYn6ZMsnlg1MWN2YVWMUkaihS4cWab3sbKV69NB/R/36WSg4UaqUJHcWeOZuLuX9zD0lBRo00GftO3fqyRkKcu7KOWp9VItslU3yC8lyMdXOHDl/hIBPA3B2dOb4uONF+v974gQ0awbp6bB4sR5OK+xXSXKnlB+wgg8+0Ik9IqLwxA6weO9iMrMz6d6wuyR2O9SociMimkSQkZXBf3f8t0jr1q2r/54Ann0Wzp61QIDCLkhyt7Dk5BtlXE0dDfrdn/qK2aCWgywTlLC559vpAeszts/gWua1Iq07ciSEhenrOE8/LaNnRP4kuVvY22/DtWt6KFtISOHtEy4ksClxE27ObnLjkh3rVK8TITVDOHvlLPP/nF+kdR0cdPkKT0/48UcZPSPyJ8ndgg4fhq+/1h/Gt982bZ3v//wegAfvehAPFw8LRidsyWAwMK7dOAA+jv24yP2q9evr4ZGg674nFjyLnyiHJLlb0Jtv6sqPQ4cWXkMG9N2+3/6pp2PLqQEu7NeA5gOo6VGTP8/8yfqE9UVef+hQfR3n4kU9RDI72/wxirJLkruF7N4N8+eDi4tpd6MCxJ2O48C5A1Rzq0b3ht0tG6CwORdHF569+1kA3t/0fpHXNxh0naKqVWHdOl09UogcktwtJKco2KhRptcCmbdb3+3bv3l/nB2dLRSZKE2evftZPF08WXt0LVuTthZ5/Ro14Msv9fNXXtEFxoQASe4WsXEj/PwzuLvDa6+Zts5V49Xc5D40aKjlghOlik9FH55r+xwA72x8p1jbePBBfRKRkQEDBugx8EJIcjez7Gx48UX9/KWXoHp109Zbsm8JF65doHWt1rSuXUhtAmFXnm/3PG7Obvwc/zO7Tu0q1jY++kjf3HToEIwda+YARZkkyd3M5s/XtWNq1YLx401fL6eQVE5hKVF+VHOvxjN3PwPApI2TirWNihVhwQJwdYVZs/RzUb5Jcjejq1dvTJs3aZLuljHF3jN72ZS4CU8XTwa2GGi5AEWp9WL7F3F1cuWnAz/xZ8qfxdpG8+bw8cf6+YgRcPCgGQMUZY4kdzOaNk2PN27ZEv6ZfdAkOWftT7R8Qsa2l1M1PWryVGs9z+qk34t39g767tXHHtP97v366RrwonySwmFmcuYM+PvDpUvwyy/Q3cSRjFeMV6j9UW0uXr9Y5BKwwr4kpyXTcHpDMrIy2PnUTlrVMqEQUT4uXYK2beHAAXjiCfjmGz1sUpQ9UjisFHjjDf2h6tnT9MQOum77xesXaevbVhJ7OedbyZcxoWMAGL92fLE/1J6e8MMPulvwu+/gv0WrTSbshCR3M9i+Hb76Cpyc9KiFopALqeJmr937Gj6uPqxPWM/qw6uLvZ3AwBvj38eOhc2bzRSgKDMkuZdQdrYuvaqUnpm+aVPT192SuIXYpFi8KnjRv1l/ywUpygyfij680ekNAF7+9WWysrOKva2BA2HMGDAa4aGHdC14UX5Ici+hr7/WZ+61a+vZ6Yvig826MPczdz+Du4uJQ2uE3Xv27mep712fv878xdzdc0u0rY8+gvvu09eEIiLkBqfyRJJ7Cfz9942hj1On6r5OUx04d4BlB5ZRwbFCbj+rEAAVnCrwn67/AeD/NvwflzOKP+TFyQkWLYKAAF3vaMgQKTBWXhSa3FevXs1dd91FQEAA779/e3GjmJgYvLy8CAkJISQkhEmTij+Mq6x5/XWd4Lt21cPPiuKjLR+hUAwJHiKzLYnbDGg+gNa1WnPy0kn+8/t/SrQtHx9Yvhy8vHT995y6R8LOFTTBamZmpmrUqJFKSEhQGRkZKigoSO3bty9Pmw0bNqg+hcz6XMhuyqSNG/VExU5OSu3dW7R1T6adVC7vuChDpEEdPHfQMgGKMm/zic3KEGlQTm87qb9S/irx9tasUcrRUf/dfv65GQIUFleS3Fngmfu2bdvw9/enfv36ODs7M2DAAJYtW5bfF4SFvnpKp6tX4d//1s8nTNAjE4pi2tZpZGRl8HDTh2lcpbH5AxR2ob1fe0a2GUlmdiYjV44kW5WsP+X++2+MoHnuOfjpJzMEKUotp4IWJicn4+fnl/u6Tp06bN2atyypwWBg8+bNBAUF4evry5QpUwjMJ9tF3jSBaFhYGGFhYSWL3IYiIyE+Xif1ov7ETbuexhc7vgBgfIciFJ8R5dJ73d7jp/0/sSlxE7PiZjG81fASbW/YMD2v7//9Hzz+OPz6K9xzj5mCFSUWExNDTEyMeTZW0Gn9kiVL1PDhw3Nff/PNN+q5557L0yYtLU1dvnxZKaVUdHS0CggIMOtPi9Jm+3alHByUMhiU2rKl6OtPjJmoiER1nt3Z7LEJ+zT/z/mKSJTPZB+Vkp5S4u1lZys1cqTunvH2VmrXLjMEKSyiJLmzwG4ZX19fEm+anDExMZE6derkaePp6YmbmxsAPXv2xGg0cv78efN885QyGRm6OyY7G8aNg3btirb+mctn+HDzhwBMDJtogQiFPerfrD89GvXgwrULjF41usTdoAYDfP45PPwwpKbqO6r/+stMwYpSo8Dk3qZNG+Lj4zl27BgZGRksXLiQiIiIPG1SUlJy/9i2bduGUorKlStbLmIbevNN2LMHGjaEd4oxr8Lbv71NekY6vQJ60bl+Z/MHKOySwWBgRq8ZeLh4sGjvotxJXUrC0VGXp+7VS4/4uu8+qSJpbwpM7k5OTnz22Wf06NGDwMBA+vfvT9OmTYmKiiIqSt82v2TJElq0aEFwcDDjxo1jgZ0Wkl63Dj74ABwcYO5c08v55oj/O56onVE4GByYfN9kywQp7FZDn4Z81vMzAJ6NfpbD5w+XeJsuLrBkiU7sKSl6SO+hQyXerCglpCqkCc6dg6AgOHlST3Z907Vhk/Vf0p9FexcxLHgYs/rOMnuMwv4ppRj4w0AW7l3I3bXvZtOTm8wy1+7lyxAerqeHrFFDVzVt2dIMAYsSk6qQFqQUDB+uE/s99+jqj0W1PXk7i/YuwtXJVfraRbEZDAb+2/u/1PWqy/aT23kz5k2zbNfdXc/5262bPoMPC4Nt28yyaWFDktwLMWMGLFum7+777jt9O3dRZKtsxq0ZB8CY0DH4efkVsoYQd+bt6s13D3+nu/f+N5mlB5aaZbseHrByJfTtCxcu6ES/YYNZNi1sRJJ7AX7/XY+KAZg5E+rVK/o2Pt/2OZsTN1PToyYT7p1g3gBFuXRv3Xt5t+u7AAz6cRB/nP7DLNt1dYXFi/X49/R06NED5pX82q2wEUnud5CYCI88ApmZ8MILRa8dA3As9RgT1umE/kWvL/B29TZzlKK8evmel/lXy39x2XiZiPkRnE4/bZbtOjvrmZvGjdOlgocM0deZyvAls3JLkns+rlyBBx/UZVK7d4d86qUVSinFUyue4rLxMo8GPsqDdz1o/kBFuWUwGJjZZybt67QnMS2RhxY+xLXMa2bZtoODnmj7s8/087ff1tP1Xblils0LK5Hkfgul9Mzxu3bp8ewLFhS9nx1gzh9zWHt0LZUrVubTnp+aP1BR7rk6ufJT/5+o61WX2KRYHl38KBlZGWbb/rPPwooVuj9+/nxo316X3RBlgyT3mygF48fD99/rEQRLl0Jx7sdKuJDAC7+8AMC0B6ZRw6OGmSMVQqvhUYOVA1dSuWJlVh5aSf8l/TFmGc22/fBw2LJF14PfswfatNGfC1H6SXK/yeTJeuYaZ2d9c0eLFkXfxhXjFR5a+BCp11Lp07gPT7R4wvyBCnGTFjVa8Ou/fsXb1ZulB5by+I+Pk5mdabbtN28OO3bocgVpaXrKvjFjpJumtJPk/o+oKHjtNV1345tv4IEHir4NpRTDlw9nd8puAioHMO+heRgMBvMHK8QtQmqFsPZfa/Gq4MWSfUvov6Q/V41Xzbb9SpX0Cc+UKbqb8tNPoVUrPcWkKJ0kuQNz5sCoUfr5jBnQv5hzVU/dMpX5f83Hw8WDpQOWyugYYVVtardhzaA1VKpQiR/3/0iXuV1ISU8x2/YNBnjxRYiN1RPBHzyo++HffBOumedarjCjcp/cp07VNa6VgnffhaefLt52ouOjefnXlwGY++BcAqsVcQYPIcwgtE4om57cRD2vemxN3kroV6HsPbPXrPto3Rp27tTDJbOydBG95s1hzRqz7kaUULlN7krpbpgXX9Svp027Mdl1UUXHR/PQwofIVtm83vF1Hm76sPkCFaKImldvztbhWwn1DeX4xeO0/7o93+35zqz7qFhRD5fcuBGaNYMjR3RX5qOPwtGjZt2VKKZyWTjs6lV45hndHePoqP8dNKh421p5aCX9FvUjIyuD0W1HM+2BadLPLkqFq8arDFs2jIV7FwIwsPlAZvSaYfbuQqMRPvlEF9S7ckUPSHj6aV2HqXp1s+6q3ClJ7ix3yf3oUX3naVycvt160SLo06d421p6YCmPLX4MY7aRsaFj+bjHx5LYRamilOLruK8Zu3osV4xX8KvkR1TvKHoG9DT7vk6c0P3v8+bpX8YeHnqu1rFjoWZNs++uXJDkbqIVK2DwYD37TMOG8MMPEBxc9O1kZWcR+Vsk/9n4HxSK59s9z0f3fySJXZRa8X/H88SPT7D9pB7eEh4QztT7p9KkahOz7+vPP3UX588/69cVKsDQoboLNCDA7Luza5LcC3HmjP7D+vZb/bpvX90V412MX6cp6Sk8/uPjrE9Yj4PBgXe6vMOEeydIYhelnjHLyLSt03j7t7e5lHEJJwcnRrUZxfgO4y1SrTQ2VpfuuPmmp/vug5Ej9WfQueSl6O2eJPc7yM6G2bP1XacXLuhumHfe0Ym+qLk4W2Xz/Z/fM37teE6nn6a6e3Xm95tP1wZdLRO8EBaSkp7CGxve4OtdX6NQODk48XiLxxnfYTzNqzc3+/7279fj47///saQyRo19JDj/v31XMQO5XZoR8Ekud8iK0t3uUyapH8igi4ANmMG+PsXfXubEzczbvW43J+0nep1Yn6/+dT2rG3GqIWwrj0pe5j8v8ks3LuQbJUNQPs67RkcNJj+zfrjU9HHrPu7cEHfIBgVBfv23Xjfz0/f9RoeDp0765MwoUly/8elS7oe9ZQp+mwBoHZt+PBDGDiwaGfrmdmZrDy0ki92fMEvR34BoKZHTd7r9h6DgwbjYCjeqUZMTAxhYWHFWtfeyLG4wZbHIuFCAh9t+Yi5u+eSnpEOgIujC90bdic8IJye/j1p4NPAbPtTSpczWLhQP5KSbiyrWBFatIjhoYfC6NhR17KpUMFsuy5zSpQ7VSFWrVqlmjRpovz9/dXkyZPzbTN69Gjl7++vWrZsqXbt2nXbchN2U2xXrii1apVSTzyhVMWKSuk/HaXq1lXqiy+UunbN9G0Zs4zq9+O/q1d/fVXV/qi2IhJFJKripIrqjfVvqEvXL5U43rfeeqvE27AXcixuKA3HIv16uvpm9zeq+7zuyhBpyP37JxLlP91fDf5psPpi+xfqj1N/qOuZ182yz6wspTZtUuq115QKCcn5/L6V+zmuUEGptm2VevpppWbOVCo2VqnUVLPsukwoSe4ssJhtVlYWzz33HL/++iu+vr7cfffdRERE0LRp09w20dHRHD58mPj4eLZu3cqoUaOIjY0t3jdNoV9EcOoU/PWXnuNx/XrYvBmuX7/RplMnfcfp44/r2d3vvC1F8qVkdp/eze6U3ew4uYP1Ceu5eP1ibpsmVZrwdJunGRw0mMoVi1EeUogyxN3FnUEtBzGo5SBOp59mVfwqVh1exS9HfuHw+cMcPn+Yebv11EyOBkcCqgTQrFozGldpTAPvBtT3rk9dr7rU8KiBVwUvkwYZODhAhw768Z//6LmKR4/W4+N//x327tWf9VvndK1VC5o0gQYNoH59PUuar69+v1Yt8PEp+nU1e1Ngct+2bRv+/v7Ur18fgAEDBrBs2bI8yX358uUMGTIEgNDQUFJTU0lJSaFGjbxlbr9ZtzP3ec73MwqyFWRl6hmPMjN1or58Rc/IfumS4u+/4exZOHNWceI4XErP+YmiwJANNbIJCMjmno6ZdLnPSNXqmVzPvM6Sg1e5YrzCFeMVLly9wIVr+nHy0klOXDxB4sVErmbeXlipSZUm9AzoSd8mfelcr7OMghHlUk2PmgwLGcawkGEYs4zsSdnDlqQtbEnawrbkbRw5f4QD5w5w4NyBfNev4FiBGh41qFyxMt6u3vi4+lCpQiU8XDzwcPHA3dkdVyfX3IeLowvOjs44OzjjXOsAvcasJOIFR65ecSLhiAOH4x2IP+TIsWMGEk8YOHXdwKkjBmIOGwADqH/+/YejI1TyNODlrYueebgb8PDQpbwrVtT9+q6uusungosBJ2d9MujkdOPh4ACODuDgqP81OICDQf9rMPzz4J8vEcONL5OCUoZVs0lBp/WLFy9Ww4cPz339zTffqOeeey5Pm969e6tNmzblvu7WrZvasWPHbT8t5CEPechDHkV/WKRbxtSzVnVLh/+t6926XAghhGUVOOTD19eXxMTE3NeJiYnUqVOnwDZJSUn4+vqaOUwhhBBFUWByb9OmDfHx8Rw7doyMjAwWLlxIREREnjYRERHMm6cvssTGxuLt7X1bf7sQQgjrKrBbxsnJic8++4wePXqQlZXFv//9b5o2bUpUVBQAI0eOJDw8nOjoaPz9/XF3d2f27NlWCVwIIUQBit1bnw9zjIm3F4Udi2+//Va1bNlStWjRQnXo0EHt3r3bBlFahyl/F0optW3bNuXo6Kh++OEHK0ZnXaYciw0bNqjg4GDVrFkz1blzZ+sGaEWFHYuzZ8+qHj16qKCgINWsWTM1e/Zs6wdpBcOGDVPVq1dXzZs3v2Ob4uRNsyX3zMxM1ahRI5WQkKAyMjJUUFCQ2rdvX542P//8s+rZs6dSSqnY2FgVGhpqrt2XKqYci82bN6vUf+7GWLVqVbk+FjntunTponr16qWWLPq/SiYAAAQNSURBVFlig0gtz5RjceHCBRUYGKgSExOVUjrB2SNTjsVbb72lXn31VaWUPg6VK1dWRqPRFuFa1MaNG9WuXbvumNyLmzfNVq7n5jHxzs7OuWPib3anMfH2xpRj0b59e7y8vAB9LJJuvgfbjphyLAA+/fRTHnnkEapVq2aDKK3DlGPx/fff069fv9yBC1WrVrVFqBZnyrGoVasWaWlpAKSlpVGlShWcnArsSS6TOnbsiI/Pnev4FDdvmi25Jycn4+d3o2xonTp1SE5OLrSNPSY1U47Fzb7++mvCw8OtEZrVmfp3sWzZMkb9M0u5vd44ZsqxiI+P5/z583Tp0oU2bdrwzTffWDtMqzDlWIwYMYK9e/dSu3ZtgoKCmDZtmrXDLBWKmzfN9jVorjHx9qAo/00bNmxg1qxZbNq0yYIR2Y4px2LcuHFMnjw5t0jSrX8j9sKUY2E0Gtm1axfr1q3jypUrtG/fnnbt2hFgZ7NcmHIs3n33XYKDg4mJieHIkSN0796d3bt34+npaYUIS5fi5E2zJXcZE3+DKccCYM+ePYwYMYLVq1cX+LOsLDPlWOzcuZMBAwYAcO7cOVatWoWzs/Ntw27LOlOOhZ+fH1WrVqVixYpUrFiRTp06sXv3brtL7qYci82bN/P6668D0KhRIxo0aMDBgwdp06aNVWO1tWLnTbNcEVBKGY1G1bBhQ5WQkKCuX79e6AXVLVu22O1FRFOOxfHjx1WjRo3Uli1bbBSldZhyLG42dOhQux0tY8qx2L9/v+rWrZvKzMxUly9fVs2bN1d79+61UcSWY8qxeP7551VkZKRSSqnTp08rX19f9ffff9siXItLSEgw6YJqUfKm2c7cZUz8DaYci7fffpsLFy7k9jM7Ozuz7dbSd3bAlGNRXphyLO666y4eeOABWrZsiYODAyNGjCAwMNDGkZufKcfitddeY9iwYQQFBZGdnc0HH3xA5cr2V5114MCB/Pbbb5w7dw4/Pz8mTpyI0WgESpY3rTJZhxBCCOuSmQuFEMIOSXIXQgg7JMldCCHskCR3IYSwQ5LchRDCDklyF0IIO2R/VXiEMMGuXbtYtmwZfn5+1KxZk4MHD/Liiy/aOiwhzEbO3EW5dPXqVTw9Palduza9e/cmOjra1iEJYVaS3EW5dM8997B161Y6deqEUorTp0/bOiQhzEqSuyi3/v77bzw8PFi/fj19+/a1dThCmJX0uYty6ciRI2RmZrJixQq2b9/OxIkTbR2SEGYltWVEufTNN99gMBgYNGiQrUMRwiKkW0aUO6dOneKrr74qcHYsIco6OXMXQgg7JGfuQghhhyS5CyGEHZLkLoQQdkiSuxBC2CFJ7kIIYYckuQshhB2S5C6EEHZIkrsQQtih/wcmnRRJ9rwqywAAAABJRU5ErkJggg==">

用pymc模拟

假设我们的知识库中没有Beta分布,如何通过模拟实验找出$p$的概率分布呢?pymc是一个用于统计估计的库,它可以通过 先验概率观测值 模拟出 后验概率 的分布。下面先解释一下这两个概率:

  • 先验概率:在贝叶斯统计中,某一不确定量p的先验概率分布是在考虑"观测数据"前,能表达p不确定性的概率分布。

  • 后验概率:在考虑相关证据或数据后所得到的不确定量p的概率分布。

拿前面抛硬币的实验来说,如果在做实验之前能确信硬币出现正面的概率大概在0.5附近的话,那么它的先验概率就是一个以0.5为中心的山峰波形。而如果是某种特殊的硬币,我们对其出现正面的概率完全不了解,那么它的先验概率就是一个从0到1的平均分布。为了估计这个特殊硬币出现正面的概率,我们做了20次实验,其中出现了8次正面。通过这个实验,硬币出现正面的可能性的后验概率就如上图中的绿色曲线所示。

pymc库可以通过先验概率和观测值模拟出后验概率的分布,这对于一些复杂的系统的估计是很有用的。下面我们看看如何用pymc来对这个特殊硬币出现正面的可能性进行估计:

  • 首先pcoin是这个特殊硬币出现正面的概率,由于我们没有任何先验知识,因此它的先验概率是一个从0到1的平均分布(Uniform)。

  • 假设我们做了20次实验,其中8次为正面。根据前面的介绍可知,出现正面的次数符合二项分布(Binomial),并且这个二项分布的概率$p$为pcoin。这个通过value参数指定了实验的结果。因此experiment虽然是一个二项分布,但是它已经不能取其它值了。

In [32]:

import pymc
pcoin = pymc.Uniform("pcoin", 0, 1)
experiment = pymc.Binomial("experiment", 20, pcoin, value=8)

接下来通过MCMC对象模拟pcoin的后验概率。MCMC是Markov chain Monte Carlo(马尔科夫蒙特卡洛)的缩写,它是一种用马尔可夫链从随机分布取样的算法。通过调用MCMC对象的sample(),可以对pcoin的后验概率分布进行取样。这里30000为取样次数,5000表示不保存头5000次取样值。这时因为MCMC算法通常有一个收敛过程,我们希望只保留收敛之后的取样值。

In [33]:

mc = pymc.MCMC([pcoin])
mc.sample(30000, 5000) [****************100%******************] 30000 of 30000 complete

通过MCMC对象trace()可以获得某个不确定量的取样值。下面的程序获得pcoin的25000次取样值,并用hist()显示其分布情况。由结果可知pcoin的分布与前面介绍的Beta分布一致。

In [31]:

pcoin_trace = mc.trace("pcoin")[:]
hist(pcoin_trace, normed=True, bins=30);
plot(p, pbeta, "r", label="n=20", lw=2)
Out[31]: [<matplotlib.lines.Line2D at 0x5182190>]

aaarticlea/png;base64,iVBORw0KGgoAAAANSUhEUgAAAXcAAAD9CAYAAABHnDf0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FGWawPFf5UBoyCkSpAmDkKgJSBJEAw5CEBkEBFF0BC9AxCweHDu6qzPuCg4yOMPiMSIDjqJ4jCAoBEjiIqZhNYQoREA8CHcSMBBISCBAjq79o5ImnaPTSbq7+ni+n0/b3bxvVz0WzZPKU2+9r6KqqooQQgiv4qd3AEIIIRxPkrsQQnghSe5CCOGFJLkLIYQXkuQuhBBeSJK7EEJ4IbuSe3V1NQkJCYwdO7bR9pkzZxIdHU1cXBw5OTkODVAIIUTL2ZXcX3/9dWJjY1EUpUFbamoqBw4cIDc3l+XLlzNjxgyHBymEEKJlmk3u+fn5pKam8thjj9HY/U4pKSlMnjwZgMTEREpKSigsLHR8pEIIIewW0FyHOXPm8Le//Y3S0tJG2wsKCoiMjLS87969O/n5+URERFj+rLEzfiGEEM1r7SQCNs/cN27cSJcuXUhISLC5g/ptjSVzVVXloaq8+OKLusfgLg85FnIs5FjYfrSFzeSemZlJSkoK11xzDZMmTeKrr77ikUcesepjNBrJy8uzvM/Pz8doNLYpKOHbgoPDURSlyUdwcLjeIQrh9mwm9wULFpCXl8fhw4f55JNPuO2221i5cqVVn3Hjxln+LCsri9DQUKuSjBAtVVZWDKhNPrR2IYQtzdbc66ottyxbtgyA5ORkRo8eTWpqKlFRUXTs2JEVK1Y4PkovkpSUpHcIbkOOxWVyLC6TY+EYitrWwo49O1GUNtePhO/QTiJsfV/k+yR8Q1typ9yhKoQQXkiSuxBCeCFJ7kII4YUkuQuXa26ooxCi7eSCqnA5ey6YygVVIeSCqhBCiHpaNM5diDYzm5kAxPISQZQRTCmXuIJ/8hh76ad3dEJ4DSnLCNfJyoKnn4bvvmu0eSNj+AvPk8lgpCwjRNtypyR34XzFxTBnDrz/PgAFwEqeo5gwSgmmD/uYxjsYuADAO8DjVGHGv4kNyvdJ+AZJ7sJ9lZfD7bfD9u3Qrh088wydFizgfL0z886cYiZv8O8spiPlfMQDTOZ9qhutHMr3SfgGSe7CPVVVwT33wIYNEBkJW7ZAdLTN0TKD+T9SGUIQsJr7eJCPqCKwXi/5PgnfIMlduB9Vhccfh3/+E8LC4JtvICYGaH4o5EAU0gkmhFLWcg/38Smq1cAu+T4J3yBDIYX7mT9fS+zt28PGjZbEbo8sYDhbKCaUCXzGU7zpvDiF8FJy5i4cbnDHYLaWlwFwN7Ch0V7N38Q0jvWsZzwXaE9/dvEzMZZ2+T4JXyBn7sJ9VFfzP+Vl+AN/ZxYbGl1wwz4p3MUKptCBi3zAwwRQ6ZyYhfBCktyFY739NolAAd34b15q8+Zm8TpH+A0D2MmfeLnt8QnhI6QsIxynsBCuvx5KSriXT1nLvU10bNncMkMx8RW3YcaPRHawiwHyfRI+Qcoywj088wyUlJAKrGWCwza7lSReZxYBVPM3nnXYdoXwZnLmLhxj+3a45RZo355eFy9yuI2zPtZvD6GEQ/QinGJGAJvl+yR8gNPO3C9evEhiYiLx8fHExsby/PPPN+hjMpkICQkhISGBhIQE5s+f36pAhIdbsEB7/vd/57ATNn+WUF7hPwH4C4DZ7IS9COE9mj1zLy8vx2AwUFVVxeDBg1m0aBGDBw+2tJtMJhYvXkxKSkrTO5Ezd++2ezfEx0OHDnD0KEqXLrR1vvbG2jtQzgGi6MYJWL0a7ruvbXEL4eacWnM3GAwAVFRUUF1dTXh4eIM+krh93MKF2vPjj8NVVzltNxcwMI8XtTd/+hNUytBIIZrS7HzuZrOZ/v37c/DgQWbMmEFsbKxVu6IoZGZmEhcXh9FoZNGiRQ36AMydO9fyOikpiaSkpDYHL9zAgQPaWXRgIPzhD07f3bs8yjP8G9G5ufDeezB9utP3KYSrmEwmTCaTQ7Zl9wXVs2fPMnLkSBYuXGiVmMvKyvD398dgMJCWlsasWbPYv3+/9U6kLOO9Hn8c3n4bHn0U3nkHcMwyerba70fhE4AePeDgQQiQNWeEd3LJUMiQkBDGjBnDd/UWWggKCrKUbkaNGkVlZSVnzpxpVTDCwxQUaGfPigL/+Z8u2+1qgKgoOHYM1q1z2X6F8CQ2k3tRURElJSUAXLhwgc2bN5OQkGDVp7Cw0PKTJTs7G1VVG63LCy/02mta3fu+++Daa122WxVg1iztzauvumy/QngSm2WZvXv3MnnyZMxmM2azmYcffphnn32WZcuWAZCcnMySJUtYunQpAQEBGAwGFi9ezMCBA613ImUZ73PpEhiNcPo07NgBN99saXJ2WQYU1LIy6N4dzp5tsH8hvIXM5y5cb9UqmDgR4uIgJ0crzdRwSXJXVXj2WVi0CCZNgo8/blH4QngCmX5AuN7bb2vP06dbJXaXevpp8PeHTz+F/Hx9YhDCTUlyFy136JC2ZF779vDAA/rF0aMHTJigLee3ZIl+cQjhhiS5ixYJDg7n5d69AVh58SJKeDiKolg9XGr2bO152TJtMW4hBCDJXbRQeVkxU7kagLfZRsOFOFx8bWXQIO1ianExrFnj2n0L4cbkgqpokXGKQgrwM9cRw09oFz/rc/YF1UCgyvJuGvBPwAQMA4KCwigtlXsthOeTC6rCZWpv9v8nj9F4YneFKur+prCas5zHQBLQm1zKyop1iksI9yHJXdjv5ElGA5UEsJJH9I7GooxgPkWbIXIK7+kbjBBuQpK7sN+aNfgDXzCSU3TROxor7/IooCV3+VILYceskMK3BAeHN1nWMAFDgVXc78qQ7PJ/3EouUURzgBF6ByOEG5CTHGFFS+wNR8B0I59bUbgIrOcuPUNsgsIKpgLUnMML4dskuQu73Men+KGSilbjdkcreYRq/BgPdK439r7uIzhYJrYT3k+Su7DL/awCqPmveyqgO18wknbAg7xG42PwVRlNI3yCJHfRrN9whEFkcR4DG/UOphnvMQWAB5CJxIRvk+QumvV7bXkMNjAWd7/BfyN3UgYkkk0vDuodjhC6keQumnW5JON+o2Tqu4CB9TWvJ2qL8Qnhk2T6AWGl/lzsUeSSy7WUEkQXTnKJDjQ/f4zz53O31T4ahU3AXvrSj72Nfl6+j8ITyPQDwmkmsBbQhj9eor3O0dhnM3CGMG7gB/rwg97hCKELSe7CpnGkAPA5d+scif0qgTXcC0hpRvguKcsIK3XLMl0o5ARXU0E7OlPEeTrRfMkEO/o4v30YW/iK4RykF1EcwHqSM/k+Cs/gtLLMxYsXSUxMJD4+ntjYWJ5//vlG+82cOZPo6Gji4uLIyclpVSDC/dzJRvxQ2cLwmsTuObYylONcTW8OcRPf6h2OEC5nM7m3b9+ejIwMvv/+e/bs2UNGRgZff/21VZ/U1FQOHDhAbm4uy5cvZ8aMGU4NWLjOXTXjTtxzugHbzPizmt8DUpoRvqnZicMMBgMAFRUVVFdXEx5ufet2SkoKkydPBiAxMZGSkhIKCwuJiIiw6jd37lzL66SkJJKSktoYunCmDpQzgs2ANnbcE33CRGbzOvezimdYhCqXmISbM5lMmEwmh2yr2eRuNpvp378/Bw8eZMaMGcTGxlq1FxQUEBkZaXnfvXt38vPzbSZ34f5u50s6cJEd3MwJuukdTqvsIJGj9OA3HCORHWQxSO+QhLCp/onvvHnzWr2tZk9l/Pz8+P7778nPz2fbtm2N/lSpX/B3+SLJwuFqSzIpjNM5krZQ+Ix7gMtDOoXwFXb/nhoSEsKYMWP47rvvrP7caDSSl5dneZ+fn4/RaHRchMLl/KhmLBsAT0/usJYJANzDZ7h88W4hdGQzuRcVFVFSUgLAhQsX2Lx5MwkJCVZ9xo0bx8qVKwHIysoiNDS0QUlGeJZEdtCFUxziGn6gr97htMl2BvErEfTiMHHs1jscIVzGZs39xIkTTJ48GbPZjNls5uGHH2b48OEsW7YMgOTkZEaPHk1qaipRUVF07NiRFStWuCRw4Ty1Ny5pZ+2eXWIz4886xvNvLGMCa9lNvN4hCeESchOTsKIoCj8QSx9+5Da2kMFt9XvgCTcx1W2/nc1s5nfsI5a+7ENuYhKeoi25U5K7j7G1RipAJHAMKCWIzhRRSbt6PTwvuQdQSSERhFPM9fzEL8TI91F4BJk4TNitqTVSax931PTbwvBGErtnqiLQcmFYu7AqhPeT5C6sjKp5TrO88g61QyIluQtfIWUZH1N/vva6AqmgiCsIBnpwlDx6NLaFJj9vfx/Xt1/BRYroTCfO0xM4It9H4QGkLCMc4hYyCQb2EdtEYvdcl2jPJsYAMF7nWIRwBUnuwmIUaYD3lWRq1U6A5tm3ZQlhH0nuwuIO0gHvTe5pjKKSAIYAFDc9YkgIbyDJXQDQjQLi2MM54GsG6x2OU5QQxjaGaHfupabqHY4QTiXJXQCXz9q/Aiq4Qt9gnMgyN31Kir6BCOFkktwFULck490sE6GlpcGlS/oGI4QTSXIX+FNlWZgjXedYnO0oPbXpw8rKYOtWvcMRwmkkuQsS2UEoZ/mFazmidzAuYCnIrF+vZxhCOJUkd2E5a/9ffqdzJK5hSe4pKSA3MwkvJcldWJL7ZkboHIlr7ATo1g3y8yEnR+9whHAKSe4+LpizJLKDKvwxkaR3OC6hAowdq72RUTPCS0ly93FJmAigmiwGUkaw3uG4zl0yJFJ4N0nuPs7XSjIWw4aBwaCVZQoK9I5GCIeT5O7jfDa5t28Pt9+uvd60Sd9YhHACm8k9Ly+PYcOG0adPH/r27csbb7zRoI/JZCIkJISEhAQSEhKYP3++04IVjhXJMa5jP2cJJpub9Q7H9Wrr7hs36huHEE5gc4HswMBAXn31VeLj4zl37hw33ngjI0aMICYmxqrf0KFDSZHapcepPWvPYBjVtr8K3mn0aO35yy/hwgXo0EHfeIRwIJtn7l27diU+XlstvlOnTsTExHD8+PEG/WQhDs/ksyWZWt26wY03aok9I0PvaIRwKLtP144cOUJOTg6JiYlWf64oCpmZmcTFxWE0Glm0aBGxsbENPj937lzL66SkJJKSklodtGg7BTPD2QL4cHIHuPNO2LlTK83UnskLoROTyYTJZHLItuxaZu/cuXMkJSXxwgsvMH689To2ZWVl+Pv7YzAYSEtLY9asWezfv996J7LMntuoXWYvgV3s4kaO0oOeHEFbmg7avsSdI7bh/HbL9/G77+CmmyAyEo4eBUWx8TkhXMupy+xVVlYyYcIEHnrooQaJHSAoKAiDwQDAqFGjqKys5MyZM60KRriOdUnGhxNa//7QtSvk5cHevXpHI4TD2Ezuqqoybdo0YmNjmT17dqN9CgsLLT9ZsrOzUVWV8PBwx0cqHKq2JPMlt+scic78/GCMtrYqGzboG4sQDmSzLPP1118zZMgQ+vXrV/PrPCxYsIBjx44BkJyczJIlS1i6dCkBAQEYDAYWL17MwIEDrXciZRm3oSgK7bhIMWEYuEAXCjlFl7o98KmyDMC6dXD33TBwIGzfbuNzQrhWW3KnXTX3tpLk7j4URWEoGZgYxh5uII499Xvgc8n93Dno3BkqKuDXX6FLl6Y/KoQLObXmLrxPbUlmC8N1jsRNdOoESUna9L/p3r5cifAVktx9kCT3RtTW3eVuVeElpCzjY4IVhTP4AxDOmUZmgvTBsgzAoUPQuzcEB0NREQQG2vi8EK4hZRlhtyFAANV8y02+NcVvc3r1gpgYKC2Fb77ROxoh2kySu4+pLcRISaYRUpoRXkSSu4+R5G7DnXdqzzIFsPACUnP3JSdPQkQEF2hPGMVcon0jnXy05g5QWQlXXQVnz8KBA1oNXggdSc1d2Kdm5sOvGdxEYvdxgYEwcqT2Ws7ehYeT5O5LtsgQyGbV1t0luQsPJ2UZX9K7Nxw6xE1k8x03NdHJh8syAKdOQUSEdhZ/+rR2g5MQOpGyjGjekSNw6BAlwC766x2NzgJQFKXxR5cuZCt+2lQEX36pd6BCtJokd19RU283Aeaam5h8VxXamX3jj/Xmaq2blGaEB5Pk7itq6u1f6RyGJ7Ck9E2btPlmhPBAktx9garCV1pal+TevN0ARiOcOAE5OXqHI0SrSHL3Bb/8oiWqLl3Yp3csnkJGzQgPJ8ndF9SctXPbbfrG4Ulq71aVqQiEh5Lk7gtq6u2S3Fvgttvgiivg22+hsFDvaIRoMUnu3s5stoyUkeTeAh07asdLVSEtTe9ohGgxSe7ebvduKC6GHj20aW2F/aTuLjyYzeSel5fHsGHD6NOnD3379uWNN95otN/MmTOJjo4mLi6OHBld4F5q6+3Dh0PNIufCTrXJ/YsvtJuahPAgNpN7YGAgr776Kvv27SMrK4slS5bw008/WfVJTU3lwIED5Obmsnz5cmbM*gEULSb299Xr2hL59oawMvv5a72iEaBGbyb1r167Ex8cD0KlTJ2JiYjh+/LhVn5SUFCZPngxAYmIiJSUlFMoFKN0EB4dbbqUPVBTO1dSLjQ8/jCJn7i0nC3gIDxVgb8cjR46Qk5NDYmKi1Z8XFBQQGRlped+9e3fy8/OJiIiw6jd37lzL66SkJJKSkloXsbCprKyY2kmzbiKTTvyWn7mO4/xc00MSfIvceSe88oqW3Bcv1jsa4eVMJhMmk8kh27IruZ87d457772X119/nU6NzJJXf9ayxs4Q6yZ34RrDqZ1yQEoyrTZwIISHQ26udjPYddfpHZHwYvVPfOfNm9fqbTU7WqayspIJEybw0EMPMX78+AbtRqORvLw8y/v8/HyMRmOrAxKOU5vcZf72NggIgNGjtdcbNugbixAtYDO5q6rKtGnTiI2NZfbs2Y32GTduHCtXrgQgKyuL0NDQBiUZ4XodKGcQ2zGjkMEwvcPxbGPHas+S3IUHsblYx9dff82QIUPo16+fpdSyYMECjh07BkBycjIATz31FOnp6XTs2JEVK1bQv7/1fOGyWIfraH9PKrezmc38jp30ZwA76/bAuQthuGIfLl7M4+xZ6NxZu6Hp5EmtTCOEC7Qld9qsuQ8ePBiz2dzsRt58881W7Vw4j5RkHCgkBIYO1YaVpqXBgw/qHZEQzZI7VL2UJHcHk9KM8DCyhqqXURSFUM5wmiupIoAwiimnY90eSFmmufZAtNWaLrsGOASUAL06hXKmrNjG54VwDFlDVVgZylb8UMliYL3ELuzTcBm+w6jsI5ZQIP5cia7RCWEPSe5eSEoyzrEBrTQzVuc4hLCHJHcvJMndOaySu5QZhZuTmruX6aYoHAfO0ZFwzlBJu3o9pObe2nY/qikkgs6chn37IDbWxjaEaDupuQuL2tuV/o9bG0nsoi3M+LOJmonEUlL0DUaIZkhy9zK1hRgpyTjHOmqm4Fi3Tt9AhGiGlGW8iapyzM+PHkACu/iehEY6SVmmLe0GzlNEJzoAHD8OV19tYztCtI2UZYRm/356ACe5it3E6R2NVyqnI5tr30hpRrgxSe7eZLOWdrYwHFX+ap1mveXFelvdhNCVZABv8uWXAGxmhM6BeLcNoK1Hu2WLtgSfEG5Ikru3qKqCjAwAvuR2nYPxbqcAbrlFWzQ7PV3v*lCR3D1N3jdS6j0GBgVBayi9AHj30DtP7jZdRM8K9SXL3MJfXSLV+jEBbjmtz0x8VjnTXXdpzaipUVuobixCNkOTuJUbUpHVJ7i4SHa3doVpSAtu26R2NEA1IcvcCnShjIFlU4Y9J72B8Se3Z++ef6xuHEI2Q5O4FhrKVQKrI5mZK9Q7Gl9x9t/b8+edgx4plQriSJHcvcLkkI0MgXWrAAOjRQ7tTdccOvaMRwkqzyf3RRx8lIiKCG264odF2k8lESEgICQkJJCQkMH/+fIcHKWy7HW18uwyBdDFFgXvu0V6vWaNvLELU02xynzp1KunNjOUdOnQoOTk55OTk8MILLzgsONG8bhTQhx8poxNZDNQ7HN9z773a89q1Mse7cCvNJvdbb72VsLAwm31kUjD9jOQLAL7iNqoI1DkaHzRokDZ52NGjsGuX3tEIYRHQ1g0oikJmZiZxcXEYjUYWLVpEbCOLGMydO9fyOikpiaSkpLbuWgB3oP1Wlc4dOkfio/z8tAurb72lnb3feKPeEQkPZjKZMJlMDtmWXVP+HjlyhLFjx7J3794GbWVlZfj7+2MwGEhLS2PWrFns37/feicy5a/DKMrl6Wj9qeIUVxFGCb04yGF6of90ua7Yh/7tVt/nr76C4cO1se+//KLV4oVwAF2n/A0KCsJgMAAwatQoKisrOXPmTFs3K+xwM9mEUcJ+omsSu9DFkCHQuTPk5sIPP+gdjRCAA5J7YWGh5SdLdnY2qqoSHh7e5sBE86Qko5cA67l9AgN5u6gIgLn9+hEcLN9/ob9ma+6TJk1i69atFBUVERkZybx586ismUsjOTmZNWvWsHTpUgICAjAYDHzyySdOD1poai+mSnJ3tSrql23Wks50RnEvfZhXtk+fsISoQ5bZ8zC1NfcrKeIkXagkkHDOUE7H2h5Izd317YFUUEgEYZQQC/wo33fhALLMng8awWb8UNnGkDqJXeilknZ8hnZD0/06xyIESHL3WFJvdz//YhKA9l85cxc6k+TugRTMUm93QyaSKKQL1wLk5OgdjvBxktw9UD/20JVC8jHyIw1vGBP6qCaAT7lPeyMDC4TOJLl7oFGkAfAFI9Eu7gl38QkTtRerVsk0wEJXktw90Fg2ALCRO3WORNSXyS0cAzh2DLKy9A5H+DBJ7h6mMzCQLC7RTuZvd0MqfqyqfSOlGaEjSe4eZjTgh0oGwzhPJ73DEY2wpPTVq6G6Ws9QhA+T5O5hagsxUpJxX7sAoqKgsBAcNMOfEC0lyd2TVFQwsualJHc3N0kb884HH+gbh/BZktw9ybZtBAN76ctReuodjbDl4Ye15zVr4Nw5fWMRPkmSuyfZuBGADYzVORDRrOhouOUWOH8ePv9c72iED5Lk7ilUFTbIEEiPMnmy9vz++/rGIXySzArpKX76CWJjOQV0pQoz/k101HvGRHeIQf92VVWhpAS6doWKCm2N1chIG58RoiGZFdIX1Jy1p4KNxC7cSmgo3HWX9lvXhx/qHY3wMZLcPYWlJCM8yiOPaM8rV8pMkcKlpCzjCX79Fbp1g8BAgisqKHPrkoQ7xKB/u+X7XlUFRiOcPAk7dsDNN9v4nBDWpCzj7dat0876RoygTO9YRMsEBMCDD2qv5cKqcCFJ7p5g7VrtecIEfeMQrTNlivb88cdQXq5rKMJ32Ezujz76KBEREdxwww1N9pk5cybR0dHExcWRIwsUON6ZM5CRAf7+MG6c3tGI1ujXTyvHlJTAp5/qHY3wETaT+9SpU0lPT2+yPTU1lQMHDpCbm8vy5cuZMWOGwwP0eSkp2uRTw4bBlVfqHY1orccf156XL9c3DuEzbCb3W2+9lbCwsCbbU1JSmFxzo0ZiYiIlJSUUFhY6NkIfExwcjqIolseGqVMB+Lcvv0RRZGEOzxBg9XeoKAodH3uMUoDMTBI7hugdoPABAW35cEFBAZF1bszo3r07+fn5RERENOg7d+5cy+ukpCSSkpLasmuvVVZWTO1IjCBK+R1XYaaSdRwHuiIrL3mCKuqPpikHPuQJnmApD5aX6hKVcH8mkwmTg2YSbVNyBxoM02nq7LJuchf2GcMmrqCCbdxKIV31Dke00TKSeYKlPALahVWDQe+QhJupf+I7b968Vm+rTaNljEYjeXl5lvf5+fkYjca2bFLUMQFtlMxaZJSMN9hDHDu4mVCQC6vC6dqU3MeNG8fKlSsByMrKIjQ0tNGSjGi5DpRbFsL+jHt0jkY4yjKSa14s0zcQ4fVs3qE6adIktm7dSlFREREREcybN4/KykoAkpO1L+lTTz1Feno6HTt2ZMWKFfTv37/hTuQOVbtpZS2V+1jNau4nm5tIJLtuD/S++1LuUG19u4HznKATwQA5ORAfb2M7wte1JXfK9ANupja5r+Mu7iKFWbzGG8yq2wN3TVzuE4N7t7+Gov2NTpkCK1bY2I7wdZLcvYiiKFzJKU5wNX6YMVJQ72Kqeycu94jBvdt7oXBQUSAwUJsKuKtcLBeNk7llvMx9fEogVWxmhIyS8UKHQLvbuKIC/vEPvcMRXkqSuxt6kI8A+JCHdI5EOM2cOdrzW2/BxYv6xiK8kiR3N9MTGMw3nMfAOsbrHY5wliFDICEBTp2Cf/1L72iEF5Lk7mYeqHlex3jO00nXWIQTKQrMnq29fvVVWchDOJwkd3eiqpZCzEc8qGsowgUmTtQupu7dC199pXc0wstIcncnOTnEACe5is2M0Dsa4Wzt2sGTT2qvX3lF31iE15Hk7k5q7vZdxf1UEahzMMIlnnwSgoJg82ZtGT4hHESSu7u4cMGS3N9jir6xCNcJC4OnntJe//nP+sYivIokd3exZg0UF/MdsIsb9Y5GuNKcOdoMkZs2wa5dekcjvIQkd3dRM5GUTCflg666Cp54Qns9f76+sQivIdMPuIN9+6BvXwgKolNZGefd+NZ5mX7AMe0N/j38+itcc412Q9OePWBj3WLhO2T6AQ9Sfxk9RVF4o29fAJaWlXFe5/iETrp2vbzOqpy9CweQM3cXq531sVYHyjlON0I5Szw57CYBvc8q5czd2e2BaEvxWeuGNu/MFQA7d0Ij02cL3yJn7h7s96wmlLPs4GZ2I3N7+4baNVatH8dReYNntC7PPit3rYo2keSus2RqL6Qm6xyJcAcL+CNnQLtjNT1d73CEB5OyjIvVLcsMZDvbuYUSQjBSQDkd0b9kIGUZvdvnoLAYtIvs338P/v42tiW8mZRlPNR/8FcA3uKJmsQuBCwB6NkTfvgB3n9f52iEp5IzdxerPXO/jp/5kVgqaEdPjtRZlMO9zyrlzN017epHH8GDD0K3brB/P3SUH/6+yKln7unp6Vx//fVER0fzSiOTG5mnZf6SAAAM+UlEQVRMJkJCQkhISCAhIYH5MozLLn/gf/BD5X0my2pLoqGJE+HGG+H4cXj5Zb2jEZ5ItaGqqkrt3bu3evjwYbWiokKNi4tTf/zxR6s+GRkZ6tixY21tRm1mNz4FULtyXL1IO7UaRY3mF1UbFlH7oN77+g93b3eHGDy/XVVVVc3MVFVFUdWAAFX94Qd9v7hCF23JnTbP3LOzs4mKiqJnz54EBgYyceJE1q9f39gPCIf8oPEVs3idK6jgM+4hl2v1Dke4q0GDIDkZqqq0Z7NZ74iEBwmw1VhQUEBkZKTlfffu3dlRb1pSRVHIzMwkLi4Oo9HIokWLiI2NbbCtuXPnWl4nJSWRlJTUtsg9VBAwg6UA/I1n9Q1GuL+//AU+/xy++QbefRcee0zviIQTmUwmTCaTQ7ZlM7lrF/9s69+/P3l5eRgMBtLS0hg/fjz79+9v0K9ucvdls4EQSjExlGwS9Q5HuKUAq3979wOfAGemT+f66dO5GBRGaekZ3aITzlP/xHfevHmt3pbNsozRaCQvL8/yPi8vj+7du1v1CQoKwmAwADBq1CgqKys5c0a+eI06edJyrv4irf9LE97O+g7WVZhJZyThwN/5PWVlxfqGJzyCzeQ+YMAAcnNzOXLkCBUVFaxatYpx48ZZ9SksLLTU3LOzs1FVlfDwcOdF7MleeokgYCNj2MZQvaMRHkPhCd6ijE7cz2oe0Tsc4RFslmUCAgJ48803GTlyJNXV1UybNo2YmBiW1cw9npyczJo1a1i6dCkBAQEYDAY++eQTlwTucXJzYdkyqoHnWKh3NMLDHKYXT/Em7zNFu8npwAGIitI7LOHG5CYmV7n/fli9mneBaW5+A03b2t0hBm9tV/kXk5jIKrjpJu0ia6CstevN2pI7Jbm7wrffws03Q/v2RF68SL5bJg5HtbtDDN7bHkIJuwnjNwDPPaeNphFeS+aWcWdmM8yerb2eOZN8faMRHu4soTwI4OcHCxfCunV6hyTclCR3Z1uyBDIztZV2nn9e72iEF/gGYMEC7c1DD2kzRwpRj5RlnOnIEW3a1vPntRtRxo9vsBJTQ57e7g4xeH+7ajbD5MnwwQcQGQnZ2doJhPAqUpZxR6qqrYl5/jzcdx+MH693RMKbKAosX65NUZCXB3ffrS2uLUQNOXN3sODgcMrKipkCrABOA7HASate7n1WKGfu7t9u+fdUWKhdrD92DO68E9auhXbtbHxWeBI5c3cjZWXF9OQQiwkFYBYfcLLO3YZCtJ02PYGiKChdu3LDsWOcBti4kc+vuILwoDC9AxRuQJK7g3UAPuduwighhbF8pI1tEMKBrKcn+AGV29lFMaHcDSw7V6LNJCl8miR3R1JV/gnEs5v9RPMIK9F+xRbCub4ngRFspoQQ7gPtprkLF/QOS+hIkrsjLV7MA0AZnRjPOs7WlGaEcIWdDGAkX3AW4LPPYNgwrSYvfJIkd0dJTYX/+A8AJvM+P9FwTnshnC2bRH4L8JvfwI4dkJgI+/bpHZbQgSR3R0hN1Yaimc3MBz7nHr0jEj5sH1xO7EePasMlP/pI77CEi0lyb6uNG7XEXlEBTz/Nf+kdjxAAERGQkaHV3svKtDtZH3gASkr0jky4iIxzb6HacewAdwGrgXbAa8AcSy/3HiMt49y9v93y701V4Z13YNYsKC/X7mZdtgxGjbLxeeEuZJy7C5WVFeNHFS/xAp+h0A5YzBzmYEbGsQv3UGccvJ8fyvTpRJeXkw3a3ayjR8OYMfDLL3oHKpxIknsLdQH+l9/xX8xHReGPvMwf+B9kyKNwH9bj4EHlACq/pYJn+BuloF0n6tsXZs7UEr7wOlKWsZfZDB9/zImHH+ZqoJAuTOJfZHBbvY76/0ouZRlpt9XeBYXCxx7TyjWqCgEBWj3+2We1hC/chizW4WyZmdqc7N9+C8BWhjCJf3GCbo10du9/2JLcpR0CgSpuAJ4D7gf8a1oygdVXGHjtRD6EyTQGepOauwMFB4ejKAoBisJ4ReELRYHf/ha+/ZYTwBRgGBlNJHbRMia9A3AjJhfuSyvb7EXlQVSiOMSbPEkZnbgFeO1SuTZ98J13wltvweHDLowNTCaTS/fnrZpN7unp6Vx//fVER0fzyiuvNNpn5syZREdHExcXR05OjsODdJmqKvqVFbOA5zhKN9YBI4FyOvBnXiCaMt4HVPmZ6CAmvQNwIybd9nyEa3iaN+nKrzzEB/wvQGUlbNoETz4JvXpBdLQ2f/w//gG7d2tDf51EkruDqDZUVVWpvXv3Vg8fPqxWVFSocXFx6o8//mjVZ9OmTeqoUaNUVVXVrKwsNTExscF2mtmNPsxmVc3LU9WNG1X15ZdV9e67VTUkRFW1KqSqgvoT16mzeFUN43SdP6Zul0Yevt7ekm286Kb/D3q0v6jz/us+AtQIUKeAugrU4kY6VYL6I6ifgvoyqI+Bejuo14IaQkDNPhp/BAWF2fyn+eKLL7omB3iAtuTOAFuJPzs7m6ioKHr27AnAxIkTWb9+PTExMZY+KSkpTJ48GYDExERKSkooLCwkIiLCemM7d7bohw5Q8z2s97rus9l8+VFVpZ1tVFXBpUvapEnl5dqjuNjyqM7L4/zPP2M4fZqARs4+fgbSmM167mIrQ5FRMML3VFGIynvAe0AAlfRjD4PYziC2czMf0xuFGFRimvj8Ra6gkAjOEE4JoRQTRinBnKMT58qWwJ//DO3bX360aweBgdrj55+1mwP9/bWLvX5+2sPfX1ukxNajrrrv67c11c+L2EzuBQUFREZGWt53796dHTt2NNsnPz+/QXJXBgxwRLwu8lrNoynNfRl8vb0l25jnpBg8sX1eM+3O3n/j7VXArprHEsufqg0+Ye0ScKzm0dDz//3fNj89b9WqZrYvmmMzuSt2/kRTVeu/6Pqfq98uhBDCuWxeGTQajeTVucEhLy+P7t272+yTn5+P0Wh0cJhCCCFawmZyHzBgALm5uRw5coSKigpWrVrFuHHjrPqMGzeOlStXApCVlUVoaGjDersQQgiXslmWCQgI4M0332TkyJFUV1czbdo0YmJiWLZsGQDJycmMHj2a1NRUoqKi6NixIytWrHBJ4EIIIWxw2JgdVVXT0tLU6667To2KilIXLlzYaJ+nn35ajYqKUvv166fu2rXLkbt3K80diw8//FDt16+fesMNN6i33HKLunv3bh2idA17vheqqqrZ2dmqv7+/unbtWhdG51r2HIuMjAw1Pj5e7dOnjzp06FDXBuhCzR2LU6dOqSNHjlTj4uLUPn36qCtWrHB9kC4wdepUtUuXLmrfvn2b7NOavOmw5O6oMfHewJ5jkZmZqZaUlKiqqn3JfflY1PYbNmyYOmbMGHXNmjU6ROp89hyL4uJiNTY2Vs3Ly1NVVUtw3sieY/Hiiy+qzz33nKqq2nEIDw9XKysr9QjXqbZt26bu2rWryeTe2rzpsFst646JDwwMtIyJr6upMfHexp5jMWjQIEJCQgDtWOTn5+sRqtPZcywA/v73v3Pvvfdy1VVX6RCla9hzLD7++GMmTJhgGbjQuXNnPUJ1OnuOxdVXX01paSkApaWlXHnllQQE2Kwke6Rbb72VMBvz+LQ2bzosuTc23r2goKDZPt6Y1Ow5FnW98847jB492hWhuZy934v169czY8YMwP4huJ7GnmORm5vLmTNnGDZsGAMGDOCDDz5wdZguYc+xmD59Ovv27aNbt27ExcXx+uuvuzpMt9DavOmwH4OOGhPvDVry/5SRkcG7777LN99848SI9GPPsZg9ezYLFy60zIBX/zviLew5FpWVlezatYstW7ZQXl7OoEGDGDhwINHR0S6I0HXsORYLFiwgPj4ek8nEwYMHGTFiBLt37yYoKMgFEbqX1uRNhyV3GRN/mT3HAmDPnj1Mnz6d9PR0m7+WeTJ7jsXOnTuZOHEiAEVFRaSlpREYGNhg2K2ns+dYREZG0rlzZzp06ECHDh0YMmQIu3fv9rrkbs+xyMzM5E9/+hMAvXv35pprruGXX35hgEfd7d52rc6bDrkioKpqZWWl2qtXL/Xw4cPqpUuXmr2gun37dq+9iGjPsTh69Kjau3dvdfv27TpF6Rr2HIu6pkyZ4rWjZew5Fj/99JM6fPhwtaqqSj1//rzat29fdd++fTpF7Dz2HIs5c+aoc+fOVVVVVX/99VfVaDSqp0+f1iNcpzt8+LBdF1RbkjcdduYuY+Ivs+dYvPTSSxQXF1vqzIGBgWRnZ+sZtlPYcyx8hT3H4vrrr+eOO+6gX79++Pn5MX36dGJjY3WO3PHsORZ//OMfmTp1KnFxcZjNZv76178SHh6uc+SON2nSJLZu3UpRURGRkZHMmzePyspKoG150yUrMQkhhHAtWXVCCCG8kCR3IYTwQpLchRDCC0lyF0IILyTJXQghvJAkdyGE8EL/D1dyLp3Ej5AeAAAAAElFTkSuQmCC">

In [34]:

pcoin_trace.shape
Out[34]: (25000,)