{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import gibbs\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from sklearn.metrics import normalized_mutual_info_score\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Mixture parameters\n", "\n", "N=1000\n", "\n", "K=4\n", "p = np.array([0.25,0.25,0.25,0.25])\n", "\n", "mu= np.array([1,3,5,7])\n", "\n", "sigma = np.array([0.5,0.5,0.5,0.5])\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "#sample the mixture to create artificial data\n", "\n", "np.random.seed(0)\n", "\n", "Zt = np.random.randint(0,K,N)\n", "\n", "X=np.zeros(N)\n", "\n", "for k in range(K):\n", " pts=np.where(Zt==k)[0]\n", " npts=len(pts)\n", " X[pts]=np.random.normal(mu[k],sigma[k],npts)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Prior parameters\n", "mu0 = np.array([0.,0.,0.,0.])\n", "sigma0 = np.array([10.,10.,10.,10.])\n", "\n", "c = np.array([1.,1.,1.,1.])\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#Gibbs sampling\n", "Niter=1000\n", "\n", "burn_in=0.\n", "\n", "replica=1\n", "\n", "sampling_rate=1;\n", "Nsamp=np.floor((Niter-np.ceil(burn_in*Niter))/sampling_rate).astype(int)\n", "Npar=N+2*K\n", "\n", "sampling=np.ones((Nsamp,Npar));\n", "\n", "\n", "gibbs.GibbsSampling(Niter,K,sampling_rate,burn_in,replica,X,mu0,sigma0,sigma,c,sampling)\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABDtklEQVR4nO3deZhcV3ng/+9bt/at903drX2xJXmRLRsbG2PANuuAIUBCIGGNMwwhhExmJiSThAy/ZIYkP0IYmAQCBJIhCQlhi9nMbmPjRbZsS9ZirS2p97329Z7541S3WlJLVquXKqvez/PUI3XVrXvfuvfc855zl3PFGINSSilVazzVDkAppZSajyYopZRSNUkTlFJKqZqkCUoppVRN0gSllFKqJnmrHcBCtLa2mrVr11Y7DKWUUkvo8ccfHzPGtJ39/vMqQa1du5Zdu3ZVOwyllFJLSET65ntfD/EppZSqSZqglFJK1SRNUEoppWqSJiillFI1SROUUkqpmqQJSimlVE3SBKWUUqomVS1BicgWEXlyzishIr9VrXiUUs9v+uigy0/VbtQ1xhwErgUQEQfoB75WrXguRaaYYffIbtLFNNlSllwpB4DH48ErXopukYncBBO5CRKFBKsiq9jSvIUrmq+gN9aLR6rXgc0UM3z72Le5/9T9tIXaWN+4nnXxdTQFmxARBKFkSoxnxxnJjDCaGaXgFma/b4yh6BYpuSWKbhHXuLjGpWzKhLwhemO99MR6WBNbw/rG9Uv6W8tumf5UP0enj9KX6MMjHuL+ODF/jIgvQsAJ4HN8+D1+vB4vXvHi8XiI+qLE/XFEBADXuAykBjg0eQiDYVV0FV2RrjOmmfmd2VKWTDFDrpyb/d0lt0TIGyLujxMPxAk6wdnvzV3P47lxXOPS4G8g5o/heJx5f1fRLTKUGiIeiJ8RQzW5xiVbyiIIYV/4jM+KbpFnxp6hP9VP0Bsk5A0R9oZpCbbQGm4l5A2RLqbZNbSLhwcf5uDkQdrD7ayJr2FtfC072nfQGek8Y57GGI5MHaE51ExzsHneePpT/RydOsrhqcP0Jfo4lTrFqeQpRjOj9MZ72d6ynW2t29jctJnuaDft4Xa8nour6tLFNABhb/iM9e8al1wph0c8BJzAoreNMYZcOUcin2AsN8ZYZoyJ3ARbW7aypXnLJc2z5JbIlXL4HT9+x7/g77vGZSw7xqnkKabz00R8EaL+KFFfdHb+JVMi4AToinRd0jIWSmqh1SEidwF/ZIy55ULT7dy509TCSBLGGL597Nt87PGPMZIZec7pY/4YMV+M4cwwZVMGIOgEWdewjo2NG+mN95ItZZnOTzOZmyRTzJAt24RXckvE/DEaAg00+BtoD7ezKrqK7mg3MX+MkcwIg+lBRjIj+B0/TYEmGoONGGPoS/TRl+hjIDVAU7CJnlgPvbFejk0f496j95IupumOdpMoJEgWks/5O87eyX0eH16PF5/HhyMOHvHgiEOqmCJRSMxO1xxs5rae27i953Y2Nm0kU8yQKWWYyk9xaPIQBycOcnDyIMYYemI9dEe7aQm1MJGbmE2O6WKaglugUC6QKqTOSJYLEXACtIXaiPljnEiemK2Qzv5dYLdz2ZQxXNw+4ohD0Bsk4ATwO36m89NkS9kzphGExkAjm5o2sbVlK1c0X8FUfoqHBx7mseHHZuPxeXy0hFpoD7XTFm6jPdxOY6CRsinPJsep/BTj2XHGc+OU3BKdkU46I520h9tnfwNAU6CJrS1b2di4EZ9j38+X8wymBkkUEuRKOXLlHNP5aY5MHeHw1GEOTx22ZbGUmZ1Pd7SbzU2bWduwlqNTR3ls6LEzPj9bzBcjW8pSMiX8Hj+bmzYznhtnMD04uy5u6rqJ1218HVe3Xc0P+n7ANw5/gyPTR/CKl1u6b+E161/DFc1X8OjQozzY/yCPDj1KqpiaXUZrqNU2hqI9tIZaOTZ9jL3jexnLjp2xXVpDrXg9XgTBIx5C3pDdpwINeD1e+lP9nEycZDI/CYBHbIPG7/hny+tcfo+fsC9MY6CR5mAzjYFGmoJNs/+W3BLDmWGG08OM5cbIlXLky3lypRzZUpZ0MT1bF5xtR/sOfmnLL3F77+1M5CYYSg8xlBniZOIkfck+TiROMJ4dn20cFt0iubKtK2Y0BhppD7fTHGy2n1eWKwh+x0/ACeB4HPKlPLlyjlwpx2h2lHw5f97tOZcgtIXa6I5184HrPsD1Hddf1PfOOz+Rx40xO895v0YS1OeBJ4wxn5zns3uAewBWr159fV/fvCNiLJuZij5TylAoF0gUEnxuz+d4YuQJtrZs5X3Xvo+OcAdhb5iQLwTYFr5rXByPQ1Og6YxK4fDUYZ6deJZDU4dmK4ORzAg+j4/GQCMNgQbi/jhBb5CgE8TxOCQLSabz00znpxnJjpxREGf4PD7Kxi53hkc8dEe7WRVdxWRukpPJk2RLWQJOgJevfTlv2vwmrmm7BoCJ3ARHp4/axGLAYPCIh9ZQK+3hdlqCLbO/42JM56c5mTzJ4anDPNj/IA/2P0iyeG4SFITV8dVsbtqMV7ycSp2iP9XPZG6SpmATHeEO2sJtRHwR/B7bMoz6o6yLr2Ndg30BJAoJEoUEmaLdToVygYJboOyWZyv1VDHFSGaE4cwwiXyC1fHVbGraxKbGTXg9XgbTgwykBhjPjoMwm3SDju0dhLwhgt4gPo/PJmWPQ6aUIZG3y04X07OVfaFcIO6P0xpqpSXUgiOO3YaFaUYzoxyYOMCzk89SdIsA9MZ6uanrJra3bidZSDKeG5/tvc68ZirmmUZBQ6CB1qCdv0c8DGeGGUoPMZGbmHeb+D1+VsdXM5mbZDw3Pu80Xo/XNpwaNtIabiXiixD2himUCxyeOszByYP0JfrojfXygs4X8IKuF7CxcSO58umKdzw7zlh2jJHMCDF/jBu7buTatmsJeoMA5Eo5+hJ9/PDED/nG4W8wkB6YXf61bdfy6vWvZiA9wLeOfuuMBuCqyCpe2P1CtrdsZ0PjBtY3rifuj5/zG4wxDGeGOTZ9jIHUAP2pfkYyI7ONjZle4cw+lS/n6Yn20Bu3ic4Rh0QhQaqYolAuEPaFZ9eDa1zy5Tz5cp50Mc1Ufoqp3BQT+Qmmc9NM5Cdm98+YP0ZHuIPWUCthb3i24RL2hYn6okR8EWL+GK2hVtpCbcQDcX568qd8+eCXOZE8Me++0hnpZHV8Ne2hdvyOf7YszNQXQW+QXCk3W2YmchP4Hf9s2QUolAvky3lKbomAN0DQsY2q9nA7PdEeemI9NAYayZQyJAvJ2UaT1+PF6/GSKWZm1+tAeoAPXvdBrmq7at7ydLFqNkGJiB8YALYZY4YvNO1K96D2ju3lzx77M3aP7D7j/eZgMx+47gPcvfHuJTl0VSwXbevuIg4blN0yo9lRBlIDJAtJ2sPtdEY6aQw0YjAkC0kmc5MYDD3RnjOSijGG8dw4QSdI1B9ddNwLUXSL7B7ezXBmeHaHj/ljrIuvO+fQEdjDDdU8BLoSim6Ro1NHifgi9MR6nnP6i10nxXIRF9tQMcYwkhlh3/g+nhl/huOJ47QEW+iMdNIV6aIp2GQrr0qZ6In1nNH7mk/ZLZ/3MOVCucZl19Au9k/s58U9L2Ztw9ozlvPY8GOcSJzghs4bWBtfWxOHPS/EGEO6mMYjnnnL9cVwjcvPB37O3rG9tIfb6Yh00BnppDvaTcAJLHHEtaGWE9TrgPcZY+56rmlXKkGNZEb4+OMf59+P/jvNwWbec9V76I52z3aNr2y+csUreKWUulydL0HVwmjmbwH+qdpBzPXbP/lt9o/v593b3817rnqPJiOllKqCqiYoEYkAdwK/Xs045koVUuwZ28M9V9/D+659X7XDUUqpulXVBGWMSQMt1YzhbE+PPo1rXK5rv67aoSilVF27vM9CX4InRp7AIx6ubru62qEopVRd0wR1lidHnmRL0xYivki1Q1FKqbqmCWqOolvk6bGn2dG+o9qhKKVU3dMENcezE8+SLWXZ0aEJSimlqk0T1BxPjDwBwI42TVBKKVVtmqDm2D2ym+5oNx2RjmqHopRSdU8TVIUxht0ju7m2/dpqh6KUUgpNULNOJU8xlh3T+5+UUqpGaIKq2D1qB4TVHpRSStUGTVAVTww/QcwXY2PjxmqHopRSCk1Qs54ceZJr26+97B/xoJRSzxdaG2Mfrndk+ojeoKuUUjVEExRwYOIAANtbt1c5EqWUUjM0QQFj2TEAvf9JKaVqiCYoTieo1lBrlSNRSik1QxMUMJ4bx+fxEfPFqh2KUkqpCk1QwHh2nNZQKyJS7VCUUkpVaILCJqiWYE092FcppeqeJijsOSg9/6SUUrVFExT2HFRLSHtQSilVS+o+QZXdMhO5CU1QSilVY+o+QU3lp3CNq+eglFKqxtR9gpq5B0p7UEopVVvqPkGN58YBvUlXKaVqjSaorE1QeohPKaVqS1UTlIg0ishXROSAiOwXkZtXOoaZBKU9KKWUqi3eKi//r4DvGmPeKCJ+ILzSAYxlxwg4ASK+yEovWiml1AVULUGJSANwG/AOAGNMASisdBzjOR3mSCmlalE1D/GtA0aBvxOR3SLyWRE5pxsjIveIyC4R2TU6OrrkQegwR0opVZuqmaC8wHXAXxtjdgBp4HfPnsgY8xljzE5jzM62trYlD2IsN6aXmCulVA2qZoI6BZwyxjxS+fsr2IS1osazOsyRUkrVoqolKGPMEHBSRLZU3noZsG8lYyi5JSZzk3oFn1JK1aBqX8X3fuBLlSv4jgLvXMmFT+WnMBg9B6WUUjWoqgnKGPMksLNay9dHvSulVO2q65EkZkeR0HNQSilVc+o6Qc32oILag1JKqVpT1wlqZqBY7UEppVTtqesENZYdI+QNEfat+AhLSimlnkNdJygdRUIppWqXJig9vKeUUjWpvhNUZaBYpZRStaeuE9RYdkwP8SmlVI2q2wRVdItM5ae0B6WUUjWqbhPURHYC0EvMlVKqVtVtgtJ7oJRSqrbVbYKaGUWiJdgCxsDJx8B1qxyVUkqpGXWboGbG4WsNtcLhH8Dn7oAHP17doJRSSs2q3wQ19xDfgXvtmz/+Uxh8qopRKaWUmlG/CSo7TsQXIeQJwMHvwvrbIdwCX70Hitlqh6eUUnWvbhPU7D1Qg7shNQTXvAXu/hSMHoAf/HG1w1NKqbpXtwlqMjdJc7DZ9p7EA5vugo13wI33wCN/DUd+VO0QlVKqrtVtgprKT9EYbISD34HVN0O42X5wxx9D6xb413fA4NPVDFEppepa3SaoyfwkTeKD4T2w+RWnP/CH4a3/Cv4o/MPdMHqwajEqpVQ9q8sEZYxhKjdFY9qOJsGWV505QdMa+NVvgjjwxdfCxNGVD1Ippeqct9oBVEO2lKXgFmicOA4tm6B147kTtW6EX/0GfOHV8LmXQ89OiLRBtB18YXveSjwQarI9sGjbiv8OpZS6nNVlgprMTwLQNHoErnr7+Sfs2GqT1I8+AlMn4NQuyIyBOWvECXHsZepXvRG2vQF8weULXiml6kRdJqip3BQAjeXiuYf3ztZ1tT0nNcN1wS2CW7aJavI47P032PsV+Pp74bHPwlu+rD0qpZRapLo8BzXbg/KGoefGhX3Z4wFvwF5MEYhC53a444/gA0/Dm74Aw/vgsy+DsUNLH7hSStWR+kxQOZugGps2grNEnUgR2PZ6eMe3oJiBz94B+++F1IgdjFYptXxKeTj1OBRz1Y5ELaGqHuITkeNAEigDJWPMzpVY7nR+GoCmSMfSz7znenjPD+BLb4Ivv9W+5wtD01rY8TbY+e7L5xxVdhKG9sDwM5BP2UOepgzBBuh9AXRdA45v+ZafGoGjP7G3BLRtsevY4yzNvF0Xpk/C+CFo3QyNq8/8fOwQPPoZOyxWsAGCjRBphZaN9hXrtI2WBS2zDPkkhBqX5jcslXLJbldv4NzPpk7C9CnwhcAfsWXdGwSv3/6bT9rzt1MnIDsBvTfZbbXQdXNOTEVb/vqfgGe+Bge/DfmEvZDphl+DG95tt8fzhetCIQWB2OLXzWIUs5AZBwTiq6obC7VxDuolxpixlVzgZH4SjzHEoquWZwFNa+HXfgzHfwZTfTDZBwNPwPd+Dx76JLz4v9pktZDKu1SwlaU/Ao1rThccY2D8CPQ/bntuHsdetOENnK44gw32nFiw0X7PLcPIfjjxc7uDT5+E5CAkBmxMjWvspfaxVVAu2NZpKWsLbyFtl5MahekTF47ZF4aeG+wwUtvfcG4FZ4zdGcYOwfhh+9t6boCGHhtndgqO/NCO6uGWbaUf67JxHPwOnHwEmNM7dQLQthk6roLOq+z/U6N23uOH7W+JddrfFW62y04M2FchbedlKhXF+BH7OwEQe6XmDe+x2/aBv4CnvwyO367TfGLOtHN+e6gZgnEIxM/6NwbekG2oOAFbwQ/stgMVF9P2xvGr3wxb7z59A/mM9Dg8+mnY9Xc2KbRfaSv8ju3Qe+OZZWOhMhO2HJ18FPp32bjSYzYRiNhE3Xm1Xeb4ETj+gC3fC9WwGjbdaRsw4Rb7cnx2mLHhfTC635aFlk022QditrwO77X/pkag0sgE7Da48rWw9habrH7yp/Czj9lGkjdo5+0LQUMvNK+DpnW2jMU67ftgy2I+aeednYTctF2GMXbahl47fTELqWG7v6SG7fpJj9kyEO+226J1k90OZzeWjLHzH3zK1gf9T9jfnJuCXAIwdpCAG3/N7jOBqP2e60LilK1Hpk/aRoHjtSPfdF59ensnh+02SY+dLm+BWOWKY7H/lotQytnfkZuy+97Ys/bf1PCZ5TjWBatvsusx2nG6AeL4bJ1QLkI5D907IbYMjX1ATBUPP1V6UDsvNkHt3LnT7Nq1a9HL/ciDf8gPDn6Fn27+Nbj1g4ue30U7dj/88CNw6lFbeDq22VfzetvCHNlnd0C3bAt7fJUtpKMHbUF2S3Y+/pitJEKN9srCytOBn5MTsAVtZucD+3fTWrus2Cpb4Cb77MUfqWFbCftCdkef20oONdok0Hm1fYUabWL0eCA5ZJNf389tghk/DJF2W8F3XlXZOR+HgSfnjz3WZeMZeNK23ENNtpeUHLIXqIBd5hWvgc0vtzvK6AH7GtkHQ3shPXJ6fuLY3+gL2WQ0d5nhVoh32e0xsyN7g7ZibNtit82x++HxL56epzdkW+i3/Nbpi2HKRVtpjR+xv3fimK3o8onK+k5W/p+w/y/nT8fgDdrfs2qHbUzs+waMHQSPz57jbNloK+vMGDzxD7axsPkVdjuMHrQNl3Lh9Pbs3mm3h8drKxO3VFl+yibfQtpWUMWM/XemwjLlyvryQPs2aFlv10+k1VauQ3ts5ZocsMl3zQth7YtshVzKQSFjE2wpX6nA8nZdNa2xPVB/1K7LQ9+3Pd9i+txt7w3ZhkUhA5PHTpd5xG6L9ivtvhFuscm7eb2Nwes/PY+RA3a4suF9dr2Ui/Z3J/rnzK8i2Gi3fWbs3EbG2cRz7hW8M7H5I3YZMxy/LXPNG+ytKRNHbdnMjJ/+Tutmu33DrXa7+4Kw75sw+KSNae2ttl4YP2zX73yinbDmZltvjB64cPzz8QZP32ozu15b7PY7+QiceNgmxwt5yz/DllcufNlziMjj8x1Bq3aCOgZMYpvBnzbGfGaeae4B7gFYvXr19X19l9BiO8tv3/frHO77Kd+84Q/h2l9e9PwWxBi7gx76nq1Ih5+BQhJ8EbvztV9hC3diwO5QuWnbqurcblvJ+aQt6MP7bEW76jrbcu690e5spmx3wlLefjc3bXsi6ZFKy2/Y7girb7avxtXL2403xvaAHvkbOHSffU8cewn/qh3QdoXdUVs22FhPPganHrM75poX2oq4Z6dtjbqu/c1u+blbbKkR2yqMdtgKcm5vtZiz8wk1X/zh1lIBDvy7TTw7fmXxLUbXtRV4MWsro7nnQo2Boaft1aFDe2DssG05exy4+hfhlg/Y5DmjXLK9jpOP2N7PwG5bwbtFWxbEsS3pQNQmCX/EJmtfuNL4qPTmgg12m6y67nTrfT7ZKRuzZxGnsMtFWx4z4/ZVyttyMPcwbbloG0v5hP3sQjFd1DJLtrKdOGb3r+SgbfTMHBqMdthXuPl0bxexPcmpPvudQMw2oGKdNvGEW+30Hsf2QMcO2cbF+GGblMaP2sGoZ5Jr+1bbSOu6xs7rbMbY8v/Ip21joHm9bQC0bLTrpnG1TST5hH2O3aH74MQjtt5Y92JY/2LbQ80nKq+ZQ++Vl+M7vb0DMYj3PPd2TA7ZxtZMA6RctEdDnICdX/M6W3YWoVYTVLcxpl9E2oHvA+83xtx/vumXqgf1rm+8kfLQU3zxZX9tu8nVNHOYK9S8uB3++WD8iD380HmVvQpSXbxi1vYGFlkRKFWLzpegqlojGmP6K/+OAF8DFnjN96WZzE/RWHZta6naROwhlMs9OYHtJa1+gSanS+ELaXJSdadqtaKIREQkNvN/4C5g70ose6qQpMmtkQSllFJqXtW8iq8D+JrY8x9e4B+NMd9d7oUaY5gq52wPKtyy3ItTSil1iaqWoIwxR4FrVnq56WKaEq4dRWKp7plRSim15OrgxMeZZoY5avTHqxyJUkqpC6m7BDUzUGxTsKm6gSillLqguktQMz2ohpCONq6UUrWs7hLUVNYmKL+/lcMjqeeYWimlVLXUXYKaTA0A8PBxL2/4Pw8yktDRj5VSqhbVXYKaSg3gGMOxdJxErsR///peqjmahlJKqfnVX4LKjNBYdjmYDNPVEOS+fcN8Z+9QtcNSSil1lvpLUNkJmtwyw6aB375zM1d1N/CH39jLZLpQ7dCUUkrNUXcJamYcvlHTyKaOGB/9hauZyhT5yLf2VTs0pZRSc9RdgpoqpWlwIU2Qda0Rtq6K897bN/DVJ/r58cGR556BUkqpFVF3CWqynCWMn9ZogIaQfUbQb7x0Ixvbo/z+V/eQypeeYw5KKaVWQl0lKNe4TJsS3nKAda2R2fcDXoeP/sLVDCZy/Nl3L+GplEoppZZcXSWoZCFJGTCl0BkJCuD6NU2884Xr+Puf9/HosYt8hLpSSqllU1cJaio/BUC+EGZd67mPj/6dl2+mpynE7/7b0+SK5RWOTiml1Fz1laAyowAUSg2sb4uc83nY7+V/veFqjo6l+asfHlrp8JRSSs1RXwlq6jgA2XIT61vPTVAAt25q5c07e/jM/UfZ2z+9gtEppZSaq64S1GTiFADpcjOrW8Lnne73X7WVprCf3/3q05TK7kqFp5RSao66SlBT6UEAvOHVBLznf5puQ9jH/3jdNvb2J/j8g8dWKjyllFJz1FWCmsyM4jWGhqa1zzntK7d3ctfWDj72/WfpG08vf3BKKaXOUF8JKjtBU7lMS0fPc04rIvyP123H5/Hwoa/u0RHPlVJqhdVVghrPTxMrw9r2houavrMhyH95xRYeOjLOrr7JZY5OKaXUXJecoERkx1IGshImi2lCrnPOTboX8obrevB7PXxnjz6SQymlVtJielBfEpH42W+KyNZFzHNZTbt5AiXfghJUNODltk1tfGfvIK6rh/mUUmqlLCZBfRr40tw3ROQO4CeLCWg5TVPG6wZZ1RBa0PdedVUng9M5njo1tTyBKaWUOsclJyhjzF8BeRH5EICIvAv4R+CtC5mPiDgisltE7r3UWC5GuVwi6YGAJ4bHIwv67suu7MDniD55VymlVtCCEpSI/LKIbBWRme+9C3iHiHwJ+EPgJcaY7y8whg8A+xf4nQVLJvtxRYj4mhb83YaQj1s3tvLtPYN6NZ9SSq2Qhfag3g88CqREZBfwl8DDwEuAu4wxzyxkZiLSA7wa+OwC41iw8cljeIwhHmq/pO+/8qouTk1m2dufWOLIlFJKzWdBCcoYczMQA3YAfw6MAO2AAPtE5JCI/MsCZvlx4L8C5x1PSETuEZFdIrJrdHR0IeGewRO7gfSBj7Bly29e0vfv2tqB1yN8e+/gJceglFLq4i34HJSxDhpjvmyM+ZAx5pXGmC6gG9vDevxi5iMirwFGjDEXnN4Y8xljzE5jzM62traFhjuruzHE19//Ul62bfUlfb8x7OfmDS18Rw/zKaXUinjOBCUi/7+IvFVErhSRc64uEJFuEQkZY4aNMd81xnz0Ipd9C/BaETkO/DPwUhH5vwuKfgH8Xg/buxtoiQYueR6v3N7F8fEM+weTSxiZUkqp+VxMD+qDwN8De4GkiDwkIp8UkXeLyHXAbwH3L3TBld5XjzFmLfBLwI+MMW9b6HxW0l3bOvAIfHuPHuZTSqnl5r2IaZqB6856vRd73mnmWFdmWaKrMa3RADdvaOHepwf4z3dtZp4OpVJKqSXynD0oY8yUMeZHxpi/MMb8sjHmCqAReAvwLDAB/NpigjDG/MQY85rFzGOlvPaaVRwfz7BHH2aolFLL6pJu1DXGJI0xXwZuANJA65JGtQwKp/o5fOddJO67b1HzecW2LnyO8M0nB5YoMqWUUvNZ1GjmxpgU9vzUB5cmnOXjCfgpnjxJaWxsUfNpCPt48eY27n1ax+ZTSqnltBSP2xgFupZgPsvKE7fj2rqJxV+B9x+uWcVQIsdjxycWPS+llFLzu5jLzEdF5D4R+V8i8mYR2TjnM8GOBLFrOYNcCp5AAPH7KScXPxLEnVs7CPkcvvmUHuZTSqnlcjFX8T0IXAvcUfnbiEgSOAC0VV5vFxGfMaa4LFEuEU9DfEl6UGG/l5dd2c639wzy4dduw+fU1XMflVJqRVzMVXx3V+5VagZeCvwO8A0gBPQCUeAr2PH59ojIP4rIf1u+kC+dE4tTTi7NTbavvWYVk5kiDx5e3DktpZRS87uYHhRgLzfHPuvpJzPviYgf2I7tYV2LHaPv1cAvAhc7osSKcWIx3MTSXB7+4i1txINevvnUALdvubQBaJVSSp3fRSeo+RhjCsATldcsEdmwmPkuF088TnlycknmFfA6vGJ7J9/eM0SuWCboc5ZkvkoppaxlOXlijDmyHPNdLCceX5KLJGa8+upVpPIlfnZID/MppdRSq6uz+554bEkukphx8/oW4kGvPmlXKaWWQV0lKCcWp5xILNnjMvxeD3ds7eAH+4cpls/7SCullFKXoL4SVDwG5TIms3Rj275yexfT2SI/PzK+ZPNUSilVZwnKE7OjSSzVpeYAL9rUSsTv8B190q5SSi2pukpQTkMlQSWW7kKJoM/hJVe0c98zw5R1bD6llFoydZWgPLEYAO4SJiiwh/nG0wUePaZj8yml1FKpqwTlxGd6UEv7yPbbt7QR8Hr4rh7mU0qpJVNfCWqmB7WE90IBRAJeXry5je8+M6SP4FBKqSVSVwnK09AALH0PCuCVV3UynMiz++TUks9bKaXqUV0lKCcaBVjS0SRmvOzKDvxeD195/NSSz1sppepRXSUo8fmQcBh3eukTVDzo4xeu6+arT5xiLJVf8vkrpVS9qasEBfY81FLeBzXXe160nnzJ5R9+3rcs81dKqXpSfwkqHl/yiyRmbGiLcseVHfz9z4+TLZSXZRlKKVUv6i5BeeLxZblIYsY9t61nMlPkK0/ouSillFqMuktQTiy2pCNJnO2GtU1c09vI5x44qiNLKKXUItRdgrKP3Fi+BCUi/Ppt6zk+nuH7+/QxHEopdamqlqBEJCgij4rIUyLyjIj88Uos14k3LNtFEjNevq2T3uYQn/jhYVL50rIuSymlLlfV7EHlgZcaY64BrgVeISI3LfdCnXgMN5nEuMv3/CbHI/zeK6/k4HCSN/3Nzxmazi3bspRS6nJVtQRlrFTlT1/ltewnbTyxOBiDm0o998SL8Mqruvj8O27g5ESGuz/1IPsGlu+wolJKXY681Vy4iDjA48BG4FPGmEfmmeYe4B6A1atXL3qZTtyOx1dOJGcHj10uL97cxr/+x5t51xce401/8xDve+lGfvnG1TSG/cu63MtNrljG8Qg+Z/HtKWMMk5kiw4kcqXyJVL5EOl9iaDrH0bE0R0dTJLIlXryljddc3cXWrjgiMm9MDx0ZYyJdpFh2KZVdYkEfL1jfTFdDaN5ll8ouJyYyDE7n6GwI0tsUxu9d+G8qlV2Oj6fpbQ4T8DoL/v6MwyMpRhI5wgEvEb9DQ8hHWyww7++dYYzBNZDKlZjOFpnOFmkM++htDl9yHACFkstEusBYKs94uoAxhoDXIejzEAl4aYsGaAz7ZmNzXUMyV8JxhGigqtXYiiiUXHyOXHDbXCpjDPmSy8BUluPjaY6PZciXXG5c18TVPY1Lst9dKlmqx58vKgiRRuBrwPuNMXvPN93OnTvNrl27FrWsxPe/T//7f5N1X/sqwSuvXNS8LtbQdI7/8pWneODQGCGfw5t39vDybZ0EfA5+x0PI72FtSwTvBQrCcCLH432TCLCqMUR3U4jmsJ+i65IruhRKLn6vh4jfwet4KLvmdIEbz3BqMsOpySynJrO4rqGnKURPU4jOhhCOgGvANYamsJ/1bRHWt0WJBbycnMywbyDB/sEEQ4kcE+kik5kCuWKZ5oiftmiA1liANS1hNnfE2NQeJeR3ODaW5tBwipOTGda3RrlhbRMt0cDs+vjZ4TF2n5gkWyhTdA3FksvG9ii/eEPvbGU3nS3y2QeO8vmfHaNYNmxoj3JFZ4zOhiBjyTwjyTyjyTwtUT/rW23MLVE/6XyJVL5MMldkMl1gLF1gPJVnOJFncDpLrjj/4d2msI/1bVH8jodHj09Qdg3rWyPcvKGFK7riXNlpGzdf3d3PvU8NkMjNf35xXWuEG9Y24XiERK5EMldisLItiuXT+5tH7LZc0xJmdXOY3uYwrdEAqVyJqWyRRLaI1yPEQz4aQj4yhTKPHhtn1/FJkvkSXQ1B/tPtG3jzDb3nTVQPHx3nxESGrV1xNrZH8XqE7+8b5gsPHeeReR4P0xLxs727ge3dcbIFl8OjKY6MpBhK5HCNYb7qQgTuvLKDX3/xBq5f03TO5yOJHD/YP8Kjx8YpuQYRQYBErsjQdI7hRI7JTHHe+OfyOx5aon7yJZepTIGZi2S7G0Nc2RVjQ1uU6WyR/qksA1NZwn4vO9c2cePaZq5f0/ScyXcuYwxPn5rmwFCCkN9LNOAQ8XvxeARj7OdT2SJHRlMcGUnTP5VhXWuUHb2NXLu6kZ6mEIIgAiXXMJrMMzSdYySZo+waIgEvEb+XcMAhWEnEAZ/DZLrAwFSW/qksJyYyHBlNc2QkRf9UFq9HaAz7aAr7Wdsa4eXbOrnjynYaw36yhTIPHBrl+/uGSeSKrG2NsK4lQk9TmFKljsiXypyazHJ4JMWhkSSnJrNkC2XypfOf7gj7HW5Y28zVPQ1s6YxxRWeMkN/LoeEkh4btfN7zovVs7ohd1Ho9HxF53Biz85z3ayFBAYjIHwIZY8xfnG+apUhQ6Ycf4cQ73sHqL36RyAtuXNS8FmrfQILP/ewY33yq/4yKCiDkc7i2t3F2R5pp3Q8ncuw6PsmJiYt/TH3A68E15oxl+B0P3ZWkJCL0VxLWhQqn3+uhUPncI9AWC9AU9tMS9RPwOoynC4wl84ym8rPTga2w5itWG9oiiAiHR+zh1XjQS0PYh8/x4IhwZDSFAV60qY3tq+L834f7SORKvPqqLnqbwxwYSnBgMMloKk9r1E9HPEhLxM9YqsDR0RTpeW6Obgj5aIn6aY0EaIsF6GoIsqoxRGdDkFjQO1tRtMcCNEVO92wn0gW+u3eIb+8Z5KlTUyTnJKOgz8MrtnXy+ut6WN8awed48DrCcCLHz4+M8/DRcZ44MYXjEWJBL7Ggj7ZogI3tUTa0RehuDDGczHFsLMPxsTQnJjKcnMgwni6csQ5jAS8l15CZ87vWt0W4aX0LW7vifH13P7v6JumMB3nXrWt5+bZO1rREADg2luZPvrWPH+wfmf3uTDxTmSLdjSF+5eY1XNPTSLZYIp0vM57K88xAgj390xwaSeF3PGxoj7ChLcqqxhA+j23Be0SIBBwaw34aQj72nJriiz/vYzpbZMfqRta1RHA8guMRDg4n2X1iCoCOeIBIwIupNIZiQS+d8SDt8SAdsSCtMT8tkQAtUT8eEfJFW3km8yVGk3lGkjnGkgWCPg9NYT+NYR/5ksvBoSQHhhIcHU3TGPazqjFIV0OQ6WyR3SemZsu43+uhIx6gMx4k5PcyU/d5PcKqxhA9TWFWNQZ58uQU39s7xMBFnjvuiAfoaghxZDR1RjlZrJDPmV3/a1silF3DRKbARKrA06emGJjO4fUI21bFOTicJFd0iQW9tMUCnJzInFPHzFjVEGRjR4w1zWHCfoeAzyHg9dAZD7K2NczaFrufPnx0fLY8HxlNMd9dM61RP3/+pmt4yZb2Rf3WmktQItIGFI0xUyISAu4DPmqMufd831mKBJXbt49jb/gFej75v4ndccei5nWpRpN5Dg0nZ3sOyXyRp05Os6tvgv2Dydn7p7weoSniZ0dvIzeua+aGtc14HWFgKkf/ZIaJTJGA1zP7ypdc0vky6UIJjwhrW8KsaYmwtjVMRyyIx3Nm69EYw3S2iDHgEQGBsVSeo6NpjoymGEvm2dgeZeuqOJs7YgR987fSXdcwMJ2dbVGl8mU2tkfZ1B6luynEs0NJHj0+wWPHJnAN3LKxhVs3tnFFZ+yMmAamsvzLrpN8+bGTDE7neNkV7fz2XZvZtqrhnLjPbgkbY1upk5ki0aCXaOWw1YV6pRfLGMPAdI4DgwnShTIv2dJGLOhb9HzPls6XGEvliQd9xEM+nMq6KZRckrkiIkLznCRqjOGhI+N8/AfP8tjxSQA2VbbXt/cM4nc8/MZLN3Hn1g6eHU6ybyDBwHS20vLumJ3/fAolF69HzikzF4r9X3ad5F93nSKVL1F2DSXXpTMe5M6tHdy5tZPNHdFlOUQ1Y75yUSi57Omf5qmTUwwncgwlcgxN52aTloidpn8qy1SlF+f3erhtUysv39bJjeuaKZRcUvkSmUIZ15jZnlE04GV9W2S2LLiu4ehYmidPTjGWytueFgZHhLZYoJKMA3g9ntlDy5lCmXypPNvDaQj56G60ibI54j/v+prp4X1n7xCPHBvn6u4G7tzayQvWN+ObcwRlYCqL1/EQ9HkI+hzaY4FLKru5YpnDIykODCXJFctsao+yqSN2RnlcjFpMUFcDXwQc7MUa/2KM+R8X+s5SJKjCqVMcueNOuv7kT2j8hTcsal7LIVMokS2UiQS8BLyeZd2ha1XZNYyn8rTHg9UO5XmjbzzND/eP8MMDw+w+McWrr+riv7xiC+0xXYcXK5krMjido6shuCwNEHV+50tQVTu7aIx5Gtix0sudfaruMo3Ht1hhv5ew//I/6Xshjkc0OS3QmpYI77p1He+6dV21Q3neigV9mphqTP2NJBGNggjuMo7Hp5RSavHqLkGJx4MnGl320SSUUkotTt0lKLADxrqJ6WqHoZRS6gLqMkF5GhqW9ZEbSimlFq8uE5R9qm5tXiShlFLKqssEZR+5oT0opZSqZXWZoJxYfFkfWqiUUmrx6jNBxePL+tBCpZRSi1eXCcoTj+FmMpiSPkxQKaVqVV0mKCc2M5qEnodSSqlaVZcJylN5JpSrCUoppWpWXSYoJ25Hxy5P63kopZSqVXWaoGZ6UJqglFKqVtVlgvLMnIPSe6GUUqpm1WWCmulB6WgSSilVu+oyQc30oPReKKWUql31maAiYXAcPcSnlFI1rC4TlIjYR27oIT6llKpZdZmgADzxuPaglFKqhtVtgnJiMcr60EKllKpZ9ZugGuL6yA2llKphdZugvF1d5I8d0wFjlVKqRtVtgoreeivu9DTZp/dUOxSllFLzqNsEFXnhC8FxSN3/02qHopRSah51m6CchgZC115L+v4Hqh2KUkqpeVQtQYlIr4j8WET2icgzIvKBlY4hettt5PbtozgystKLVkop9Ryq2YMqAf/ZGLMVuAl4n4hsXckAoi++DYD0Az9bycUqpZS6CFVLUMaYQWPME5X/J4H9QPdKxhDYsgVvezupB/Qwn1JK1ZqaOAclImuBHcAj83x2j4jsEpFdo6OjS71cIre9iPSDD+rl5kopVWOqnqBEJAr8G/BbxphzBsczxnzGGLPTGLOzra1tyZcfve023GSS7JNPLvm8lVJKXbqqJigR8WGT05eMMV+tRgyRF74QvF5SP72/GotXSil1HtW8ik+AzwH7jTEfq1YcTjRK+LrrSN2vCUoppWpJNXtQtwC/ArxURJ6svF5VjUCiL76N/MGDFE6dqsbilVJKzaOaV/H9zBgjxpirjTHXVl7frkYssbvuQoJBTr3vNyhPTVUjBKWUUmep+kUStcDf20vPJz9J4ehRTtzz65RT6WqHpJRSdU8TVEX01lvo/vhfknvmGU6997242Wy1Q1JKqbqmCWqO2MtexqqPfpTMrl2cet/7cNPak1JKqWrRBHWWhte8mq7/+aekH36Evne9S89JKaVUlWiCmkfj3XfT84m/Ir9vP32/8qsUh3UwWaWUWmmaoM4jdscd9P7tZyj299P31reSP3Kk2iEppVRd0QR1AZGbbmL1F7+Am81y/Bd/idRP9eGGSim1UjRBPYfQVVex7l//BV9vLyf/43sZ/9znMcZUOyyllLrsaYK6CL5Vq1j7pf9L7OUvZ+TP/5zRT3yi2iEppdRlTxPURfKEw3T/5cdouPtuxj/9GbLPPFPtkJRS6rKmCWoBRISO3/sQTkszg//9DzDFYrVDUkqpy5YmqAVy4nE6/+APyO/fz/gXvlDtcJRS6rKlCeoSxO+6i9iddzD2yU9ROH682uEopdRlSRPUJer473+A+P0M/uEf6VV9Sim1DDRBXSJfRzttH/wtMo8+Svqhh6odjlJKXXY0QS1C4xvfiLe9nfHP/G21Q1FKqcuOJqhF8Pj9NL/znWQeeYTsU09VOxyllLqsaIJapMY3vQlPQwNj2otSSqklpQlqkZxohOa3vpXUD39I/vDhaoejlFKXDU1QS6DpV96GhEKM/+1nqx2KUkpdNjRBLQFvUxNNb34T09/6FsX+/mqHo5RSlwVNUEuk+Z3vREQY/rM/1/uilFJqCWiCWiK+zk5af/P9JL/3PRLf/Ga1w1FKqec9TVBLqOVd7yJ0/fUMfeT/00N9Sim1SFVLUCLyeREZEZG91YphqYnjsOqj/wtcl4Hf/RDGdasdklJKPW9Vswf1BeAVVVz+svD39NDx+79P5rHHGP/M3+r5KKWUukRVS1DGmPuBiWotfzk1vOH1xO68k9GPf5zjb3ozifvu096UUkotkFSzhS8ia4F7jTHbLzDNPcA9AKtXr76+r69vhaJbHFMoMPX1rzP+2c9RPHEC/7p1hHbswNfZgbezk8CGDYSuugrx+6sd6oKZchk3m8MT8CM+34osD48HEVmS+ZUmJ8kfOIAEgzgNjTgNcdxUisKJkxROnsBNpoi88GaC27YhnqVtwxljMIUCnkBgSee7GMYYyuPjlMbG8K9ZgycUmn+6YpHCiRMUjh0DQEIhPKEw3rZWfD09F719yqk0xf5+vK0teFtazomlcPQoplDAaWzEaWw8bzxq6bmFArm9z1A40Yc4XsTnQwJ+Qtdcg7e5edmWKyKPG2N2nvN+rSeouXbu3Gl27dq1vEEtMVMqkfje95j65y9T6OujNDoKlXUuwSDh63YQuv56fN3deFvb8La1gutSmpigPDFBeWoaN5fF5AuYfM7O1HEQj4PTECf8ghcQ2Lz5nIrUTafJHzpE7tlnKRw5ipvNYsolKJUxxSJuIY8pFKDs4m1rw9vVia+rC19HB05rK97WVsTvJ7f3GbJPPUX26acoDQ5SGp+gPDk5+xvwevEEgwS3biVyyy1EbrkF36oucs/sI/fMM+SffRYJBvE2N+E0NWPcMsWBAUoDg5TGxhC/H08kgiccxtvRQWDDBgKbNuLE46QfeZT0Aw+QfvRRxO8neOWVBK+8El9XJ8WBQYr9/RQHB/EEgzhtrXhb2/BEI1As4hYKUCoh/gASCuIJhSkNDZJ5bBf5Q4cuatt529uJvuylBLdcgdPUhNPUCK5Lbu9esnv2ktu/H1MsIt7Kjux4gEolLYInFsOJx3HicdxCnmLfCQonT+Imk3jb2ghs2kRg00a8HZ14wmE8kTDi8+PmsriZDCaTwc3nMaUSlEoY17WVhteL+Lx2vcXiOA1xfN3dthxUkoQxhuwTTzDxhS9SHB4msG4d/g0b8HWvwk0kKI2OURodpXD8OPnDh+02BfB48K9bR3DLFiQYxE2lcFMpSqOj5I8fh/M8RdppbCR49VUEt23DEwyBcTHlMiabpTQxSXligtL4OMX+fsoTE7PrKHj1VcRe8hICmzeTfujnpH70I4oDA2fM2xMO41u9Gv/q1fh6esB1KScTuMkU4nXwr12Hf8N6/L29lMbGKBw7bp/TJmLX8ebNBDasR4IhxO9DvF4QsWXYdTHFoo2vss9JIIi3tQWnpQVTLJJ57DEyjz5G9sknQQQnHscTj+FEY7Nl1xON4u1ox9e1Ct+qLtxk0u43Tz5ly5uILSM+H962Vvzr1uNftw5vWxulkRGKAwMUhwbx+AN429vxdnTgCQZsOR8YoDg4SHlykvL0NOVEAvF6CWzZTPCKKwlsWE95eppi/wDFgQFMuYy3uRmnpRknFsdNpyhPJyhPT2Pyedvgc8sYY/AEgkgoiPh8FA4dJrtnDyafP3cDezyEb7iB2MvvIrRtG+VkCjeZoDydIPqiW/F1d1/UPnU+mqBqhCkWKY2MkNu/n/Sjj5J55FHyBw9e1HfF5wMRW8DK5dn3nZYWIi94AWBmC3RpZOT090IhPJEI4ji2cvN6kUDA9t48Hkqjo5SGh+F8hyFFCGzciH/tWpzmZrwtzXgiUUwhj5vN4aZSZHbvJr9//zlf9fX0zFYAplK5OQ0NeFetwtvWapNlOoObTlMcHMRkMmd+f81qorfcgimWyO3fT/7ZZzGFAhIM4uvuxtfZicnnKY2NURobw02nEb/t2YnXi1soYLJZMAZPOEzouusI79xJ6JqrMaUS5alpyolpPMEQ/jWr8fWuRvw+0vffT/KHPyL1wAP2+/P8ruC2bXjCYUyphCkVoXR6mxjXxU2lKCcSlKenEJ8Pf+9q/Kt7cVpbKZ44Sf7wYfKHD2NyuefY8GIrVY/HbvtSad7JnLZWorfcSnDbNhL33kv2qadwGhoIXHEFhWPHzigTiOA0N+Pv6SGweROBTZvwtraSP3yE3IED5A8cwLgunmgEJxLFaWoisHED/o0bCWzYgDgObjaLm8lSHBggu+dpck89bYf7mlOniM+H09xsy01Tk91mvb34uldR6Osj9eOfkNuzx04bCBB54QuJvuR2nMZGylNTlKenKY2MUjxxgsKJExRPnUJ8vkryj+HmCxRPnTqn7DrNzVAuU56evvC6vUieSITQjh2I3085MY2bSFJOJWfL7nm3SUsLwSuvRBxntmFYGh6xMZ9V93oaGmwCOas8SDCIr6vLJpx4A05DA242S/7AAQp9fafn4/Ph6+xEvF7buJ3z2yUctok1GLQNXMcBwOTzuLkcbi6Hv7eX8PXXE9p5PcFNmzCuqeyfKVIPPEDye/dROHr0nN/Y/b8/QfzOOxezejVB1TI3nZ6tYEujYyBiW3BNzThNjXiCQZtQ5vSSjDGURkZIP/Rz0g8+SOaxx5BAAN+qVbYn1NtDcMsWAps34+vufs5DVaZUojQyYl/j45RGx3CzWdtr2b4dJxp5zt9RGh0l/dBDlMYnCG7dSnDbVpxYbDZeN51GRPBE5p+XcV1Kg4PkDx+mND5B+Prr8K9Zc06c5UQCp6npog8pGWMw+fxscl4IUyzaHsDUJOXJSUy5THDrVrxNTQuaz3nn77q4mYyt6DJpTKGIJxyyrfJQyG73SmUy9/dQKuGm0zYBJpLkn32W9M9+RvrBBylPT+NbvZrmd7ydxrvvxhMOA1BOpSgODOA0NOJtaV7wurio31MoYIyx5a3yeq7tVBweoXDsGKFrrr6kw3luPk/heB/FUyfxtrfjX7MGJx63+8joKPlnD1E4ftxW/qWSbSgZAx6pJH+f7eE3t+A0NWLyBUpjY5THxzCuIbzzeptkLrC+3FyO0tAQxUHbQJRgkNA11+LrXjXv73fzeQp9fZTHxvB2dODr6rKNHWNsj3V4GDeXx7eq64Jl3c1kKJw4gdPUZI96zCkrpliknErhRCJLdiohf/gwhRMncRrilZ5kA96mxkXPv+YSlIj8E3A70AoMA39kjPnchb5zuSYopZaKKZcpnDiBf/XqcxKbUrXqfAlq6ZtQF8kY85ZqLVupy5U4DoF166odhlJLQkeSUEopVZM0QSmllKpJmqCUUkrVJE1QSimlapImKKWUUjVJE5RSSqmapAlKKaVUTdIEpZRSqiZVdaijhRKRUWCxw5m3AmNLEM7lQtfHuXSdnEvXybl0nZzrUtfJGmNM29lvPq8S1FIQkV3zDalRr3R9nEvXybl0nZxL18m5lnqd6CE+pZRSNUkTlFJKqZpUjwnqM9UOoMbo+jiXrpNz6To5l66Tcy3pOqm7c1BKKaWeH+qxB6WUUup5QBOUUkqpmlQ3CUpEXiEiB0XksIj8brXjqQYR6RWRH4vIPhF5RkQ+UHm/WUS+LyKHKv8uzfPMn0dExBGR3SJyb+XvdSLySKW8fFlEluaZ2c8TItIoIl8RkQMisl9Ebq7nciIiH6zsM3tF5J9EJFiPZUREPi8iIyKyd85785YLsT5RWT9Pi8h1C11eXSQoEXGATwGvBLYCbxGRrdWNqipKwH82xmwFbgLeV1kPvwv80BizCfhh5e968wFg/5y/Pwr8pTFmIzAJvLsqUVXPXwHfNcZcAVyDXTd1WU5EpBv4TWCnMWY74AC/RH2WkS8ArzjrvfOVi1cCmyqve4C/XujC6iJBATcCh40xR40xBeCfgddVOaYVZ4wZNMY8Ufl/ElvpdGPXxRcrk30RuLsqAVaJiPQArwY+W/lbgJcCX6lMUlfrREQagNuAzwEYYwrGmCnqu5x4gZCIeIEwMEgdlhFjzP3AxFlvn69cvA74e2M9DDSKSNdCllcvCaobODnn71OV9+qWiKwFdgCPAB3GmMHKR0NAR7XiqpKPA/8VcCt/twBTxphS5e96Ky/rgFHg7yqHPT8rIhHqtJwYY/qBvwBOYBPTNPA49V1G5jpfuVh0vVsvCUrNISJR4N+A3zLGJOZ+Zux9B3Vz74GIvAYYMcY8Xu1YaogXuA74a2PMDiDNWYfz6qmcVM6pvA6buFcBEc49zKVY+nJRLwmqH+id83dP5b26IyI+bHL6kjHmq5W3h2e63pV/R6oVXxXcArxWRI5jD/2+FHv+pbFyOAfqr7ycAk4ZYx6p/P0VbMKq13JyB3DMGDNqjCkCX8WWm3ouI3Odr1wsut6tlwT1GLCpctWNH3uC85tVjmnFVc6tfA7Yb4z52JyPvgm8vfL/twPfWOnYqsUY8yFjTI8xZi22XPzIGPNW4MfAGyuT1ds6GQJOisiWylsvA/ZRv+XkBHCTiIQr+9DM+qjbMnKW85WLbwK/Wrma7yZges6hwItSNyNJiMirsOcaHODzxpg/qW5EK09EbgUeAPZw+nzL72HPQ/0LsBr7OJM3G2POPhF62ROR24HfMca8RkTWY3tUzcBu4G3GmHwVw1tRInIt9qIRP3AUeCe2QVuX5URE/hj4ReyVsLuB92DPp9RVGRGRfwJuxz5WYxj4I+DrzFMuKsn8k9jDoRngncaYXQtaXr0kKKWUUs8v9XKITyml1POMJiillFI1SROUUkqpmqQJSimlVE3SBKWUUqomaYJSCnuJuYgYEXlHtWO5FCLyjkr8t1c7FqWWiiYopeYhImtF5MOV+4FqQiWJflhEGqsdi1IrQe+DUgoQEQ/2ptSiMaZc6Yn8GHtz4ReqGNosEfkw9sbIdcaY42d95gA+oGCMcc/9tlLPP97nnkSpy1+lUs+t1PJEJFZ55MmSMMaUgfJSzU+pWqCH+JTizHNQlfNQP6589HeV942I/GTO9CIi7xWRx0UkIyIpsU8rfslZ811b+e6HReQXK9Nngf9d+fwKEfk/lae1JivzelxE3nPWfL6A7T0BHJsT04crn897DkpEWkXkUyJyUkQKlX8/JSItZ0038/2XisjviMgREcmLyLMi8naUqgLtQSl1rvuBP8WOU/gZ7PiFYMcem/EPwFuwI33/HRAA3gp8X0TeYIw5ezDiu7FPZf1r4G+Amcec3I59OOC9wDHsoxzeBPytiLQZY/5nZbpPA3Hg9cAHgbHK+0+f70dUHjz4ELAR+DzwBPYZYO8FXioiN87Ti/tTIFRZXr4y7RdE5LAx5sHzLUupZWGM0Ze+6v6FTRQGeMd8f5817esrn91z1vteYBc20cyc311bmbYIXDnPvCLzvOcBfoJ9MJ5vzvsfrsxr7TzfeUfls9vnvPcnlff+01nTvq/y/kfm+f5uwD/n/W5sovqnam8jfdXfSw/xKbVwbwOSwNcrh9BaRaQVaAT+HZuUNp31nW8ZY/afPSNjTHrm/yISrBx6awbuw/aYrlhEnK/HPhn3M2e9/+nK+6+f5zv/xxhTmBNfP/As5/4epZadHuJTauGuBGKcecjvbB3Yin3Gs/NNVHm68YeBN3Pmw91mNF1aiIB9Auwuc/qx5AAYY0oi8iz2IYRnOzrPe+PAmkXEodQl0QSl1MIJtgfyyxeYZu9Zf2fOM90/Aq/B9nLuxyaDMvAq7LmmlT7Kcb4rAWVFo1AKTVBKnc+FbhA8BGwGHjbGpC51AZUbbl8D/IMx5j+e9dkdC4xpPkeBLSLinduLqjymfDPz95aUqhl6Dkqp+c0knuZ5Pvt77L7zP+f5DBHpuMhlzPRWzuidiEgX9omtC4lpPl8H2uaZ169V3v/aRc5HqarQHpRS89uHvRDiP4lIBpgCRowxPzLGfEVE/g74DRG5DnuJ+BjQA9yMvax7/XMtwBiTFJH7gLdV7o16DHuu59exVwK2nPWVhyv/flREvoS9sXivMebsw4kz/gx7yfqnKnHuxl5m/m7gYOVzpWqW9qCUmocxJgv8EvZ+pY8D/wT84ZzP3wX8KuACH8LeePt2bC/nQwtY1Nuw9yj9B+CT2Pulfh/41DwxPQj8N2AD8LeVmN54gd8wDdyCvWrvVcAnKv/+DXCrWcKRLJRaDjoWn1JKqZqkPSillFI1SROUUkqpmqQJSimlVE3SBKWUUqomaYJSSilVkzRBKaWUqkmaoJRSStUkTVBKKaVqkiYopZRSNen/Af95JK6QLddyAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.figure()\n", "plt.plot(sampling[:100,0])\n", "plt.plot(sampling[:100,1])\n", "plt.plot(sampling[:100,2])\n", "plt.plot(sampling[:100,3])\n", "plt.xlabel('iteration',size=18)\n", "plt.ylabel(r'$\\mu_k$',size=18)\n", "plt.tight_layout()\n", "plt.savefig('convergence_means.png')\n", "\n" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvG0lEQVR4nO3deXhTZd4+8PvbhlYoWJaWWlqgLBVBBNEKKogsOoMLII4LiCBuqMCo8PoKjg4yOO8wwqvj4KA/FfcZQWUcRMVBX0ERRqBBUAtVloJQQC2LgGj3+/dH0pCmKaRbTqD357p60Zw8OfkmnObO85znnGMkISIiEmminC5AREQkGAWUiIhEJAWUiIhEJAWUiIhEJAWUiIhEJJfTBdS2hIQEpqWlOV2GiIiEaO3atXtJJgYuP+kCKi0tDW632+kyREQkRGb2bbDlGuITEZGIpIASEZGIpIASEZGIpIASEZGIpIASEZGI5GhAmdkgM/vGzLaY2ZQg97c1s4/M7Esz+9jMUp2oU0REws+xgDKzaABzAFwGoAuAEWbWJaDZ/wJ4hWQ3ANMBzAhvlSIi4hQne1A9AWwhmUOyEMB8AEMD2nQBsNT7+7Ig94uIyEnKyYBKAbDT73aud5m/LwBc7f19GIAmZtYicEVmNtbM3GbmzsvLq5Ni5eT34osvIj09HWvXrnW6FBFB5E+SuA/AxWa2DsDFAHYBKAlsRPJZkhkkMxITK5wtQwLMmzcPO3bscLqMiDN37lxs2bIF/fr1wwcffOB0OSL1npMBtQtAa7/bqd5lPiR3k7yaZA8AD3qX/Ri2Ck9Cq1atwg033IA//OEPTpcSUfbt24dVq1Zh7NixaN++Pa644gr8/e9/D9q2pKQEO3fuDHqfiNQeJwMqE0C6mbUzsxgAwwEs8m9gZglmVlbjAwBeCHONJ53p06cDAN577z2UlpYes+3ChQuxcePGcJTluCVLlqC0tBS33norli9fjosuugijRo3CkiVLKrSdNWsW0tPT8cMPPzhQaXgUFBRgxYoVTpchtWzt2rV4//33sWHDBhw+fNjpco6PpGM/AC4HsAnAVgAPepdNBzDE+/s1ADZ728wFEHu8dZ577rmU4NasWUMAPPfccwmAq1evrrTt6tWrCYDNmjXj+vXrw1ilM2644QYmJiaypKSEJJmfn8/k5GReeeWV5dqVlJQwLS2NAPjqq686Uepxff311+zduzeXLVtW7XVMmDCBALhx48baKyxMsrKyeNttt/HZZ5/lli1bWFpa6nRJtW7z5s185plnWFhYGPJjduzYwVNOOYUAfD9dunRhdnZ2HVYaGgBuBsuIYAtP5B8FVOUGDx7MZs2acfv27YyOjuZDDz0UtF1JSQl79erFpKQktm7dmgkJCdywYUPY6lyzZg1Xrlx53HY7duzgZ599VqV1r1y5kmeeeSYzMzN9y4qLi9m8eXOOHj26XNvf/e53jIqK4s6dO33LPvroI98f94033lil5/Z34MAB7t69u9qPr8z+/fuZnp5OAExMTGRubm6V15GZmUkzIwA++eSTtV5jVRw+fJj5+flB7ysoKKiwbO/evUxLS/PVD4Bt27blBx98UKFtYWEhly5dyuLi4gr3rV69mvPnzw97uBUUFPCbb76p9P5NmzZx9OjRjIqKIgD+/ve/D3ndo0aNYmxsLN977z3OmzePM2bMYMuWLdm0adMafZmpDQqoeua7777j5s2bfbc///xzAuD06dNJkn379mX37t2DPvbll18mAL700kvctGkTTzvtNCYnJ5dbn7+yXkdt+Pnnn5mcnMzk5OSgHxxljhw54vsg7t+/P1esWHHcde/Zs4etWrXyPabMihUrCICvv/56ufZbtmwp956R5MiRIxkfH89hw4axZcuW1XrtX331FVNSUnjaaafx8OHDVX58ZYqKinjppZeyQYMGfP7559m4cWNeeOGFVfqWXVxczHPOOYfJyclMSUnh1VdfXWv1VUV+fj5nzZrF+Pj4CgGTn5/PadOmMSYmhjfddBOPHDlC0vP6Bw4cyJiYGK5atYrZ2dmcM2cOO3XqxKSkJO7du7fcc5T1Evv27cvt27eT9GzLM2bMYHR0NAHw9ttvr9L7F4rt27fz0ksv5Z///Gdf7aTni0HXrl2DboulpaWcMmUKo6Ki2LBhQ06aNInDhw9nVFQUP/300+M+59q1awmAkydPLrc8JyeHnTt3ZoMGDfjSSy/VzgusBgXUSeiDDz5gnz59gn7InX/++QTAoUOHcs2aNbzqqqsYHx/PH3/8kSQ5a9YsAuC3335b7nGHDh3iaaedxp49e/o+fDds2MCEhASecsop7NWrF8eNG8fZs2dzwoQJ7NWrF2NjYzlhwoQKNRQVFXH27NncsWNHyK/pr3/9q++b70cffVRpu4kTJxIAJ02axKSkJALgoEGDKu2VFBUV8eKLL2bDhg15xx13EAD/7//+jyT5wAMPMDo6mgcOHKjwuIEDB7Jt27YsKSnhgQMHeMopp/Cuu+7yhfjatWtDfm0k+fHHHzM+Pp4JCQkEwGnTplXp8cdyzz33EADnzp1Lknz99dcJgPfee2/I65g9e7bvA3LMmDFs3rx5rX4BOZ7S0lK+9dZb7NChg+//tFOnTgTA2267jUuWLGGXLl0IgH369KGZsVu3bty8eTMnTZpEAHzhhRfKrXP9+vV0uVwcOXKkb9nbb79NAPz1r3/NJk2a8NRTT+UzzzzDQYMGEQCvu+46Tp482fdlZt++fSHVX1RUdMz7d+7cyfbt2zMmJoYA2KpVKz7zzDP83e9+x+joaLZq1Yo9evRgw4YNy21b06ZNIwCOGTOG3333HUnP32r79u3Ztm1b39816Qlw/+ArLS1lv379mJCQUK5dmQMHDnDAgAEEwPfffz+k11nbFFAnoQsvvLDcB1KZrKwsAuDAgQPZrFkz3wf+ww8/7GuTnZ1NAJwzZ065x5b9Ua5atarc8uzsbE6cOJF9+/ZlkyZNCIBxcXHs27cvzzvvPMbGxjIvL6/cY1555RXfEMuWLVuO+3p++eUXtmrVihdccAGbNGnCW265JWi75cuX08w4fvx4kp7e1KxZsxgXF8e2bdsG3W9y3333+fYb/fLLL0xNTWWvXr1YWlrK7t278+KLLw76XPPnzycALlmyhE899RQB0O12c8+ePQTAP/3pT8d9XWUWLFjA2NhYdu7cmdu3b+c111zDuLg47tmzJ+R1BFNQUMCHHnqIAHjPPfeUu+/uu+/29YaPZ9euXWzSpAl//etfs7S01BfCVdkHWZMhsYMHD3LEiBEEwDPPPJNLliwh6dkuJk+e7BvWatOmDd977z2S5Pvvv8/mzZszLi6OAPjb3/426LoffvhhAuDbb7/N3NxcNm/enD169GB+fj63bdvGPn36EABjY2P59NNP+17HK6+8wpiYGHbo0IF//etfmZOT41vn3r17+e6773LatGkcNmwY27dvTwBMSUnhZZddxsmTJ/Pdd9/1DVHu3r2bp59+Ops0acLVq1fzk08+4QUXXOD7+7z55pt54MABfvfdd2zTpg1TU1O5e/du35eGMWPGVHh/P/vsM0ZHR3PkyJHcunUr//u//5stWrRgo0aN+NBDD/HgwYNctGgRAfBvf/tbpe99QUEBW7duXW5kIZwUUCeQdevW8d577+W8efMq7RFkZmYSAKOiotirV69y99133310uVz84YcfeOjQIc6aNYtDhgzh/v37fW1KS0vZsWNHDho0yLds/fr1vmGTYykpKWFubq5vCK4sEGfMmFFu/WeddRbbt2/P5s2bs1WrVvz666+Pud4nn3ySALh06VKOHj2a8fHx/OWXX8q1+emnn9ihQwe2a9euQs/R7XYzKSmJzZo146effsqSkhJmZmby/vvvJwBfoJHkc889RwC+0Hn00UeD1pSfn88WLVrwmmuuYUZGBrt16+b7kDj77LMrBNs777zDefPmVVjPv//9b0ZFRfHCCy/0fRvftGkTXS4X77zzzkrfk+eff56XXHIJZ8+eHXR/0urVq33DQjfeeGOFb/AFBQXs27cvAXDUqFHltoFA1113HWNjY31fJnbs2EEAfPzxxyt9jL/c3Fx27NiR99xzT6VBVVJSwg8//JCjR4/mnXfeyYULF/LQoUNcs2YN27dvz6ioKE6fPj1oTyQzM5MzZ86s8P++fft29u7dm5dddlmlw3EFBQXs1q0bTzvtNF500UWMi4srt6+nuLiYL7zwAr/44osKj125cqWv11YWnqeffrrvtpmxU6dOvPbaa/nggw9y1KhR7N69Oxs0aEAAbNKkCa+//np26dKFcXFx5YajS0tL+cEHH1QYplu3bh0bNWrk60kOHTq00t7ZH/7wB18t0dHR/M1vfsPrr7+eAJiQkMDU1FR26tTpuEOVjz76KAFw3bp1x2xXFxRQESiwx0F6dtyWfeCU/XTq1KnCjLvRo0ezcePGvq5/2R9WYWEhk5KSeNVVVx33+SdOnMiYmBgePnyYWVlZTExMZKtWrar1jX7gwIFs3bq1749oyZIlvm/uX3zxBVu2bMmkpCTOnj2bU6ZM4fDhw3ndddf56s7Pz2dKSgovuugilpaW+h7/z3/+s9zzlPUIPv7446B15OTk8PTTT2dsbKxv6M/MeNVVV5XbqV5YWMiOHTvS5XIRALOysip9bZMmTfJ9e3/iiSd8yydPnkyXy8VDhw75nrthw4YEwL/85S++dtnZ2YyPj2f37t0rfLhOmDCB0dHRQWdS5ebmMi4uztdjBcDzzjuPV1xxBYcNG8Yrr7ySUVFRTElJ4dtvv11p/QUFBZw6dSqjo6OZnJzMd999t0KbxYsXEwAfeeSRcss7duzIwYMHl1sWLHyOHDnCc8891zc5Yfbs2eXu/+mnn/jII4/4ZkA2bdqUjRs3JgC6XC66XC62adMmpH2J1fX555/79i29+OKLVX785s2b+fjjj/PSSy/lkCFDOGPGDH788cf86aefgrbPz8/n+++/z9tvv50tW7Zkw4YNqzQZ4Z///CcBsF+/fhW+qPkrKiriuHHjOHXq1HJfYjIzM9m/f38C4DvvvHPc59u/fz8bNWp03C+odUEBFUFKS0v5+9//Pui30yeeeIIA+OabbzIzM5OPPfYYW7duzbS0NB48eJCkZwJETEwMJ0yYwL179/p+J+nrzh/rA6vM0qVLfcNUSUlJTE5OPuYMomNZuHAhAXDBggUkyUsuuYStWrXyhUJ2drZvgoLL5WKHDh3YrFkzRkdH8+677+aMGTPK7RcqKipiUlJSuZ30Za+tsmGcMnl5eRwxYgSvv/56vvrqq0G/CJDkP/7xD98Q5LGGpjZu3EgAjImJKbejfdmyZQTAhQsXsrS0lJdffjnj4uJ4xRVX+MJs37597NixI1u2bOnbEe/v+++/Z5MmTYJ+oRg5ciRjY2O5detWZmdn85FHHuFFF13Ec845h127dmWnTp04fvz4oPsVglm7di3POussmlm5nfBHjhxhWloazzjjjAoz5m6//XbGx8eXm7AyZswYdujQgYsWLSLp6RVde+21NDO+/fbbvOqqqxgVFeXbn+F2u309joEDB3LevHn85ZdfWFBQwGXLlnHKlCm89957j9m7qy1z587lww8/HPbZecXFxSH/P/n76quv+PPPP1f7eUtLS/n999+H3H78+PGMiYkp9yV18eLFHDlyZNDJGAUFBdy1a1e16yujgIoQhYWFvPnmm31j1dHR0Vy6dClJT/Cceuqp/NWvflXuD+g///kPo6OjOWrUKJLk9OnTCcA3ZHbDDTewadOmPHLkiG92WSgzjwoLCxkfH08ATEpKqtHxEMXFxUxLS2Pfvn25bt06AuCf//zncm2OHDnCb7/91vdht2/fPt51112+b929e/cu97rvvvtuxsbG8sCBA8zJyWHTpk3Zo0ePY36brIqSkhIOGDCgQq8hmKFDh/KOO+4ot6ygoICNGzfmXXfdxQULFhAAH3vsMRYWFvLqq68mAHbs2JExMTHHnDb/xz/+sUKvq2xm4YMPPljt1xfMzz//zN69ezMmJsb3bX7KlCmV9kpfe+01AvBNy//kk098PSAAHDx4sG9yxsyZM0l6poZ3796dp556KqdMmcIGDRowJSXFt51L5Nq0aRPNzDd9fcGCBXS5XOX+Rt944w3Onj2bV155JePi4srtJqguBVQEOHz4sG+W0MMPP8yDBw+yc+fOTEhI4Lfffsubb76ZDRo0CLqvpmwn78svv8zk5ORyG0XZN/nHHnuMLpeLkyZNCrmmsuGH2jjOaebMmQTAXr16sXHjxkFnxQXjdrt53XXX0e12l1tedrDwU089xXPOOYdNmzbl1q1ba1xnbRo8eDDbtGnDlJQUdu/e3TfEWVhYyGHDhvn+z44lPz/fF2j33Xcfi4qKePbZZzM1NbXS4aOa2LdvHzt37sz4+Hi+9tprdLlcHDNmTNC2u3fv9oVPcXExe/TowdatW/PHH3/kzJkz2ahRIwLgTTfdVO7LxY4dO3xDrFdffXXIs+DEeYMHD2ZCQgJffPFFRkdH88ILL+SePXv45JNPsk2bNr7h5o4dO3LcuHG+CSs1oYCKAPfffz+joqL43HPP+ZZ9/fXXPPXUU33H9Nx///1BH1tUVMTzzz/fty9k8eLFvvtKS0uZnp7u2yn75ZdfhlxTQUFBpQdCVtW+fft8+2AmTpxY4/WVlpayQ4cOvtcVyrBluM2ZM8e3nytw5mNRUVHIQ6bFxcUcN24cAfCss84iAM6fP78uSiZJfvvtt74h1xYtWlQ6DEqSnTt35qBBg/j8888TAF977TXffTt27OBTTz0VdBvKzs7mv/71r5PyTA4ns7Khf8BzjFjZPlbS88Xrww8/rPUvigqoCDBgwIAKM+7Io8dkJCcnl9sYAm3dupWNGzdmenp6hWNTynovTr/+sWPH0uVyBd3fUh1Tp04lUPEAw0ixdetWmtkxZ+KFqrS01Dfcd/HFF9f5B/sXX3zBjh07VjgoNNC4ceMYFxfHpKQkXnDBBQqck1xpaSn79+/PK664ok568MEooCJAy5YtKz2256233qowxBVMVlZW0DM6fP/992zWrJmjR4OTnmHMqvTgQlnfyy+/fNwDIJ2UmZkZ9LQ71bVixYpj9mjC7c033/R9o16zZo3T5UgYhPtLSGUBZZ77Th4ZGRl0u91Ol1FBXl4eWrZsiccffxwTJ06sk+coLi6Gy+Wqk3VL/bV3714kJyfjhhtuwMsvv+x0OXISMrO1JDMCl+vTLEw2bNgAADjzzDPr7DkUTlIXEhISsGbNGpxxxhlOlyL1jD7RwiQrKwsA0LVrV4crEam6Hj16OF2C1EORfsn3k8aGDRvQrFkzJCcnO12KiMgJQQEVJllZWTjzzDNhZk6XIiJyQlBAhQFJbNiwQcN7IiJVoIAKgz179uDAgQN1OkFCRORko4AKg7IZfOpBiYiEztGAMrNBZvaNmW0xsylB7m9jZsvMbJ2ZfWlmlztRZ02VzeBTD0pEJHSOBZSZRQOYA+AyAF0AjDCzLgHNHgLwBskeAIYDeCq8VdaOrKwstGzZEomJiU6XIiJywnCyB9UTwBaSOSQLAcwHMDSgDQGc6v09HsDuMNZXazRBQkSk6pwMqBQAO/1u53qX+ZsG4EYzywWwGMBvg63IzMaamdvM3Hl5eXVRa7WVzeDT8J6ISNVE+iSJEQBeIpkK4HIAr5pZhZpJPksyg2RGpA2j7dixAz/99JN6UCIiVeRkQO0C0Nrvdqp3mb9bAbwBACQ/A3AKgISwVFdLNEFCRKR6nAyoTADpZtbOzGLgmQSxKKDNDgADAcDMOsMTUJE1hncc4ThJrIjIycixgCJZDGACgCUAsuGZrbfBzKab2RBvs/8CcLuZfQFgHoAxPMGuD5KVlYXU1FQ0bdrU6VJERE4ojp7NnORieCY/+C+b6vf7RgC9w11XbdIECRGR6on0SRIntIKCAmzYsAHdunVzuhQRkROOAqoOud1uFBQUoHfvE7oTKCLiCAVUHVqxYgUA4MILL3S4EhGRE48Cqg6tWLECZ5xxhk5xJCJSDQqoOlJaWoqVK1eiT58+TpciInJCUkDVkezsbBw4cED7n0REqkkBVUfK9j+pByUiUj0KqDqycuVKJCUloUOHDk6XIiJyQlJA1ZEVK1agT58+MDOnSxEROSEpoOrArl27sG3bNg3viYjUgAKqDqxcuRKA9j+JiNSEAqoOrFixAnFxcTj77LOdLkVE5ISlgKoDK1aswPnnnw+Xy9Fz8YqInNAUULXs0KFD+OKLL3T8k4hIDSmgatnq1atRWlqqgBIRqSEFVC3LzMwEAPTs2dPhSkRETmwKqFrmdruRnp6uK+iKiNSQAqqWZWZmIiMjw+kyREROeI4GlJkNMrNvzGyLmU0Jcv9fzGy992eTmf3oQJkh++6775Cbm4vzzjvP6VJERE54js2DNrNoAHMAXAogF0CmmS0iubGsDcmJfu1/C6BH2AutgrVr1wKAelAiIrXAyR5UTwBbSOaQLAQwH8DQY7QfAWBeWCqrpszMTERFRaFHj4jOURGRE4KTAZUCYKff7VzvsgrMrC2AdgCWhqGuanO73ejcuTMaN27sdCkiIie8E2WSxHAAC0iWBLvTzMaamdvM3Hl5eWEuzYMk3G63hvdERGqJkwG1C0Brv9up3mXBDMcxhvdIPksyg2RGYmJiLZYYutzcXHz//fcKKBGRWuJkQGUCSDezdmYWA08ILQpsZGZnAGgG4LMw11clbrcbADSDT0SkljgWUCSLAUwAsARANoA3SG4ws+lmNsSv6XAA80nSiTpD5Xa74XK50K1bN6dLERE5KTh6um2SiwEsDlg2NeD2tHDWVF1utxtdu3ZFw4YNnS5FROSkcKJMkohomiAhIlL7FFC1YNu2bdi/f7/2P4mI1CIFVC0omyChHpSISO3RJV9rYNeuXfj000/x9NNPIyYmBl27dnW6JBGRk4YCqhqKiopwySWXYPny5QCAJk2aYPz48YiJiXG4MhGRk4cCqhr+9a9/Yfny5ZgyZQquvfZadOvWDS6X3koRkdpkEX54UZVlZGSwbJ9QXbnooouwe/dubNq0CdHR0XX6XCIiJzszW0uywk58TZKoovXr12PFihUYP368wklEpA4poKroySefRKNGjXDLLbc4XYqIyElNAVUFe/fuxWuvvYZRo0ahadOmTpcjInJSU0BVwdy5c5Gfn48JEyY4XYqIyElPARWi4uJiPPXUU+jfv7+OdxIRCQMFVIg+/PBD7Ny5U70nEZEwUUCFKCsrCwAwcOBAhysREakfFFAhysnJQYsWLRAfH+90KSIi9YICKkQ5OTlo376902WIiNQbCqgQKaBERMJLARWCkpISbN++XQElIhJGCqgQ7Ny5E8XFxQooEZEwcjSgzGyQmX1jZlvMbEolba4zs41mtsHMXgt3jYBneA+AAkpEJIwcu0aEmUUDmAPgUgC5ADLNbBHJjX5t0gE8AKA3yQNm1tKJWssCqkOHDk48vYhIveRkD6ongC0kc0gWApgPYGhAm9sBzCF5AABI/hDmGgF4AsrlciE1NdWJpxcRqZecDKgUADv9bud6l/k7HcDpZrbSzFaZ2aBgKzKzsWbmNjN3Xl5erReak5ODtLQ0XV5DRCSMIn2ShAtAOoB+AEYAeM7MmgY2IvksyQySGYmJibVehKaYi4iEn5MBtQtAa7/bqd5l/nIBLCJZRHIbgE3wBFZYKaBERMLPyYDKBJBuZu3MLAbAcACLAtoshKf3BDNLgGfILyeMNeLgwYPYt2+fAkpEJMwcCyiSxQAmAFgCIBvAGyQ3mNl0MxvibbYEwD4z2whgGYD/JrkvnHVu27YNgKaYi4iEm2PTzAGA5GIAiwOWTfX7nQAmeX8coWOgREScEemTJByngBIRcYYC6ji2bt2qy2yIiDhAAXUcmsEnIuKMkPZBmVlV9wGR5F+qUU/EycnJwbnnnut0GSIi9U6okyT+t4rrJYATPqDKLrNx7bXXOl2KiEi9E2pA9a/TKiJUbm6uLrMhIuKQkAKK5Cd1XUgk0gw+ERHnaJLEMSigREScE+okib5VXTHJ5VUvJ7KUXWajdevWx28sIiK1KtR9UB/DM/EhFOZte8JfmyInJwdt27bVZTZERBxQlVMd5QN4C8D6uikl8uzfvx91cfkOERE5vlAD6kkAN3h/zgTwAoB/lF3p9mSVn5+PU045xekyRETqpZAmSZC8B0AreC6JsQeeY5x2m9k8M/tVHdbnKAWUiIhzQp7F571o4JskrwDQFsAjAM4B8G8z2+G9TEZqXRXqBAWUiIhzqjXNnORukn8i2QlAXwDfAHgQwC21WZzTFFAiIs6p9vWgzCwWwNUAbgYwAJ5JFGG92m1dU0CJiDinyj0oMzvPzJ6GZ1/UPwDEAxgHIJnk32u5PkcpoEREnBPqgbotAYyCp7fUBcAP8Mzke4Hkxrorz1kKKBER54Q6xJcLz8G378Ozr+ldkiV1VlWEUECJiDgn1CG+siAbAOBVAAfM7NAxfg6GslIzG2Rm35jZFjObEuT+MWaWZ2brvT+3hVhvjZFEYWGhAkpExCGh9qCWI/RTHYXEzKIBzAFwKTw9tEwzWxRkyPB1khNq87lDUVBQAAAKKBERh4R6uY1+dfDcPQFsIZkDAGY2H8BQABGxTys/Px8AEBsb63AlIiL1k5OX20gBsNPvdq53WaDfmNmXZrbAzIKeVtzMxpqZ28zceXl5tVJcWUCpByUi4oxQZ/EtquJ6SXJoNeoJ9A6AeSQLzOwOAC/Dsx8s8MmeBfAsAGRkZNTKUKQCSkTEWaHug7qyiusNJSR2AfDvEaV6lx1dCbnP7+ZcADOrWEe1KaBERJwV6slio473A6A/gEzvQ/aEsNpMAOlm1s7MYuA5EW25npqZJfvdHAIgO5R6a4MCSkTEWTXeB2VmXc3sPQBLAXQC8HsA6cd7HMliABMALIEneN4gucF70tkh3mZ3m9kGM/sCwN0AxtS03lApoEREnFWTc/G1hueM5iMBlACYDeCPAcNyx0RyMYDFAcum+v3+AIAHqltjTSigREScVeWAMrNm8JxNYhyAWADzADxEcnvtluYsBZSIiLNCDijv2cvvBTAZQFMAHwKYTHJ9XRTmNAWUiIizQtoHZWa3AtgC4E8AtgK4lOSvT9ZwAhRQIiJOC7UH9Rw8U8fdAN4A0N3Muh+jPUn+pabFOUmnOhIRcVZV9kEZgPO8P8dDACd0QOlURyIizgo1oPrXaRURSEN8IiLOCvVksZ/UdSGRRgElIuIsJ08WG9E0xCci4iwFVCXy8/MRExODqCi9RSIiTtCnbyV0uXcREWcpoCqhgBIRcZYCqhIKKBERZymgKqGAEhFxlgKqEgooERFnKaAqoYASEXGWAqoSBQUFOgZKRMRBCqhKqAclIuIsBVQlFFAiIs5SQFVCASUi4ixHA8rMBpnZN2a2xcymHKPdb8yMZpYRrtoUUCIiznIsoMwsGsAcAJcB6AJghJl1CdKuCYB7AKwOZ30KKBERZznZg+oJYAvJHJKFAOYDGBqk3SMAHgWQH87iFFAiIs5yMqBSAOz0u53rXeZjZucAaE3yvWOtyMzGmpnbzNx5eXm1UpwCSkTEWRE7ScLMogA8DuC/jteW5LMkM0hmJCYm1vi5SSqgREQc5mRA7QLQ2u92qndZmSYAugL42My2AzgfwKJwTJQoKioCSQWUiIiDnAyoTADpZtbOzGIADAewqOxOkgdJJpBMI5kGYBWAISTddV2YLvcuIuI8xwKKZDGACQCWAMgG8AbJDWY23cyGOFUXoMu9i4hEApeTT05yMYDFAcumVtK2XzhqAjzn4QPUgxIRcVLETpJwkob4REScp4AKQgElIuI8BVQQCigREecpoIJQQImIOE8BFYQCSkTEeQqoIBRQIiLOU0AFoYASEXGeAioIBZSIiPMUUEEooEREnKeACkKnOhIRcZ4CKgid6khExHkKqCA0xCci4jwFVBD5+fmIjo6Gy+XouXRFROo1BVQQupquiIjzFFBBKKBERJyngApCASUi4jwFVBAKKBER5ymgglBAiYg4TwEVhAJKRMR5jgaUmQ0ys2/MbIuZTQly/51m9pWZrTezFWbWJRx1KaBERJznWECZWTSAOQAuA9AFwIggAfQaybNIng1gJoDHw1Fbfn6+TnMkIuIwJ3tQPQFsIZlDshDAfABD/RuQPOR3Mw4Aw1GYelAiIs5z8lQJKQB2+t3OBdArsJGZjQcwCUAMgAHBVmRmYwGMBYA2bdrUuLCCggIFlIiIwyJ+kgTJOSQ7AJgM4KFK2jxLMoNkRmJiYo2fUz0oERHnORlQuwC09rud6l1WmfkArqrLgsoooEREnOdkQGUCSDezdmYWA2A4gEX+Dcws3e/mFQA2h6MwBZSIiPMc2wdFstjMJgBYAiAawAskN5jZdABukosATDCzSwAUATgA4KZw1KaAEhFxnqPXkyC5GMDigGVT/X6/J+xFQQElIhIJIn6SRLgVFxejuLhYASUi4jAFVABd7l1EJDIooALocu8iIpFBARWgLKB0qiMREWcpoAKoByUiEhkUUAEUUCIikUEBFUCTJEREIoMCKoB6UCIikUEBFUABJSISGRRQARRQIiKRQQEVQAElIhIZFFABFFAiIpFBARVAASUiEhkUUAEUUCIikUEBFUCnOhIRiQwKqADqQYmIRAYFVICygIqJiXG4EhGR+k0BFaDsarpm5nQpIiL1mgIqQEFBgYb3REQigKMBZWaDzOwbM9tiZlOC3D/JzDaa2Zdm9pGZta3rmsp6UCIi4izHAsrMogHMAXAZgC4ARphZl4Bm6wBkkOwGYAGAmXVdlwJKRCQyONmD6glgC8kckoUA5gMY6t+A5DKSP3tvrgKQWtdFKaBERCKDkwGVAmCn3+1c77LK3Arg/WB3mNlYM3ObmTsvL69GRSmgREQiwwkxScLMbgSQAWBWsPtJPksyg2RGYmJijZ5LASUiEhlcDj73LgCt/W6nepeVY2aXAHgQwMUkC+q6KAWUiEhkcLIHlQkg3czamVkMgOEAFvk3MLMeAJ4BMITkD+EoKj8/X6c5EhGJAI4FFMliABMALAGQDeANkhvMbLqZDfE2mwWgMYA3zWy9mS2qZHW1Rj0oEZHI4OQQH0guBrA4YNlUv98vCXdNCigRkchwQkySCCcFlIhIZFBABdCpjkREIoMCKoB6UCIikUEBFUABJSISGRRQfkhqiE9EJEIooPwUFHiOA1ZAiYg4TwHlR5d7FxGJHI4eBxVpGjdujDVr1iA1tc5Pmi4iIsehgPLjcrlw3nnnOV2GiIhAQ3wiIhKhFFAiIhKRFFAiIhKRFFAiIhKRFFAiIhKRFFAiIhKRFFAiIhKRFFAiIhKRjKTTNdQqM8sD8G0NV5MAYG8tlHMy0HtxlN6Lo/ReHKX34qjqvhdtSSYGLjzpAqo2mJmbZIbTdUQCvRdH6b04Su/FUXovjqrt90JDfCIiEpEUUCIiEpEUUME963QBEUTvxVF6L47Se3GU3oujavW90D4oERGJSOpBiYhIRFJAiYhIRFJA+TGzQWb2jZltMbMpTtcTTmbW2syWmdlGM9tgZvd4lzc3sw/NbLP332ZO1xouZhZtZuvM7F3v7XZmttq7fbxuZjFO1xguZtbUzBaY2ddmlm1mF9TXbcPMJnr/RrLMbJ6ZnVJftg0ze8HMfjCzLL9lQbcD85jtfU++NLNzqvp8CigvM4sGMAfAZQC6ABhhZl2crSqsigH8F8kuAM4HMN77+qcA+IhkOoCPvLfri3sAZPvdfhTAX0h2BHAAwK2OVOWMvwL4N8kzAHSH532pd9uGmaUAuBtABsmuAKIBDEf92TZeAjAoYFll28FlANK9P2MBPF3VJ1NAHdUTwBaSOSQLAcwHMNThmsKG5B6Sn3t/PwzPB1AKPO/By95mLwO4ypECw8zMUgFcAWCu97YBGABggbdJfXov4gH0BfA8AJAsJPkj6um2AcAFoKGZuQA0ArAH9WTbILkcwP6AxZVtB0MBvEKPVQCamllyVZ5PAXVUCoCdfrdzvcvqHTNLA9ADwGoASST3eO/6DkCSU3WF2RMA7gdQ6r3dAsCPJIu9t+vT9tEOQB6AF71DnnPNLA71cNsguQvA/wLYAU8wHQSwFvV32wAq3w5q/JmqgJJyzKwxgH8CuJfkIf/76Dkm4aQ/LsHMrgTwA8m1TtcSIVwAzgHwNMkeAI4gYDivHm0bzeDpGbQD0ApAHCoOedVbtb0dKKCO2gWgtd/tVO+yesPMGsATTv8g+ZZ38fdl3XLvvz84VV8Y9QYwxMy2wzPUOwCefTBNvcM6QP3aPnIB5JJc7b29AJ7Aqo/bxiUAtpHMI1kE4C14tpf6um0AlW8HNf5MVUAdlQkg3TsbJwaeHZ+LHK4pbLz7WJ4HkE3ycb+7FgG4yfv7TQDeDndt4UbyAZKpJNPg2Q6WkhwJYBmAa7zN6sV7AQAkvwOw08w6eRcNBLAR9XDbgGdo73wza+T9myl7L+rltuFV2XawCMBo72y+8wEc9BsKDInOJOHHzC6HZ99DNIAXSP6PsxWFj5n1AfApgK9wdL/L7+DZD/UGgDbwXMbkOpKBO0lPWmbWD8B9JK80s/bw9KiaA1gH4EaSBQ6WFzZmdjY8E0ZiAOQAuBmeL7j1btswsz8AuB6ema/rANwGz76Vk37bMLN5APrBc1mN7wE8DGAhgmwH3gD/GzxDoD8DuJmku0rPp4ASEZFIpCE+ERGJSAooERGJSAooERGJSAooERGJSAooERGJSAookSDMrJ+Z0czGOF1LdZjZGG/9/ZyuRaS6FFAiITCzNDOb5j0eKCJ4Q3SamTV1uhaRuqDjoESCMLMoeA5KLSJZ4u2JLIPnYMOXHCzNx8ymwXOgZDuS2wPuiwbQAEAhydKKjxaJfK7jNxGpf7wf6vnhej4za+K9zEmtIFkCoKS21ifiBA3xiQThvw/Kux9qmfeuF73LaWYf+7U3M7vLzNaa2c9m9pN5rlDcP2C9ad7HTjOz673tfwHwpPf+M8zsKe8VWw9717XWzG4LWM9L8PSeAGCbX03TvPcH3QdlZglmNsfMdppZofffOWbWIqBd2eMHmNl9ZrbVzArMbJOZ3QSRMFAPSuT4lgP4EzznJnwWnnMWAp5zkZV5FcAIeM70/SKAWAAjAXxoZleTDDzx8FXwXJn1aQD/D0DZpU36wXNxwHcBbIPncg7XAnjOzBJJzvC2ewbAqQCGAZgIYK93+ZeVvQjvhQf/A6AjgBcAfA7Pdb/uAjDAzHoG6cX9CUBD7/MVeNu+ZGZbSK6s7LlEagVJ/ehHPwE/8AQFAYwJdjug7TDvfWMDlrsAuOEJmrL9vWnetkUAOgdZV1yQZVEAPobn4ngN/JZP864rLchjxnjv6+e37H+8y8YFtB3vXf5IkMevAxDjtzwFnqCa5/T/kX5O/h8N8YnU3I0ADgNY6B1CSzCzBABNAbwDTyilBzzmPZLZgSsieaTsdzM7xTv01hzAB/D0mM6oQZ3D4Lky7rMBy5/xLh8W5DFPkSz0q28XgE2o+HpEap2G+ERqrjOAJig/5BcoCZ4P9jKbgjXyXtF4GoDrUP5ib2WaVa9EAJ6rwLp59NLkAACSxWa2CZ6LEAbKCbJsH4C2NahDJCQKKJGaM3h6IDcco01WwO2fK2n3GoAr4enlLIcnDEoAXA7PvqZwj3pUNhPQwlqF1EsKKJHQHOuAwc0ATgewiuRP1X0C7wG3VwJ4leSdAfddUsWagskB0MnMXP69KO+lyk9H8N6SiGO0D0okNGXB0zzIfa/A87c0I8h9MLOkEJ+jrLdSrndiZsnwXLW1KjUFsxBAYpB13e5d/q8Q1yMSFupBiYRmIzwTIcaZ2c8AfgTwA8mlJBeY2YsAJpjZOfBMEd8LIBXABfBM625/vCcgedjMPgBwo/fYqEx49vXcAc9MwBYBD1nl/fdRM/sHPAcWZ5EMHE4sMxOeKetzvHWug2ea+a0AvvHeLxIx1IMSCQHJXwAMh+d4pScAzAMw1e/+WwCMBlAK4AF4Dry9CZ5ezgNVeKob4TlGaTCAv8FzvNSDAOYEqWklgMkAOgB4zlvTNcd4DQcB9IZn1t7lAGZ7//1/APqwFs9kIVIbdC4+ERGJSOpBiYhIRFJAiYhIRFJAiYhIRFJAiYhIRFJAiYhIRFJAiYhIRFJAiYhIRFJAiYhIRFJAiYhIRPr/7KPkZOC2pAkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "NMI=np.zeros(Niter)\n", "\n", "\n", "for n in range(Niter):\n", " NMI[n] = normalized_mutual_info_score(Zt,sampling[n,8:])\n", "\n", "\n", "plt.figure()\n", "plt.plot(NMI[:100],c='k')\n", "plt.xlabel('iteration',size=18)\n", "plt.ylabel('NMI',size=18)\n", "plt.tight_layout()\n", "plt.savefig('convergence_nmi.png')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.7.9" } }, "nbformat": 4, "nbformat_minor": 4 }