commit em algorithm examples

EM/EM.ipynb

@@ -0,0 +1,324 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 91,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import copy\n",
+    "import math\n",
+    "import matplotlib.pyplot as plt\n",
+    "# mathplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 创建数据集"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "一个样本中有两个高斯分布，均值不同，方差相同；\n",
+    "实例：从一个学校里随机抽取学生，抽取n个，其中包含a个女学生，b个男学生。"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 92,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "# n 样本数量 ，k 几个高斯分布，sigma：方差，real_mu 包含两个分布的真实期望\n",
+    "\n",
+    "def create_data(n,k,real_mu,sigma):\n",
+    "    X = np.zeros((n,1))\n",
+    "    for i in range(0,n):\n",
+    "#   随机抽取 男女生\n",
+    "        if np.random.random()>0.5:\n",
+    "            X[i][0] = np.random.normal()*sigma + real_mu[0]\n",
+    "        else:\n",
+    "            X[i][0] = np.random.normal()*sigma + real_mu[1]\n",
+    "    return X"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 高斯分布"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 93,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "\n",
+    "def gaussian_dist(x,mu,sigma):\n",
+    "    exponent = np.exp(-np.power((x-mu),2)/(2*np.power(sigma,2)))\n",
+    "    result = 1/np.power(2*np.pi*sigma**2 , 1/2) *exponent\n",
+    "    return result"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Expectation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "固定$\\theta$ ,优化条件期望"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "简单而言：已知 $\\theta$ ,求条件期望，即在某个样本结果出现的条件下，是男生或者女生的概率"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 94,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "def e_step(X,mu,sigma,expectations):\n",
+    "    length = X.shape[0]\n",
+    "    k = mu.shape[0]\n",
+    "\n",
+    "    for i in range(0,length):\n",
+    "        denominator = 0;\n",
+    "        for j in range(0,k):\n",
+    "            denominator += gaussian_dist(X[i][0],mu[j],sigma)\n",
+    "        \n",
+    "        for j in range(0,k):      \n",
+    "            numberator = gaussian_dist(X[i][0],mu[j],sigma)         \n",
+    "            expectations[i][j]= numberator / denominator\n",
+    "            \n",
+    "    return expectations"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Maximization"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "已知条件期望，求解$\\theta$在此次迭代中的最大值"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 95,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "#  这种方法更直观，\n",
+    "# def m_step(X,mu,sigma,expectations):\n",
+    "#     n = X.shape[0]\n",
+    "#     k = mu.shape[0] \n",
+    "#     for i in range(0,k):   \n",
+    "#         numberator = 0\n",
+    "#         denominator = 0\n",
+    "#         for j in range(0,n):\n",
+    "#             numberator += expectations[j,i]*X[j,0]\n",
+    "#             denominator += expectations[j,i]\n",
+    "        \n",
+    "#         mu[i] = numberator/denominator\n",
+    "    \n",
+    "#     return mu"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 96,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "# 向量法，简单高效\n",
+    "def m_step(X,mu,sigma,expectations):\n",
+    "    numbers = np.transpose(X).dot(expectations)\n",
+    "    denominator= np.sum(expectations,axis=0)\n",
+    "\n",
+    "    result = numbers/denominator\n",
+    "#     result是1*2的矩阵，\n",
+    "# 我们需要返回一个一维向量，包含两个分量而已。\n",
+    "    mu[0] = result[0][0]\n",
+    "    mu[1] = result[0][1]\n",
+    "    return mu"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# 训练"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 97,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "def train(X,sigma,expectations,iteration,epsilon):\n",
+    "    mu = np.random.random(2)\n",
+    "    for i in range(0,iteration):\n",
+    "        # deepcopy 使得两个值相互不影响\n",
+    "        old_mu = copy.deepcopy(mu)\n",
+    "        # E步：求出条件期望\n",
+    "        expectations = e_step(X,mu,sigma,expectations)\n",
+    "        \n",
+    "        print(i,mu)\n",
+    "        # M步，利用E步的期望，求最大的theta（）\n",
+    "        mu = m_step(X,mu,sigma,expectations)\n",
+    "        # 如果两次迭代的差值小于阈值，则训练结束，也可以认为找到局部最优值，\n",
+    "        if np.abs(np.sum(old_mu) - np.sum(mu)) < epsilon:\n",
+    "            break"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 98,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0 [0.71085143 0.97453594]\n",
+      "1 [168.60550916 169.62356611]\n",
+      "2 [167.44608098 171.32745846]\n",
+      "3 [163.31669663 175.570302  ]\n",
+      "4 [160.18984499 179.48883559]\n",
+      "5 [160.04240354 180.32256705]\n",
+      "6 [160.14223833 180.53002579]\n",
+      "7 [160.19102044 180.60482685]\n",
+      "8 [160.21080304 180.63380363]\n",
+      "9 [160.21861444 180.64513455]\n",
+      "10 [160.22168114 180.64957042]\n",
+      "11 [160.22288306 180.65130723]\n",
+      "12 [160.22335384 180.65198726]\n",
+      "13 [160.22353819 180.65225352]\n",
+      "14 [160.22361038 180.65235778]\n",
+      "15 [160.22363865 180.6523986 ]\n",
+      "16 [160.22364971 180.65241458]\n",
+      "17 [160.22365405 180.65242084]\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADYhJREFUeJzt3X+MZXdZx/H3hy5FKyCtO8WGdtxK\noEljIpSRVBHQFgSKofgD0yYgCmYjkdo2IilBgcR/KiKKidEstFItKQg0gopKrVRiAgvdpdCWpZQf\nBbYtLYQI/EXBPv5xz2aHzcx2554zs3Mf3q9kcs895+zc59nvvZ8599wf31QVkqReHna8C5AkTc9w\nl6SGDHdJashwl6SGDHdJashwl6SGDHdJashwl6SGDHdJamjHVt7Yzp07a9euXVt5k5K08Pbt2/f1\nqlrayL/Z0nDftWsXN99881bepCQtvCRf2ui/8bSMJDVkuEtSQ4a7JDVkuEtSQ4a7JDX0kOGe5Ook\n9ye5bdW6U5LckOTO4fLkzS1TkrQRx3Lk/nbguUesuwK4saqeANw4XJckbRMPGe5V9WHgG0esvhC4\nZli+BnjhxHVJkkaY95z7Y6vqXoDh8tTpSpIkjbXpn1BNshvYDbC8vLzZN6cFt+uKf11z/V1XPn+L\nK5EW27xH7vclOQ1guLx/vR2rak9VrVTVytLShr4aQZI0p3nD/f3AS4fllwLvm6YcSdIUjuWtkNcB\nHwHOSnIwycuBK4FnJ7kTePZwXZK0TTzkOfequnidTedPXIskaSJ+QlWSGjLcJakhw12SGjLcJakh\nw12SGjLcJakhw12SGjLcJakhw12SGjLcJakhw12SGjLcJakhw12SGjLcJakhw12SGtr0OVS1eJzH\ndHubanwc5948cpekhgx3SWrIcJekhgx3SWrIcJekhgx3SWrIcJekhgx3SWrIcJekhgx3SWrIcJek\nhgx3SWrIcJekhgx3SWrIcJekhgx3SWpoVLgnuTzJ7UluS3Jdkh+aqjBJ0vzmDvckjwN+H1ipqp8C\nTgAumqowSdL8xp6W2QH8cJIdwEnAPeNLkiSNNXe4V9XdwJuALwP3At+sqg9OVZgkaX5zT5Cd5GTg\nQuBM4H+Bdyd5cVVde8R+u4HdAMvLyyNK1UPpPOGxk0JLGzPmtMyzgC9W1deq6rvA9cDPHblTVe2p\nqpWqWllaWhpxc5KkYzUm3L8MnJvkpCQBzgcOTFOWJGmMMefc9wLvAfYDtw6/a89EdUmSRpj7nDtA\nVb0eeP1EtUiSJuInVCWpIcNdkhoy3CWpIcNdkhoy3CWpIcNdkhoy3CWpIcNdkhoy3CWpIcNdkhoy\n3CWpIcNdkhoy3CWpIcNdkhoy3CWpIcNdkhoaNVmHFsN6k0Ifz9t2Qur+jna/c/w3n0fuktSQ4S5J\nDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnuktSQ4S5JDRnu\nktSQ4S5JDY0K9ySPSfKeJJ9JciDJz05VmCRpfmNnYnoL8O9V9etJTgROmqAmSdJIc4d7kkcDzwB+\nC6CqHgAemKYsSdIYY07L/CTwNeDvknwiyduS/MhEdUmSRhhzWmYHcA5wSVXtTfIW4Argj1fvlGQ3\nsBtgeXl5xM394DmeE1uvZcoJj7dbb1uha89d+1p0Y47cDwIHq2rvcP09zML++1TVnqpaqaqVpaWl\nETcnSTpWc4d7VX0V+EqSs4ZV5wOfnqQqSdIoY98tcwnwjuGdMl8Afnt8SZKksUaFe1XdAqxMVIsk\naSJ+QlWSGjLcJakhw12SGjLcJakhw12SGjLcJakhw12SGjLcJakhw12SGjLcJakhw12SGjLcJakh\nw12SGjLcJakhw12SGho7WYcEOI/mdrZIY7NerRudo1ceuUtSS4a7JDVkuEtSQ4a7JDVkuEtSQ4a7\nJDVkuEtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ4a7JDVkuEtSQ4a7JDU0OtyTnJDkE0n+\nZYqCJEnjTXHkfilwYILfI0mayKhwT3I68HzgbdOUI0mawtgj978EXg08OEEtkqSJzD1BdpJfBu6v\nqn1JfuEo++0GdgMsLy/Pe3OtLdIExtuNEyofttn3o0W6n3q/GHfk/jTgBUnuAt4JnJfk2iN3qqo9\nVbVSVStLS0sjbk6SdKzmDveqek1VnV5Vu4CLgP+qqhdPVpkkaW6+z12SGpr7nPtqVXUTcNMUv0uS\nNJ5H7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z\n7pLUkOEuSQ0Z7pLUkOEuSQ1NMlmHvt8iTSTc1UbHYKMTKjvG4/j/t/k8cpekhgx3SWrIcJekhgx3\nSWrIcJekhgx3SWrIcJekhgx3SWrIcJekhgx3SWrIcJekhgx3SWrIcJekhgx3SWrIcJekhgx3SWpo\n7nBPckaSDyU5kOT2JJdOWZgkaX5jZmL6HvAHVbU/yaOAfUluqKpPT1SbJGlOcx+5V9W9VbV/WP42\ncAB43FSFSZLmN8kcqkl2AU8G9q6xbTewG2B5eXmKm5O0ibbj/KYbneNWE7ygmuSRwHuBy6rqW0du\nr6o9VbVSVStLS0tjb06SdAxGhXuShzML9ndU1fXTlCRJGmvMu2UCXAUcqKo3T1eSJGmsMUfuTwNe\nApyX5Jbh54KJ6pIkjTD3C6pV9T9AJqxFkjQRP6EqSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEu\nSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ0Z7pLUkOEuSQ1NMkH2VjheE+QebbJgJ+eVtNp2\nmsjbI3dJashwl6SGDHdJashwl6SGDHdJashwl6SGDHdJashwl6SGDHdJashwl6SGDHdJashwl6SG\nDHdJashwl6SGDHdJashwl6SGRoV7kucmuSPJ55JcMVVRkqRx5g73JCcAfw08DzgbuDjJ2VMVJkma\n35gj96cCn6uqL1TVA8A7gQunKUuSNMaYcH8c8JVV1w8O6yRJx1mqar5/mLwIeE5V/c5w/SXAU6vq\nkiP22w3sHq6eBdyxavNO4OtzFbD9de2ta1/Qt7eufcEPTm8/UVVLG/nHO0bc8EHgjFXXTwfuOXKn\nqtoD7FnrFyS5uapWRtSwbXXtrWtf0Le3rn2BvR3NmNMyHweekOTMJCcCFwHvH/H7JEkTmfvIvaq+\nl+SVwH8AJwBXV9Xtk1UmSZrbmNMyVNUHgA+M+BVrnq5pomtvXfuCvr117QvsbV1zv6AqSdq+/PoB\nSWpoU8M9ydVJ7k9y2xrbXpWkkuwcrifJXw1fZfCpJOdsZm1jrNVXkjckuTvJLcPPBau2vWbo644k\nzzk+VR+b9cYsySVD/bcneeOq9QvR2zpj9q5V43VXkltWbVuIvmDd3p6U5KNDbzcneeqwfmEeZ7Bu\nbz+d5CNJbk3yz0kevWrbQoxbkjOSfCjJgeExdemw/pQkNyS5c7g8eVi/8XGrqk37AZ4BnAPcdsT6\nM5i9EPslYOew7gLg34AA5wJ7N7O2qfsC3gC8ao19zwY+CTwCOBP4PHDC8e5hg739IvCfwCOG66cu\nWm/r3RdXbf9z4HWL1tdRxuyDwPOG5QuAm1YtL8Tj7Ci9fRx45rD8MuBPFm3cgNOAc4blRwGfHep/\nI3DFsP4K4E/nHbdNPXKvqg8D31hj018ArwZWn/C/EPj7mvko8Jgkp21mffM6Sl9ruRB4Z1V9p6q+\nCHyO2Vc3bEvr9PYK4Mqq+s6wz/3D+oXp7WhjliTAbwDXDasWpi9Yt7cCDh3R/iiHP4OyMI8zWLe3\ns4APD8s3AL82LC/MuFXVvVW1f1j+NnCA2Sf8LwSuGXa7BnjhsLzhcdvyc+5JXgDcXVWfPGJTh68z\neOXwlOnqQ0+n6NHXE4GnJ9mb5L+T/MywvkNvAE8H7quqO4frHfq6DPizJF8B3gS8ZljfobfbgBcM\nyy/i8IcpF7K3JLuAJwN7gcdW1b0w+wMAnDrstuHetjTck5wEvBZ43Vqb11i3SG/l+Rvg8cCTgHuZ\nPc2Hxe8LZm+ZPZnZ08E/BP5xONrt0BvAxRw+aocefb0CuLyqzgAuB64a1nfo7WXA7yXZx+yUxgPD\n+oXrLckjgfcCl1XVt4626xrrjtrbVh+5P57ZubBPJrmL2VcW7E/y4xzj1xlsV1V1X1X9X1U9CLyV\nw08HF7qvwUHg+uEp4ceAB5l978XC95ZkB/CrwLtWrV74voCXAtcPy++m0f2xqj5TVb9UVU9h9kf5\n88OmheotycOZBfs7qurQWN136HTLcHnoFOiGe9vScK+qW6vq1KraVVW7mBV8TlV9ldlXF/zm8Krw\nucA3Dz09WQRHnP/6FWZPHWHW10VJHpHkTOAJwMe2ur6R/gk4DyDJE4ETmX2hUYfengV8pqoOrlrX\noa97gGcOy+cBh045LfTjDCDJqcPlw4A/Av522LQw4zY8870KOFBVb1616f3M/jAzXL5v1fqNjdsm\nvyJ8HbNTFN9lFuQvP2L7XRx+t0yYTf7xeeBWYGWrX8Ee0xfwD0PdnxoG4rRV+7926OsOhncwbNef\ndXo7EbiW2R+s/cB5i9bbevdF4O3A766x/0L0dZQx+3lgH7N3j+wFnjLsuzCPs6P0dimzd5d8FriS\n4cOYizRuw/jUkBe3DD8XAD8G3Mjsj/GNwCnzjpufUJWkhvyEqiQ1ZLhLUkOGuyQ1ZLhLUkOGuyQ1\nZLhLUkOGuyQ1ZLhLUkP/DzazO+/xdYuZAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<matplotlib.figure.Figure at 0x1b011947b38>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# 200个抽样样本，（男女生集合）\n",
+    "n=200\n",
+    "k=2\n",
+    "# 真实期望\n",
+    "real_mu = np.array([180,160])\n",
+    "sigma=6\n",
+    "iteration=1000\n",
+    "epsilon=0.00001\n",
+    "expectations = np.zeros((n,k))\n",
+    "X = create_data(n,k,real_mu,sigma)\n",
+    "train(X,sigma,expectations,iteration, epsilon)\n",
+    "plt.hist(X,50)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.3"
+  },
+  "toc": {
+   "nav_menu": {},
+   "number_sections": true,
+   "sideBar": true,
+   "skip_h1_title": false,
+   "title_cell": "Table of Contents",
+   "title_sidebar": "Contents",
+   "toc_cell": false,
+   "toc_position": {
+    "height": "calc(100% - 180px)",
+    "left": "10px",
+    "top": "150px",
+    "width": "249px"
+   },
+   "toc_section_display": true,
+   "toc_window_display": true
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}