{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### 一.简介\n", "这一章开始EM算法的探索，EM即Expectation Maximization，称为期望极大算法；如果直接介绍EM算法会有些云里雾里，所以这一节会先介绍GMM（Gaussian mixture model,高斯混合模型），然后由此引出EM算法的求解，我们首先看看单个的高斯模型求解" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 二.单个高斯分布\n", "\n", "假如，我们有如下的一堆一维数据点，并且知道它iid采样于某一一维高斯分布" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.1, 1)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEYxJREFUeJzt3X+sZGddx/H3p3e3UAu21F4UdhdayQKWH9pyUyo1WgOk2ypdUJRdRZEQG5X6I5omRUnRotHQP1SSIjZI+BFpqfyoG7K4EsVggJbeUii0dXVZfuxlTbpSqD+otlu+/jGzZfbu3J0ze+fO3D68X8nNnvOcZ57z3SfPfnrmnJnbVBWSpLacNOsCJEmTZ7hLUoMMd0lqkOEuSQ0y3CWpQYa7JDVoZLgneXuSe5N8foXjSfLmJPuS3JnkvMmXKUkaR5cr93cA245z/BJga//ncuAvVl+WJGk1RoZ7VX0MuO84XbYD76qeW4DTkzxpUgVKksa3YQJjbAIODOwv9dv+fXnHJJfTu7rn1FNPfd4zn/nMCZxekr5z3H777f9RVfOj+k0i3DOkbejvNKiq64HrARYWFmpxcXECp5ek7xxJvtyl3yQ+LbMEbBnY3wwcnMC4kqQTNIlw3wX8Yv9TMxcA91fVMbdkJEnTM/K2TJIbgIuAM5MsAW8ANgJU1VuB3cClwD7gm8Cr16pYSVI3I8O9qnaOOF7AaydWkSRp1fyGqiQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDeoU7km2JdmbZF+Sq4Ycf0qSjya5I8mdSS6dfKmSpK5GhnuSOeA64BLgHGBnknOWdXs9cFNVnQvsAN4y6UIlSd11uXI/H9hXVfur6kHgRmD7sj4FfHd/+zTg4ORKlCSNq0u4bwIODOwv9dsG/T7wyiRLwG7g14cNlOTyJItJFg8dOnQC5UqSuugS7hnSVsv2dwLvqKrNwKXAu5McM3ZVXV9VC1W1MD8/P361kqROuoT7ErBlYH8zx952eQ1wE0BVfRJ4LHDmJAqUJI2vS7jfBmxNcnaSk+k9MN21rM9XgBcCJPkBeuHufRdJmpGR4V5Vh4ErgD3APfQ+FXNXkmuSXNbv9jvALyf5LHAD8EtVtfzWjSRpSjZ06VRVu+k9KB1su3pg+27gwsmWJkk6UX5DVZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDWoU7gn2ZZkb5J9Sa5aoc/PJrk7yV1J3jPZMiVJ49gwqkOSOeA64MXAEnBbkl1VdfdAn63A64ALq+rrSZ64VgVLkkbrcuV+PrCvqvZX1YPAjcD2ZX1+Gbiuqr4OUFX3TrZMSdI4uoT7JuDAwP5Sv23Q04GnJ/l4kluSbBs2UJLLkywmWTx06NCJVSxJGqlLuGdIWy3b3wBsBS4CdgJvS3L6MS+qur6qFqpqYX5+ftxaJUkddQn3JWDLwP5m4OCQPn9bVQ9V1ReBvfTCXpI0A13C/TZga5Kzk5wM7AB2LetzM/DjAEnOpHebZv8kC5UkdTcy3KvqMHAFsAe4B7ipqu5Kck2Sy/rd9gBfS3I38FHgyqr62loVLUk6vlQtv30+HQsLC7W4uDiTc0vSo1WS26tqYVQ/v6EqSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDOoV7km1J9ibZl+Sq4/R7eZJKsjC5EiVJ4xoZ7knmgOuAS4BzgJ1JzhnS7/HAbwC3TrpISdJ4uly5nw/sq6r9VfUgcCOwfUi/NwJvAv53gvVJkk5Al3DfBBwY2F/qtz0iybnAlqr60PEGSnJ5ksUki4cOHRq7WElSN13CPUPa6pGDyUnAnwK/M2qgqrq+qhaqamF+fr57lZKksXQJ9yVgy8D+ZuDgwP7jgWcD/5TkS8AFwC4fqkrS7HQJ99uArUnOTnIysAPYdeRgVd1fVWdW1VlVdRZwC3BZVS2uScWSpJFGhntVHQauAPYA9wA3VdVdSa5JctlaFyhJGt+GLp2qajewe1nb1Sv0vWj1ZUmSVsNvqEpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGdQr3JNuS7E2yL8lVQ47/dpK7k9yZ5B+SPHXypUqSuhoZ7knmgOuAS4BzgJ1JzlnW7Q5goaqeC7wPeNOkC5Ukddflyv18YF9V7a+qB4Ebge2DHarqo1X1zf7uLcDmyZYpSRpHl3DfBBwY2F/qt63kNcCHhx1IcnmSxSSLhw4d6l6lJGksXcI9Q9pqaMfklcACcO2w41V1fVUtVNXC/Px89yolSWPZ0KHPErBlYH8zcHB5pyQvAn4P+LGq+r/JlCdJOhFdrtxvA7YmOTvJycAOYNdghyTnAn8JXFZV906+TEnSOEaGe1UdBq4A9gD3ADdV1V1JrklyWb/btcDjgL9J8pkku1YYTpI0BV1uy1BVu4Hdy9quHth+0YTrkiStgt9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQRu6dEqyDfhzYA54W1X9ybLjjwHeBTwP+Brwiqr60mRL/bab7/gq1+7Zy8FvPMCTTz+FKy9+Bi89d9O6G3M91DqsLzD09a+/+XPccOsBHq5iLmHn87ew8NQzhvbtOu7il+87Zsw/fOlzxqprWN9h465VravV6tqa9bjWenypquN3SOaAfwVeDCwBtwE7q+rugT6/Bjy3qn4lyQ7gZVX1iuONu7CwUIuLi2MXfPMdX+V1H/gcDzz08CNtp2yc449/6jknPFlrMeZ6qHVY340nBQIPPVxHvf68p5zGx79w3zHnmzspPPyto/v+9PM28f7bvzpy3JOAbw35O1z4tDP49Ffu71TXsHOtNO5a1PrKC56yqoBvdW3Netzv5FqT3F5VC6P6dbktcz6wr6r2V9WDwI3A9mV9tgPv7G+/D3hhkoxTcFfX7tl71CQBPPDQw1y7Z++6GnOtxh1nzGF9H/pWHRVqR14/LNiBo8LySN8bbj3QadxhYQnw8S/c17muYedaady1qPWGWw+scKSbVtfWrMe11tG6hPsmYHCFL/XbhvapqsPA/cD3LB8oyeVJFpMsHjp06IQKPviNB8Zqn9WYazXuOGOutv6VPDzi3d56OtesX9/q2pr1uNY6WpdwH3YFvnzFd+lDVV1fVQtVtTA/P9+lvmM8+fRTxmqf1ZhrNe44Y662/pXMrc2bsjU516xf3+ramvW41jpal3BfArYM7G8GDq7UJ8kG4DRg+Pv8Vbry4mdwysa5o9pO2Tj3yEOy9TLmWo07zpjD+m48KWycOzqwTtk4x4VPO2Po+eZOOrbvzudv6TTuSovrwqed0bmuYedaady1qHXn87escKSbVtfWrMe11tG6fFrmNmBrkrOBrwI7gJ9b1mcX8Crgk8DLgX+sUU9qT9CRBxCTfPK8FmOuh1pX6rvS68f5tMyw9mHjTuLTMsPONc6nZVZb62q0urZmPa61jjby0zIASS4F/ozeRyHfXlV/lOQaYLGqdiV5LPBu4Fx6V+w7qmr/8cY80U/LSNJ3sq6flun0Ofeq2g3sXtZ29cD2/wI/M26RkqS14TdUJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDOv0+9zU5cXII+PIqhzkT+I8JlDNJ67EmsK5xrMeawLrGtR7rmkRNT62qkf+f0pmF+yQkWezyS+unaT3WBNY1jvVYE1jXuNZjXdOsydsyktQgw12SGvRoD/frZ13AEOuxJrCucazHmsC6xrUe65paTY/qe+6SpOEe7VfukqQhDHdJatC6DPckb09yb5LPr3A8Sd6cZF+SO5OcN3DsVUn+rf/zqinW9PP9Wu5M8okkPzhw7EtJPpfkM0kWJ1VTx7ouSnJ//9yfSXL1wLFtSfb25/GqKdd15UBNn0/ycJIz+sfWZL6SbEny0ST3JLkryW8O6TOLtdWlrqmvr451TXV9daxpFmvrsUk+leSz/br+YEifxyR5b38+bk1y1sCx1/Xb9ya5eCJFVdW6+wF+FDgP+PwKxy8FPgwEuAC4td9+BrC//+cT+ttPmFJNLzhyLuCSIzX1978EnDmjuboI+NCQ9jngC8D3AycDnwXOmVZdy/q+BPjHtZ4v4EnAef3txwP/uvzvPKO11aWuqa+vjnVNdX11qWlGayvA4/rbG4FbgQuW9fk14K397R3Ae/vb5/Tn5zHA2f15m1ttTevyyr2qPgbcd5wu24F3Vc8twOlJngRcDHykqu6rqq8DHwG2TaOmqvpE/5wAtwCbJ3He1dZ1HOcD+6pqf1U9CNxIb15nUddO4IZJnXslVfXvVfXp/vZ/AfcAm5Z1m8XaGlnXLNZXx/layZqsrxOoaVprq6rqv/u7G/s/yz+tsh14Z3/7fcALk6TffmNV/V9VfRHYR2/+VmVdhnsHm4ADA/tL/baV2qftNfSu/o4o4O+T3J7k8hnU88P9t4sfTvKsftu6mKsk30UvJN8/0Lzm89V/S3wuvSusQTNdW8epa9DU19eIumayvkbN1bTXVpK5JJ8B7qV3IbDi2qqqw8D9wPewRnO1YbUDzEiGtNVx2qcmyY/T+8f3IwPNF1bVwSRPBD6S5F/6V7bT8Gl6v4viv5NcCtwMbGUdzFXfS4CPV9XgVf6azleSx9H7B/9bVfWfyw8PeclU1taIuo70mfr6GlHXTNZXl7liymurqh4GfijJ6cAHkzy7qgafOU11bT1ar9yXgC0D+5uBg8dpn4okzwXeBmyvqq8daa+qg/0/7wU+yATecnVVVf955O1iVe0GNiY5kxnP1YAdLHvbvJbzlWQjvVD466r6wJAuM1lbHeqayfoaVdcs1leXueqb6toaOMc3gH/i2Nt2j8xJkg3AafRuXa7NXE3qgcKkf4CzWPkh4U9w9EOvT/XbzwC+SO+B1xP622dMqaan0LtX9oJl7acCjx/Y/gSwbYpz9X18+8tq5wNf6c/bBnoPBc/m2w+8njWtuvrHjyzuU6cxX/2/97uAPztOn6mvrY51TX19daxrquurS00zWlvzwOn97VOAfwZ+clmf13L0A9Wb+tvP4ugHqvuZwAPVdXlbJskN9J7Cn5lkCXgDvQcUVNVbgd30PtWwD/gm8Or+sfuSvBG4rT/UNXX0W7K1rOlqevfP3tJ7RsLh6v32t++l9xYNegv+PVX1d5OoqWNdLwd+Nclh4AFgR/VW1OEkVwB76H2y4e1VddcU6wJ4GfD3VfU/Ay9dy/m6EPgF4HP9e6MAv0svOGe2tjrWNYv11aWuaa+vLjXB9NfWk4B3Jpmjd0fkpqr6UJJrgMWq2gX8FfDuJPvo/YdnR7/mu5LcBNwNHAZeW71bPKvirx+QpAY9Wu+5S5KOw3CXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDfp/oxsCIhpLT8cAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data=np.sort(np.linspace(1,3,15).tolist()+np.linspace(1.5,2,15).tolist())\n", "plt.scatter(data,[0]*len(data))\n", "plt.ylim(-0.1,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "我们知道一维高斯分布的概率密度函数如下： \n", "\n", "$$\n", "N(x\\mid u,\\sigma)=\\frac{1}{\\sqrt{2\\pi}\\cdot\\sigma}exp(-\\frac{(x-u)^2}{2\\sigma^2})\n", "$$\n", "\n", "那么该如何求该高斯分布参数呢？即估计高斯参数$u,\\sigma$，极大似然估计便是一种常用的方法\n", "#### 极大似然估计\n", "\n", "极大似然估计是一种点估计，它目标函数即似然函数，如下： \n", "\n", "$$\n", "L(u,\\sigma)=\\prod_{i=1}^MN(x_i\\mid u,\\sigma)\n", "$$\n", "\n", "$M$表示样本量，我们要求解： \n", "\n", "$$\n", "u^*,\\sigma^*=arg\\max_{u,\\sigma}L(u,\\sigma)\n", "$$ \n", "\n", "可以简单的理解为让样本出现的概率尽可能的大，由于累乘不好计算，通常会对似然函数取对数，将乘号变为加号，这时的目标函数表示如下： \n", "\n", "$$\n", "L(u,\\sigma)=log(\\prod_{i=1}^Mp(x_i\\mid u,\\sigma))=-Mlog(\\sqrt{2\\pi}\\cdot\\sigma)-\\sum_{i=1}^M\\frac{(x-u)^2}{2\\sigma^2}\n", "$$ \n", "\n", "显然取对数后的最优解和取对数前一致，由于极值点必然是导数为0的点，所以，我们分别让$L(u,\\sigma)$对$u,\\sigma$求偏导，并令其为0，即可得到$u^*,\\sigma^*$： \n", "\n", "$$\n", "\\frac{\\partial L(u,\\sigma)}{\\partial u}=\\frac{1}{\\sigma^2}\\sum_{i=1}^M(x_i-u)=0\\Rightarrow u^*=\\frac{\\sum_{i=1}^Mx_i}{M}\n", "$$\n", "\n", "$$\n", "\\frac{\\partial L(u,\\sigma)}{\\partial \\sigma}=-M\\frac{1}{\\sigma}+\\sum_{i=1}^M(x-u)^2\\sigma^{-3}=0\\Rightarrow\\sigma^*=\\sqrt{\\frac{\\sum_{i=1}^M(x-u^*)^2}{M}}\n", "$$\n", "\n", "结果非常简单且好理解，对于上面的伪数据，我们可以很方便的求解出其高斯分布的参数" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "#定义一维高斯概率分布函数\n", "def gaussian_1d(x,u,sigma):\n", " return 1/(np.sqrt(2*np.pi)*sigma)*np.exp(-1*np.power(x-u,2)/(2*sigma**2))\n", "u=np.mean(data)\n", "sigma=np.std(data)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.1, 1)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FPX9x/HXJ5uTEAiBcIVw32eAcKvFqyBVYhEhIKDoTxRvrVqxv6qFemK9fmIVEamIHAIqWhSt9agoR8JNQiDcIWACARKOkOv7+2MXG0MgG7K7s8fn+XjwMDs7O/tmHsPbycx8Z8QYg1JKKf8SZHUApZRSrqflrpRSfkjLXSml/JCWu1JK+SEtd6WU8kNa7kop5YeqLHcRmS0iOSKy5Tzvi4i8JiKZIrJJRHq5PqZSSqnqcGbPfQ4w9ALvXwO0c/yZBPy95rGUUkrVRJXlboz5Hsi7wCxJwHvGbhUQLSJNXBVQKaVU9QW7YBlxwP5yr7Mc0w5WnFFEJmHfuycyMrJ3x44dXfD1SikVOFJTUw8bY2Krms8V5S6VTKv0ngbGmJnATIDExESTkpLigq9XSqnAISJ7nZnPFVfLZAHx5V43A7JdsFyllFIXyRXlvgyY4Lhqpj9w3BhzziEZpZRSnlPlYRkRmQ8MBhqISBbwJBACYIx5E1gODAMygVPARHeFVUop5Zwqy90YM6aK9w1wt8sSKaWUqjEdoaqUUn5Iy10ppfyQlrtSSvkhLXellPJDWu5KKeWHtNyVUsoPabkrpZQf0nJXSik/pOWulFJ+SMtdKaX8kJa7Ukr5IS13pZTyQ1ruSinlh7TclVLKD2m5K6WUH9JyV0opP6TlrpRSfkjLXSml/JCWu1JK+SEtd6WU8kNa7kop5Ye03JVSyg8FWx1Aqeo6fqqY9EP5bDuYz/HTJee8HxYSRPtGtenUpA6N64QjIhakVMpaWu7K66UfzOfzLYfYeuA46QfzyT5e6PRno2uF0KlxHTo1qcPlHWMZ2KYBtiAte+X/tNyVVzp0vJBPNhzgo/UH2HaoAFuQ0CY2kj6tYujYuA6dmkTRuUkdGtQOO+ezBWdKyDhUwLZD+aQfzCftYAEfrNnL7JW7aVQnjKSEOH7fM45OTepY8DdTyjPEGGPJFycmJpqUlBRLvlt5J2MM/0rP4R8/7mHlzsMYAwnx0YzoFce13ZsSExl60csuLC7l6/QcPlqfxbcZuZSUGTo2jmJM3+Yk940nLNjmwr+JUu4jIqnGmMQq59NyV95gze48nv9iG6l7j9KsXgQjejXj+oSmtI6t7fLvyjtZxGebslmy7gAb9x8jLjqCh65uz/U94/SQjfJ6Wu7KJ2w7lM8LX2Tw7205NIwK44Gr2jMqsRnBNs9cyPWfHbk8/8U2thzIp2PjKB4Z0oErOjbUk7DKa2m5K692/FQx0/6ZxpJ1WdQOC2by4DZMHNiKiFDPHx4pKzP8c/NBXvwyg71HTtG3ZQwvJycQFx3h8SxKVUXLXXmtlZmHefjDjeQWnOHWS1px1+A2RNe6+OPprlJcWsaCtft54fNt1A4PZu5tfWnbMMrqWEr9irPlroOYlMcUFpcy7bM0bpq1mohQG0vvGsjjwzp5RbEDhNiCGN+/BQvvGEBxqWHkmz+xbt9Rq2MpdVGcKncRGSoiGSKSKSKPVfJ+cxH5RkTWi8gmERnm+qjKl6Vl55P0+kre+WE34/u34J/3Xkr3ZtFWx6pU56Z1WDp5IHUjQrjp7dV8m5FjdSSlqq3KchcRGzADuAboDIwRkc4VZvtfYJExpieQDLzh6qDKdy1cu4/rZ6wk71QR707sw7Tru1pybL06mtevxeI7B9KqQST/848UPtlwwOpISlWLM3vufYFMY8wuY0wRsABIqjCPAc6OCKkLZLsuovJVZWWGZz9P549LNtOvdQxf3H8pl3doaHUsp8VGhbHgjv4ktqzH/Qs28O7K3VZHUsppzpR7HLC/3Ossx7TyngLGiUgWsBy4t7IFicgkEUkRkZTc3NyLiKt8xemiUu6at463vtvFuP7NefeWPtSvZDSpt6sTHsKciX0Z0qURf/k0jekrtmHVRQhKVYcz5V7ZBb8Vt+4xwBxjTDNgGDBXRM5ZtjFmpjEm0RiTGBsbW/20yifk5BcyeuZPrEg7xBPXdmZaUlePXbfuDuEhNt64qTdj+sYz45udTFm6mZLSMqtjKXVBztxbJguIL/e6GecedrkNGApgjPlJRMKBBoCeiQow6QfzuXXOWo6fLubt8Ylc1bmR1ZFcwhYkPPP7btSPDOP1bzI5eqqIV5N7Eh7i3ecOVOByZndqLdBORFqJSCj2E6bLKsyzD7gSQEQ6AeGAHncJMBv3H2P0Wz9hDHx45wC/KfazRISHh3Tgyes6s2Lrz9w8ew35hcVWx1KqUlWWuzGmBLgHWAGkY78qZquITBWR4Y7Z/gDcLiIbgfnALUYPTAaU9fuOMu6d1dStFcLiyQPo0rSu1ZHcZuKgVryanEDq3qOMfmsVOQXO34JYKU/REaqqxlL3HuWW2WuoFxnK/En9A2bY/rcZOUx+fx2xUWHMva0vLepHWh1JBQAdoao8InVvHjfPXkNM7VAW3hE4xQ4wuENDPri9H/mFxdzw95/Ymn3c6khK/ULLXV20tXvymPDOGmKjwlg4aQBN6gZOsZ/Vs3k9Ft85gBCbkPzWKlbtOmJ1JKUALXd1kTbsP8bNs9fQqG44Cyb1p3HdcKsjWaZtwyiWTB5IwzphTJi9hhVbD1kdSSktd1V9uw+f5NY5a6lfO5QFt/enUZ3ALfazmkZHsPjOgXRuUofJ76eycO0+qyOpAKflrqolp6CQCbNXA/Derf1oqMX+i3qRoXxwez8uaRfLH5dsZsY3mTqaVVlGy105raCwmInvruVwQRHv3tKHVg306pCKaoUGM2tCIkkJTZm+IoNpn6VTVqYFrzzPmRGqSnGmpJQ7308l41ABs25OpEe8d96u1xuEBgfx8qgE6tUKZfbK3eSdPMP0G3sQ4sO3YFC+R8tdVamszPDwh5tYmXmEl0b1YLAP3dnRKkFBwpPXdSY2KozpKzI4drqYN27qRa1Q/SenPEN3JVSVnl+xjU83ZvPYNR0Z0auZ1XF8hohw9+VteW5EN77fnsttc1L0hmPKY7Tc1QUtXZf1y21777istdVxfFJy3+Y8f0N3ftp1hJf/td3qOCpAaLmr81q/7yiPLd3MgNb1efK6LohUdvdn5YwbE+MZnWi/ZbA+tk95gpa7qtSh44XcMTeVRnXCeOOmXnoy0AX+ktSFjo2jeHDhBg4eP211HOXn9F+sOkdhcSmT5qZw8kwJsyb0oV5kqNWR/EJ4iI0ZN/WiqKSMez9YT7Eef1dupOWufsUYw6OLN7H5wHFeSe5Jh8ZRVkfyK21ia/PMiG6k7D3Ki19mWB1H+TEtd/Urf/9uJ8s2ZvPwbztwtZ89bMNbJCXEMbZfc976bhdfp/9sdRzlp7Tc1S/+syOX6SsyuK5HU+4a3MbqOH7tiWs707lJHR5atJGso6esjqP8kJa7AiD72GnuX7CBdg1r8/wN3fTKGDezP3S7F6Vlhns+WE9RiR5/V66l5a4oKinjng/Wcaa4lL+P662jKD2kZYNIXhjZnQ37j/H8F9usjqP8jJa74tnP01m37xgvjOxBm9jaVscJKMO6NeHmAS1454fdeh945VJa7gHus03ZvLtyDxMHteR33ZtYHScgPf67TnRvVpeHP9zI/jw9/q5cQ8s9gO3MPcEfF2+iV/NoplzTyeo4ASss2MaMsb0AuPuDdZwpKbU4kfIHWu4B6lRRCZPfTyXMMbAmNFg3BSvFx9Ri+sgebMo6zrPL9fi7qjn9Fx2gnlq2lR05J3g1OSEgH2ztjYZ2bcytg1ox58c9LN980Oo4ysdpuQegZRuzWZSSxV2D23Bpu1ir46hyHrumIwnx0Ty6eBN7Dp+0Oo7yYVruAWbfkVM8vnQzvZpH88BV7a2OoyoIDQ7i9bE9sQUJd81bR2GxHn9XF0fLPYAUl5Zx74L1iMCryT31To9eqlm9Wrw0qgdpB/OZ9lma1XGUj9J/3QHkxS8z2Lj/GM/f0J34mFpWx1EXcGWnRtxxWWvmrd7HJxsOWB1H+SAt9wDx/fZc3vpuF2P7NWdYN72e3Rc8PKQDvVvU4/Glm9mZe8LqOMrHaLkHgJyCQh5atIH2jWrzxLWdrY6jnBRisx9/Dw0O4m49/q6qScvdzxljePjDTRQUlvD62F6Eh9isjqSqoUndCF4ancC2QwU8tWyr1XGUD3Gq3EVkqIhkiEimiDx2nnlGiUiaiGwVkQ9cG1NdrPd+2sv323P53991on0jffCGL7q8Q0PuvrwNC9bu56P1WVbHUT6iytv/iYgNmAFcDWQBa0VkmTEmrdw87YApwCBjzFERaeiuwMp5mTkFPLM8ncEdYhnXv4XVcVQNPHhVe1L2HOXxpVvo2rQu7fR/1KoKzuy59wUyjTG7jDFFwAIgqcI8twMzjDFHAYwx+nh3ixWVlPHAwg3UCrXxwg3d9f7sPi7YFsRrY3pSK9TGXfPWcaqoxOpIyss5U+5xwP5yr7Mc08prD7QXkZUiskpEhla2IBGZJCIpIpKSm5t7cYmVU179ejtbDuTz7IjuNKwTbnUc5QKN6oTzanJPMnNP8OeP9fi7ujBnyr2yXT5T4XUw0A4YDIwBZolI9DkfMmamMSbRGJMYG6vD3t0lZU8ef/92Jzf2bsbQro2tjqNc6JJ2Dbj3inYsWZfFopT9VX9ABSxnyj0LiC/3uhmQXck8nxhjio0xu4EM7GWvPKygsJgHF20grl4ETw7vYnUc5Qb3X9mOgW3q88QnW8g4VGB1HOWlnCn3tUA7EWklIqFAMrCswjwfA5cDiEgD7IdpdrkyqHLO1E/TOHD0NC+PSqB2mD4uzx/ZgoRXkhOoHRbC5HmpnDyjx9/Vuaosd2NMCXAPsAJIBxYZY7aKyFQRGe6YbQVwRETSgG+AR4wxR9wVWlXuiy0H+TA1i8mD25DYMsbqOMqNGkaF89qYBHYfPslf/6n3n1HnEmMqHj73jMTERJOSkmLJd/ujnPxChrzyPXH1Ilg6eZA+fCNAPPf5Nt78biezJiRyVedGVsdRHiAiqcaYxKrm0wbwA8YYHl2yiVNFpbwyOkGLPYA8eHU7OjWpw2NLN3HkxBmr4ygvoi3gB95ftZdvM3J5fFgn2jbUwS2BJCzYxiujE8g/XcKUpZux6jdx5X203H3cztwTPL08ncvaxzJhgI5CDUQdGkfxyJAOfJn2M4tT9fYEyk7L3YcVl5bx4MINhIfYmD5SR6EGstsuaUW/VjH85dM09uedsjqO8gJa7j7sta93sCnrOM/+vhuNdBRqQAsKEv42qgcAf1i0kdIyPTwT6LTcfVTq3jxmfJPJDb2acY0+fENhfzzfU8O7sGZPHrP+o8NMAp2Wuw86caaEBxdupGl0BE8N14dvqP+6oVccQ7s05m9fbif9YL7VcZSFtNx90NRPt5J19BQvj04gKjzE6jjKi4gIz4zoRp2IEB5cuIEzJfr0pkCl5e5jVmw9xKKULO78TRv66ChUVYmYyFBeGNmNbYcKeOmr7VbHURbRcvchOQWFTFm6ma5xdXjgqvZWx1Fe7IqOjRjbrzkzv9/F6l16J5BApOXuI4wxPLp4EyfPlOgoVOWUPw3rRPOYWvzhw40UFBZbHUd5mDaEj9BRqKq6IsOCeWlUAtnHTjP1U725WKDRcvcBOgpVXazeLepx1+C2fJiaxYqth6yOozxIy93L6ShUVVP3XdmOrnF1mLJ0M7kFenOxQKHl7uXOjkJ9boSOQlUXJzQ4iJdHJXDiTAlTlm7Sm4sFCC13L3Z2FOrI3s0Y2lVHoaqL165RFI8N7ci/0nNYuFafvRoItNy9VPlRqE9ep6NQVc3dMrAlg9rWZ+pnaew9ctLqOMrNtNy9lI5CVa4WFCRMH9kDW5DwkN5czO9puXuhL7bYR6FOHqyjUJVrNY2OYFpSV1L3HuXN73ZaHUe5kZa7l8nJL2TK0k10javD/VfqKFTlekkJTfld9ya88q/tbDlw3Oo4yk203L2IPgtVeYKI8PT1XalXK5QHF26gsFhvLuaPtD28iI5CVZ4SXSuU6Tf2YEfOCV5ckWF1HOUGWu5eIjPHPgr1NzoKVXnIb9rHMr5/C2b9sJsfdx62Oo5yMS13L1BUYh+FGqGjUJWHTRnWkdYNInl40Uby9eZifkXL3Qu89vUONh84zrMjutFQR6EqD6oVGsxLoxP4ueAMT32y1eo4yoW03C2WujePN77N5EYdhaoskhAfzT2Xt2Xp+gMs33zQ6jjKRbTcLXTiTAkPLNxAXL0Inhzexeo4KoDdc0Vbujery+MfbSYnv9DqOMoFtNwt9JdlWzlw9DQvj0qgdliw1XFUAAuxBfHy6AQKi0t5dIneXMwfaLlb5IstB/kw1T4KNVFHoSov0Ca2NlOu6cS3GbnMW73P6jiqhrTcLWAfhbqZbnF1dRSq8irj+7fg0nYN+Os/00g/mG91HFUDWu4eZozhkcWbOF1cyss6ClV5maAg4aVRCdQJD+Gueev02as+zKlmEZGhIpIhIpki8tgF5hspIkZEEl0X0b/MXbWX77afHYVa2+o4Sp0jNiqM18f2Yl/eKR5drMfffVWV5S4iNmAGcA3QGRgjIufcYFxEooD7gNWuDukvMnNO8PQ/038ZGaiUt+rbKoY/Du3A51sOMXvlHqvjqIvgzJ57XyDTGLPLGFMELACSKplvGvACoNdRVeLsKNRaoToKVfmG2y9tzW87N+LZ5emk7s2zOo6qJmfKPQ4o/1yuLMe0X4hITyDeGPPZhRYkIpNEJEVEUnJzc6sd1pf9dxRqdx2FqnyCiDD9xh7E1Yvg7nnrOXJCH67tS5wp98p2MX85CCciQcDLwB+qWpAxZqYxJtEYkxgbG+t8Sh+Xsqf8KNTGVsdRyml1I0J446Ze5J0q4v4FG/TpTT7EmXLPAuLLvW4GZJd7HQV0Bb4VkT1Af2CZnlS1O3aqiAcX6ShU5bu6NK3LtKQu/JB5mFe/3mF1HOUkZ8p9LdBORFqJSCiQDCw7+6Yx5rgxpoExpqUxpiWwChhujElxS2IfUlpmuHf+eg4dL+TV5J46ClX5rNF9mnNj72b837938G1GjtVxlBOqLHdjTAlwD7ACSAcWGWO2ishUERnu7oC+7IUV2/jPjsNMS+pKr+b1rI6jVI1MTepKh0ZRPLBwAweOnbY6jqqCU9e5G2OWG2PaG2PaGGOedkx7whizrJJ5B+teOyzbmM1b3+1iXP/mJPdtbnUcpWosItTG38f1prTUcNe8dRSVlFkdSV2ADo90g63Zx3l08Ub6tKzHE9fqcXblP1o1iGT6jd3ZuP8YzyxPtzqOugAtdxfLO1nEpPdSiY4I5Y2beuvtBZTfGdq1Cf9zSSvm/LiHTzdmV/0BZQltHhcqKS3j7nnryD1xhrfG9yY2KszqSEq5xR+v6UjvFvV4bMkmMnNOWB1HVULL3YWeWb6Nn3Yd4dnfd6NHfLTVcZRymxBbEDPG9iI8xMbk91M5VVRidSRVgZa7iyxJzWL2yt1MHNSSG3o3szqOUm7XuG44r43pSWbuCR5fullvMOZltNxdYOP+Y0z5aDMD29TnT8M6WR1HKY8Z1LYBD13Vno83ZPPBGn3AhzfRcq+h3IIz3DE3ldja9tukBtt0larAcvflbRncIZa/LEtjU9Yxq+MoB22iGigqKeOueakcO13EzAm9iYkMtTqSUh4XFCS8PCqB2KgwJr+/jmOniqyOpNByr5Gpn21l7Z6jvDCyB12a1rU6jlKWqRcZyoybepFTUMgfFm2kTG8wZjkt94s0f80+3l+1jzt+05rhPZpaHUcpyyXER/Pnazvz9bYc3vx+p9VxAp6W+0VI3ZvHE59s4dJ2DXh0SEer4yjlNcb3b8F1PZry4ooMftx52Oo4AU3LvZp+zi/kzvfX0TQ6gv8b0xNbkD5RSamzRITnRnSjVYNI7pu/np/z9cFsVtFyr4bC4lLumJvKyTMlzByfSHQtPYGqVEWRYcG8Oa43J8+Ucu8H6ykp1RuMWUHL3UnGGJ74ZAsb9h/jpVE96NA4yupISnmtdo2ieO6GbqzZk8f0LzOsjhOQtNydNHfVXhalZHHfFW0Z2rWJ1XGU8npJCXGM69+ct77bxZdbD1kdJ+BouTth1a4jTP00jSs7NuSBq9pbHUcpn/HnazvTvVld/vDhRvYeOWl1nICi5V6FA8dOc/e8dTSvX4uXkxMI0hOoSjktLNjGjLG9CBJh8vvrKCwutTpSwNByvwD7CdQUikrKeHtCInXCQ6yOpJTPiY+pxcuje5B2MJ97PlivT3DyEC338zDGMGXpZrZm5/NKcgJtYmtbHUkpn3VFx0ZMS+rCv9J/5r756ynWK2jcTsv9PN75YTcfrT/AQ1e158pOjayOo5TPGz+gJU9e15kvth7igYUb9BJJNwu2OoA3+mHHYZ5Zns7QLo25+/K2VsdRym9MHNSKklLD08vTCQ4SXhqVoAMB3UTLvYJ9R05xz/x1tG1YmxdH9dATqEq52O2XtaakzPD8F9uwBQnTR/bQgncDLfdyThWVMGluCmVlhrcnJFI7TFePUu4weXAbSkrL+NtX2wkOEp4b0V13pFxM28vhVFEJd8xNZfvPBbw7sS8t6kdaHUkpv3bvle0oLjO89vUObEFBPH19Vy14F9JyB46fLubWOWtZv+8oz93Qnd+0j7U6klIB4cGr2lFSWsYb3+4kOEiYmtQFES14Vwj4cj984gwT3lnDjpwCXh/bi2Hd9NYCSnmKiPDIkA6UlBlmfr+LYJvwxLWdteBdIKDLPfvYacbNWk328dPMurmP7rErZQERYco1HSkpNcxeuZsQWxBTrumoBV9DAVvuu3JPMP6dNeSfLmbubf3o0zLG6khKBSwR4c/XdqK0rIyZ3+/CFiQ8OqSDFnwNBGS5p2XnM2H2aoyB+ZP60zVOn3+qlNVEhKeGd6G4zPD3b3cSEiQ89NsOVsfyWQFX7ql785j47loiw4KZe1s/2jbU2woo5S1EhL8mdaW01PDavzMJtgVx35XtrI7lk5y6/YCIDBWRDBHJFJHHKnn/IRFJE5FNIvK1iLRwfdSa+8+OXMbNWkNMZCgf3jlAi10pLxQUJDw7ohsjesXx0lfbmfFNptWRfFKVe+4iYgNmAFcDWcBaEVlmjEkrN9t6INEYc0pEJgMvAKPdEfhifbHlEPfNX0/r2Ejeu60vDaPCrY6klDqPIMfI1dIyw/QVGYTYhEmXtbE6lk9x5rBMXyDTGLMLQEQWAEnAL+VujPmm3PyrgHGuDFlTS1KzeHTJJro3q8ucW/pSt5beulcpb2cLEv52o73gn1m+DVtQELdd0srqWD7DmXKPA/aXe50F9LvA/LcBn1f2hohMAiYBNG/e3MmINTNn5W6e+jSNQW3rM3N8IpF6SwGlfEawLYiXRydQWmaY9lkaITZhwoCWVsfyCc4cc6/sWiRT6Ywi44BEYHpl7xtjZhpjEo0xibGx7r2m3BjD/329g6c+TeO3nRvxzs19tNiV8kEhtiBeTe7JVZ0a8cQnW5m3eq/VkXyCM+WeBcSXe90MyK44k4hcBfwJGG6MOeOaeBfHGMMzy9P521fbGdEzjjdu6kV4iM3KSEqpGggNDmLGTT25omND/vTRFhau3Wd1JK/nTLmvBdqJSCsRCQWSgWXlZxCRnsBb2Is9x/UxnVdaZn+C0tv/2c3NA1rw4o09CLbpM0mU8nVhwTbeuKkXl7WP5bGlm1mcmmV1JK9WZesZY0qAe4AVQDqwyBizVUSmishwx2zTgdrAhyKyQUSWnWdxblVUUsZ9C9azYO1+7r2iLU8N76J3mVPKj4SH2Jg5vjeD2jTgkcUb+Xj9AasjeS0xptLD526XmJhoUlJSXLa800WlTJ6XyrcZuTw+rKNeNqWUHztdVMrEOWtYszuPV5J7MrxHU6sjeYyIpBpjEquazy+OVxQUFnPz7DV8tz2XZ0d002JXys9FhNp45+Y+JLaI4b7565n6aRqFxaVWx/IqPl/ueSeLGPv2atbtO8pryT0Z09czl1gqpawVGRbMP27ty4QBLZi9cjfDX/+BtOx8q2N5DZ8u90PHCxn11k9s/7mAtyckcl0A/WqmlLLvwU9N6sq7E/tw9FQx189YyVvf7aS0zJrDzd7EZ8t9z+GTjHzzRw4dL+S9W/tyeceGVkdSSlnk8g4NWfHAZVzeMZZnP9/G2LdXceDYaatjWcony33P4ZPc+NZPnDxTwvzb+9OvdX2rIymlLBYTGcqb43rzwsjubDlwnKGvfM/H6w9g1UUjVvPJcm9cN5wBreuz6I4BdGum92JXStmJCKMS4/n8/sto3yiKBxZu4N756zl+qtjqaB7nN5dCKqVUeSWlZbz53U5e+dcOYqPCePHGHgxq28DqWDUWUJdCKqVURcG2IO65oh1L7xpIRKiNm2atZtpngXPJpJa7UsqvdW8WzT/vvZTx/Vvwzg+7SXp9JekH/f+SSS13pZTfiwi1Me36rrx7Sx+OnCwi6fWVzPx+J2V+fMmklrtSKmBc3rEhKx64lMEdYnlm+TbGzvLfSya13JVSAaV+7TDeGt+bF27ozuYs+yWTn2zwvxuQabkrpQKOiDCqTzzL77+Udg1rc/8C/7tkUstdKRWwWtSPZNEdA/jD1e35fPNBhr76PT9mHrY6lktouSulAlqwLYh7r2zHkskDiQixMXbWav7qB5dMarkrpRTQIz6az+67hHH9mzPrh90MfO7fvPDFNrJ99ISrjlBVSqkK1uzO450fdvFV2s+ICEO6NOKWga3o07IeItY+3c3ZEarBngijlFK+pG+rGPq2imF/3ineX72XBWv2s3zzITo1qcPEgS0ZntCU8BCb1TEvSPfclVKqCqeLSvlkwwHm/LiHbYcKqFcrhOS+zRnXvwVx0REezeLsnruWu1JKOckYw+rdecxZuYcv0w4BMKRLY24e2JJ+rWI8cshGD8sopZSLiQj9W9enf+v6ZB09xfur9rFg7T4+33KIjo2juGVgS5KcXE7ZAAAJhElEQVQS4ogItf6Qje65K6VUDRQW2w/ZvLvSfsgmulYIo/vEM75/C5rVq+Xy79PDMkop5UHGGNbszmPOj3v4Mu1njDFc3dl+lU3/1q47ZKOHZZRSyoNEhH6t69OvdX0OHDvN+6v2smDNPlZs/ZkOjaK4ZVBLrvfgIRvdc1dKKTcpLC5l2YZs5vy4h7SD+dSNCGH5/ZfW6Aob3XNXSimLhYfYGNUnnhsTm5Gy9yj/SvuZpnXDPfLdWu5KKeVmIkKfljH0aRnjse/Ue8sopZQf0nJXSik/pOWulFJ+SMtdKaX8kFMnVEVkKPAqYANmGWOeq/B+GPAe0Bs4Aow2xuxxbdT/+nj9AaavyCD72GmaRkfwyJAOXN8zzuuW6Q1ZK5sXqPTz//vxZuav3k+pMdhEGNMvnsQWMZXO6+xyU/bmnbPMv17frVq5Kpu3suW6K2tN+eu2ZfVyNeuFVXmdu4jYgO3A1UAWsBYYY4xJKzfPXUB3Y8ydIpIM/N4YM/pCy73Y69w/Xn+AKUs3c7rcU1IiQmw8O6LbRa8sdyzTG7JWNm9IkIBAcan51ed7Na/Lyp1553yfLUgoLfv1vDf0jmNJ6oEqlxsElFXydxjUJoZ1+447lauy7zrfct2RdVz/5jUqeH/dtqxebiBndfY6d2cOy/QFMo0xu4wxRcACIKnCPEnAPxw/LwauFDfdHm36ioxfrSSA08WlTF+R4VXLdNdyq7PMyuYtLjO/KrWzn6+s2IFfleXZeeev3u/UcisrS4CVO/OczlXZd51vue7IOn/1/vO84xx/3basXq5mrZoz5R4HlN/CsxzTKp3HGFMCHAfqV1yQiEwSkRQRScnNzb2owOd75FVNHoXljmW6a7nVWaa7Hg9W6sFRzTX9Lqs/76/bltXL1axVc6bcK9sDr7jFOzMPxpiZxphEY0xibGysM/nO0fQ8w3bPN92qZbprudVZZk3zn4/Ng48Zq+l3Wf15f922rF6uZq2aM+WeBcSXe90MyD7fPCISDNQFKv89v4YeGdKBiAqPt4oIsf1yksxblumu5VZnmZXNGxIkhNh+XVgRITYGtal85Jwt6Nx5x/SLd2q559u4BrWJcTpXZd91vuW6I+uYfvHnecc5/rptWb1czVo1Z66WWQu0E5FWwAEgGRhbYZ5lwM3AT8BI4N/GTXckO3sCwpVnnt2xTG/Ier55z/f56lwtU9n0ypbriqtlKvuu6lwtU9OsNeGv25bVy9WsVXPqrpAiMgx4BfulkLONMU+LyFQgxRizTETCgblAT+x77MnGmF0XWqbeFVIpparPpXeFNMYsB5ZXmPZEuZ8LgRurG1IppZR76AhVpZTyQ1ruSinlh7TclVLKD2m5K6WUH9JyV0opP6TlrpRSfkjLXSml/JCWu1JK+SEtd6WU8kNa7kop5Ye03JVSyg9puSullB/ScldKKT+k5a6UUn7Iqfu5u+WLRXKBvTVcTAPgsAviuJI3ZgLNVR3emAk0V3V5Yy5XZGphjKnyOaWWlbsriEiKMzet9yRvzASaqzq8MRNoruryxlyezKSHZZRSyg9puSullB/y9XKfaXWASnhjJtBc1eGNmUBzVZc35vJYJp8+5q6UUqpyvr7nrpRSqhJa7kop5Ye8stxFZLaI5IjIlvO8LyLymohkisgmEelV7r2bRWSH48/NHsx0kyPLJhH5UUR6lHtvj4hsFpENIpLiqkxO5hosIscd371BRJ4o995QEclwrMfHPJzrkXKZtohIqYjEON5zy/oSkXgR+UZE0kVkq4jcX8k8VmxbzuTy+PblZC6Pbl9OZrJi2woXkTUistGR6y+VzBMmIgsd62O1iLQs994Ux/QMERniklDGGK/7A1wG9AK2nOf9YcDngAD9gdWO6THALsd/6zl+ruehTAPPfhdwzdlMjtd7gAYWravBwGeVTLcBO4HWQCiwEejsqVwV5r0O+Le71xfQBOjl+DkK2F7x72zRtuVMLo9vX07m8uj25Uwmi7YtAWo7fg4BVgP9K8xzF/Cm4+dkYKHj586O9RMGtHKsN1tNM3nlnrsx5nsg7wKzJAHvGbtVQLSINAGGAF8ZY/KMMUeBr4ChnshkjPnR8Z0Aq4Bmrvjemua6gL5ApjFmlzGmCFiAfb1akWsMMN9V330+xpiDxph1jp8LgHQgrsJsVmxbVeayYvtycn2dj1u2r4vI5KltyxhjTjhehjj+VLxaJQn4h+PnxcCVIiKO6QuMMWeMMbuBTOzrr0a8stydEAfsL/c6yzHtfNM97Tbse39nGeBLEUkVkUkW5Bng+HXxcxHp4pjmFetKRGphL8kl5Sa7fX05fiXuiX0PqzxLt60L5CrP49tXFbks2b6qWlee3rZExCYiG4Ac7DsC5922jDElwHGgPm5aV8E1XYBFpJJp5gLTPUZELsf+j++ScpMHGWOyRaQh8JWIbHPs2XrCOuz3ojghIsOAj4F2eMG6crgOWGmMKb+X79b1JSK1sf+Df8AYk1/x7Uo+4pFtq4pcZ+fx+PZVRS5Lti9n1hUe3raMMaVAgohEAx+JSFdjTPlzTh7dtnx1zz0LiC/3uhmQfYHpHiEi3YFZQJIx5sjZ6caYbMd/c4CPcMGvXM4yxuSf/XXRGLMcCBGRBli8rspJpsKvze5cXyISgr0U5hljllYyiyXblhO5LNm+qsplxfblzLpy8Oi2Ve47jgHfcu5hu1/WiYgEA3WxH7p0z7py1QkFV/8BWnL+k4S/49cnvdY4pscAu7Gf8Krn+DnGQ5maYz9WNrDC9EggqtzPPwJDPbiuGvPfwWp9gX2O9RaM/aRgK/57wquLp3I53j+7cUd6Yn05/t7vAa9cYB6Pb1tO5vL49uVkLo9uX85ksmjbigWiHT9HAP8Brq0wz938+oTqIsfPXfj1CdVduOCEqlcelhGR+djPwjcQkSzgSewnKDDGvAksx35VQyZwCpjoeC9PRKYBax2Lmmp+/SuZOzM9gf342Rv2cySUGPvd3xph/xUN7Bv8B8aYL1yRyclcI4HJIlICnAaSjX2LKhGRe4AV2K9smG2M2erBXAC/B740xpws91F3rq9BwHhgs+PYKMDj2IvTsm3LyVxWbF/O5PL09uVMJvD8ttUE+IeI2LAfEVlkjPlMRKYCKcaYZcA7wFwRycT+P55kR+atIrIISANKgLuN/RBPjejtB5RSyg/56jF3pZRSF6DlrpRSfkjLXSml/JCWu1JK+SEtd6WU8kNa7kop5Ye03JVSyg/9P2ISQi59VZ+aAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.scatter(data,[0]*len(data))\n", "plt.plot(data,gaussian_1d(data,u,sigma))\n", "plt.ylim(-0.1,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 二.多个高斯分布\n", "可以发现结果比较make sense，数据密集的地方，对应的概率密度也很高，接着我们再造一些数据看看，比如下面，这部分数据明显可以分为两部分，中间有隔断，如果再用一个高斯分布去拟合，其结果就有些牵强了,概率分布高的地方零星只有几个点的数据" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.1, 1)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGBhJREFUeJzt3X1wJPdd5/HPVzOjp9GsVg+jtS1pn9d2fCEpxyIkuADngbKTo2zggNhXuQtUChcFAe5ImXKACilfXSVFru7g6gycK5cLcFyME1Jhi1pwDsdUqFQcLOPEiXfteL27trRP0q72Qc+rh+/90TPSaHakaY1mNNqf3q8qlfrh191ftVqf6flNd8vcXQCAsDTUuwAAQPUR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASob7mb2eTMbMbPvrzLfzOy/m9lxM3vJzN5R/TIBAOsR58z9C5LuW2P+ByQdyn09LOlPNl4WAGAjyoa7u39D0tgaTR6Q9OceeU7STjO7uVoFAgDWL1mFdfRKGioYH85NO1vc0MweVnR2r3Q6fdftt99ehc0DwPbxwgsvXHD3bLl21Qh3KzGt5DMN3P0JSU9I0sDAgA8ODlZh8wCwfZjZG3HaVeNqmWFJ/QXjfZLOVGG9AIAKVSPcD0v697mrZt4l6Yq7X9clAwDYPGW7Zczsi5LukdRtZsOSfl9SSpLc/U8lHZH0QUnHJU1J+qVaFQsAiKdsuLv7Q2Xmu6Rfq1pFAIAN4w5VAAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACFCvczew+M3vVzI6b2aMl5u82s2fN7EUze8nMPlj9UgEAcZUNdzNLSHpc0gck3SHpITO7o6jZ70l6yt3vlPSgpD+udqEAgPjinLm/U9Jxdz/h7tckPSnpgaI2LmlHbrhd0pnqlQgAWK844d4raahgfDg3rdCnJH3YzIYlHZH066VWZGYPm9mgmQ2Ojo5WUC4AII444W4lpnnR+EOSvuDufZI+KOkvzOy6dbv7E+4+4O4D2Wx2/dUCAGKJE+7DkvoLxvt0fbfLRyU9JUnu/i1JzZK6q1EgAGD94oT785IOmdk+M2tU9IHp4aI2b0p6nySZ2VsUhTv9LgBQJ2XD3d3nJX1M0tOSjim6KuZlM3vMzO7PNfu4pF82s+9K+qKkX3T34q4bAMAmScZp5O5HFH1QWjjtkwXDRyXdXd3SAACV4g5VAAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AEKBY4W5m95nZq2Z23MweXaXNL5jZUTN72cz+b3XLBACsR7JcAzNLSHpc0k9KGpb0vJkddvejBW0OSfqEpLvd/ZKZ9dSqYABAeXHO3N8p6bi7n3D3a5KelPRAUZtflvS4u1+SJHcfqW6ZAID1iBPuvZKGCsaHc9MK3SrpVjP7ppk9Z2b3lVqRmT1sZoNmNjg6OlpZxQCAsuKEu5WY5kXjSUmHJN0j6SFJnzOzndct5P6Euw+4+0A2m11vrQCAmOKE+7Ck/oLxPklnSrT5G3efc/eTkl5VFPYAgDqIE+7PSzpkZvvMrFHSg5IOF7X5qqT3SJKZdSvqpjlRzUIBAPGVDXd3n5f0MUlPSzom6Sl3f9nMHjOz+3PNnpZ00cyOSnpW0iPufrFWRQMA1mbuxd3nm2NgYMAHBwfrsm0AuFGZ2QvuPlCuHXeoAkCACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQoFjhbmb3mdmrZnbczB5do93PmZmb2UD1SgQArFfZcDezhKTHJX1A0h2SHjKzO0q0y0j6DUnfrnaRAID1iXPm/k5Jx939hLtfk/SkpAdKtPtPkv5A0kwV6wMAVCBOuPdKGioYH85NW2Jmd0rqd/e/XWtFZvawmQ2a2eDo6Oi6iwUAxBMn3K3ENF+aadYg6b9J+ni5Fbn7E+4+4O4D2Ww2fpUAgHWJE+7DkvoLxvsknSkYz0h6q6R/NLNTkt4l6TAfqgJA/cQJ9+clHTKzfWbWKOlBSYfzM939irt3u/ted98r6TlJ97v7YE0qBgCUVTbc3X1e0sckPS3pmKSn3P1lM3vMzO6vdYEAgPVLxmnk7kckHSma9slV2t6z8bIAABvBHaoAECDCHQACRLgDQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgAAR7gAQIMIdAAJEuANAgAh3AAgQ4Q4AASLcASBAhDsABIhwB4AAEe4AECDCHQACRLgDQIAIdwAIEOEOAAEi3AEgQLHC3czuM7NXzey4mT1aYv5vmdlRM3vJzJ4xsz3VLxUAEFfZcDezhKTHJX1A0h2SHjKzO4qavShpwN3fJunLkv6g2oUCAOKLc+b+TknH3f2Eu1+T9KSkBwobuPuz7j6VG31OUl91ywQArEeccO+VNFQwPpybtpqPSvq7UjPM7GEzGzSzwdHR0fhVAgDWJU64W4lpXrKh2YclDUj6bKn57v6Euw+4+0A2m41fJQBgXZIx2gxL6i8Y75N0priRmb1f0u9K+gl3n61OeQCASsQ5c39e0iEz22dmjZIelHS4sIGZ3Snpf0q6391Hql8mAGA9yoa7u89L+pikpyUdk/SUu79sZo+Z2f25Zp+V1CbpS2b2HTM7vMrqAACbIE63jNz9iKQjRdM+WTD8/irXBQDYAO5QBYAAEe4AECDCHQACRLgDQIAIdwAIEOEOAAGKdSkksFW5u2bnF3V1Zk4TM/OanF1Qc6pB6aak0k1JtTUllWgo9QQNIGyEO+pmfmFRk7MLUTDPzmt8Zl7jueGrM/OayI2Pz8zn5s/l2qwcn18s+aijJS2phNJNSWWak0o3JZRuzA9HX5mmZMGLQUJtTSmlmxLKNKfUk2lSNtOk5lRik/YKUB2EO9bN3TU9t6CJmVwIFwRtNO36sM6HcmFYT11bKLutRIMp0xyFcVtTSpnmpG7Z2ay2pqQyzSm15eZlmlPKNCXV2pjQ7PyiJmejbUzMzueGFwqG53Xm8owmr0Xj4zPzmp1fXLOOTHNS2UxTLuyblW1rUs+OpuXvmWi4o7VRDbxTwBZAuG8zcwuLmsiFa74rY3xmXuOzc6uGdTR/ZTAvlDlblqR0YyIXvqlcGCfVu7NlaTg/L9McnT2vCOvceHOqQWa1D8u5hUVNzS5ofHZOk7kXgqvTcxqdmNXo+PLXyPiMvjd8WSPjsyVfnJINpu62poIXguj7rvZm9XW0qq+jRb07W3gngJoj3G9QM3MLGr40rTfHJjVydXbtroyCs+iZubXPUKUooDIFwdvWFIVypjlz3Vl0puDMOR/amaYopG+kvu5UokHtrQ1qb03FXmZydj4X+MvBv/wiMKuzV2b03eErujg5Ky96LcxmmtTX0bIU+IXDhD+qgXDfwq7OzOnUhUmdvDCpNy9O6c2xKb0xNqWhsSmduzpzXWBIUlvuQ8R86La3NqqvszV3JlwUwkUBnh9uSm7O2fKNLt9Pv7c7vWa7+YVFjYzP6vTlaQ1fmtLw2LSGL01r+PKUXhq+rL///lnNLaz8ZWYzTervaNGerrT2dLVqb1dau3PfO1pT/H5QlnmphNgEAwMDPjg4WJdtbyVT1+Z1YjQK8FMXJnXyYvT91MUpjU1eW9G2J9OkPV2t6u9s1Z7OtHZ3tWh3Z6t27WjWjpaU0o031tkyIguLrpHxmSjwC8L/jbHoRf1s0Qt5pjmpfd3pFV/7u9u0t7tVmeb47zxwYzKzF9x9oFw7ztw3wcKi68zlab0+OqETo5M6cWFiKdDPXplZ0famHc3a292qe//VLu3tSmtvdzo6a+tsVUsjb9VDlGgw3dzeopvbW/TDezuvmx91wU3pjYtTOnlhUm9cnNKpi5N64Y1LOvzdMyuCP5tpyoV9WvuzaR3Itml/tk39HS1KJritZTsh3Ksofxb++uiEXs9/H5nQyQuTK67G2NGc1P5sm969v0v7s2ntz7ZpX3f09ru1kV8JVmpOJXSwJ6ODPZnr5s3MLejNsamlk4WTF6Lj7WtHz69455dKmPZ2LQf+gWybDvS06UA2zdl+oEiSClyavKbjoxN67fyEjo9M6HguxE9fnl5q02BSf2erDmbb9OO3ZrW/O60DPVGId6Ub6TNFVTSnErp1V0a37ro++C9PXVs+yci9a3xtZELPHBtZcW/ATTuadbCnTQd7osA/mI2Gu9s4Tm9khPsq3F3nrs5E4T0yoddy318fmdDFgjOi5lSDDmTbNLC3Qx/K9kd/INk27elq5YoH1NXO1kbdtadRd+3pWDF9bmFRb45NRcfz6PJx/aXBIU0WXN7Z3pKKQj8X9gd3RcO9O1u4lv8GsO3DfWHRNXwpOtA//Xev6PjIhN7ev1Ovj0xoYnZ+qV3+QH//W3YtneUc7OFAx40nlWhY6popVHxCkz+p+Ydj5/VXg0NL7VpSCR3oSetgtk2HdmV0IBf+e7palaJff8vYNuE+M7egUxcnVxy4x0v0h0vRzTc/+45eHexp06GeDG9RsS2YLX+w+2OHsivmleqK/OeTY/rqd84stUklTHu68qHftvQu9kC2jYsB6iC4cL8yNbfUB55/y3l8dEJDY1PKdzOaSX0dLTqYbdOPHerWoZ5M1NfY06b2Fj5cAop1pBv1w+nO667mmZyd1+v50M99f/X8uL529NyKv7fenS1LXTwHlkI/rU4+f6qZGzLcFxZdpy9NL31QlP/Q6MTohC5MLPeHNyYatK87rbfe0q4H3n5LwUHFmQRQDemmpN7Wt1Nv69u5Yvrs/IJOXZhafpecO9F67sTFFXdJ72xNRd062baCSzfT6u+ki2ejbrhwf+XcVd3/P76pawVdKR25A+S9t/cshffBnjb1cW0vUBdNyYRuuymj225aeRXP4qLr9OXp6y4XfuaVEf3V4OxSu2SDaXdXq/Z3R2f4+7Np7cvdqJVta+JsP4YbLtz7Olr1iz+6VwcKbtDoTDfWuywAMTQ0mPo7o7us77lt5bx8l+qJ0eizsPwNf9/4waiuLSyfzKUbE9rTFd2Zu7c7eiRDNMxlxoV4/ACALS3fDbv8aI7lR3QMjU2tuGY/k3vWz56uVu3rju7szr+Y3LSjOYjHc/D4AQBBSOS6aHZ3teonbl15Fc/cwuLK4L8wqZMXp/TS8BUd+d5ZFT6ZOtlgumVni/o7W9TfEQV+/mmc/Z0twXX3EO4AblipREP0/KXutFTUzXNtflFnr0xraGxaQ5eis/yhS9MaGpvSPxw7v+LiCym6IbGvo1X9HS1LwZ9/EejvaNWOluQNFf6EO4AgNSYbco9MLv1I5qlr80tP4hwam86FfzQ8+MYljc/Mr2jfkkro5vZm7drRrJvac187ovGbc+PdbU1bpuuHcAewLbU2Jld9Lo8UfcCbP+MfvjSts1dmdP7qjM5dndE/nxzT+asz1/3/3kSDRf95Kxf4+ReCpeHc+GY8moRwB4AS2ltTam9t11t720vOX1x0XZy8pnNXosA/d3VG565M69yVWZ2/OqMfnB/XP712YcVjTPL+6bffo/7O1prWT7gDQAUaGiz6x+iZJv2QSr8ASNL4zFx0xn9ldukFoGdHU83rI9wBoIaif1+ZKvk8/lqKFe5mdp+kP5KUkPQ5d/9M0fwmSX8u6S5JFyV9yN1PVbfUyFdfPK3PPv2qzlye1i07W/TIvbfpp+/sLTsvP/9Th1/W5ek5SdEz1xddSphpocT1/vnppeabJFf0zIz33J7Vs6+MrrndjdRc6bLbTbX2RzX3ayXrWm2Z/PTTl6eXjsneCo+nUuuQtO5jMb/MWjWVq6vcz1XNfViJrXhcxVH2JiYzS0j6gaSflDQs6XlJD7n70YI2vyrpbe7+K2b2oKSfcfcPrbXeSm5i+uqLp/WJr3xP03PLz5xuSSX06Z/9IUladV7+AHnkS9/V3GLtb9oq3u5Gaq502e1mrX21nv1RrfVUuq7Vlvk3d/Xqr184vWJ68Tql9R9PeakGk0wr/lF3uWVLLRN32bVqLvdzV7oPN+t3WOua4t7EFCfc3y3pU+5+b278E5Lk7p8uaPN0rs23zCwp6ZykrK+x8krC/e7PfH3FfzvK693ZIkmrzvvmo+9dddlaKbfdjdQcZ9ntZq19tZ79Ua31VLqu1ZZZ7d1l4Tql6v8N1HLZtWrOW+3nrmQfbtbvsNY1VfMO1V5JQwXjw5J+ZLU27j5vZlckdUm6UFTUw5IelqTdu3fH2PRKZ1Y5CFabXjhvrTa1UG67G6k5zrLbTSX7qpbrqXRdq81bK9jjrnMjP0Mtlo2zztV+7kr24Wb9DmtdU1xxHplY6or84j0ep43c/Ql3H3D3gWw2W2KRtd2Se6UvNX2teWstWyvltruRmuMsu91Ua39Uc79Wsq7V5iXK3BlZq7+BWi67Vs15q/3clezDzfod1rqmuOKE+7Ck/oLxPklnVmuT65ZplzRWjQILPXLvbWopuvi/JZXQI/fetua8/LKpTbpzrHi7G6m50mW3m2rtj2ru10rWtdoyD/1I/3XTi9dZyfGUl2owpRIr/z7KLVtqmbjLrlVzYZtSP3el+3Czfoe1rimuON0yz0s6ZGb7JJ2W9KCkf1vU5rCkj0j6lqSfk/T1tfrbK5X/4GGtT5xXm5f/vtlXy1Sj5kqW3W7i7KvNXE+l61prmYE9nbGuTIlzPK33apnV6lprfes9jldbR/7nrsY+XK+teFzFFeuRv2b2QUl/qOhSyM+7+382s8ckDbr7YTNrlvQXku5UdMb+oLufWGudPPIXANavqo/8dfcjko4UTftkwfCMpJ9fb5EAgNrgf9ABQIAIdwAIEOEOAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CACHcACBDhDgABItwBIECEOwAEiHAHgADFep57TTZsNirpjbpsvHq6VfR/Yrc59scy9sVK7I9lG90Xe9y97P8prVu4h8DMBuM8NH+7YH8sY1+sxP5Ytln7gm4ZAAgQ4Q4AASLcN+aJehewxbA/lrEvVmJ/LNuUfUGfOwAEiDN3AAgQ4Q4AASLcK2Bm/Wb2rJkdM7OXzew3611TvZlZwsxeNLO/rXct9WZmO83sy2b2Su4YeXe9a6oXM/uPub+R75vZF82sud41bSYz+7yZjZjZ9wumdZrZ/zOz13LfO2qxbcK9MvOSPu7ub5H0Lkm/ZmZ31LmmevtNScfqXcQW8UeS/t7db5f0dm3T/WJmvZJ+Q9KAu79VUkLSg/WtatN9QdJ9RdMelfSMux+S9ExuvOoI9wq4+1l3/5fc8LiiP97e+lZVP2bWJ+lfS/pcvWupNzPbIenHJf0vSXL3a+5+ub5V1VVSUouZJSW1SjpT53o2lbt/Q9JY0eQHJP1ZbvjPJP10LbZNuG+Qme2VdKekb9e3krr6Q0m/LWmx3oVsAfsljUr637luqs+ZWbreRdWDu5+W9F8kvSnprKQr7v61+la1Jexy97NSdKIoqacWGyHcN8DM2iT9taT/4O5X611PPZjZT0kacfcX6l3LFpGU9A5Jf+Lud0qaVI3edm91ub7kByTtk3SLpLSZfbi+VW0fhHuFzCylKNj/0t2/Uu966uhuSfeb2SlJT0p6r5n9n/qWVFfDkobdPf9O7suKwn47er+kk+4+6u5zkr4i6UfrXNNWcN7Mbpak3PeRWmyEcK+AmZmiPtVj7v5f611PPbn7J9y9z933Kvqw7Ovuvm3Pztz9nKQhM7stN+l9ko7WsaR6elPSu8ysNfc38z5t0w+XixyW9JHc8Eck/U0tNpKsxUq3gbsl/TtJ3zOz7+Sm/Y67H6ljTdg6fl3SX5pZo6QTkn6pzvXUhbt/28y+LOlfFF1h9qK22WMIzOyLku6R1G1mw5J+X9JnJD1lZh9V9AL48zXZNo8fAIDw0C0DAAEi3AEgQIQ7AASIcAeAABHuABAgwh0AAkS4A0CA/j+1cdTZbGihCwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data1=sorted(np.linspace(1,3,10).tolist()+np.linspace(1.5,2,10).tolist())\n", "data2=sorted(np.linspace(5,10,10).tolist()+np.linspace(6.5,8.5,10).tolist())\n", "u=np.mean(data1+data2)\n", "sigma=np.std(data1+data2)\n", "plt.scatter(data1+data2,[0]*len(data1+data2))\n", "plt.plot(data1+data2,gaussian_1d(data1+data2,u,sigma))\n", "plt.ylim(-0.1,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "既然一个高斯分布不能很好的拟合这些数据，那就多个去拟合？这便是**高斯混合模型**，我们分别构建两个高斯分布去拟合左右两侧的数据看看效果" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.1, 1)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW9//HXJ5OFrIRANkIggEkAUQSCslyt+0J73Soq7q2KWLdar/en7f319uft7aK9t9Ve1OLeUhdcarmWSutWFxaJoIiQBcKSBLIQSMieTOb7+yMzGGOWSTIzZ+bM5/l4+DBz5iyfHCbvc+Z7vud7xBiDUkope4mwugCllFK+p+GulFI2pOGulFI2pOGulFI2pOGulFI2pOGulFI2NGi4i8jTIlIjItv7eV9E5BER2SUi20Rkju/LVEopNRTenLk/C5w/wPsXALnu/5YBj428LKWUUiMxaLgbY94HDg8wy0XA7023jUCyiGT6qkCllFJDF+mDdWQB5T1eV7inHew9o4gso/vsnvj4+LnTpk3zweaVUip8fPLJJ4eMMamDzeeLcJc+pvU5poExZiWwEqCgoMAUFhb6YPNKKRU+RGSfN/P5ordMBZDd4/UE4IAP1quUUmqYfBHua4Dr3L1m5gMNxpivNckopZQKnEGbZUTkBeB0YJyIVAD/DkQBGGMeB9YCi4FdQAvwHX8Vq5RSyjuDhrsxZukg7xvgNp9VpJRSasT0DlWllLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXellLIhDXcv1TW182HpIdo6u6wuRSmlBhVpdQHBrqXDyfJVW/igtBZjYPbEZJ64roBxCTFWl6aUUv3SM/dB/GztTj4oreW204/jpxfPZOfBo1z66HoONbVbXZpSSvVLw30A7xbVsGrjfm76p8n8y3n5XDN/En+8aT4VR1p4/L3dVpenlFL98ircReR8ESkWkV0icl8f708UkXdFZKuIbBORxb4vNbCcXS5++KfPyU9P5J5z849NnztpDBfPzmLVpn3UNLZZWKFSSvVv0HAXEQewArgAmAEsFZEZvWb7N2C1MWY2cCXwqK8LDbR3imo42NDGPefmMSrK8ZX37jwzl84uw+PvlVlUnVJKDcybM/eTgV3GmDJjTAfwInBRr3kMkOT+eTRwwHclWuOFj/eTlhjDmdPSvvZezrh4LpmdxR837aOhpdOC6pRSamDehHsWUN7jdYV7Wk8/Aa4RkQpgLXBHXysSkWUiUigihbW1tcMoNzAq61v5R0ktlxdkE+noexctPXki7U4X75cG7++hlApf3oS79DHN9Hq9FHjWGDMBWAz8QUS+tm5jzEpjTIExpiA1NXXo1QbI6s3luAxcMS+733lOyk4mOS6K94o13JVSwcebcK8AeqbcBL7e7HIjsBrAGLMBGAWM80WBVljz2QEWHTeW7JS4fudxRAin5abyj5IaXK7exzqllLKWN+G+GcgVkckiEk33BdM1vebZD5wFICLT6Q73kDyl3VfXzJ5DzZwzPX3QeU/PT+VQUwdfHDgagMqUUsp7g4a7McYJ3A6sA3bS3SvmCxF5QEQudM92D3CziHwGvADcYIwJydNZTzPL6flfv5Da22l5qe5lavxak1JKDZVXww8YY9bSfaG057Qf9/h5B7DIt6VZ473iGnLGxpEzLn7QecclxDBrwmjeK6nljrNyA1CdUkp5R+9Q7aGts4sNZXVenbV7fCM/ja37j1Df0uHHypRSamg03HvYtOcwbZ0uvpHvfU+e0/NTcRl4v/SQHytTSqmh0XDv4R/FtcRERrBgylivl5k1IZkxcVHa7q6UCioa7j1s2lPH3EljvjbcwEAcEcJpeam8X1KrXSKVUkFDw92tud3JzoNHKZg0ZsjLerpEbj/Q4IfKlFJq6DTc3T4tr8dlYG5OypCXPS03FRH0blWlVNDQcHf7ZN8RRLqftDRUYxNiODFrNP8o0XBXSgUHDXe3wn1HyE9PJGlU1LCWn5eTwvbKBpxdLh9XppRSQ6fhDrhchq37jjB3GO3tHtMzk2h3uthb1+zDypRSang03IGSmkYa250jDneAHQcbfVWWUkoNm4Y7ULj3CAAFk4Z+MdVjalo8kRHCzoM6iJhSynoa7sBn5fWMjY8mOyV22OuIiXRwXFoCRRruSqkgoOEOfF7ZwMys0Yj09VwS703PTGKnNssopYJA2Id7W2cXpTVNnJA1esTrmp6ZSNXRNo406yBiSilrhX24F1U10uUyzMxKGnzmQXguqmq7u1LKamEf7tsru4cMmOmDM/dpGe5wr9KmGaWUtTTcKxtIjosiK3n4F1M9UhNjGJcQo2fuSinLabgfaOAEH1xM9ZiemajhrpSyXFiHe7uzi+KqRo4fP/ImGY8ZmUmUVjfRqcMQKKUsFNbh3h3Cxic9ZTymZSbS0eVizyEdhkApZZ2wDvcvL6aOvKeMh/aYUUoFg7AO96KqRuKjHWSPifPZOqemJhDtiGCHhrtSykJhHe47Dx4lPyORiAjfXEwFiHJEcFxagt6pqpSyVNiGuzGGoqpGpmX6rknGY1pmoo4xo5SyVNiGe9XRNhpaO5mekejzdc/ITKKmsZ26pnafr1sppbwRtuFe5G42yc/w/Zn7lxdVtWlGKWWNsA33nVXdzSb5fjhzn+Zep/aYUUpZJWzDvbiqkazkWEbHDu+ZqQMZmxBDWmLMsQOIUkoFWtiGe9HBxmNn2P6gY7srpawUluHe7uxid22TX5pkPKZnJrGrppEOpw5DoJQKPK/CXUTOF5FiEdklIvf1M8/lIrJDRL4Qked9W6Zv7a5pxukyfukG6TE9M5HOLsPu2ia/bUMppfoTOdgMIuIAVgDnABXAZhFZY4zZ0WOeXOB+YJEx5oiIpPmrYF8oru5uC/d3swxAUdXRYz8rpVSgeHPmfjKwyxhTZozpAF4ELuo1z83ACmPMEQBjTI1vy/StkuomIiOEnLHxftvGlHHxREdGaLu7UsoS3oR7FlDe43WFe1pPeUCeiHwkIhtF5Py+ViQiy0SkUEQKa2trh1exD5RWNzLZHb7+EumIIC89QbtDKqUs4U269TXwiun1OhLIBU4HlgJPikjy1xYyZqUxpsAYU5CamjrUWn2mpLqJPD82yXhMz0jScFdKWcKbcK8Asnu8ngAc6GOePxtjOo0xe4BiusM+6LR0ONl/uIW8NP+H+7TMJA41dVDT2Ob3bSmlVE/ehPtmIFdEJotINHAlsKbXPK8DZwCIyDi6m2nKfFmor+yq6e69kpee4PdtTc/sPoAUabu7UirABg13Y4wTuB1YB+wEVhtjvhCRB0TkQvds64A6EdkBvAvca4yp81fRI1FS3R3uuen+P3OfoQ/uUEpZZNCukADGmLXA2l7TftzjZwP8wP1fUCutbiTaEUHOWN89oKM/yXHRZI4epeGulAq4sLtDtaS6kSmp8UQ6AvOrT8tI1O6QSqmAC8NwbyIvAE0yHtMzk9hd20S7sytg21RKqbAK96Z2J5X1rQG5mOoxPTMJp8scu5CrlFKBEFbhXlrd3TwSiIupHvrgDqWUFcIs3D3dIAMX7pPHxRMTGaEXVZVSARVW4V5S3UhMZAQTU/zfU8bDESHkZyRSpA/uUEoFUHiFe00Tx6Ul4Ijoa0QF/+kehqCR7h6jSinlf+EV7lWNAW2S8Ziemcjh5g5qGtsDvm2lVHgKm3BvaO2k6mgbuQHsKePhuai6Q9vdlVIBEjbhvqumu7dKIAYM621aRne4l1RpjxmlVGCETbh7xpTx53NT+zM6Lor0pBiKqzXclVKBEUbh3khslIOs5FhLtp+XnkiJhrtSKkDCKtxz0xOICHBPGY/89ERKq5vocmmPGaWU/4VRuDeRa0F7u0deRiLtThflh1ssq0EpFT7CItzrWzqobWwP6JgyveW7u2Bqu7tSKhDCItw9F1MD8dzU/ni6YGqPGaVUIIRJuLu7QVpwA5NHXHQkE1Pi9MxdKRUQYRPuCTGRjB89ytI68tITtMeMUiogwibcj0tLQMSanjIeeemJlNU20+F0WVqHUsr+wiLcS6ubjl3QtFJ+RiJOl2HPoWarS1FK2Zztw72uqZ265g5LxpTpLU97zCilAsT24V5iwQM6+jMlNR5HhBx7IpRSSvlLGIS79T1lPGIiHUweF0+xdodUSvlZWIR74qhI0pNirC4F6L6ZSXvMKKX8zfbh7rmYanVPGY+89ET2HW6htaPL6lKUUjZm63A3xlBS00huEDTJeORnJGAM7KppsroUpZSN2Trca5vaqW/ptHRMmd5ytceMUioAbB3uJVXB01PGY1JKHNGREdrurpTyK3uHuztAg6GPu0ekI4LjUhO0x4xSyq9sHe6lNY2MiYsiNSE4esp45Gckal93pZRf2TrcS6qbyA2injIeeemJHGho42hbp9WlKKVsyqtwF5HzRaRYRHaJyH0DzHeZiBgRKfBdicNjjKGkujGoLqZ65Gd016Rn70opfxk03EXEAawALgBmAEtFZEYf8yUCdwKbfF3kcFQdbaOxzRlUF1M9jo0xU6XdIZVS/uHNmfvJwC5jTJkxpgN4Ebioj/n+A3gQaPNhfcPmGVPGyuem9icrOZb4aIf2mFFK+Y034Z4FlPd4XeGedoyIzAayjTFvDLQiEVkmIoUiUlhbWzvkYoei9NiYMsHXLCMi5GUkao8ZpZTfeBPufV2NNMfeFIkAfg3cM9iKjDErjTEFxpiC1NRU76schpLqRsYlRDM2yHrKeOSl6RgzSin/8SbcK4DsHq8nAAd6vE4EZgLvicheYD6wxuqLqiXVTUHZJOORl5FIXXMHh5rarS5FKWVD3oT7ZiBXRCaLSDRwJbDG86YxpsEYM84Yk2OMyQE2AhcaYwr9UrEXjDGUBmlPGQ/Pk6FKtGlGKeUHg4a7McYJ3A6sA3YCq40xX4jIAyJyob8LHI7K+laaO7rIywjmM/fuA482zSil/CHSm5mMMWuBtb2m/bifeU8feVkjUxpET1/qT2pCDGPioiiu1u6QSinfs+UdqseevhTEbe4iQp4+uEMp5Sc2Dfcm0hJjGB0XZXUpA8rPSKSkqhFjzOAzK6XUENg03BuDuknGIy89kcZ2JwcbguK+L6WUjdgu3F0uw66appAI9/wMfXCHUso/bBfuFUdaae3sCupukB6eawLaHVIp5Wu2C/cvH9AR/Gfuo+OiSE+K0TN3pZTP2S7ci4Pw6UsDyUtPPNZ1UymlfMWrfu6hpKS6kczRo0gaFdw9ZTzy0xNZtWkfXS6DIyK4HiqiwovLZSiqamR3bf8nGyLdI63mpiUQoZ/XoGa7cP+0vJ4TJ4y2ugyv5WUk0tbpovxwCznj4q0uR4URT5hvLKtjY1kdH+89TH2Ld08HS4mP5pTJKcyfMpb5U8Zq2AchW4V7XVM7++pauOrkiVaX4jXPGDPF1Y0a7iogiqqO8uLH5az57ACHmzsAmJgSx7kz0pk/ZSwzxicR2U9Qd3YZtlc2sLHsMBvL6vjr9ioAxiXEcPFJ47ny5IkclxYaTaJ2Z6tw37q/HoDZE8dYXIn3PNcGSqoaOe/4DIurUXa291AzD/2tmL9sO0i0I4Jzjk/nrGlpnDJlLFnJsV6vZ3pmEksKugeKLT/cwsayOt7aWc2z6/fy1Ed7uOSkLO4+J4/slDh//SrKC7YK9y37jxAZIZyQFTrNMnHRkUxMidMeM8pvahvbeeTtUl74eD9RjgjuPPM4vrNoMmPio0e87uyUOLJT4lhSkM2hpnae+KCMZz7ayxvbDnLtgkncdsZxpPhgO2robBXuW/fXMz0zidhoh9WlDEleeiJF2tdd+VhTu5OV75fx5AdltDtdLD05mzvPyiUtcZRftjcuIYb7L5jO9Qty+M1bJTzz0R5Wby5n+elT+c6iHOKibRU3Qc82e9vZ5eKzinqWzJ1gdSlDdlL2aN7aWU1DS2fQj4ejgl+H08Xzm/bx23d2UdfcwTdPyOSec/OYkhqYtvDxybE8eNksbjp1Cg++WcxD64p5bv1evn92HpcXTCDSYbse2EHJNuFeXN1IS0cXcyaFTnu7R0FOCgCf7D/MmdPSLa5GhbJ9dc3c/vxWPq9sYP6UFJ66YDonZSdbUkteeiJPXl/A5r2H+cVfi/jhnz7n1S0V/HbpbMYPoY1fDY9tDqHHLqZmh164z5qQTJRD2Lz3iNWlqBD2l20H+dYjH7KvrpnHrp7DCzfPtyzYe5qXk8IryxfwmytOoujgURY/8gFv7ai2uizbs024b9l/hHEJ0WSnhN4ZQWy0g5lZoynce9jqUlQIauvs4t9e/5zbnt/CcekJrL3rVC44IROR4Ol3LiJcPDuLN+48lfGjY7np94X89I0ddDhdVpdmW7YJ90/31zN74pig+kAPxbycFD4rb6Cts8vqUlQIKatt4pJH17Nq435uOW0Kq29ZwIQxwdsFcfK4eF773kKuWzCJJz/cw5LfbaD8cIvVZdmSLcL9SHMHZYeamT3R+q+gw1UwaQwdXS62VzZYXYoKEa9vreRbv/2QqoZWnr6hgPsXTycqBC5Wjopy8MBFM3ns6jmU1Tax+JEPeHP7QavLsp3g/yR4YWt5d1v1nBC6eam3ue4LwdrurgbT7uzivle38f2XPmXm+NGsvevUkLwQf8EJmfzljlOZMi6e5au28JM1X9DZpc00vmKPcN9fjyNCQmpMmd7GJsQwNTVe293VgOqa2rn6iU28uLmc286YyvM3n0Lm6NC7zuQxcWwcLy9fyI3/NJln1+/lhmc+psHL8W3UwGwR7lv2H2FaRmLI3yRRMCmFT/YfweXSZ6qqr9tzqJlLHl3P55UNrLhqDveeN80WfcajIyP4v9+awa+WzOLjPYe59LGPqDii7fAjFfKfjC6X4bPyhpBukvEoyBlDfUvngEOuqvD0aXk9335sPU3tTl5cNp9vnphpdUk+d9ncCfzhxlOobWzn0kfXs+PAUatLCmkhH+6lNY00tTtD+mKqxzz3zUza7q56ereohqUrNxIf4+DVWxeG1MB4QzV/ylheXr6QCBGu+N0G1u8+ZHVJISvkw33Lvu6bl+xw5j5pbBzjEmK03V0d83JhOTf9vpCpafG8dusiJofBsND5GYm89r2FZIwexQ1Pb+Z/PztgdUkhKeTDfev+I6TERzNpbPD27fWWiDAvZwyb92m4hztjDP/zTin3vrKNhVPH8uKyBaQmxlhdVsCMT47lleULOSk7mTte2MpTH+6xuqSQE/LhvmX/EWZnJ4fszUu9FeSkUH64laqGNqtLURZxuQz/73938Ku/lXDJ7Cyeun4eCTGh3VlgOEbHRfH7G0/m/OMz+I83dvDLN4swRjsbeCukw72hpZPdtc0hOVhYf+bldP8uhXr2Hpa6XIb/8+o2nl2/lxv/aTL/tWQW0ZEh/Wc6IqOiHKy4eg5LT57IY+/t5idrvtDeZF4K6dMBz81LdriY6jEjM4m4aAeFe4/wrRPHW12OCqDOLhd3v/Qpb2w7yF1n5fL9s3Nt8410JBwRws8umUl8tIMnP9xDc0cXv/z2ifpA+UGEdLhv2V9PhHSPqmgXkY4IZk9MZrNeVA0rbZ1d3P78Ft7aWcMPF09j2WlTrS4pqIgIP/rmdOJjInn47VJaO7v49eUnhfW3msF4tWdE5HwRKRaRXSJyXx/v/0BEdojINhF5W0Qm+b7Ur9u6/wj5GUnE26w9cu6kFHYePEpjm96pFw6a253c+Nxm3i6q4acXz9Rg74eIcPc5efxo8XT+su0gt676RAfaG8Cg4S4iDmAFcAEwA1gqIjN6zbYVKDDGnAi8Ajzo60J7c7kMn5bXM8dGTTIe83LG4DJfjlGv7KuhtZNrn9rEht11/NeSWVwzPyDnRSHt5tOm8NOLZ/JOcQ3ffXYzze1Oq0sKSt6cuZ8M7DLGlBljOoAXgYt6zmCMedcY47lfeCPg92fd7a5torHNacsbOmZPHEOEQOE+vZnJzuqa2rnqiY18XtnAo1fP4dI5ofeISKtcM38S/7VkFhvL6rj2qU00tOq33N68CfcsoLzH6wr3tP7cCPy1rzdEZJmIFIpIYW1trfdV9mHLfs9IkPY7c0+IiWTG+CS9mcnGqo+2ccXKjeyqaeKJ6wo4f6b9hhPwt0vnTODRq+fweWUDS1dupK6p3eqSgoo34d7XJek++yKJyDVAAfBQX+8bY1YaYwqMMQWpqaneV9mHLfvqSY6Lsu0dewWTUti6v16HQLWh8sMtLHl8AwfrW3nuuydzen6a1SWFrPNnZvLEdQXsrm3iipUb9f6QHrwJ9wogu8frCcDX7gcWkbOBHwEXGmP8fgjdWm6vm5d6m5eTQmtnlw6eZDO7a5tY8vgG6ls6WHXTKcyfMtbqkkLe6flpPPfdkzlY38rl+mSnY7wJ981ArohMFpFo4EpgTc8ZRGQ28Du6g73G92V+1dG2Tkprmmwxnkx/CnI8D+/Qphm72HnwKFf8bgOdXS5eXLbAlteLrDJ/ylhW3XQK9S0dLHl8g46sihfhboxxArcD64CdwGpjzBci8oCIXOie7SEgAXhZRD4VkTX9rM4nPiuvxxhs/ceRnjSKiSlxFOoIkbbwaXk9V67cSJQjgtXLFzBjfJLVJdnO7IljeOmWBThdLq743QZ2Hgzvb71e9XM3xqw1xuQZY6YaY/7TPe3Hxpg17p/PNsakG2NOcv934cBrHJkt++oRgVnZofvkJW8U5IyhcN9hHU8jxK3ffYirn9jI6NgoVt+ygKmpCVaXZFvTM5N46ZYFRDkiuOJ3G451vAhHIXl715b9R8hPTyRxVJTVpfhVwaQUDjV1sLdO2xBD1erCcq5/+mMyk2N5efkCslNCf/TSYDc1NYHVtyxgTHw0S1duZE2YDhkccuHuuXnJTuPJ9GeetruHrC6X4Wdrd/Kvr2zjlMljeXX5QtKTRlldVtjITonjtVsXcuKE0dz5wlZ+/feSsPsGHHLhXnaomYbWTlu3t3tMTU0gOS5K+7uHmKZ2J7f8oZCV75dx7fxJPPOdeYyOs/e3zGA0NiGGVTedwrfnTODht0u5/YWtYTVcQcgNypKeFMOjV89hro2G+e1PRIRQMGmM3qkaQiqOtHDTc4WU1jTxwEXHc92CHKtLCmsxkQ5+teREctMT+OWbRVQcbuGJ6wpIC4NvUSF35p44KorFJ2SGzVfcgpwUymqb9e67ELBl/xEuXvERlUdaeeaGeRrsQUJEWP6NqTx+zVxKqpu4aMVHbK9ssLosvwu5cA83nnb3j3bXWVyJGsifP63kypUbiYuO5E+3LeS0vJHdga1877zjM3jl1gUALHl8A+u+qLK4Iv/ScA9ysyYkk5Ucy/Ob9lldiuqDy2X4778Vc9eLn3JSdjKv37aI49ISrS5L9eP48aP5822LyMtIZPmqT3jsvd22vdCq4R7kIh0RXLtgEhvLDlNUFd43ZQSb1o4u7nhhK4+8s4slcyew6sZTSImPtrosNYi0pFG8tGw+3zwhk1++WcS/vLyNdqf9LrRquIeAKwqyiYmM4Ln1e60uRbl1j+q4gbXbD/LDxdN48LIT9alAIWRUlIPfLp3N98/O5dUtFVzz5CbbXdfST2MIGBMfzSWzs/jT1krqWzqsLifsba9s4ML/+ZBdNU2svLaAZadNte0AdnYmInz/7Dx+u3Q22yoauGjFR5RUN1pdls9ouIeI6xfm0Nbp4qXN5YPPrPzmze0Huezx9ThEeGX5Qs6ZkW51SWqE/nnWeF66ZQHtTheXPrqed4v9PvZhQGi4h4jpmUmcMjmF32/YR5fLnheAgpkxhhXv7mL5qi1Mz0zi9dsX6eBfNnJSdjJ/vm0RE1PiuPHZzTz94Z6Qv9Cq4R5CbliYQ2V9K2/trLa6lLDS7uzintWf8dC6Yi6cNZ4Xbp5PWmJ43GcRTsa7x/85e3o6D7yxgx+9vj2kH5aj4R5CzpmRzvjRo/TCagAdamrnqic28drWSn5wTh4PX3kSo6IcVpel/CQ+JpLHr5nL8m9M5flN+7nhmY9paAnN57NquIeQSEcE1yyYxPrddRRX2efCT7Aqrmrk4hUf8cWBBlZcNYc7z8rVC6dhICJCuO+CafxqySw+3nOYSx79iD2Hmq0ua8g03EPMlfMmdneL3LDX6lJsq7PLxaqN+/j2Y+vpcLpYfcsCvnmiPsA63Fw2dwLP3zyf+tZOLl7xEasLy0PqepeGe4hJiY/mopPG86ctlSH7dTFYGWNY90UV5/3mff7t9e3MGJ/En29fxIkT7D+8tOrbvJwUXv/eIiaPi+dfX9nG4oc/4N3impC42KrhHoKuX5hDa2cXqwu1W6SvfLLvMJc9voFb/vAJAjxxXQEvLZtP5uhYq0tTFps4No4/fW8hK66aQ5uzi+88s5mrntjEtop6q0sbkFh1BCooKDCFhYWWbNsOLn98AwcaWvnHvWfgiNB24OHaXdvEg28Wse6LalITY7j77DwuL5hApEPPe9TXdThdvPDxfh5+u5TDzR3886zx3HtuPhPHBu4JWyLyiTGmYLD5Qm48d9XthkU5fO+PW3inqEZvpBmGmsY2Hn6rlBc3lzMqMoIfnJPHTadOJi5a/yRU/6IjI7h+YQ6Xzsli5ftlPPFBGW9uP8i183O4/czjgmpsIT1zD1HOLhenPvguU1Lj+eNN860uJ2Q0tzuP/VF2OF1cdcpE7jwrl3EJMVaXpkJQ9dE2fv33ElYXlhMfHcmtZ0zlu4sm+7W7rLdn7hruIWzFu7t4aF0xf7/7NHLTdZjZgXR2uXhxczkPv1XKoaZ2Fp+Qwb3nTWPyuHirS1M2UFrdyC/fLOKtnTVkJI3iB+fm8e05E/zSZKrhHgYON3cw/+dvc3JOCmdMSyMu2uH+L5K4aAex7tfx0ZHHfo6NcoRVX+3uHjDVPPhmEWWHmjk5J4X7Fk9jThg8g1cF3qaya4W+AAAIsklEQVSyOn721yI+K68nPz2R+y6Yxun5qT79m9NwDxMP/O8Onv5oz5CWiY1yEB/jDv+o7uCPj3EQGxV57AAR2+ug0D0tkvhjB43IYweL+Jjun2MiI4LqwFG49zA//2sRn+w7wnFpCfyf86dx9vS0oKpR2Y8xhrWfV/HQuiL21rWwYMpY7l88zWddajXcw0iH00VrRxctnU5aOrpoae+ipcNJS2dX9/QO92v3z60dTpo7PO853dO63NN6vnYylHs2IqT7wBEX82Xw9/wm4TlA9PcNo695PQeTaIf3B45dNd09YP62o5q0xBjuPiePJXO1B4wKLE/PmkfeLqXOhz1rNNzViBljaD924HAfFNrdB4jOLw8WLe0DHUi++tpzcGnp7GIoHz1HhBAX5f5GERN57MDRu+mpqd3JG9sOEhvl4JbTpnCj9oBRFmts62Tl+2U8+cEenC4X18yfxF1n5ZIcN7yeNdoVUo2YiDAqysGoKAe+bqE2xtDW6frym0NnF83tzi8PED0OJq2d/X/DaGxzUn207dh7nV0urjllIndoDxgVJBJHRXHPuflcM38Sv3mrhFcKK7j1G1P9vl09c1dKqQBqaOlkdFzUsJf39sxdGyGVUiqARhLsQ6HhrpRSNqThrpRSNuTVBVUROR94GHAATxpjftHr/Rjg98BcoA64whiz17eldnt9ayUPrSvmQH0r45Njufe8fC6enTXoe573f7LmC+pbu4fKjRBwGXCI0NXHtQfP9L7eF8AAWcmxnDEtlXeLagfc7khqHu6y4cZX+8OX+3U46+pvGc/0yvrWY5/JrGF+nvpaBzDkz6JnmYFqGqyuwX4vX+7D4QjGz5U3Br2gKiIOoAQ4B6gANgNLjTE7eszzPeBEY8xyEbkSuMQYc8VA6x3OBdXXt1Zy/2uf09rZdWxabJSDn196AkC/73k+IPe+/BmdARhsv/d2R1LzcJcNNwPtq6HsD1+tZ7jr6m+Zb8/N4tVPKr8yvfc6YeifJ4+oCAGBzi7j9bJ9LePtsgPVPNjvPdx9GKh/Q3/X5LN+7iKyAPiJMeY89+v7AYwxP+8xzzr3PBtEJBKoAlLNACsfTrgv+sU7VNa3fm16VnL3mNv9vffRfWf2u6y/DLbdkdTszbLhZqB9NZT94av1DHdd/S3T37fLnusE3/8N+HPZgWr26O/3Hs4+DNS/ob9r8mU/9yyg51MhKoBT+pvHGOMUkQZgLHCoV1HLgGUAEydO9GLTX3Wgnw9Bf9N7vjfQPP4w2HZHUrM3y4ab4ewrf65nuOvq772Bgt3bdY7kd/DHst6ss7/fezj7MFD/hv6uyVveXFDt657v3nvcm3kwxqw0xhQYYwpSU1O9qe8rxif3/VSc8cmxA7430LL+Mth2R1KzN8uGG1/tD1/u1+Gsq7/3HIMMveCvvwF/LjtQzR79/d7D2YeB+jf0d03e8ibcK4DsHq8nAAf6m8fdLDMaOOyLAnu697x8YnuNkxwb5eDe8/IHfM+zbFSAnljUe7sjqXm4y4YbX+0PX+7X4ayrv2WWnpL9tem91zmcz5NHVIQQ5fjq38dgy/a1jLfLDlRzz3n6+r2Huw8D9W/o75q85U2zzGYgV0QmA5XAlcBVveZZA1wPbAAuA94ZqL19uDwXHga64tzfe57/B7q3jC9qHs6y4cabfRXI9Qx3XQMtUzApxaueKd58nobaW6a/ugZa31A/x/2tw/N7+2IfDlUwfq685dXwAyKyGPgN3V0hnzbG/KeIPAAUGmPWiMgo4A/AbLrP2K80xpQNtE4dfkAppYbOpwOHGWPWAmt7Tftxj5/bgCVDLVIppZR/6B2qSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQxruSillQ16N5+6XDYvUAvss2bjvjKPXc2LDnO6PL+m++CrdH18a6b6YZIwZ9DmlloW7HYhIoTeD5ocL3R9f0n3xVbo/vhSofaHNMkopZUMa7kopZUMa7iOz0uoCgozujy/pvvgq3R9fCsi+0DZ3pZSyIT1zV0opG9JwV0opG9JwHwYRyRaRd0Vkp4h8ISJ3WV2T1UTEISJbReQNq2uxmogki8grIlLk/owssLomq4jI3e6/ke0i8oKIjLK6pkASkadFpEZEtveYliIifxeRUvf/x/hj2xruw+ME7jHGTAfmA7eJyAyLa7LaXcBOq4sIEg8DbxpjpgGzCNP9IiJZwJ1AgTFmJuAArrS2qoB7Fji/17T7gLeNMbnA2+7XPqfhPgzGmIPGmC3unxvp/uPNsrYq64jIBOCbwJNW12I1EUkCTgOeAjDGdBhj6q2tylKRQKyIRAJxwAGL6wkoY8z7wOFeky8CnnP//BxwsT+2reE+QiKSA8wGNllbiaV+A/wr4LK6kCAwBagFnnE3Uz0pIvFWF2UFY0wl8CtgP3AQaDDG/M3aqoJCujHmIHSfKAJp/tiIhvsIiEgC8CrwfWPMUavrsYKIfAuoMcZ8YnUtQSISmAM8ZoyZDTTjp6/dwc7dlnwRMBkYD8SLyDXWVhU+NNyHSUSi6A72PxpjXrO6HgstAi4Ukb3Ai8CZIrLK2pIsVQFUGGM83+ReoTvsw9HZwB5jTK0xphN4DVhocU3BoFpEMgHc/6/xx0Y03IdBRITuNtWdxpj/troeKxlj7jfGTDDG5NB9sewdY0zYnp0ZY6qAchHJd086C9hhYUlW2g/MF5E499/MWYTpxeVe1gDXu3++HvizPzYS6Y+VhoFFwLXA5yLyqXvaD40xay2sSQWPO4A/ikg0UAZ8x+J6LGGM2SQirwBb6O5htpUwG4ZARF4ATgfGiUgF8O/AL4DVInIj3QfAJX7Ztg4/oJRS9qPNMkopZUMa7kopZUMa7kopZUMa7kopZUMa7kopZUMa7kopZUMa7kopZUP/H7KjrTm3KbSuAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data1=sorted(np.linspace(1,3,10).tolist()+np.linspace(1.5,2,10).tolist())\n", "data2=sorted(np.linspace(5,10,10).tolist()+np.linspace(6.5,8.5,10).tolist())\n", "u1=np.mean(data1)\n", "u2=np.mean(data2)\n", "sigma1=np.std(data1)\n", "sigma2=np.std(data2)\n", "plt.scatter(data1+data2,[0]*len(data1+data2))\n", "plt.plot(data1+data2,gaussian_1d(data1+data2,u1,sigma1)+gaussian_1d(data1+data2,u2,sigma2))\n", "plt.ylim(-0.1,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "nice,该凸的地方凸了,该凹的地方凹了，但这里还有一个问题，如果对该概率分布做积分，其结果为2，因为包含了两个高斯分布，每个高斯分布的积分都为1，所以我们还需要对每个高斯分布的权重进行约束，假设$\\alpha_1$为第一个高斯分布的权重，$\\alpha_2$为第二个高斯分布的权重，那么，需要满足如下约束： \n", "\n", "$$\n", "\\alpha_1+\\alpha_2=1,\\alpha_1\\geq 0,\\alpha_2\\geq 0\n", "$$ \n", "\n", "我们不妨假设$\\alpha_1=0.5,\\alpha_2=0.5$，所以一个比较合理的分布如下： " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-0.1, 1)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt4VfWd7/H3NzcI4U4CgSTcwx0RjMhFFLUWRMVbW7HWaqfq9GbnzHQ8R+ec0/F02mc67TzT6cxjbaltFdtqvUuRqvUOAkoABQExAQJJCBDCHUKu3/NHEowxl52wd3ay9uf1PDzsvdZv/dZ3r+x89spvXba5OyIiEixx0S5ARETCT+EuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIB1Ga4m9lvzeygmX3Ywnwzs/8ys3wz22xmM8JfpoiItEcoe+6PAAtbmX8VkF3/727goXMvS0REzkWb4e7ubwOHW2lyHbDM66wD+pvZ0HAVKCIi7ZcQhj4ygMJGz4vqp5U0bWhmd1O3d09KSsoFEyZMCMPqRURix4YNGw65e1pb7cIR7tbMtGbvaeDuS4GlADk5OZ6bmxuG1YuIxA4z2xNKu3CcLVMEZDV6ngnsC0O/IiLSQeEI9+XAV+vPmpkFHHP3zwzJiIhI52lzWMbMHgfmA6lmVgT8M5AI4O6/BFYCi4B84DTwtUgVKyIioWkz3N39ljbmO/DtsFUkIiLnTFeoiogEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJAIYW7mS00sx1mlm9m9zUzf7iZvWFmm8xss5ktCn+pIiISqjbD3czigQeBq4BJwC1mNqlJs/8DPOnu04ElwC/CXaiIiIQulD33mUC+u+9y90rgCeC6Jm0c6Fv/uB+wL3wliohIe4US7hlAYaPnRfXTGnsA+IqZFQErgXua68jM7jazXDPLLS0t7UC5IiISilDC3ZqZ5k2e3wI84u6ZwCLgMTP7TN/uvtTdc9w9Jy0trf3ViohISEIJ9yIgq9HzTD477PJ14EkAd18L9ARSw1GgiIi0Xyjhvh7INrNRZpZE3QHT5U3a7AWuADCzidSFu8ZdRESipM1wd/dq4DvAy8B26s6K2WpmPzCzxfXNvgfcZWYfAI8Dd7h706EbERHpJAmhNHL3ldQdKG087fuNHm8D5oa3NBER6ShdoSoiEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCaCQwt3MFprZDjPLN7P7WmjzJTPbZmZbzeyP4S1TRETaI6GtBmYWDzwIXAkUAevNbLm7b2vUJhu4H5jr7kfMbHCkChYRkbaFsuc+E8h3913uXgk8AVzXpM1dwIPufgTA3Q+Gt0wREWmPUMI9Ayhs9Lyoflpj44BxZvaOma0zs4XNdWRmd5tZrpnllpaWdqxiERFpUyjhbs1M8ybPE4BsYD5wC/CwmfX/zELuS909x91z0tLS2luriIiEKJRwLwKyGj3PBPY10+YFd69y993ADurCXkREoiCUcF8PZJvZKDNLApYAy5u0eR64DMDMUqkbptkVzkJFRCR0bYa7u1cD3wFeBrYDT7r7VjP7gZktrm/2MlBmZtuAN4B73b0sUkWLiEjrzL3p8HnnyMnJ8dzc3KisW0SkuzKzDe6e01Y7XaEqIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQmghGgX0B1s3XeMLUXHOHC8gkvGpTJ9+IBolyQi0qqQwt3MFgI/B+KBh939xy20+wLwFHChu+eGrcooOVNVw7+99BG/e6fg7LSfvfoxM4b35+dLppM1sFf0ihMRaUWb4W5m8cCDwJVAEbDezJa7+7Ym7foA3wXejUShne10ZTW3LF3HB0XHuGPOSP5m7ij69UrkuY1F/OzVPJYsXcfjd81i+CAFvIh0PaGMuc8E8t19l7tXAk8A1zXT7l+AnwBnwlhfVLg79z+7hc3Fx3jo1hk8sHgywwf1ol9yInfMHcUf7ryIU5XVLFm6lmOnq6JdrojIZ4QS7hlAYaPnRfXTzjKz6UCWu69orSMzu9vMcs0st7S0tN3FdpZH1hTwwvv7+N6V47hq6tDPzJ+S0Y9HvjaTkuNn+PlreVGoUESkdaGEuzUzzc/ONIsDfgZ8r62O3H2pu+e4e05aWlroVXaikmPl/NtLH3HFhMF8a/7YFtudn9WfJRdmsWxtATtLT3ZegSIiIQgl3IuArEbPM4F9jZ73AaYAb5pZATALWG5mOeEqsjP95KUd1Do8sHgycXHNfa594nufH0/PxHh+9OL2TqpORCQ0oYT7eiDbzEaZWRKwBFjeMNPdj7l7qruPdPeRwDpgcXc8W+bdXWU8t6mYu+aNCulMmNTePbjn8rG8/tFB3vq46w4ziUjsaTPc3b0a+A7wMrAdeNLdt5rZD8xscaQL7Ew/fy2PQSlJfLOV4Zim7pg7khGDevHDFduorqmNYHUiIqEL6QpVd1/p7uPcfYy7/6h+2vfdfXkzbed3x732D4uPsWZnGXddMprePUK/tqtHQjz/tGgieQdP8sf39kawQhGR0On2A/V+9fYuevdI4MsXDW/3sp+fNITZowfx36/na+9dRLoEhTtQePg0K7eU8OWLhtO3Z2K7lzczbp8zgtITFazOPxSBCkVE2kfhDjy6pgADvjZ3ZIf7uGzCYPolJ/LcpuKw1SUi0lExH+4V1TU8s7GIBZPTGdovucP99EiI55rzhvLy1v2crKgOY4UiIu0X8+H+8tYDHDldxZKZWW03bsONMzI5U1XLX7aUhKEyEZGOi/lwf+K9vWQNTGbumNRz7mvG8P6MHNSLZzdqaEZEoiumw31P2SnW7Czj5pysNq9GDYWZccP0TNbtLqP4aHkYKhQR6ZiYDvenNxQRZ/CFC859SKbBDdMzcIfndWBVRKIoZsPd3Xn+/WLmjk0lvV/PsPU7fFAvLhw5gOc2FePubS8gIhIBMRvuG/ceofBwOdefn9F243a6cUYm+QdPsqX4WNj7FhEJRcyG+3ObiumZGMeCKelh73vR1KEkJcTpwKqIRE1MhntldS0vbi7hyknp7bqPTKj6JSdy5cQhLP9gH1W6HYGIREFMhvuqvFKOnK7i+vOHRWwdN0zP4PCpSt7aoVsBi0jni8lwf3FzCf2SE5mXHblvg7p0fBoDU5J0OwIRiYqYC/czVTW8su0ACyYPISkhci8/MT6OxdOG8dftBzhWri/RFpHOFXPh/vbHpZysqObq8yI3JNPgxhkZVFbXslK3IxCRThZz4f7ilhIG9EpkzphBEV/X1Ix+jElL4dmNRRFfl4hIYzEV7meqanh12wEWTkknMT7yL93MuHFGJusLjrC37HTE1yci0iCmwv3NHQc5VVnD1VMjPyTT4PrpdRdJ6cCqiHSmmAr3P28uYVBKErNGD+y0dWb0T2b26EE8t6lItyMQkU4TM+F+urKa17cfZOGUdBI6YUimsRtmZFBQdpqNe4926npFJHbFTLi/8VEp5VU1XH3e0E5f91VT0klKiOPPH+zr9HWLSGyKmXBfsXkfqb17cNGoyJ8l01SfnolcPn4wL24poaZWQzMiEnkxEe6nKqp5/aODLJqaTnwYvpSjI66dNozSExW8u7ssKusXkdgSE+H+6vYDVFTXck0nXLjUkssnDKZXUjx//kAXNIlI5MVEuL+4uYQhfXuQM2JA1GpITorncxOH8JcPS3SnSBGJuMCH+4kzVbz5cSmLpg4Ny/eknotrpw3j6OkqVucfimodIhJ8gQ/3V7cfoLK6lmuicJZMU5eMS6VPzwRWaGhGRCIs8OG+4oMShvXryfSs6A3JNOiREM/Cyem8snU/Z6pqol2OiARYoMP9WHkVb+d1jSGZBtdMG8aJimre+lhf4iEikRNSuJvZQjPbYWb5ZnZfM/P/wcy2mdlmM3vNzEaEv9T2e2XrfqpqnGumRe8smabmjBnEwJQkXdAkIhHVZribWTzwIHAVMAm4xcwmNWm2Cchx9/OAp4GfhLvQjnhxSwkZ/ZOZltkv2qWclRgfx1VT0nlt+0FOV1ZHuxwRCahQ9txnAvnuvsvdK4EngOsaN3D3N9y94Z6264DM8JbZfkdOVbI67xDXnDcUs64xJNPg2mnDKK+q4bXtB6NdiogEVCjhngEUNnpeVD+tJV8H/tLcDDO728xyzSy3tDSyY86vbNtPda1H9cKlllw4ciCD+/TQ0IyIREwo4d7cbm+zN0gxs68AOcBPm5vv7kvdPcfdc9LSIvfl1AArNpcwfGAvpmT0jeh6OiI+zrj6vKG8uaOU42f0/aoiEn6hhHsRkNXoeSbwmV1OM/sc8L+Bxe5eEZ7yOqbsZAVrdpZ1ySGZBtdOG0ZlTS1/3Xog2qWISACFEu7rgWwzG2VmScASYHnjBmY2HfgVdcEe9YHkl7ceoKbWo3J731BNz+pPRv9k/rxZQzMiEn5thru7VwPfAV4GtgNPuvtWM/uBmS2ub/ZToDfwlJm9b2bLW+iuU6zYvI/RqSlMGtr1hmQamBnXTBvK6rxDHDlVGe1yRCRgQjrP3d1Xuvs4dx/j7j+qn/Z9d19e//hz7j7E3c+v/7e49R4jp/REBet2lXF1Fx6SaXDtecOornVe2ro/2qWISMAE7grVlz4sodbpkmfJNDV5WF9Gp6borBnpMtydiuqaVv/pu4C7h4RoFxBO7s6TuUVkD+7NuCG9o11Om+qGZobx36/ncfD4GQb37RntkiQGuDtFR8rZVnKcwsOnKTpSTtGRhv/LOVnR+sV1fXokkDEgmcwBvcgckEzmgGRGDEph4tA+ZPRP7vJ/MceKQIX7+oIjbCk+xg+vn9Jt3mDXnjeU/3otj5VbSrhj7qholyMBVXj4NKvzD7E67xBrdh7iyOlPTsHt3SOhPqR7MWv0IFJ7J7X4++PulJ6oOPuBsHbnIU5VfnITvEEpScwZm8q8salcnJ3KsP7JEX9t0rxAhftvVu+if69EbpoR9QtkQ5Y9pA8T0vvw/Pv7FO4SNsfKq1i7s4zV+aWszjtEQVndBeTpfXtyxcQhzBg+gEnD+jJqUAp9kxM6vDPk7hwrr2LXoVNs23ecDXuOsDr/0NmhxtFpKfVBn8as0QPp0zMxbK9RWheYcN9bdppXth3gm5eOITkpPtrltMvNF2bx//68jdyCw+SMHBjtcqSbKq+s4aWtJTyzoZg1Ow9R65CSFM+s0YO4fc5I5mWnMiatd1j/qjUz+vdKYsbwJGYMH8BXZo3A3fn4wElW5ZWyOv8QT+YW8ejaPSTEGRdnp3LTjEyunDSEnond6/e0uwlMuP9uzW7izfjq7JHRLqXdbr4wi/96LY9fvrWThxXu0g7uznu7D/PMxiJWbtnPyYpqMgck8835Y7h03GCmD+9PYnznnjdhZoxP78P49D7cOW80FdU1bNxzlLc+LuWF94u55/FN9O2ZwDXThnHTjExmDO/fbYZRu5NAhPvxM1U8ub6Qa6cNI71f9zso2SspgTvmjOJnr37Mjv0nGJ/eJ9olSRdXePg0z2ws4tmNxew9fJqUpHgWTR3KTRdkMnPkwC7z/QVQ9yU1s8cMYvaYQdy7YDxrd5bV117EH9/dy+jUFG66IJMbpmdojD6MLFqnNeXk5Hhubm5Y+vr127v40crtrLjnYqZkdJ3b+7bHkVOVzP2311k4JZ3/+NL50S5HuqDaWuftvFIeXVPAm/Vf9jJnzCBumpHJwinp9ErqXvtqJyuqWbmlhKc3FPHe7sPEGVwxcQh3zBnJnDGDtDffAjPb4O45bbXrXu+GZlTX1PLImgJmjhrYbYMdYEBKEksuHM6ytQX8w5XjyBzQK9olSRdxpqqG5zcV8/Dq3eQfPElq7x7cc3k2N1+YRUY33tPt3SOBL+Vk8aWcLPaWneaJ9Xt5Yn0hf912gIlD+3LnxaO4dtowkhICdzlOp+j2e+4vbi7h23/cyNLbLuDzk9PDUFn07DtaziU/eYOvzBrBA4snR7scibKykxX8ft1eHltXwKGTlUwc2pe75o3imvOCG3hnqmp44f1iHl61m7yDJxncpwe3zxnJrRcNp3+vpGiX1yWEuufe7cP9hl+8w+FTlbz+vfnEd6Fxxo76x6c+YMXmfay57woGpujNHIsKD5/mV2/v5KncIiqqa7lsfBp3zRvN7BgaqnB33vq4lN+s3s2qvEMkJ8azZGYWd80bHfPj8jExLLNx7xE27T3KA9dOCkSwA3zj0tE8vaGIR9bUDc9I7Mg7cIKH3tzJCx/sI87gxumZ3HXJKMYOjr0D7GbG/PGDmT9+MNtLjvPrVbtYtnYPv1+3hxunZ/KN+WMYlZoS7TK7tG4d7r9ZvZs+PRP4Yk5W2427ibGD+/D5SUN4dE0Bf3vJaFJ6dOsfkYTg/cKj/OKNfF7ZdoDkxHjumDOSO+eNYmi/2N5DbTBxaF/+40vn8/efG8evV+3iT+sLeWpDIVdNHcq35o9h8rDue6wtkrptchQfLeelD/dz58WjAheA35g/hle2HeCJ9YV8/WJdtRpE7s7anWU8+GY+7+SX0bdnAt+9Ips75ozUcFwLsgb24gfXTeGey7P57Tu7eWztHl7cXMJl49P41mVjuVDXiHxKt03FR9cUAHD7nJFRrSMSZgwfwEWjBvLwql3cNmtEYA+exaLaWue1jw7y4Bv5vF94lLQ+Pbj/qgncOmsEvQO2kxIpaX168L8WTuAbl47hsbUF/PadAr74y7XMHDmQb102hkvHpcXMsYnWdMt308mKah5/dy9XTUkP7MGVb84fwx2/W88L7xcHatgpVlXX1LJicwkPvbmTHQdOkDkgmR9eP4UvXJCpy/A7qF9yIt+5PJu/uXgUf1pfyNK3d3HH79YzeVhfvn3ZWBZMTg/MsbiO6Jbh/lRuIScqqgM9ZHHpuDQmDu3LL9/ayU0zMrvUFYcSuhNnqnhuU92pfXsPn2bckN787OZpXHveMBI6+bYAQdUrKYGvzR3FrReN4PlNxTz01k6+9YeNjE5N4c55o7l++rBud4FXOHS7UyFrap3L/v1N0vr04JlvzolAZV3H8g/28d3HN/Gr2y5gQTc/hz/W5B04wbK1e3h2YxGnKmuYltWfb88fw+cmDtEHdYTV1DovfbifX7yZz9Z9x+tOurggi9tmjwjEGTaBPRVy675j7Dtazn1XTYh2KRG3aEo6Px2YzENv7uTzk4ZoHLGLq66p5a/bDrBs7R7W7iojKT6Oa6YN5auzR3J+Vv9olxcz4uOMq88byqKp6WzYc4RH1+5h2doCfvvObuZlp3L77JFcNmFw4Idsut2eO0DJsXLSeveIiT9rH1u3h//7/If8+qs5XDlpSLTLkWaUnqjgiff28sf39lJy7AwZ/ZO5ddZwbs7JYlDvHtEuT4CDx8/w+HuF/PG9PRw4XkHmgGS+MmsEN+dkMaCbnZ0UM1eoBt2Zqhpu+MUaCg+f5ulvzmZCet9olyTUncq4ce9Rlq0tYOWWEqpqnHnZqdw2awRXTBwS+L3C7qqq/q+rR9cU8O7uwyQlxLF42jC+OnsE52V2j7+uFO4BUnKsnOsffIeEuDie+/YcBvfpfrc1DoryyhqWf1DMsrV76sZzeyRw0wWZ3DZ7BGPSuv739sonduw/wWPrCnh2YzGn64+L3D57BIumDu3SZzAp3APmw+JjfPGXaxk3pDdP3D27233bVHe3p+wUv1+3hydzizhWXsX4IX24bfYIbpieEbiL6GLN8TNVPLuhiGXr9rCr9BQDU5K4+cIsbr1oeJe8O6vCPYD+uu0Adz+Wy4JJ6fzi1hk66yLCamvrbl61bG3d/dPjzFg4OZ2vzh7BzFEDdYA7YNydd/LLWLa2gFe3HwDq7i9/++yRzB3bdW7apnAPqN+s3s2/rNjG314ymvsXTYx2OYFzpqqG3IIjrMov5aUP97On7DRpfXpwy8zhfHnm8G75TV/SfsVHy/nDuj08sb6Qw6cqGZ2WwsLJ6czLTuOCEQOietW4wj2g3J3vv7CVx9bt4V9vnMotM4dHu6RurfGXOb+dd4h3d5VRUV1LYrxx4ciB3DJzOAsmp+sWEDHqTFUNK7eU8Kf1hWzYc4TqWqdXUjwXjRrIvOw0LhkX/i8db4vCPcCqa2q5c1kuq/IO8cjXLmRedlq0S+pWSk9U8E7+Id7OK2V13iEOnqgAYOzg3szLTuWS7DRmjhqosXT5lBNnqli36zCr6t83uw6dAmBov57My07l4uw0Lh6bGvEbvyncA+7EmSq++Mu1FB8pZ9nXZzJiUAqJ8UZifByJ8XE6Fa+Rs0Mt9Xvn20uOAzCgVyIXZ6cxb2wqF2enBvY+RRIZhYdPszr/0NmwP36mGjOYMqwf87JTmZedxowR/emREN6THxTuMaD4aN0pkqX1e56NxRkkxMeRFB9HQn3oN36cEGckJcQ1+zgxIY7EuPoPikaP6/r75HHjD5PGjxPijaRGj5vOb+5xQryRGBcXloPE7s6OAydY9fEhVuV/eqjlghED6v6czk5j8rC+OigtYVFT62wpPsaqj0tZlXeIjXvrhnCSE+OZNTq8QzgK9xhRePg0b+eVUlVdS3WtU1lTS3WNU1VT+6nHVWf/r5tWN69ueuPHDW2qapzqmloqa5zq2lqqqj9pG0kJcdboA+GzHwQJcXUfOEmNHjf+IKqprSW34IiGWiSqGoZwVufVhX3jIZyLx6byjwvGM6Rvxw7OK9wlItydmlo/+0Fy9kOl/v+qmtpPPT77wVJdS3Vt/YdF/fTGj5v/8AnxQ6q2lqrquum17kzN7K+hFulSGg/hvLf7CG/eO7/D9+8P643DzGwh8HMgHnjY3X/cZH4PYBlwAVAG3OzuBe0tOhTPbyrmpy/vYN/Rcob1T+beBeO5fnpGm/Ma5j+wfCtHy6uAuqGLWod4M2qa+ZBrmN7cfAMcyOifzGUT0njjo9JW13suNXd02Ugwq9uzToiny13F17A9Vnywj417jpCUENeh7RHO7dqRvlpapmF68dHys+/JjA6+n5rrA2j3e7FhmdZqaquutl5XOLdhR4Sjr6yBvUhOjOeDwmMcOlnBgp+9HfHf1zb33M0sHvgYuBIoAtYDt7j7tkZtvgWc5+7fMLMlwA3ufnNr/XZkz/35TcXc/+wWyqtqzk5LToznX2+cCtDivIY3yL1PfUBVbeT/Umm63nOpuaPLxprWtlV7tke4+uloXy0tc9MFGTyzofhT05v2Ce1/PzVIjDMwqKrxkJdtbplQl22t5rZed0e3YWf9DCNdU9iGZcxsNvCAuy+of34/gLv/a6M2L9e3WWtmCcB+IM1b6bwj4T73x69TfLT8M9Mz6v/0bmneO/dd3uKykdLWes+l5lCWjTWtbav2bI9w9dPRvlpapqW/Lhv3CeH/HYjksq3V3KCl192RbdhZP8NI1xTOYZkMoLDR8yLgopbauHu1mR0DBgGHmhR1N3A3wPDh7b/4Zl8Lb4KWpjee11qbSGhrvedScyjLxpqObKtI9tPRvlqa11qwh9rnubyGSCwbSp8tve6ObMPO+hlGuqZQhXLZXXPn7TTd4qG0wd2XunuOu+ekpbX/wpuWDo4N65/c6rzWlo2UttZ7LjWHsmysCdf2COd27UhfLc2Lb+P0uUj9DkRy2dZqbtDS6+7INuysn2GkawpVKOFeBDT+huZMYF9LbeqHZfoBh8NRYGP3LhhPcpODeMmJ8dy7YHyr8xqWTeykc5qbrvdcau7osrEmXNsjnNu1I321tMwtF2V9ZnrTPjvyfmpQdzrpp38/2lq2uWVCXba1mhu3ae51d3QbdtbPMNI1hSqUYZn1QLaZjQKKgSXAl5u0WQ7cDqwFvgC83tp4e0c1HHho7ch1S/Ma/u/ss2XCUXNHlo01oWyrzuyno321tkzOiIEhnZkSyvupvWfLtFRXa/21933cUh8Nrzsc27C9uuL7KlQhneduZouA/6TuVMjfuvuPzOwHQK67LzeznsBjwHTq9tiXuPuu1vrUee4iIu0X1vPc3X0lsLLJtO83enwG+GJ7ixQRkcjQfUxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBFNL93COyYrNSYE9UVh4+qTT5ntgYp+3xCW2LT9P2+MS5bosR7t7m95RGLdyDwMxyQ7lpfqzQ9viEtsWnaXt8orO2hYZlREQCSOEuIhJACvdzszTaBXQx2h6f0Lb4NG2PT3TKttCYu4hIAGnPXUQkgBTuIiIBpHDvADPLMrM3zGy7mW01s7+Ldk3RZmbxZrbJzFZEu5ZoM7P+Zva0mX1U/x6ZHe2aosXM/r7+d+RDM3vczHpGu6bOZGa/NbODZvZho2kDzeyvZpZX//+ASKxb4d4x1cD33H0iMAv4tplNinJN0fZ3wPZoF9FF/Bx4yd0nANOI0e1iZhnAd4Ecd58CxANLoltVp3sEWNhk2n3Aa+6eDbxW/zzsFO4d4O4l7r6x/vEJ6n55M6JbVfSYWSZwNfBwtGuJNjPrC1wC/AbA3Svd/Wh0q4qqBCDZzBKAXsC+KNfTqdz9beBwk8nXAY/zfrU7AAABo0lEQVTWP34UuD4S61a4nyMzGwlMB96NbiVR9Z/A/wRqo11IFzAaKAV+Vz9M9bCZpUS7qGhw92Lg34G9QAlwzN1fiW5VXcIQdy+Buh1FYHAkVqJwPwdm1ht4Bvgf7n482vVEg5ldAxx09w3RrqWLSABmAA+5+3TgFBH6s7urqx9Lvg4YBQwDUszsK9GtKnYo3DvIzBKpC/Y/uPuz0a4niuYCi82sAHgCuNzMfh/dkqKqCChy94a/5J6mLuxj0eeA3e5e6u5VwLPAnCjX1BUcMLOhAPX/H4zEShTuHWBmRt2Y6nZ3/49o1xNN7n6/u2e6+0jqDpa97u4xu3fm7vuBQjMbXz/pCmBbFEuKpr3ALDPrVf87cwUxenC5ieXA7fWPbwdeiMRKEiLRaQyYC9wGbDGz9+un/ZO7r4xiTdJ13AP8wcySgF3A16JcT1S4+7tm9jSwkbozzDYRY7chMLPHgflAqpkVAf8M/Bh40sy+Tt0H4Bcjsm7dfkBEJHg0LCMiEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIAP1/1MJzNjXBn3IAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "alpha1=0.5\n", "alpha2=0.5\n", "plt.scatter(data1+data2,[0]*len(data1+data2))\n", "plt.plot(data1+data2,alpha1*gaussian_1d(data1+data2,u1,sigma1)+alpha2*gaussian_1d(data1+data2,u2,sigma2))\n", "plt.ylim(-0.1,1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "到这里，我们通过两个高斯分布去模拟数据的概率分布，并取得了看起来不错的效果，但整个过程都是靠“人”的先验知识去完成的，比如（1）通过观测数据的分布，人为的选择了两个高斯分布；（2）一个高斯分布只去学习左侧数据，另一个高斯分布只去学习右侧的数据；（3）人为的为每个高斯分布的权重设置为0.5，但我们更加希望这个过程是能自动完成的，即自动去学习参数$\\theta=\\{u_1,\\sigma_1,u_2,\\sigma_2,\\alpha_1\\}$，所以让我们从头开始，假设，现在的高斯混合模型为： \n", "\n", "$$\n", "P(x\\mid\\theta)=\\alpha_1N_1(x\\mid u_1,\\sigma_1)+(1-\\alpha_1)N_2(x\\mid u_2,\\sigma_2)=\\sum_{i=1}^2\\alpha_iN_i(x\\mid u_i,\\sigma_i),\\alpha_2=1-\\alpha_1\n", "$$ \n", "\n", "对这些参数的求解，我们同样可以使用最大对数似然估计的方式，所以，损失函数： \n", "\n", "$$\n", "L(\\theta)=log(\\prod_{i=1}^MP(x_i\\mid \\theta))=\\sum_{i=1}^MlogP(x_i\\mid \\theta)=\\sum_{i=1}^Mlog[\\sum_{k=1}^2\\alpha_kN_k(x_i\\mid u_k,\\sigma_k)]\n", "$$ \n", "\n", "接下来，按之前的思路，我们分别对各个参数求偏导并令其为0，不是就可以解出最优参数了？可问题是现在很难解...,不信，那下面推导一下损失函数$L(\\theta)$对$u_k$的偏导（$k=1,2$）： \n", "\n", "$$\n", "\\frac{\\partial L(\\theta)}{\\partial u_k}=\\sum_{i=1}^M\\frac{\\frac{\\alpha_k(x_i-u_k)}{\\sqrt{2\\pi}\\sigma_k^3}exp(-\\frac{(x_i-u_k)^2}{2\\sigma_k^2})}{\\sum_{j=1}^2\\alpha_jN_j(x_i\\mid u_j,\\sigma_j)}=0 \\Rightarrow 谁能直接算出来，请联系我\n", "$$ \n", "\n", "所以，需要转换一种思路去求解，EM是一种迭代逐步优化的思路，它主要分为两步，第一步求期望(E)，第二步对其求极大化(M)，重复执行这两步，直到收敛，即找到一组$\\theta^1->\\theta^2->,...,->\\theta^i$使得$L(\\theta^1)