mixture/GaussianMixture.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.sys.path.append(os.path.dirname(os.path.abspath('.')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 数据准备"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.datasets.samples_generator import make_blobs\n",
    "from model_selection.train_test_split import train_test_split\n",
    "\n",
    "X, _ = make_blobs(cluster_std=1.5, random_state=42, n_samples=1000, centers=3)\n",
    "X = np.dot(X, np.random.RandomState(0).randn(2, 2))    # 生成斜形类簇\n",
    "\n",
    "plt.clf()\n",
    "plt.scatter(X[:, 0], X[:, 1], alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test = train_test_split(X, test_size=0.2)\n",
    "\n",
    "n_samples, n_feature = X_train.shape\n",
    "n_cluster = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 模型基础\n",
    "\n",
    "### 初始化\n",
    "首先需要随机初始化$K$个高斯分布。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # 随机初始化均值，维度为(n_cluster, n_feature)\n",
    "# # 生成范围/2是为了限制初始均值的生成范围\n",
    "# mu = np.random.randint(X_train.min()/2, X_train.max()/2, size=(n_cluster, n_feature))\n",
    "\n",
    "# # 一个协方差矩阵的维度为(n_feature,n_feature)\n",
    "# # 多个分布的协方差矩阵维度为(n_cluster,n_feature,n_feature)\n",
    "# cov = np.zeros((n_cluster, n_feature, n_feature))\n",
    "# for dim in range(len(cov)):\n",
    "#     np.fill_diagonal(cov[dim], 1)\n",
    "    \n",
    "# # 初始均匀的类分布概率\n",
    "# pi = np.ones(n_cluster)/n_cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### E-step\n",
    "E-step的实质就是根据已有的概率分布计算样本属于各类的概率。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # 概率矩阵\n",
    "# P_mat = np.zeros((n_samples, n_cluster))\n",
    "\n",
    "# for k in range(n_cluster):\n",
    "#     g = multivariate_normal(mean=mu[k], cov=cov[k])    # 根据均值与方差生成多维分布\n",
    "\n",
    "#     # 计算X在各分布下出现的频率\n",
    "#     P_mat[:, k] = p*g.pdf(X_train)\n",
    "\n",
    "# # 归一化使频率变成概率\n",
    "# P_mat/=P_mat.sum(axis=1).reshape(-1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### M-step\n",
    "M-step的实质就是根据样本的归类情况，更新高斯分布的参数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # M-step，更新参数\n",
    "# for k in range(n_cluster):\n",
    "#     N_k = np.sum(P_mat[:, k], axis=0)    # 类出现的频率\n",
    "#     mu[k] = (1/N_k)*np.sum(X_train*P_mat[:, k].reshape(-1, 1), axis=0)    # 该类的新均值\n",
    "#     cov[k] = ((1/N_k)*np.dot((P_mat[:, k].reshape(-1, 1)\n",
    "#                               * (X_train-mu[k])).T, (X_train-mu[k])))\n",
    "#     pi[k] = N_k/n_samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 预测\n",
    "各高斯分布的参数训练好之后，预测的实质就是对测试样本再执行一次E-step。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # 测试集的概率矩阵\n",
    "# pred_mat=np.zeros((X_test.shape[0], n_cluster))\n",
    "\n",
    "# # 计算测试样本出现于各类的频率\n",
    "# for k in range(n_cluster):\n",
    "#     g=multivariate_normal(mean=mu[k],cov=cov[k])\n",
    "#     pred_mat[:,k]=pi[k]*g.pdf(X_test)\n",
    "\n",
    "# # 归一化得到概率\n",
    "# totol_N = pred_mat.sum(axis=1)\n",
    "# totol_N[totol_N == 0] = n_cluster\n",
    "# pred_mat /= totol_N.reshape(-1, 1)\n",
    "\n",
    "# # 概率最大者为预测值\n",
    "# Y_pred=np.argmax(pred_mat,axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 完整模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 随机初始化均值，维度为(n_cluster, n_feature)\n",
    "# 生成范围/2是为了限制初始均值的生成范围\n",
    "from scipy.stats import multivariate_normal    # 生成多维概率分布的方法\n",
    "mu = np.random.randint(X_train.min()/2, X_train.max() /\n",
    "                       2, size=(n_cluster, n_feature))\n",
    "\n",
    "# 一个协方差矩阵的维度为(n_feature,n_feature)\n",
    "# 多个分布的协方差矩阵维度为(n_cluster,n_feature,n_feature)\n",
    "cov = np.zeros((n_cluster, n_feature, n_feature))\n",
    "for dim in range(len(cov)):\n",
    "    np.fill_diagonal(cov[dim], 1)\n",
    "\n",
    "# 初始均匀的类分布概率\n",
    "pi = np.ones(n_cluster)/n_cluster\n",
    "\n",
    "# 概率矩阵\n",
    "P_mat = np.zeros((n_samples, n_cluster))\n",
    "\n",
    "\n",
    "max_iter = 20\n",
    "for i in range(max_iter):\n",
    "    # 对每一组参数进行计算\n",
    "    for k in range(n_cluster):\n",
    "        # 实时生成高斯分布，免去了存储\n",
    "        g = multivariate_normal(mean=mu[k], cov=cov[k])\n",
    "\n",
    "    # E-step，计算概率\n",
    "        # 计算X在各分布下出现的频率\n",
    "        P_mat[:, k] = pi[k]*g.pdf(X_train)\n",
    "\n",
    "    # 计算各样本出现的总频率\n",
    "    totol_N = P_mat.sum(axis=1)\n",
    "    # 如果某一样本在各类中的出现频率和为0，则使用K来代替，相当于分配等概率\n",
    "    totol_N[totol_N == 0] = n_cluster\n",
    "    P_mat /= totol_N.reshape(-1, 1)\n",
    "\n",
    "    # M-step，更新参数\n",
    "    for k in range(n_cluster):\n",
    "        N_k = np.sum(P_mat[:, k], axis=0)    # 类出现的频率\n",
    "        mu[k] = (1/N_k)*np.sum(X_train * P_mat[:, k].reshape(-1, 1),\n",
    "                               axis=0)    # 该类的新均值\n",
    "        cov[k] = ((1/N_k)*np.dot((P_mat[:, k].reshape(-1, 1) * (X_train-mu[k])).T,\n",
    "                                 (X_train-mu[k])))\n",
    "        pi[k] = N_k/n_samples\n",
    "\n",
    "\n",
    "# 迭代更新好参数之后，开始预测未知数据的类\n",
    "pred_mat = np.zeros((X_test.shape[0], n_cluster))\n",
    "for k in range(n_cluster):\n",
    "    g = multivariate_normal(mean=mu[k], cov=cov[k])\n",
    "    pred_mat[:, k] = pi[k]*g.pdf(X_test)\n",
    "\n",
    "totol_N = pred_mat.sum(axis=1)\n",
    "totol_N[totol_N == 0] = n_cluster\n",
    "pred_mat /= totol_N.reshape(-1, 1)\n",
    "Y_pred = np.argmax(pred_mat, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3VdwXFma4Pf/uSZ9JhIJbwiAoPcsmmL5qvbdsz3dM6uxoZB2FRvbEdLoQRF62JFe9KKNmBdtxEZIK0WHNKGRtNqZ0Zjtnpme7mkz1V3dZWmKRU8ABAhvEukTaa45ejjJJEGCpooAARbOL4JBInHz3gtU93dPfuc73xFSSjRN07TPP2Ozb0DTNE17NnTA1zRN2yZ0wNc0TdsmdMDXNE3bJnTA1zRN2yZ0wNc0TdsmdMDXNE3bJnTA1zRN2yZ0wNc0TdsmrM2+gXu1t7fLoaGhzb4NTdO058q5c+fSUsqOxx23pQL+0NAQZ8+e3ezb0DRNe64IIW4/yXE6paNpmrZN6ICvaZq2TeiAr2matk3ogK9pmrZN6ICvadpzy3VcyvkyTt3Z7Ft5LmypKh1N07QnIaVk4vIk1z8cxXM9hGGw5+ROdh/fiWHocezD6ICvadpzZ3Zsnos/v0JbbwrLtvBcj6vv3iQQsBk6PLDZt7dl6UehpmnPndELEyTaEli2GrOalklrd5Kb52+ht219OB3wNU177lRKFeygveo1O2BRq9R1wH8EndLRNG1Lyi7mGTl/i9xinlRXkt0ndpLsaAGga6iT+VsLJDtbmseXcmU6+tp0Dv8R9G9G07QtJzOf5Zd/+T65hTzhWIjMfI53/vJ9sgs5APae2IlhGmTms6wUK2QX8jg1lwMv7dnkO9/adMDXNG3LufbBCKFYiHgqhmVbxFMxQtEQNz4aAyDaEuWN33qJXcd3EkmE2Xl4B2/9zsu0tCc2+c63Np3S0TTtmZNSsjiZZuLyJPWaQ9/uHnbs78UOqLx8di5Hqrd11XsiiQjLc9nm1+FYmP2ndz/T+34UKSXpmQyzo/MIA3p399DW04oQYrNvrUkHfE3TnrmRC7e4+u5NYskopmVy+VfXmbu1wEvfPIlpmcTb4lTLNcKxUPM91XKVeCq6iXf9aFffu8nI+VuEYyGklIxfmmLfi7u31ENJp3Q0TXumlqbSfPT3HxOOh4gkwoSiQTr621iezbAwmQZg3+ldFJaLVFdqAFRXapQyJfad2jrB816F5SJjF8fp2NFOoi1OS3uC9r4UI2fHKOfLm317TTrga5r2xJy6Q2Y+S2G5+KnLHz3X49xPPuEn//4dJq5Mc/2DEW58NIbruAAEQgEyjZRN91AnL37jBaTnk57JID2f07/2Al2Dj93jY1PklgoIYWAYd9M3hqnCaz5d3KzbeoBO6Wia9kRmRuf45OdXcR0PKSWtXS2c/MoxIvHwE73/9tVppm7M0DnQzuJkmkRbnEKmyOzoPAMH+nEcl2gi0jy+d1c3PcNdeK6HaZlbKhd+PztgrfkAlIAV2DphVo/wNU17rEKmyPkff0K0JUJ7X4qO/jbKuRXO/fjiE4/0Jy5PkuxoIRIPk+xIUMwUicQjLE6lKWZKWJZJoj3OxNUpJq5OUc6XEUJg2daWDvYAbX0pguEA5cJK87VSrkwkFqKtp/UR73y2ts6jR9O0LWt2dB7Ttlatbm3pSJCeWaaYLZFIxR97Dt/3MQ0VcoaPDTJ5fZb0VJpyYYVQNEhbX4p3v/chAhXcJZJjbx1iYH//xvxQ6ygQtDnzT05w/iefsDS9DEAiFePEV45iWuYm391dOuBrmvZY9Wod03owISAAz/Wf6BwDB/q5/uEoHf1tmJbJjn09xFrC9Ax3c+T1/fz0/3mHZGdLsz+O67h88vZVOvrbCMeeLG1UKVVUdU88TCgSfOKfbz0kO1p463dfpZRTn0xiyeiW+2SiA76maY/VPdTJrU8miadizSBWzJZZms1y7f2btPWmGNjfuyow+756ENxpdTB0eID0TIZbn9xmYXyRykqVSDzM0JFBNTErZTPYA1i2he9LMvM5+naHkVJSKVUxTOOBYO65Hpd/dZ3Ja9Pq/iTsPrGTfad3rwq69WqdW5cmmbk5ixWwGT46QN+ennVrx2AYxhN92tksOuBrmvZY7f1t7Njfx/SNWYLhAKV8mdELEwwc7KNSrDJ64Rbjl27z2m+ewbRNrn84yszNWQzTYPBgP3tP7SIQtDn48l4mrk7Ru7eHlrY40ZYIE5cmKSwlWXMmQKggmk8XuPj2lUbFi6R7qIvDr+8nHFV1+qMfTzBxZYr2vjYMQ+B7Ptc/HCWajLJjby+gPjG8/7fnKKSLxNvjeI7HuR9/QmG5yKFX9j+z3+Vm0gFf07THMgyDF754mB37elmcTDN6YZwDL+2hc0c7AJFEmPxSgWsf3KSYLVMt1WjtSiKl5NYnkyxOLrPn5E5ufXKblvYEqa5k89xtfSmWp5cRpqBWqRMMB4BGGsk0iCYjvPu9j7Bsi/a+FFJKlqaXOfujj3ntN880rjFBqjvZLIs0TIOW9gS3Lt5uBvzFyTS5xTwdjXsmaBOMBLn1yW2Gjw4+cdroeaYDvqZtY57nYZpPNqloGAadO9rp3NHOxOVJWrtXV5/EUzFGzo0Tb42S6EgwfXOW9EyG7EKelXyZxck0c+PzmJZF5LW7OXbTNDAskyOv7+fa+yMUMyWQEtM2OfX14+QW8zg1p9knRwhBa1cL6ellcksFEm0x3LqLYRiUcmV8zyeSCGPaJrVyrXl/2cU8gVDgvp9JAIJyoaIDvqZpn08Lk0uN4Fokloyy78Xd9A53P/H7Q7EQTs1pjsYB6jUHwxBIATc+HKW6UsP3fPLpAvWaS2G5wI79fYyeG2fiylSz5UC95hAI2uzY30fv7h5yjY6Yya4kgaDNzXNja1e6GAKn5mCaJom2OB/98DzCMACBYQhau5Oc+NKR5uHx1ugDe99KKZFSPvMJ3s2i6/A1bZtZnErz/t+cw/d82vvaAMGHP7jA/MTiquNKuTLpmeVVteV37DmxsznyzqcLLE4tszCxxP4zu8nO5ShmS0QTYcq5MuFYiEg0qEbjqTiJ9hhzY/OUsiWK2RK5hTyHX9+PaZoEgjadAx10DnQQaJSAtnYlcWqrA7Xn+SAl8dYovu9Tyq/g+4AQWAELp+4xP75wN32DmngORoLk0wWklHiux/JMhr7d3cSSW7dHz3rSI3xN22ZGzt0ikgg3R+fCEBSWS/zF//Q3nPjKUQYO9DMzOsf1D0aoVxxC0QCHXz/AsTcPNatZBvb3k1sq8JP/+x1qK1WEEHQMtFNYLjJ1fZr522nirRFqFQdhCDr62giEAviez4GX9jJyfhwpBK2dCU58+SjtvamH3m9bbyu9u7qZvbVAtCWC7/msFCocfHkv4ViY7GIez/E4+dVjZOZzVAoVoskoSMgu5ujobwNU64ZXvnWaq+/fZGFiEcMy2f3CTvacHN74X/oWoQO+pm0j9ZrD1fdusFKqAtDSFqeYLYOUODWHcmGF//g//z2Lt5dI9bTiux6T1/Nc+uV18ukir377NHbARgjBSq7CgTO7CMfC2EGLarnOz//iA7oGO3A9H7fu4hYqarK1v43cYp70TIZyfoWuoU6+8LuvNNshP4phGJz48lG6x+aZGZnHCpgcf+tQc/Tuez5CCEzLJBoPE4mFiLZEKGZKeI636lyxZJQXv/4CnushDLHtdsfSAV/TtgkpJef+4SLVUg3LNgnHwkyPzFFIFxk6tINIIkI0EWFhYpGVwgptva2kpzMYloHruPz4T97GNA1e/vVTOHWX9MwybX2pZp378uw8obCNYQi6Bjtx6y6Dhwa4fek245dvU1+pU686RBIhnGqd975/lpd+/VQzdfMoaqFWHzv29T3wvZb2OLVKjfPNNg8CwzRo7Urw4q+deOj5tqOnDvhCiB3A/wV0Az7wXSnlvxVCpIA/A4aACeB3pJTZh51H07TPznVcFm4vkZ7OEE6E6Nvds6oRGUA+XWBpOs3+l/Zw9b0bVMtVpK9G9tnFPKe+egzfl5RzKwRCNnPji4SjQQzTxLRNKqUqmfkcs2PzKk0iRDPYu3WXiSvTzE8skF0ssPPIAMGQTSlbJtoaJdnZQqqrhURbgngqhmmZpKeXmb4xy/DRwaf62e/k452qgx2yEYagVq6xUqgQSXz+K28+jfUY4bvAfyulPC+EiAPnhBA/Bv458FMp5R8JIf4Q+EPgX63D9TRNu4dTd/jwBxdYns0QjARxHZeRc+O89M2Tqxp3Vcs1hDAIx4IM7O9nYXIJkJiWyeCBfuKtMXzfI9oaJb+YJ5IIYzRKNmvlGi1tcSLxMIuTywzs7yfV3UIxUyLWGmX0wjhOrY70JF2DHdQrdQzL5PBr+ynly0gfOgfaV913NBllbmz+qQN+dj5HJB7hxFePkVvM47keLW1xKuUambksfbt7nur8nydPncCSUs5JKc83/l0ErgF9wLeBP2kc9ifAbzzttTRNe9DMyBzLs5nm5hup7lbCsRAX3768qpNlJBEmt5jnws8uM3ljhtpKjVRXknhbjEqpwsf/eJkPf3gRz3FVDftKnepKjXJhBd/32XN6N07NJRxTJYzH3jqMYQomr88we2ueSCLSXBhlBSwKS3mmbsxy+PUDmJaB763uueM6LsHow8shfd+nXFihfl+Fzv0810eiSiu7hzrp291DrFW1gHDvy+Fvd+uawxdCDAEvAB8AXVLKOVAPBSFE53peS9M0Ze7WgqpKuUc4FiI9s0ylVG32qzcMg2KmRL3qEG+NIX3J7etTpGeyLNxeItmeoGuogwNn9pKZz5FPFyjny8RSMXYeGqC9t5VyttzMo8eSUd763Ve5/uEotZUafbt7CIRslmczZOfzJNriHHptP8OHByllyty+Ok1bb2sjELtUilWGvrRjzZ9p9tY8l395ndpKHcNQfXj2v7hnzdx7sqsFAxXcLVt93/N8pJSkupMPHL+drVvAF0LEgL8E/hspZeFJu8QJIb4DfAdgYGBgvW5H07aFhdtLXHnvJsuzWXYe3kHvri5MSzUdE4hmAASYn1hkx/5ePNdn4fYSkzdmyC0UyC8WCMeDlPMVpIDjXzxC50A7oUgACaRnMgTDAdy6y4u/doJEWxzP9fBcj9xSgXAsSDQeJhwPYxiC7p1ddO/sYmkqzeAB1dr4wEt7cOousyNzIASmZXDsrUONdQCrZRdynP3hxyTaEyrN5PmMfjwBsGbPm3A0xJE3DnDx7SvNB4Lrehx8eS/x1tj6/9KfY+sS8IUQNirY/3sp5V81Xl4QQvQ0Rvc9wOJa75VSfhf4LsCpU6c+3Z5pmraNXfrlNf7m3/2Ies1heSbD7Ng8O/b18uI3XqCQLtG3t2dVKwGn6iCEoHOgnUg8zOS1aXxfpUMisQhWwGT6+iwf/uC82l5QSt76nVeprtRw6y7heIh6pc65n3zCrU8muH1pimgyQvdQJ+nZLJmFPLuOD2GYBsXlIh072mnrVXMIdsDm5JePcuDMHurVOtGWyENLMieuTBGMBAkEbaQvMUyDtt4UE5en2Htq15rvGzy4g9buJIuTaaSUdO5ob7Zi0O5ajyodAfwfwDUp5b+551vfB/4Z8EeNv7/3tNfSNE2plCv88I9/RrK7hXA0TEtHC4uTi9y6eJtEa5QTXz7GoVf2NY+fG1/g5vkxLv3iOm19rYTjIXLpItLzVWAFnLqLHbSYuj7NrqND7D29C0C1HYgE8VyP9/72HNVSlcxcjkhLBMfxyCzmOfjyPm5duo1TrROKhzn82gEGDvQ9UOceiYcfuyViIVNi4fYSN8+NIX1JW28r/Xt78X0fp+Y+9EGRSMW3dGvirWA9RvivAv8ZcEkI8XHjtf8eFej/XAjxL4BJ4LfX4VqapgHztxapV+u0R1VKJNmRIJ6KsjC+SKonxamvHUdKSSFTZGl6mUu/uEaiLc6u40PM3JxlZmSOYqaIbVvEWqPkFnKYloFb93BqDoXlAjsP382vSym5fXWaxdtpUj1JquUqibYETs1hbmyeUDRIsiNB395ejrx2YM17zi3luXVxgny6RFtvK8NHBx9oaSClJDuXZfrGLF1DnQhDkJ3Pk13IsefELkKPmOTVHu+pA76U8pfAwxL2X3ra82vadrNSrDw27WEFLQQCKWWzFt40TUzbItERb+5BW8iUGLs4gWWZHHh5L4MH+0l1Jxk5P8HU9RmK+RWCVZtEe5xKqYIpIdXdSiwVI9wYiVdKFc7/5BIj528xOzbP1A2bSrGC9CWLU2mqpRqGZWIIQSgaXDPgp2eWefd7HxGMBAlFg8yMzDF9c5bX/umZVaPywnIR35e09aUo5crq04UhyC8W6BnubH5iKGZLTF6bppQt096fon9vL8Gwfhg8jl5pq2lbhFN3+OTnV5kdW0AI1dP90Cv7GDz4YCVL91AnbX0pMvM5Ut1JhBDUKlU8x+XAmT188LfnQUBHfxuT16YxTIMbH41x9I2DrBQrOLU6e07t4taFcXLpArVKjVRvK6muJLuP7yQQClDKlYm3xjj/k0sUMyU6B9oZ/XgcZ7HA/OQShmHQNdROMBKktbNFLcyay1HOl4m2rB65X33vJtGWaHMhVCAUIL9UYOzibV74wuHmcdVyDTtoc+DMHtIzGXJLBeJtMboGO4g0FpJl5rO8+72PVLO1cIDFqWXGL0/x6m+82NwQRVvb9mokoWlb2NV3bzA7tkBbbyttvSniqRgf/+Nl0rOZB461Azbf+q++Rrw1yvytBWbH5smnS3zjX36JYChApVxtpkvirepvp+6QX8ozMzKHYQqk59O9q5tkRwuhSJCWtgRH3zjInpPD3CmyK2ZLZOaytHQkWJ7N4Hs+hqW28atXasyMLCClpLZSV10nWyMsz+VW3avneuTThQdWvUaTEZam0qteiyTCzTr+nuEuDpzZw87DAwRCNolUDCkll965TjgWJtnVQiQRpr0vRa1cY+LK1Hr9p/jc0iN8TdsC6tU6k9dnSfW0NlM0lm0RjoWZvDq1ZjfJnp1d/Bf/4+8zP67y+T27uojGI8yMzq06rndPD1ffvUG1WGXyxgwj525RzJXoHuqia6gd2VgQJUyD3t3dVMs14q0x4q0xcot5hCGoV+tkFvIMHxuimC5Sq8ySaIvT1peid3c3e0/uIpqIkJ5ebu46dYfagzZEvVpfVTVUW6k/kMOPt8bYsa+XyWsztHQkMAxBPl0k1Z2krS/VnF+4v5wzmoyyMLHIgRf3fPb/CNuADviatgV4rqfSOPcFSytgUSnXH/q+QNBmYP/qhmJ3yhF9z8cwDWItUQ6+tJd3/uoDnKpDMGxjWAmEgOXZLIlUnPnbiwTDASYuTdIz3M2JrxxFCNHse1MpVUGqeYKWzhb8Rv+acDRErCWKHVQLrmrVOm19qx9OQgj2nh7mwk8vk+pOYgdtais1yvkVjr116IGf6eibB0l2tDB+eRKn6rHn5DDDRwbUzlxSNT7zXG/VIiy1I5au0HkcHfA1bQsIRUNEEmEqpSrh2N08dDlXZueRT7cgMZaMsufkMDc/GiUYCSKEILdUoG9PD0deP8DohVuc/dFFnLrL4tQywXAA0zYIRoJ07ezkC7//KpatQoNlWxx76xAf/v0FStkS9WodwzRIdSfpGe7i/E8v4dQdpm7MYhiCgQP9XHz7Cie+dGTVaH5gfz/ShxsfjVJIFwnFQ5z++vFmr/p7mabJziMDa/7cpmUyfHSQGx+N0d6XwjANnJrDSqHC8XvmArS16YCvaZugsFxk5Pwt0jMZEm0x9pwY5thbh3nvb85SKVaxQ5baCLy7hb49n7751/7Tu+noSzE7No/0YfDQDj555yrz44sUMmWqK3WWZzNUGk3Rjr51iMx8lg//7jylbJnhY0PsO7WLSGMFrWmZBCMhRs+Pqc6ZxS6yC3l2HhmgVqnT2Zci3hbHsi3S08tcee/mqslYIQRDh3YwcKAPt+5iBazP3It+z8lhPNdj/NIkAKZt8sKXDq/58NBWE/c2V9psp06dkmfPnt3s29C0DVVYLvKLv3wfO2ARSUSa6Y2XvnmSWDLC9Mg8K4UV2vtS9Ax3NUfbT+P6R6P85b/5G+xwgLmROZbncxTSBeyQTSgSRPqS3r099A520bunm3hrDDts8/Kvn+IXf/EegVCA6RuzjTLMKk7dZd/p3WTmc+w+PkSys6V5Ld+XZOezfPWff+GJet1/VvWq6q8fjoW2bX/7O4QQ56SUpx53nB7ha9ozNvrxOJatNt4GsFoimJbJ1fdu8Nbvvsq+U7vW9Xq1So2xC7fo39fLyNlb5JeLVFequI7X7Lnj1ByqpRpmwMSpuyS7WliaWVbbHFYdgmG1J+2dydLCcoF4a5QbH42STxfo2dlJ984uWrta1EbmUj7QHXO9BUKBVWkj7fF0wNe0Zywzm23WlN8RigZZns3gud5Tj+jrNYfsfK7ZLTK/VMD3JLuP72Ts4wkV2G0LwzSIxMOYpoEvfYrLRaZGZhlqrLANBG3ySwUMw3hwE3HH4+bZMaSUVMo1nLr6eujwDmJJteFJKBKkXq0zcv4WU9dnQaiul7uPD63Lpxbt09O/dU17xhLtcXJLBezA3U6O9WqdUDT41KmJxak0Z3/0Ma7TqPoxzWYAd2pqR6hYawQpwa+pfV0l4LsSYRhIH4IRtWK1XnMYPNTPyLlxAuEAAtWjHqCUKxNNRkl1JTEDJrVyFdMyufHhKIdfP8CL33gBz/P44AfnyS8VaOloQUrJzbNjFNIFTn/9BZ60o662fnTA17RnbPcLO/nlX32AaZlE4mFqKzXyS4VmKeRnVa/W+eiHHxNtiRAMq1RHveYw9vEExWyJ6x+OkJnLghAsz2WxLJNoS0T1xw9axFojdO5ow627ZOdzxBIRdh3bCQhufjRGoiPB7OicmsCNBfFdj9beFHtODFPMlSlmSlSKFV76JydpaU+wOJUmu5BfNZna0d/G/MQi+XSBZEfLQ34SbaPogK9pz1iqu5WXv3Waa+/fJD2zTDge5uTXjtG/p/epzrs8l8VzvWawB5WWqa3UWMmvEIqFCMfCCCGoVxxq1RrxVAzP9ekZ7mBg/w5mR+Zwag47jwyy58ROAkGb/ad307mjndlb8wwfHaC2UmN5Nkttpc6u44OYlkVrZwstbXGyCzkSbeqTSylXxjAfrMQRQlApVnXA3wQ64GvaJujob6Pjt17G8zwMw1hzZC+lpLBcxHM9Eo2Sx0dao+BOSsni1DKReIgTL+8jczjH2IVx+nb3MH7pNoZpMHCgn72nhnGqLr/+B19n/+ndq84hhKCtp3XV/rjlwgpv/9m7VEo1oi0mruORnc+y9+Su5kRqrCWy5sStlJJwXPe82Qw64GvaJjLNtXP25XyZj370McXlEqBW3B7/4mF6dnY99FzJRoWMU3exAxYrhQojF8aYuDRJS2cc35MMHx3g1NdfoJgp0daXYs/JnUjXJxAKMHhoxwMbjT9MNBHh1d84zZV3b7A8k8EO2Rx8Zd+qDcnb+lK0drWwPJsh2dmC70tyi3m6hzr15iSbRNfha9oW4/s+v/j/3qNWqTdLN+vVOoVMiS/+3qsPdKK81/TILB//9DKO43Hz7CjCEPTu6mZpeplA0MayTQ6/fgDX8VjJr/CV//xNnLpLIV3ECliNsspPtyDKdVwM01jzfbVKjdEL40xem0EY4nNRpSNlDfwsYIPRihCb34NS1+Fr2nOqsFykmCnRfs9kZyAUwDAEc+OL7D6+86Hv7d/TS6q7lWsf3CSfLjB4sJ9AKIDvS5ZnMxSzJWZuzhKORzj5laPcvjrFtfdGkEiQEG2NcuYbLzzyoXK/RwXvYDjIoVf2r7kX7fPId0bBOQuykaoykhB8A2E8H3vnbv6jSdO0VVzHA+PBnP5a9fBricTDtPWmaOtpJRQJYhiCXceG2P/iHlraE/Tt7uHN33mZYCTI5V9eJ9nVQntfG+39bThVh/M/vcRW+uS/VUg/A/X3QaQQZjfC7AZZQdZ+9dz8vnTA17QtpqU9jmmuDu5SSupVh44dT5ZjT3Yk1GpXXwUiwxAkUjH69/Rw/IuHSaTizIzOEQyvrv1PtMXJLeQp51fW94f6HJDuBIggQtz9RCOMVvCXQRY278Y+BZ3S0bQtxg7YHHvrEOd//AmGaWJYBvVKncFD/asqZR6lpT3B4MF+xi9Nqk1FfEm1XGPvqeFmusZz1cKr+0nuLrDaKE7dYfLaDNM357CDFkOHdtAz3PVMF2NJWQc/A1hgpB6ai5dSgiyDV2LtMbIAvA280/WjA76mbUF9u3uIp2LMjy/i1F26BtpXbY7yOEIIjrx+gO6hTqZvzmKYBv17e2m/p1d9365uJq9ME0/FmuetlKpEEuEHNiZZT57r8eEPLrA8lyWeilEpVvnw7y+w7/QuDpzZu2HXvZfv3AbnA5CNQG3EIfg6wli9NkD6WWTtA5AZ8HLgzSMDZxCG2r1LygqIEIjnY02BDviatkUlUvFVG3x/WoZh0DXYQddgx5rfb+9vY/DwDm5fmcayTXzfx7ItXvrmyc/cuvhJLE6lSc9m6LwnPRWKBhm9MM7Q4YEN35dW+nmo/0qN6o1A47UCsvYOhH6tOdKXsoas/gwwEUY3kg6QRaj/CmkfBSFBCgi+hRDPR7dOHfA1bZsyDINjbx5iYH8f2YU8gZBN50A7wXBwQ6+bW8g/0OXyzorccq688QHfmwEMhLh7D8JIIL15VW5pquoo6c6BrKnJWUAYJtI+Ce4omD1g9iDMPoSxcZ+G1psO+Jq2jQkhSHW3kup+srmB9RBpieDWH6w2klISCD+DdsfSYc1cvLg/F18BzPsOMZBGC8LejTBXby35PNABX9O0Z6p7qIPrHwQoZkrEWqNIX5JdyNM12PlUKawnJcwepHMJKf170jd11OTt3QefMFJIHKSUzTkOKX1APjJnL2UV6VwG95ZK+VhdYO5EmO0IEd7IH+2xdMDXNO2ZCoaDvPytU1z+5XUys1kwBAMH+th/Zs+GXldKD/wFpJcDIwXeLFKEgEYQD76GEPfs0GV0gDkE3jhSJNRxsgz20YcutJLSQ1bfBj8HRgxkjMwtAAAgAElEQVTca1B/F4wWpLkXGTyGsA6umnyX/grIPBBoVAttXKWSDviapj1ziVScV751urkp+ka3WpCyhqz9HLwlwFJpHREAaycYMYTZjzBWf7oQwoDgy0i3H7wJECbCOgPGI/YY9hfAX258irgCsgLWoJobECGon1cPG7MHKaU6xrl05y7V94JvIIzIw6/xFHTA1zRt0zyrLQqlcw08FYibr3mLABj2AaRfwq+9B96UCszWAYS1CyFMhD0E9tCTXccvAiaSGvhpEMnGdwTgAFFk7UOktUs9dJxPwOxtVvlIbxlZ/wAR+sJ6/eirrEvAF0L8MfBNYFFKebjxWgr4M2AImAB+R0qZXY/raZqmfSrumBo938toA+8Wvn8Iqj9WAdhoBRyov6/KNwPHV62sfRxhJJB4d3vtNNMzPkhTpXf8WTC7G6P+Ngh30JwcNlLgzSH9lQ0Z5a9Xse3/CXz9vtf+EPiplHIP8NPG15qmac+eMHhwwwAJGOBNAtXGpKoJ0gV3EVb+X+TKn+LXzyOlg5S+Wojl5x7eO8foBLMD/AIQAq8IXkbNB7g3wRsFcxisARAxde3aB0i/rG5TCPVhgI1Z6bwuI3wp5S+EEEP3vfxt4K3Gv/8EeBv4V+txPU3TtE/F2gf180ij5+6kqEyDdajR6ljV/ktZA+ciYIKIg4iAex3pLoBwoBGYMRJqkve+lblCmBB8U6WQcMC5pj5JiFZwfwVGN1h3PmnY6tr1cyA9pJkCcwcYLSA2prZ/I3P4XVLKOQAp5ZwQonOtg4QQ3wG+AzAwMLCBt6Np2nYlrD1ILwPebaQwVMrF7EPYB5DuLTUpS0tjUtdtBNw6GBHwQ1D9W7APg9WLEPHGytyfN1bmrg6jQgQRgeMQOI6UFaQ7rRqsOX2NRV0G+BWQuUZ5ZxmECd4syDIi/gcbVqmz6ZO2UsrvAt8FtQHKJt+OpmnPGSl9pDsK7nWQNTAHEfahVStghbARodeQ/kE1SjeiIFRvImn2QfUdpDOpyi6pqWOsPvAdqH8A/hx4UfDnkEY32HvBW1QTs42VuGsRIoyw9yDlLmT9onqwGOFGsLfADIIxBNaQqhrCA+yHnu9pbWR75AUhRA9A4+/FDbyWpmnblHQ+VkEZS1XFuBPI2k9VeuY+wkghrB3qbyGQ0oX6RypvLoSq0ql/BF5B5d6rPwBnCURY5edFayP4L9y5+BPdoxAGhL6uHjTejPok4WXUJ4nQ6whrsFFBZAPuev1qHrCRI/zvA/8M+KPG39/bwGtpmrYNSX8FnBtg9Nxtb2y2I/15pDuFsHc/+v3uDHizCGsXEgf8FVRYrKgg7xeBGtCOyusLkFFwZ8DqX7Uy93EMexBp/JfI+nl1zzINgRcR4s78QVWN8sXG7fe7XmWZ/wE1QdsuhJgG/gdUoP9zIcS/ACaB316Pa2ma9vkhZRW8RaR0EWY7yBrSHQG/pOrT7d3NgLj2CVZAsEYv+1Cj1/09h/oF1e2S0N0Vrf4MYKtJWX8JqIG9D9zbKj1khIAONdkqMyCD6nU8CHzzibY2lLLeSDmNAwGw9kDwLai/B94kUoRVZRA+BN/c0M6b61Wl8/sP+daX1uP8mqZ9/khvQU18SgcQSH8Z/CrYu0AEwbmE9MYh9FWEeEgHTxEGKe/ri+OrVgVib/Nr6ZwD5yYqiy3B7EAGXlOpldp7qupGlsCfV/1vjBQETkHdUuWURhKMHSrYexNgDoC3gDRaQSTUillvorF6dz/C2oMQJlK6yNrbdxdhySrU3wH7MCL4Kvi7kO4siJBK6xgbN7qHLTBpq2na9qMC4TsgoggjjMQHZwSVPjFVkzEzjPQXkO4thH1gzfMII4q094JzA2m0qclW5wLggDTx/Zwqt6//CoxesNoQIoB0b0P9f1MNzvyFxmStBW4WxDLYL6hcuywAnppkddKNVgkDEBgGP4+s/LCxYCvZSO+4UD+LlCVE4BR48+AtIszexg2DlGFwroO1V7VXfoZdN3XA1zTt2fMzgINakI8aOQsJBNX3mn1tYuDNwX0BX207mANZB+uAmvx0PlY17WYP2PsBG2r/qFIpRg+wAvUJpHVIlUD6C2Am1ei+fh7Vy6ZFBXlnUqV4zD4IfU1Ntjq3wRsD+xhCRFXw9vLg3QD7awhM9eAiCPUL+OY+9bPc9+lECEMtAZNF4Nn20tcBX9O0DSP9DLJ+WQVX0QL2IQyrD1UWc08V9p1a9gdWsFZBrC57lLKKrL2rRuD+CuBC4CWwjoIUCFMt+ZHeQiPfLlUu3oirEbpzoXH9O3MDqiZfPTwG1NdSqjRM4EQzzSKNefBjat6AxgNJqE8kyKrqo+OOqOv5RRB/C9ZhVA+dB34zqmfPM6YDvqZpG0L6OWTlxyDsu/nr2s/w5RsIqx9EGOmX79bLSx+cK4CjHgVGG0jngUobWT8L7g3wpsGdUwG49o8qDWPtvHugv6RG/kZCrWiV1cZm5Dk1Yrd2qff6y40aeBqTpy5Yu8GpqHu/Q0TU9+59TYZVSsevqnsSCVQ1D4AE7zYQRPrZRiM1X1XnmH0II8mztpF1+JqmbWPSuaZaChtJ1XXSiKog7n4MGIjA60BNTVrWfqlSKfZxwAD3E7WQKvAW4p7SRymrKtfvTKuFTyLQWPgUhepHUP0Z0p1F4nJn5A02eHk1qeplVfD35sHoUqtnsRoblC+AdwukB+6kSse400g/r3rpi5D6IxudLaUHwlftELxJkI1w6mcb2x92q58p8KKaBPYXQS6DuRsRfPlZ/qdo0iN8TdM2hp9WDcLuIURIbUCCo8oww99SFS61HAReVt/HVUHXW0aY96U9pNvI3VfUMcJQo3aZB6qqaVntV+B1AuFGr/koGF6jN05AjeyNGLhXIHBaBeTqj9Vx1gtgiMangwA4Z8GNA4Y6Lvovwb2h9r8VJgSOgbkTKn+tHkDY6vzNiViBMKKI0Bcau2oZn6r75nrTAV/TtI1htKkJ13s2C5eypkopG6FHiABSRMBob9bbCywQFtIwVT3+vW2NRRSINFI0jRG+X1IjedEN1g4VgP05tWo28Lr6RCEt9WAwAhA8DQRUCsZIgdkJ4osqPURFNar0llSKSUoInlDdLoWHMLtUGwTqqGoiS22uYu0A9ypIu/mzSVlRE7aN7RDv3TR9s+iAr2nahhD2fqQ7oSYzRQw1As9A4JVVC6WEEW+kYO4jfTBiSOk0KloCCCOGDL4I1b8D3wXjzkKoYCPoh8DeCa6tSjTNTvAHAKdR+VOjEdHB7EAEX0MIE98vg90Ghg3VXzRKLVOqNt+bVvdjtiC9BQx7p7oe4PsrsPLn4I2rOWh/rvGw6QRreMMXUn1aOuBrmrYhhJGC0FeQzsVGlU4cAm9g2IOrDzS6wGhXO1AZbaiJzWUwe5FeAZyfqQAsXaSwgACYh0GON+rnA+qPYYHRgRBhpHcb5Dx4N0FU1CcNYx9gqAeEzKm9ae8EY2sn1D4AZxGcq2pPWj8DZq+6J3++8cnk7oNJyjpU/qMqBzW6wDBVSkl0qXRP8E0M6+GN1TaDDviapm0YYXYgzC8jpXxoy18hTAi9pXL57hggVDmj2Qm1n6pdoWQWar9otDfwVDVMYJ9KofizQEkdZw4gvWwjn97eWAzVqiZanREwu1T+3z6EsA/evQdrCFl7D+o/V5uWUG9U7Ew3mqbZIGtq3qFBulPgTqsHghEGwo2Vv1Uwh9nAvcg/Mx3wNU3bcI/r7y5ECBE4CYGTgOp7IyvfU+0QjFRjUxIDjH7w0mpCuJ4Gay/Yg410ygJQBndWlVWKkCq5JKjKMGUKwv8JInhi7VYNsq5y9b4EGVA19tJTpaLWIIS/uapiCH9OrbC9k/KBxoOhrB4qm1Bn/zg64GuatqVIP4es/gO4U6h2B+MqqJq7G3vENiZD/ZIa+fshNRqP/KcIqx/pjqvmZ2YKvGX1cECC2YWw9z4k2DfSPMJUpZZGkOYWhZRUeib44ur3iIiamzAaXTVFTN2HX4DA4KbU2T+OrsPXNG1Lkc5VwFR5dXFnQxBb5dT9kgrORj+Y7Woxk7Vb7SPr3kT66UZJZL2xatdVo3xvBpxxpDOiJoHvJ0JqktjLqVE+vpqwFZZKFQUOP/gWa1gtsLL2qdSRn1WfMuxjiNCbG/kr+sz0CF/TtK3Fm1MTvAbgp9QoW3rgT6nWBndq7xFgH1ELmvwcYENdqrSKtV/V0DvX1chbRME+BN4Y0gkgAidWXVIICxk4CpWfgFFqlFMKoNZIBz1YaSOMJDL4FjgfqLSTSII1hAi+uCVKMNeiA76maVuLSICsIEQMaR9SVZTuiGpjYLSp1gTeBATOoLpiOkBI9c83O1UbA7miKmz8FTX6NnoRho2UMfVJwD6CuLdFAqgySrNLlYDe6YUvulTtvlh7ZyvD6kWa31YPIGGpLp9bmA74mqY9E1JWGpuEo2rgHxYc7UNQ+ylSWo1Ok3kIHGo0UauorQEpgRlTXS+h0Z64Tf1bJNSWh15OdeCUefCXkfYhhAgifR9VXnn/3rEW2DuAQKOZm6GOk5U1jr1LCEN9InkO6ICvadqG891JqL2rJkSlVCtpAy9jWAMPHGtYvfjyDbVC1suoXLp1TPWnaVT7SG8B3Ak1cWq0qVTKnZp6f0n1xTEPgJxSDwO/AM4I0tqpWiBzZ1tBT+2w5d5sPIzCjRROoXEzKRD9CCOy8b+kZ0AHfE3TNpT0V1SwN5LN3LaUdai/i3zISN+wB1W7AurIShuqd/49pZ3ChtBbasFT/Vc02lOq3a7cG2AOQ6BPlW76OVVm6d8GowUR/sbdB0f9rEoXGW1qgxTGwHdUusiw1SSxMFUt/ueADviapm0sfwmEv2oiU4gA0vfUqHqNUb46xgBCyMAxKP8HpMw2Ui0JsLrBHFI5etEK3jhSxAABZj+IDgQBZOC4WoTlZ9XEa+ird/vl+0W10MtofHIQIIOvqnSQO67KOs0ORODMg/n+55QO+JqmbVlSysYIPKQWTvkrQBacFah/AlQb+fse1fc+9KZ6iNTfBRIIbDD7kCIOIqRaFjdPrip97v3kIEQYaR8FcwBhH7rbq/9zQgd8TdM2ltEBUiCl0xwpq1bBpqqffwTpXIPqD1TLAr/Y6KlTUp0wjcHG1ocrquWCHwZvDmENI70Z1U8HG/BABBHBM/elhSIg5INtH2QdjM7PXbAHHfA1TdtgwoggA6+onL2QgFTpleArjyxjlH5G9dLBBKxGWqaG6nhpNconl9QGKO4yiDLULzRq4V8Ff486B2GE1fvAClthJJDmILgTSKNdXcfPghFHWL0b9vvYTDrga5q24Qx7EGl2qHw+qK6Wj6l8kc7NxqYlS6pFgggChiq3BPU9P6P6zfszUF8Gs61RRXkGwx5ancJZgwicQYqEmujFVQ8L+8iWXTj1tHTA1zTtmRBGRKVhnpRfAKO9MbKfQpVSWkBNTdQKGjtjTaqRv5lSKR6RbFQAJR/bz0YICxE4AoEjj+zo+Xmhe+lomrZppKwjvTTSzz/4TbNHpW3sg6oiRzbaFosu9TWykc+fUw3MrGEwOxvzBBbSnfxU9/J5D/agR/iapm0S3xlT/W7wAYk0uhHBV+5udWjtQnq31Ag/cLSxn20BQq8CNfBCqj5eZlTDMiNx9+TCaLRc0O614QFfCPF14N+iZl7+dynlH230NTVN29qktwT19xt72drN12TtQ0ToDaCRAgp+BemOqhYK4S8BQRB1EIlGt8qI6pt/z4SslBJkFWH1rXXpbW1DA75Qa53/F+ArwDTwkRDi+1LKqxt5XU3Ttjbpjqu6+HsXNBnt4N3G90sYRgxQQV8EjgJHH36uwMtQ+wUSgcpS18Ha87lZHbueNnqE/yIwKqW8BSCE+FPg24AO+Jq2rVVRH/oV6efUqldvDgA/cAph7Vu12fnDGFYf0vh1pDcF0kGYXaqOfhvk5D+tjQ74fcDUPV9PA2c2+Jqapm11Rn9jR6sEUpbA+QSkUAuxRDvUzyGlh1hj45G1CCOGMA5s7D1/Dmx0lc5aj1i56gAhviOEOCuEOLu0tLTBt6Np2lYgrAEwe5DeLNRvNHriSDD3IoygSse4V9fenUr7zDY64E8DO+75uh+YvfcAKeV3pZSnpJSnOjoevcxa07TPByEsRPANCL6h2hFbwxA4jTBbm98HT7U50NbNRgf8j4A9QoidQi1d+z3g+xt8TU3T1pn0i/j1S/i19/CdsUYvnKcjhKX64QdfBaN1VZsFKWuoipzQU19Hu2tDc/hSSlcI8V8DP0LN0PyxlPLKRl5T07T1Jb0lZPVn6gsRAHdcbRoS+sID/Wk+C2HtRrqjSC+tOl7KmlpkFXj17qYm2rrY8Dp8KeUPgB9s9HU0TVt/Ukpk/UMQkXu6R7Yg/XmkewthP/1EqTCiEPoa0r2u6u2NOMI+89g+ONqnp1faapr2cLICfkGVOq6SUFU26xDwoVFlEzi1LufSHk4HfE3THk5YIARS+vfVxDsq/bKBpKyDvwyIxopcHa6elv4Napr2UEIEkOYweKNIuhBCbWSCLCOslzfsur47qVovSLdxIyEIvokw2zbsmtuB7papadojicDxxv6xC0hvAWQeAhuXY5d+CWq/Uv1yzO7GdWxk7W3knQeA9pnoEb6maY8kRAARfAVpH0P1oo9t6AYh0psH5OpNz40o0iuBn1Y7XGmfiQ74mqY9EVWl8yz2efVYO/ngN/5on5VO6WiatqUIsxOEh5Re8zUp6yBMMHQO/2noEb6maVuKMFqR1lFwLiKxUO23ZGMh1tMv9NrOdMDXNG3LMQJHkFafmiTGRJi9iEaPfO2z0wFf07QtSRgphJHa7Nv4XNE5fE3TtG1CB3xN07RtQgd8TdO0bUIHfE3TtG1CB3xN07RtQgd8TdO0bUIHfE3TtG1CB3xN07RtQgd8TdO0bUIHfE3TtG1CB3xN07RtQgd8TdO0bUIHfE3TtG1CB3xN07RtQgd8TdO0bUIHfE3TtG1CB3xN07Rt4qkCvhDit4UQV4QQvhDi1H3f+++EEKNCiBtCiK893W1qmqZpT+tpR/iXgX8K/OLeF4UQB4HfAw4BXwf+nRDCfMprac8hx/Ooue5m34amaTzlnrZSymsAQoj7v/Vt4E+llDVgXAgxCrwIvPc019OeHzXX5eOFecYyy/hS0hWLcbq3j2QovNm3pmnb1kbl8PuAqXu+nm68pm0T701PMppJ0xGJ0h2NUajW+PHYGFXX2exb07Rt67EBXwjxEyHE5TX+fPtRb1vjNfmQ839HCHFWCHF2aWnpSe9b28Jy1QrThQLd0TiGEAghSIZCOJ7HVD6/2benadvWY1M6Usovf4bzTgM77vm6H5h9yPm/C3wX4NSpU2s+FLRnR0rJeC7L5cUFSk6dvliCo13dtIafPBVTdV2MNZ75lmlQcurrebuapn0KG5XS+T7we0KIoBBiJ7AH+HCDrqWtoxvpNO9M3sZA0BGOkl4p86OxEQq12hOfIxEM4SPx5ernd81zaQ9H1vuWNU17Qk9blvmbQohp4GXg74QQPwKQUl4B/hy4CvwQ+AMppfe0N6ttLNf3ubg4T1ckSti2MYSgJRSi6rp8sjD/xOeJ2DZHOruYKxUp1etUHIe5UpGuaIyeeGIDfwJN0x7laat0/hr464d8718D//ppzq89W1XXxfV9bFNV0BZqNa6nl8jVqowsp6m5Di/1DxANBB57rqNd3bSGwtxYTuP4Hid7ehluTWEZeq2fpm2Wpwr42udLyLKwDQPH8/CRXFyYwzZMQqZFbyzOcqXCz2+P8/XdezEeLMVdRQjBQDLJQDL5jO5e07TH0cMtrckyDI52drO4Uma2WMTxPDzpg4DeRIK2cIRMpcJyZWWzb1XTtM9Aj/C3oXK9TrFeI2RZDyyE2tfeTsAy+fvRESqOQ2c0ys5kiohtN4+pe3o6RtOeRzrgf47UXJfxXJbZQoFYMMjuVIrUPVUxUko+np/j6tIiIJBI+hMJXu4fIGip/ykIIRhuTfFbBw7xD2Mj9MTizZXUnu8D0KpXy2rac0kH/M+Jquvw47ExCrUasUCAdGWFkeU0bw0N05dQlTGT+RyXFhboicebOfiZYoGLC/O82Ne/6nwd0SiDyVYmclnigSC+9Ck5dU729K0a7Wua9vzQAX+T+FKSqVSoey7JUPipg+hYJkOhXqU7FgcgRoCK4/DR7DQ98QMYQnBjeZlkKLRqwrUjHGUss8wL3T3N6hwAQwhe3THAYEsLE/kcljDYlUo1z69p2vNHB/xNUK7X+fntCTKVFYQAJBzv7uFQZ9dnPudMsUAiEGx+LaXE8X3mSkUWSkUqrsu5uWlMDAZakrRHIgghMITAhwcWSQGYhsFgspXBZOtnvi9N07YOHfA3wfszU5TqNXoao2XP9zk3N0tbJPKZR9BRO0ChWiNiqxWtVxYXyVUrlGp1/tdigbBtEw8EmcjnyNeqDLYk2ZVqI1+r0h2NNXP4j+N4HhXXIWhaT/weTdO2Bv3/2GesVK8zVyzSHY01XzMNg6htcyub/cwBf297O7eyGSKuzUhmmWK9hkCwI5kkW6lQqjsMJFtZcR3y1SrX0mlMw6AlFOJU7+MbmUopuZFO8/HCHJ6UGAIOtndypKv7sTX5mqZtDTrgP2O+lIhGB8l7GcLA8T97uWNHJMobgzt5d+o2t3IZEoEg/YkEtmlSrFWxTZN0ucTx7h4ylQoTuSz72zs41dv3RCP1yXyOD2am6IrGsE0Tz/e5uDBPwDQ50NH5me9b07RnRwf8ZyweCBAPBCjV68QaLQqklBTrNU709DzVuQeTSVLhMGXHUcHeMFkolfAlCASulJjCoCMSxfV99ra1P3Fa5urSEq2hcHNi1zTUea4sLbK/vWOtTXA0Tdti9ErbZ0wIwSs7Bql6LvOlEksrZW5mlqk4DudmZ3l7YpzFUukznz8WCJAIBrm+tMRoZhlTCGxDkFlZoSsaQ0rJUqVMKhyhLfLknStXnDoBc/UulbZhUHPdNSd8NU3bevQIfxO0RyJ8a+9+ZooFbueyLK2UyVZq3MpmqLoef2cZfHPPfl4bGMT8lM3GrqeXyFYqzBQLuL6PL6ElFKS3JYEvJfPlEjtaWjjd2/+pcu/9LUnGsxk6ItHma4Vaja5Y7FPfo6Zpm0MH/E0Stm2GW1NcWlwgIAxWGm0MhBBkViq8Nz1FRyTKwc7V+fFyvc7tXI5srUrINAlZFjXPJRWOEAsEOD83y1CylYFkC5mVChXXpea6/N7ho42WxxCyPn3N/8H2DqbzeRbKJaJ2gBXXQQAnenrX6TeiadpG0wH/PqV6nXK9TjQQaObY14vjeYxlM9xcTgPQHgkzupzmdi5PMhxCovaGjAYC+NLnemZpVcDPVv7/9u49NrLrPuz499w7d+68SQ5nyOFrucvVah96Wlo9bBd2BLmN3CpRnSaA2z9iJAEMAw76b2IIaP8oDAQw0D/aNAX0R9AUSOq6aAoLTYxoZSSR3Fi2FEm2diXtLneXXL45T8577uv0jxlySS1XXIkrcan5fQBiZ+6de+fMwc6Pl+f87vm1OHd1Fi/wWarWmC0VKbVajMTjxCyLXCKBYRi9TB9jK+NntVFjuVYll0h+7Bu8krbNV0/cy5VyiXyjzvTgIMeH0iRte++DhRB3BQn4PX4Q8MbyEpdLxe7NSFpzYjjD2bHxOzJkobXmHxauc726QToSpdppc+7KFYrNBprupO1Gp82R1CCu75OybRxvZ9bOmyvLmMoAA0rtFoEOiIRCxMNhbNPkvXyefKtBpdXmnvQwuWQ39fP6RoVCo0k2HkdrzT0f83NtFjaBj3+DmBDi4EjA77lYKHCxWGAskdwK+BcLeVLh8B1JOyy2WlyvVhhPpNBo5tYq5BIJPN+j6bmYhqLWcSg2G4RMk7gdZmYovXW86/us9qpGXSisAxov0CTtMMVWEwNF3XGImCHKrQbvFwI6vovWsFqv889mThAJhQi05lIhT8Ky9nVnrxDi8JGA3/N+IU82Ft+ayDSUIhOL8V4xf1sBv9RqslKr0XI9XO2DhtFEgqnUQDcXvncjFEDL7Qb5oUiE4XicccPkcrlI3XEJdMBjE5Pk4knu3xaQTcPANA38zYyY3j++1jieh2WYRC2LoUgUrbovOL+2TswKcXZskpbn0nRdkrZNNhbnvUJeAr4QfUYCfo8T+CTUzjH7kDJuGlbZzfn1Nd5aWablebxfyONrzanhDEnbJh2N8fSxGWLbJkpNQ+H5Pk3Hpe25nBgZ5uGxca6WS6Rsm6eOzTCeTO1IgzSU4vRwll+urzEaT7Baq3WHgtoOUcvECwLiVpiYZfHI+DgGivVmg7bnMlsqEGjQaEzD4HQm033euwlMCNEfJOD3HB0YZG6jsiPtsNRucWzbsMpuNtpt3lpdYSQe563VFYajUcKmyUqtxvTAIKVWk6vlEiczWTKxGIvVDRarVS4VC6w3G9iGScv1yMbiDMdifPXEvbdcXuFMdoR8o8G7hTwtz2O90cT1fRTdXyL3jSQ5lk4zYEcItCbqdCg0Gji+Tzbe/VyuH/D60hL/4t6TEuyF6DMS8HseGM2x2qizWq8TNk0c3ycRDnP/HsM5640Giu7QSsNxGIp2i4P4BJTaLQbsCPMbFU5nR3h8YpLv/eRVLhTWaTouAJZhslKv0XA7bHRazJZK2GZo6zybqp02r87Pcz6/xmyxyFqjTi4eJx2LMRSJcqVcYjBicyQ1QMfzyLeajCUSuH6ArnczfCzTpNrusFavM1sokosnOZ5OyyJoQvQJ+ab3RC2Lx8YnWK7V8LVmJB5nIpnaMxiaxuaYv4FSCs8PyLcaLGxUaTouQ5EoD+e6Sya8l8/jaJ9Hxsa5VimTtG3KzRZN18Ht3ST1+tICc5UyX5o+ypGBbsbO/EaF//3uBRqOQ93pMBCJ0PF9UIpoKETYNHnu5Gl+sbbCSmrhpykAABMLSURBVL1Gwg7zxMQkMctivdHg0bFxCs0mc5Uybc8jE4+BoXhzdZlrlTJfmTkuQV+IPiDfcqDSbvF3c3PUnQ4KsEyTqdTewR4gl0hiKgM/CMglkt2xfNclblmMxuPkm01WG3XqnQ5XK2UiIQtDKYzeFK6poNBskPK7wzAaGI7F+dnSIiPxBD+5PsflYpF8s0HTcSi320RDJlErhOP7aA1Vp4OhFMfTwzx74l4GeiUIW667NQmdjceZLRWZTKWouw6j8QR1p8Or1+e4vlHhiclJTmVGpJqVEJ9hfX9PvB8E/O3cVdCasUSSXCJJwgrz99fnqDvOnsfHLIsvHz1Kw3UIGwZ1xyFAk45GaXk+Z7IjhE2TuUqZsGGStMK9yVIIAs16s0HHD0iGw1iGyWDE5mKhQKnZLVG4XKsxHIsRtUKEe4G+7rq0Pa+3IFoAWtNwOqzX6/zDwnV+cn2efKNB1LJ4cmKKYqvJQrVCpdOi7jocH0qTb9a5kM8TNgxansv7hQIvX52l43mfQq8LIQ5C31/hF1tNmo67Y6LUDoWgA8u1KvcOZ246xg8CfK23smjGkym+duoMl0tFXB0wkUjS8X2WazXmNiq0XZeNTpuHR3NU2i2W6zUSYZv1Wo31RhN0QKHVIhEOs1Kr03Ad5iplAh2QDEdI2mEUilQ4QrndIWKGQEPDdSAWJ9Dw5soKo70x+7V6jWuVMl+cOsJMOk0mHudquUil3eb40DARy+SnCwukIxFqjkMybHczf+p1FqtVjqc/fKJaCHE49X3A94IAdklWUShcf2dKpuv7nF9f42KxgK81o/E4j45NMBTtLht8fCjN26sr2CGLq5UypWaLgYiN1hrbDLFUqzE9ONQreBLG8T1Sje4SDpl4nJbrca1cJmFbHEkN0nI9rpSXeXA0R8IKs9KoETYU6+0OI7E4KdsGQzGRSuIHmjPbJphjvs8/Li9zZGCQlG3zcG68d74Svg6h6aaiaq3JJbp35EZCIdYbdQn4QnxG7SvgK6W+B/wa4ABXgN/RWld6+74D/B7gA/9Wa/03+2zrJ2I4GsNQ3eC+udZ7oDVe4DP6gfTIN5aXmC0VGYknMJWi2u5w7uosz957iphlYYdCPDw6xivz11iqVhmOxdlod4iHwxwdHKTYanF0YIDHJyZoui4/W1xkOBZDB5rleo1Ku4Vthuj4Pg/lciTDNj9ZmKfhuIzEExi6m9XzT6ameWJiiolUkolUip8tLVJutna0tZtp5NFwHAYiEQDOjk8QMgwurK9TbbeJJBI8lBsjZnXvP2j77k3ZQUKIz479juGfA+7XWj8IXAK+A6CUOgN8HbgPeAb4E6WUecuzHCA7FOKJ8SkKrSbrze769CuNGveNjDK8LfjVHYcr5RJjiSQho5uRMxCJ4AUBc+Xy1utOZ7M8OTmFHQqh0RwdHOThXA7LMLtj/K5LOhpjMjWAoRRnhkcwDINMLMZgNEI6FiUTjTOWTLFYqzIRT2KZBgGaWNgmG0/w1RP38vjkBKOJJLYZIhkO0/Z3jr1vTgBvv3nLMk0em5jk3zz4EL9x+j6mkgPELQutNeVWi7ARYnpg8BPvcyHEwdjXFb7W+qVtT18DfrP3+Dng+1rrDnBNKTULPA78dD/v90mZSacZjsVYqG7gBwHjyRSZWGzHjUltz8PXAWuNOtVOm2gozEg8jm2G2Oi0t16nlOJkJsuZ7AiZaAzTMKg5HdbadQrNJvdlR9Fa03BdMvEYy77H2bEJ5je6Q0CDkQjZeAJDKUqtFgPRKKcyWQYiNqYyaPkuL12Z5Rdr3WGYVNjmdDZLEAQ0XZeYZRFozVqjzon0MNFdsm7CpskXj0wzFInwbiFPudVkKBrjiYlJIpKeKcRn1p38dv8u8D97jyfo/gLYtNjbdtcaiES2hj52EzYMLhbyKBR2KITn15nfqDCWTG7l2W+KWhYPjOb4+eIClXabtXqdTuCRtGzOXZ3lf114h7rnYmiNoQxOZbLMDKVpez5LtQ2mBwdx/YCG4zAST5CNxwn1VrZ8a22ZjuvzwMgohlI0HIc3lpd4YnKK8+vrrDbqKODkcOamdm0XMgzuGc6wUNvA8X28IOD/LVxnoVrlC1NHtt5PCPHZsWfAV0q9DOR22fW81vqHvdc8D3jAn28etsvrd62Dp5T6JvBNgCNHjtxGkw/GfHWDwUiMSrtNSBnY4RD5RgPLNJhMpXa8tuN51NptKu02P124Tjxs8cjYGOPJQX548V3ansfncmNoFEvVDZZrG4ylkvzKsRlCCpZqNRzf46HRMZShtoJv1emwWqvz+ckjW/n18XCYeu/GrV87eYqW6xIyjNu6h+CtlWWqHYeJZLf9WmvmK2VyicSu2UlCiMNtz6igtf7Kh+1XSn0DeBZ4Wuut4qaLwNS2l00Cy7c4/wvACwBnz569a4ujzlfKnMlmqHdcFmtVOr7Hyczw1kSv1pqNThvXD7iwvsbVSolio4FpmAQa/nFlhdV6HS8ISIZt3CAgZdscGxxiuV7jkbFxMr11fB7tvWfbc3llfo6Veg1Ft6TgeCK5tS7OJsswaDgOhlLEb7Noi+v7zFcqO86llGIwEmW2VJSAL8Rn0H6zdJ4B/gD4sta6uW3Xi8BfKKX+IzAOnAB+vp/3Omi2GaLhOmTj8a0g6Qfd9XI6vs8r83OsN+p0fJ83lhZpex4bThvPD0hFEnQ8j1+srTEUiYLSW38Chczu1Xul1d4K+JsiIYt/OnMPhVYTx/OxTZOXrs4SaL2jHm3b82+54NpHpYDNX9u1Tod8s4GhFKPxxK7zAUKIw2O/Y/h/DNjAud4E52ta629prS8opX4AvEt3qOfbWuu91xm+i53KZPjbuWvEQhamYaB19y7Z05ksry8vUmw1ySWSrDcarDbqlFstjgwMsu40WK7WSEejWIZBudViPHWj1KDnB0C3sPlulFI7VvD8XG6M15eXiFthQoZBrdNhPJnYyqW/XZZpcmRggMVa9aYVQp+cnOJiscAbS4tb43CmMvjS9DQTqYGP9D5CiLvHfrN07vmQfd8Fvruf899NJlMDfG5snHfWVtF0h3CODgxybHCIH81eYjTeDbiFVhODbiaMF/iMJxOsNRqs1GqMp5JstNvYZgjH83EDl3yrwVNHZxi8zfz309kRhqJRZkslHN/j/pFRjgwMfKwyjI+MT7BxrcNKrdabddFMDwySicX468uXyMZuTBZ3PI9Xr8/ztVNnZKE1IQ4p+ebeJqUUD4yMciKdptZxiIRCJG2bSru1tR+6yxgfHUrz1soyTdcl0SuC0rZcHh4d4+mZ47y9usJsqUgkZPGbp+/nS9NHP1Jbcr01f/YrZlk8c88J1hsNWl63GlYmGuNSsYBS7MjUsUMh/E6LQqu5NckrhDhcJOB/RJGQRWRb9apk2CZqWbRcl6hlkbDCBIFmemCIlucQBEF33J4IZycmeWA0x0O5MVzfxzSMHWPxB8E0DMaSO395GMqgVydxJ612Tb8SQhwOkmy9h7rjMF+psFjt5qt/kGkYfH5yiqrTYb3RIGnbrDXq3DOc5guTRwibIcqtFo+NT/D4+MRWgLdM86Zg3/ZcNtrt7vo+B2g0kQD0jrWE2p6LZRo3TSwLIQ4PucL/EO/n87y5ssRm+A2bJk8dm9kxyQkwlkzx7L2nuL5Roe15PDSa48fXrnCpVGQwGuWh3BgBmteWFvjS9LGbAr3r+7y5sszlUhGlFGHD5LGJCY4ODn1Kn3SnlG3z+ckjvLa0gO7Vwg0ZBl+ePrpjqQYhxOEiAf8WSq0mr68sMbJt4rLpOrwyd43nTp2h6brkG3WUUuQSSVK2zf0jo1vH5vIJHhgZJWSaGL2BkIXqBoVmg5H4zoyat1dXuFwqMtpbUqHje7w6P0fcCt+Uc/9pmUmnySWTFHppmdlYXCZrhTjk5Bt8C0u1GiFl7Ji4jFlhqp0ary8tMFsub41nG0rx5emjWymLtY6DqUzC5s7uNVC8s7pGzb1O03OZTKY4lcnuCPbQzfmPhiwuFQsHFvChO6l7RBZTE+IzQwL+rWi96wRly/N4e3WVmaH0TSmLz957kiulEq8vL3F+fY0z2SxTAwOEjW43L9VqFJstjqfTZKNh1uoN5isV3CC4aZgnbJo0XPeT/pRCiD4ik7a3MJ5M4ekAf9sEast1aXkuiXD4ppRFLwj40eXLvLO+xngiSS6e4L18gV+uruIGPquNGuV2k3vSaSKhEIZSpKNRDGXQcDo0PxDca47DVErSH4UQd45c4d/CcCzG53JjvLW6srUtZBg8Nj7JlVLpptfXHYeNoMWpTLfq1IO5HNfKJS4ViiRtm+NDaU6kMzeNgyfCYcKmQbXTpuk62KEQdcdhKBJlZkgqTwkh7hwJ+B/ivpFRpgYGKDSbhAyD0XgCx/e5WiruqJDVcl00kAjbW8faZohTmRFS4QiPTkxwT3qYv3zvAl4Q7PjroOm6nMlmmR4c5EqpRM3pcCY7wvTAoEySCiHuKIkoe0jZEVL2jXXy7VCIJ7dSFrtVpSzT5KnpY7yxsnTT8cpQDEYihE2T+0dGeXN5meFYlLAZotJuYRqK4+lhEuEwj47f1SUDhBCHnAT8j2EmnWYsmaTQbGIoRSYWI2yarDZqLFarZKLdalmFZpNMLLa1zs592RFilsX59TWKvSUKHhzNkbjNJY2FEGI/JOB/TFHLYmpg58qRX5ya5v1CnkvFAr7WnMl2Sx1uLmymlGJmKC1j80KIAyEB/w6yTJMHRnM8MLpbgTAhhDhYkpYphBB9QgK+EEL0CQn4QgjRJyTgCyFEn5CAL4QQfUICvhBC9AkJ+EII0Sck4AshRJ849DdeOb7PYnWDjU6boUiU8WRKyvAJIcQuDnXAbzgOL1+9Qt3pYJkmju8zYEd4euY4Mcs66OYJIcRd5VAP6byztkrLc8klkgxHY4wlktScDu/m1w+6aUIIcdc5tAFfa83VjTLD0diO7cPRGFfLNxcoEUKIfndoA75SCkuZO0oQAvhBsFWYRAghxA37CvhKqf+glPqlUuptpdRLSqnx3nallPpPSqnZ3v5H7kxzdzqVzZJvNtFaA92r/kKryelM9pN4OyGEONT2e4X/Pa31g1rrh4H/C/y73vavAid6P98E/us+32dXpzNZjqfTrDbqWz8nhzOcSA9/Em8nhBCH2r6ydLTW1W1P44DuPX4O+O+6e+n9mlJqUCk1prVeuekk+xAyDL4wdYT7R0Zpug5xK0zStvc+UAgh+tC+0zKVUt8FfhvYAJ7qbZ4AFra9bLG37Y4G/E0p2yYlgV4IIT7UnkM6SqmXlVLnd/l5DkBr/bzWegr4c+D3Nw/b5VR6l20opb6plHpDKfVGPp//uJ9DCCHEHva8wtdaf+U2z/UXwF8B/57uFf3Utn2TwPItzv8C8ALA2bNnd/2lIIQQYv/2m6VzYtvTXwfe7z1+EfjtXrbOk8DGnR6/F0II8dHsdwz/j5RSJ4EAmAe+1dv+18A/B2aBJvA7+3wfIYQQ+7TfLJ1/dYvtGvj2fs4thBDizlKbNy3dDZRSebp/KfSDDFA46Ebc5aSPPpz0z976pY+mtdZ73nF6VwX8fqKUekNrffag23E3kz76cNI/e5M+2unQrqUjhBDio5GAL4QQfUIC/sF54aAbcAhIH3046Z+9SR9tI2P4QgjRJ+QKXwgh+oQE/E+ZUup7Sqn3e3UC/o9SanDbvu/0aghcVEr96kG286AopX5LKXVBKRUopc5+YF/f988mpdQzvX6YVUr94UG356Appf5UKbWulDq/bVtaKXVOKXW59+/QQbbxbiAB/9N3Drhfa/0gcAn4DoBS6gzwdeA+4BngT5RS/Vi66zzwG8Ar2zdK/9zQ+9z/hW7diTPAv+71Tz/7b3T/X2z3h8CPtdYngB/3nvc1CfifMq31S1prr/f0NboLy0G3hsD3tdYdrfU1ustSPH4QbTxIWuv3tNYXd9kl/XPD48Cs1vqq1toBvk+3f/qW1voV4IPFrJ8D/qz3+M+Af/mpNuouJAH/YP0u8KPe41vVEBBd0j83SF/cntHNRRt7/44ccHsO3L4LoIibKaVeBnK77Hpea/3D3mueBzy6dQTgI9QQOOxup392O2yXbZ/J/rkN0hfiY5GA/wnYq4aAUuobwLPA0/pGXuxt1xA47D5CjYXt+qZ/boP0xe1Z2yytqpQaA9YPukEHTYZ0PmVKqWeAPwB+XWvd3LbrReDrSilbKXWMbgH4nx9EG+9S0j83vA6cUEodU0qF6U5mv3jAbbobvQh8o/f4G8Ct/nrsG3KF/+n7Y8AGzimlAF7TWn9La31BKfUD4F26Qz3f1lr7B9jOA6GU+hrwn4Es8FdKqbe11r8q/XOD1tpTSv0+8DeACfyp1vrCATfrQCml/gfwK0BGKbVIt/LeHwE/UEr9HnAd+K2Da+HdQe60FUKIPiFDOkII0Sck4AshRJ+QgC+EEH1CAr4QQvQJCfhCCNEnJOALIUSfkIAvhBB9QgK+EEL0if8PtkIKg7IR0L0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.scatter(X_test[:, 0], X_test[:, 1], c=Y_pred, alpha=0.3)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 注意\n",
    "在GMM模型中存在一个奇异矩阵的问题。回顾多维高斯分布的情况，我们在计算协方差矩阵时的公式为：\n",
    "$$\n",
    "\\Sigma_{k}=\\Sigma_{i}P(c_{k}|x_{i})(x_{i}-\\mu_{k})^{T}(x_{i}-\\mu_{k})\n",
    "$$\n",
    "奇异矩阵问题只有在用户设置的参数$K_{user}$大于数据实际的类数$K$时才有可能发生。在这种情况下，随着参数的不断迭代优化，某一个多余的component可能只包含了单个样本，此时有$x_{i}=\\mu_{k}$，那么由上式得到的协方差矩阵为：\n",
    "$$\n",
    "\\Sigma_{k}=0\n",
    "$$\n",
    "再来看生成多维高斯分布的式子：\n",
    "$$\n",
    "N(x_{i}|\\mu_{k},\\Sigma_{k})=\\frac{1}{(2\\pi)^{n/2}\\Sigma_{k}^{1/2}}exp(-\\frac{1}{2}(x_{i}-\\mu_{k})^{T}\\Sigma_{k}^{-1}(x_{i}-\\mu_{k}))\n",
    "$$\n",
    "如果协方差是奇异的，那么就无法生成多维高斯分布，程序会报错。为了防止出现奇异协方差矩阵，每次在对协方差矩阵操作时人为加上一个微小值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_covar=1e-06"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}