linear_discriminant_analysis(BiClass).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 线性判别分析"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "线性判别分析（Linear Discriminant Analysis或者Fisher’s Linear Discriminant）简称LDA，是一种监督学习算法。\n",
    "\n",
    "LDA的原理是，将数据通过线性变换(投影)的方法，映射到维度更低的空间中，使得投影后的点满足同类型标签的样本在映射后的空间比较近，不同类型标签的样本在映射后的空间比较远。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 一、线性判别分析(二类情形)\n",
    "\n",
    "\n",
    "在讲解算法理论之前，先补充一下协方差矩阵的定义。\n",
    "\n",
    "\n",
    "### 1. 协方差矩阵定义\n",
    "\n",
    "矩阵$X_{mxn}$协方差的计算公式：\n",
    "\n",
    "设$x, y$分别是两个列向量，则$x, y$的协方差为\n",
    "$$cov(x, y) = E(x-\\bar{x})(y-\\bar{y})$$\n",
    "\n",
    "若将$x, y$合并成一个矩阵$X_{mxn}$，则求矩阵$X_{mxn}$，则矩阵$X_{mxn}$的协方差矩阵为\n",
    "\n",
    "$$A = \\sum_{i=1...m,j=1...n}a_{ij} = cov(X_i, X_j) = E(X_i-\\bar{X_i})(X_j-\\bar{X_j})$$\n",
    "\n",
    "\n",
    "### 2. 模型原理"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/jpeg": "R0lGODlhSAE+AfcAAAkBAg8HHA4UGxIKCBEKFhwTCRoVFR44GRQhISEaFyUbJiQhGjEoFCYqJiEw\nMTUlMjU4NTUvQipGJUMkK0c5Om49O0U2RlhJOkdESU9JUElQUFdDSF1OXFVXTFtbW1xbZ1thV2Rb\nWmleanNdYX5hXGNnZ2tteGVxamhwcHlpaHZ6a3ZwcAp05wF39Bd26RR99HZqmHl4hWh7zwKD7wGA\n8BmH6RWG8iaG7iiL8iWT7iyR8TaJ6TaZ6zmT8T2j6nWAglmW3USM5UWW6kSX9lWf70um6kCi8Uuw\n51el6lil8Faz6Wud1H2h3Weh6GKn82u16mey8nKn6HW47HG48XrH55NUVIRwcKRsbcolJuAVCPYB\nAvIKEvMVC/QYF/AmF/MiJPI2KvQ3NvI5RPJPP9VIQMtvb8Zyd+dPUfNIR/JIVvJXSvJWU/JHavJo\nWu96afNnafNpcvF5a/J5d4B5hrJ4gM9vh8Z7geB7iPN6ifF8kIaMf5iBeMWOePaJfYiEiYWQgpSI\nhpSOkpabjZaWmJiboZailJ+ipp20rbmFhbmQnb6Vq7WbtaOknKenq6y4tbuotrG1qbO1uIqr54W6\n6pG76JO/96iuxqe5w4nI8o7R6JjH6ZXH9pnV7LzMzqfN6qTH8qrV6avV9LLN7LfY6LnX8seHicaK\nlcqWhsmWmN6FhtiLmNSQitSbmMaDosqYo8aasNWbpsqmndeoncumqMipscy2tNinp9aqtNuxqde1\nuO+Fj/SCh/OOlvGVhfWZmOWbofWZqeummPaim+qopeqrtuCwr+O4tvWlpfaqsfK0qvS4tcm9wM27\n1d27xNG+2Oe9w/K6xcPButnBuujGufXDvMfFx8zD1cjW2t3GydLG1NrRytjY2MPP6sXW7sPf9d/P\n4Nrb5N3g3Mvj8tnm59fq9ebJxuXL0e/XyOfW1vbIyPfP0/LXzPbW1ufb5PXb5eji3O73wvXi3Orp\n6OXs9Orz9vfo5/bt8fjy7f7+/v/M////AP//M///Zv//mf//zAAAACH5BAEAAP8ALAAAAABIAT4B\nAAj/APEJHEiwoMGDCBMqXMiwocOHECNKnEixosWLGDNq3Mixo8ePIEOKHEmypMmTKFOqXMmypcuX\nMGPKnEmzps2bOHPq3Mmzp8+fQIMKHUq0qNGjSJMqXcq0qdOnUKNKnUq1qtWrWLNq3cq1q9evYMOK\nHUu2rNmzaNOqXcv2q7x2MfzIa0v3bCQPDwAMUEGvrt+wjRoAEEABwAJIfxNzbaQAQAFC1QQAKDFX\nsWWrBwZH2IbPhOFol0NLJaQXkKFx+MARAHCir+jXS+U1ptCuUDiBKwAk4Ay791FAAwjowTeIt6MF\nABr5Xj5UHggADVAXF7gNA4AQzLP/jCS4g0BC4AbO/wFgQbt5nfRKAAhwTeB0fNsKAzh0vr7Ndn/k\n1XPPe1sHAxOgZt+AN72HTzjWCEjggjMZyOCDNDkI4YQv2UbhhS9JiOGGJmnI4YcheQjiiByJSOKJ\nF1mI4ooamcjiiw65COOMCQ1yG404xnhjjjzWyFuPQBYkY5AzDknki4LseCSQRi6JYpNOkghllCCq\nSGWOU165oZVaFvljlzBmCeaEXI65ophmMohmmgSuyaZ9Nr55ZnhynuhmndoF8iWeHBKiIJ8c3gmo\nb3EOymGZhk4oKEv1KANNOum4k2hRi67UixaYaiHHpEQh2hMwmWoRB6dDVarSpZmOSmpQpqbkS6iq\nrv/6U6souXMMMscck46sQNHKK1u+/qpWksIuGGyxZxWKbH3HLkuWn84yu2e0yzVLLVjWXuuVINNm\nlI4cccDxhjDaZnWsMF1kuka5WHmKkTBcZIoGu1cdC8wW6g7FDqTpsENvjN1ilAwaYRS8aVD1pJHv\nvwwpy5E97ETMzn4Io7EwwwrpeVY9Y1yMMULZFlWPGB5/bFDIRcVTTjnnUGzyQQ6//FfMMtc1CJ01\nzxxwzmqhzHNT0P5s885Cm0Vz0Wj5jDRS7i5NltJOF9V01GFBXVE97rDjjstUl2Q1Rb6kEUYacMTT\ndYdEy/RGpmGYfTZJX08kB9v+vj2SxjutjWkYddv/HRKxOx2jBhtvxMG13x7FjXhOii9+E+COf9V4\n5BGmPWI3mmQ+ypKT+zVEC6An4VqQR69oBOg0OMG55SCeDrrqR3ZelxE11DCDFKMDCTnlXMnOe4as\n/w4Vt8L3HnzxTfmOvEp4L3/V7g3Vow6/hzvPU5b1bLFFF11QY32vxwtUT6jef+8T9AyNn6ky5vuU\n5Tqh7oo4Os0kJWYyyiiTjNt+18NKOUhRHr3aYYt2HKV0QVkHHnYhh10cI03lKMUBw5eTZIRKDWyK\nxS2MIsCLWDBTb3iTKTZIFPQN5YOYSsOb4IGK6rmPgjhRBve41wY5sSMVvejbrGCIE3fEY2su7NIx\n/9K1hiDmhHjta4jCMKXDnnSQYcr4ghYwyCoeOo8d09AGKoLSvCQ6pAzoAIoJvagQVADwhWSESDlY\noQxfGLFBVixe9YYBhi684Y0xGWMaEQKHTK1jJ0+U2dwwBQ1AxtF86UBDF85wvUNeBQk7iGQN5nEi\ndiQjF7YApJLSwoMW0AB0fyKRGcLIOEdaxQeoa0EoR1QPV+AiDrwokCmr0knQqfJF7PACpvJgEz2S\nxRNUCCYVclfJTIHBJoH8mT3w0IUvJAOZs0xiPYRxBzyuBIl7rIgtjAHNbFqEHrCQHxy9eZEqWBMl\nXSSnRM7BinScsyS+XJIUIsmDKFByK9QIQxfk8P/OkSTzNUmwpRDIwZU1YKoL6ojJP0WThBmAbqBc\nCQamsiBOl8TzSEgQKEG5Qgw5qIIYCo3miThBBSQ8QRNhQcU0YLLQqKnjFAZ0SUujpg1T8G8lCFSn\nRDApU5Hq1Bi5aEmzjoGHBT7TfOhgxU1RctGJ1EMNJXNeOVqIU58O5Kny8iIsuKmSY9ljiVqYlzRR\nkQ2VYJMj93gDGL4AhhCqhB7kiOs9l4UOO/SzRVYdyL6U0a+VTKIHgG0CtdBB1ZPM1CJSsCUPrjWL\n+t3VIulcSkZBh4Rr1UMVFqPiSJo6lMQ+VFuvwpQv4JZXmWjiCajlhLaiqIUvIIO0OhXJMdCgC3j/\nbjK2HanHKeyxjsdG5LCO2wYWthCGQoIEuItjR6ba4FuA4VYk7lBXcxuS0+diZBdc6MJoj1ta6w5E\nGb+AxXQb1l3vDgQXK/0IZ81bETNUpiPI5V0zUHEM424kvpSzBxla+8CNnJW9HkkHvrRQxPuWF8D4\nWGIfSnRgBLvhDpnEx3gJElkEe+QZEc7IeqUCig6DYnPIkocrSIkR/IZkErZsQQ2chQ4JZsTEIJns\nJ1fsrFxsESMbhsoTWuBQFVNrGje2CIw/ooQU0zhapXDkkD3CiSII4QhHAHG02MGKJv62wRYuSDky\nPJHqZpkj9chFUCmyZOTZoxTvvfKXUbKNUzQ3/8drvkgpQCqRMjtPpXXGcpwJcuaYPuS/ex5JOWaR\nZvKuJBRSmIQUPgHgbUIkbvSYbAuGAGB5FLYhFR4JPZJAgxnMwAgIrkcprIwQOGOEHlD4JA0ojeBb\n4OIe1NWzQaTQAx30AApZRsWYF6I4cohDHHFd2oQtImISJ8TUgc5IOy4NMlkn2yC1IKGPns0SW+wa\nZs6mNkEIu1SCePkppKhEJTaxCXEsDR2rmLZWdmxLlFJp2BvRYLO3wu5PuntpfJDGQZB9FCqk+N5I\nc4cp7HGybM/kCS+wgQ1e4AmqTfVwdr7IOCa+SqTRgqsDAbS2U1IPW2DD2wbfuEEiCPK3gcMZjv+Q\nFFekEeRMC+0RjQGABHCmlTkLhN/sckQDDDAAAGyg4lWJByrCGPEZbSMSKgAAADQQCT9jJRstLDqO\nflAA3XCgPVopRy6kjqNGfMAAAIBAoWOyjbKb/exoT7va1Y6ObdBhBJ1QOzPMPveyZyMbdsd72eu+\nDb7rfRt3p7vg7W72wO998IAv/N+tgfi/rx3tfzf8Nhh/eMLn3eyU73vhFY/5xnO+7Jn3++cn7/nL\ng770gI884tGODnCUfRzWEEHPge4SQqDACii4fe5XsILc4x4FvPe98IHfe90TP/cnSEEFSKD74OM+\n+M3v/e+df/zoG5/60H++9He//ewff/rdL/7/9Ief+/Jff/vAT3/zhY/98F9f/dq3fvzHX/zei//7\n+Kd++u/ffvnv/v38d38lkALcdwGzxxV7YAqREwkJAAC0NxWBgAgftziQAHYhJxOFoA2lMIF+w4AA\ncIEhhQ6uUA/EFDXbIAJK94BSYSHmAAM3gAQlWDQrAHZKB4IsxRnkwAKg8wRRUw0MoHQD4ACE0BXT\nQQ82ADpSsDSNEAIIoBslYAhjhxXv4QlFsAQqiDHgkAFKhwGMUDU/Mg/NAAv4AFc/Qw9+gAFc92J7\n4grMoAQ4sAk/Uw9RmFs7ZBBBADo3YG4ihxIaUmQtsAMbtYdocxCQFAqCyIfdMg+PYGyH6E8B/4MO\nswBvjRhrCbFNowAK2SSJJ+FLkgA6kzCJIbIz81BLOKCHoAhfRAMFDgUEc3WKBqYQ84AJUZAIrpg4\nljMPxlA/tehft3UQ6KAK3baLFBE0DKENtnANSNBwwlgRUGIJLtACL2CIy9hlx8MJR9gCUzCN1PgQ\nkdYCEKWNEbFh9MANiwCOeSYR0yCG5vgQ31YQ8sAKCbWO1NWLDGEOrDCH8lhyE5ELxjAOn+AN+Thv\nTkULQvACPUAKAVlwFaEJqAM7CTkQU7MQ45ADnoQJD8kR4FAETHCRHsGPHAlmqOB0H3kRLGcNMTiS\nEdEOMuACR4CSF7EJM/BJmXBVLhkRoXCNlP9QkwtJBJNwC7CmkxNBSYpwCcAGlBMhBC6AAzkgCkYJ\nEX4IOjvQlA+BYp7UAkUglQ8RBTzQA0BAc1i5EPRAD+8AC/j4lQeRCwpolg8RDNKmlgxhCtfmlgmR\nVGUplwIhDVxmlwdhD7mQXnqJEPWACJo4ks/gCn+pELXAB4eZEKVgDot5EO/ACoz4mPhgDiNImQWR\nC8WAmQTRDq8wmB+5DnbAmQRxC8vwk6RpC3lJmfdgCpP5mP6jcqRpY6QpEPbgcbWJD+UgC3Vpl+UQ\nC7mJD7oWnKhgX6R5BcGJYQRXm2KWm/UgCzZ4is1gmLkpC8EQnHVwRqQZD63wDrn5DW6Wm7kIgAvB\n2ZykEhAAOw==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import Image, display\n",
    "\n",
    "display(Image('lda.jpg'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上图中红色的正方形形状的点为0类的原样本点、蓝色的正方形形状的点为1类的原样本点，过原点的那条直线就是投影的直线，从图上可以清楚的看到，红色的点和蓝色的点被原点明显的分开了。下面具体看看二分类LDA二类问题的情形：\n",
    "\n",
    "现在我们觉得原始特征数太多，想将$n$维特征降到只有一维(LDA映射到的低维空间维度小于等于$n_{labels}-1$)，而又要保证类别能够“清晰”地反映在低维数据上，也就是这一维就能决定每个样例的类别。\n",
    "\n",
    "假设用来区分二分类的直线（投影函数)为:\n",
    "\n",
    "$$y=w^Tx$$\n",
    "\n",
    "注意这里得到的$y$值不是0/1值，而是$x$投影到直线上的点到原点的距离。\n",
    "\n",
    "已知数据集\n",
    "\n",
    "$$D=\\{x^{(1)}, x^{(1)}, …, x^{(m)}\\},$$\n",
    "\n",
    "将$D$按照类别标签划分为两类$D_1, D_2$, 其中$D_1\\bigcup D_2=D, D_1\\bigcap D_2=\\emptyset$, \n",
    "定义两个子集的中心：\n",
    "\n",
    "$$\\mu_1 = \\frac{1}{n_1} \\sum_{x^{(i)} \\in D_1} {x^{(i)}},$$\n",
    "$$\\mu_2 = \\frac{1}{n_2} \\sum_{x^{(i)} \\in D_2} {x^{(i)}},$$\n",
    "\n",
    "则两个子集投影后的中心为\n",
    "$$\\tilde{\\mu_1} = \\frac{1}{n_1} \\sum_{y^{(i)} \\in \\tilde{D_1}} w^T {x^{(i)}},$$\n",
    "$$\\tilde{\\mu_2} = \\frac{1}{n_2} \\sum_{x^{(i)} \\in D_2} w^T {x^{(i)}},$$\n",
    "\n",
    "则两个子集投影后的方差分别为\n",
    "$$ \\tilde{ \\sigma}^2_1 = \\frac{1}{n_1} \\sum_{y^{(i)} \\in \\tilde{D_1}} ({y^{(i)}-\\tilde{\\mu_1}})^2 = \n",
    "\\frac{1}{n_1} \\sum_{x^{(i)} \\in {D_1}} ({w^T x^{(i)}-w^T \\mu_1})^2 = \n",
    "\\frac{1}{n_1} \\sum_{x^{(i)} \\in {D_1}} w^T  ({x^{(i)}- \\mu_1})  ({x^{(i)}- \\mu_1})^T w  ,$$\n",
    "\n",
    "同理可得\n",
    "$$ \\tilde{ \\sigma}^2_2 = \\frac{1}{n_2} \\sum_{y^{(i)} \\in \\tilde{D_2}} ({y^{(i)}-\\tilde{\\mu_2}})^2 =\n",
    "\\frac{1}{n_2} \\sum_{x^{(i)} \\in {D_2}} w^T  ({x^{(i)}- \\mu_2})  ({x^{(i)}- \\mu_2})^T w  , $$\n",
    "\n",
    "令\n",
    "$$  S_1 = \\frac{1}{n_1}  \\sum_{x^{(i)} \\in {D_1}} ({x^{(i)}- \\mu_1})  ({x^{(i)}- \\mu_1})^T, $$\n",
    "\n",
    "$$  S_2 = \\frac{1}{n_2}  \\sum_{x^{(i)} \\in {D_2}} ({x^{(i)}- \\mu_2})  ({x^{(i)}- \\mu_2})^T, $$\n",
    "\n",
    "则有\n",
    "$$\\tilde{ \\sigma}^2_1 = w^T S_1 w, \\,  \\tilde{\\sigma}^2_2 = w^T S_2 w.  $$\n",
    "\n",
    "令\n",
    "$$  \\tilde{S_1} = \\frac{1}{n_1} \\sum_{y^{(i)} \\in \\tilde{D_1}} ({y^{(i)}- \\tilde{\\mu_1}})  ({y^{(i)}- \\tilde{\\mu_1}})^T, $$\n",
    "\n",
    "$$  \\tilde{S_2} = \\frac{1}{n_2}  \\sum_{y^{(i)} \\in \\tilde{D_2}} ({y^{(i)}- \\tilde{\\mu_2}})  ({y^{(i)}- \\tilde{\\mu_2}})^T, $$\n",
    "\n",
    "则有\n",
    "$$  \\tilde{S_1} = w^T S_1 w, \\,\\,\\,  \\tilde{S_2} = w^T S_2 w,  $$\n",
    "\n",
    "现在我们就可以定义损失函数：\n",
    "\n",
    "$$ J(w) = \\frac{|\\tilde{\\mu}_1 - \\tilde{\\mu}_2|^2 }{\\tilde{S}^2_1 + \\tilde{S}^2_2 } $$\n",
    "\n",
    "我们分类的目标是，使得类别内的点距离越近(集中)，类别间的点越远越好。分母表示数据被映射到低维空间之后每一个类别内的方差之和，方差越大表示在低维空间(映射后的空间)一个类别内的点越分散，欲使类别内的点距离越近(集中)，分母应该越小越好。分子为在映射后的空间两个类别各自的中心点的距离的平方，欲使类别间的点越远，分子越大越好。故我们最大化$J(w)$，求出的w就是最优的了。\n",
    "\n",
    "因为\n",
    "\n",
    "$$ |\\tilde{\\mu}_1 - \\tilde{\\mu}_2|^2 = w^T (\\mu_1 - \\mu_2) (\\mu_1 - \\mu_2)^T w = w^T S_B w,$$\n",
    "\n",
    "其中，\n",
    "\n",
    "$$ S_B = (\\mu_1 - \\mu_2) (\\mu_1 - \\mu_2)^T .$$\n",
    "\n",
    "设$S_w = S_1 + S_2,$ 则\n",
    "\n",
    "$$ J(w) = \\frac{w^T S_B w}{w^T S_w w} $$\n",
    "\n",
    "这样就可以用最喜欢的拉格朗日乘数法了，但是有一个问题，如果分子、分母是都可以取任意值，就会导致有无穷解，我们将分母限制为长度为1（这是用拉格朗日乘子法一个很重要的技巧，在下面将说的PCA里面也会用到，如果忘记了，请复习一下高数），并作为拉格朗日乘子法的限制条件，带入得到：\n",
    "\n",
    "$$ loss(w) = w^T S_B w - (\\lambda w^T S_w w -1) $$\n",
    "\n",
    "令$$\\frac{dloss}{dw}=2 S_B w - 2 \\lambda S_w  w = 0, $$ \n",
    "\n",
    "则有$$ S_B w = \\lambda S_w w. $$ \n",
    "\n",
    "很显然，$S_B w$和$\\mu_1 - \\mu_2$是平行的，又因为对$w$扩大缩小任何倍(平移$w$)不影响结果，因此，只要找到的$w$满足条件$S_B w$与$\\mu_1 - \\mu_2$平行即可。如果$S_w$是非奇异的，则有\n",
    "\n",
    "$$w = S_w^{-1}{(\\mu_1 - \\mu_2)}. $$\n",
    "\n",
    "下面看看具体的数学推导，\n",
    "$$ S_B w = (\\mu_1 - \\mu_2) (\\mu_1 - \\mu_2)^T w = (\\mu_1 - \\mu_2) \\lambda_w .$$\n",
    "\n",
    "将上式代入特征值公式中可得$$ S_w^{-1} S_B w = S_w^{-1} (\\mu_1 - \\mu_2) \\lambda_w = \\lambda w , $$\n",
    "\n",
    "因为$w$的平移不影响结果，故可以扔掉$\\lambda_w ,  \\lambda $，因此可得\n",
    "\n",
    "$$w = S_w^{-1}{(\\mu_1 - \\mu_2)}. $$\n",
    "\n",
    "\n",
    "得到$w$之后，就可以对测试数据进行分类了。\n",
    "\n",
    "一个常见的LDA分类基本思想是假设各个类别的样本数据符合高斯分布，这样利用LDA进行投影后，可以利用极大似然估计计算各个类别投影数据的均值和方差，进而得到该类别高斯分布的概率密度函数。当一个新的样本到来后，我们可以将它投影，然后将投影后的样本特征分别带入各个类别的高斯分布概率密度函数，计算它属于这个类别的概率，最大的概率对应的类别即为预测类别。\n",
    "\n",
    "\n",
    "但是这里还有另外一种分类的思想，就以LDA二值分类为例，我们可以将测试数据投影到低维空间(直线，因为二分类问题是投影到一维空间)，得到$y$，然后看看$y$是否在超过某个阈值$y_0$，超过是某一类，否则是另一类。但是又该怎么去寻找这个$y_0$呢？\n",
    "\n",
    "因为\n",
    "\n",
    "$$ y = w^T x, $$\n",
    "\n",
    "根据中心极限定理，独立同分布的随机变量和服从高斯分布，然后利用极大似然估计求\n",
    "\n",
    "$$ p(y|label_i), $$\n",
    "\n",
    "然后用决策理论里的公式来寻找最佳的$y_0$，详情请参阅PRML。这是一种可行但比较繁琐的选取方法。\n",
    "\n",
    "其实，还有另外一种非常巧妙的方法可以确定$y_0=0$，投影之前的数据集的标签$y_{label}$是用0和1来表示，这里我们将其做一个简单的变换，\n",
    "\n",
    "$$ \\begin{cases}\n",
    "\\tilde y_{label}=\\frac{m}{n_1} & \\text{ if } x\\in D_1 \\\\ \n",
    "\\tilde y_{label}=-\\frac{m}{n_2} & \\text{ if } x\\in D_2\n",
    "\\end{cases} $$\n",
    "\n",
    "从变换后的$\\tilde y_{label}$的定义可以看出，对于样本$x^{(i)}$， 若$\\tilde y_{label}^{(i)}>0$，则$ x^{(i)} \\in D_1$，即$ y^{(i)}_{label}=0$，若$\\tilde y_{label}^{(i)}<0$，则$ x^{(i)} \\in D_2$，即$ y^{(i)}_{label}=1$.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 算法流程\n",
    "\n",
    "\n",
    "输入：数据集$D=\\{x^{(1)}, x^{(1)}, …, x^{(m)}\\}$；\n",
    "\n",
    "输出：投影后的样本集$D′$；\n",
    "\n",
    "* 计算类内散度矩阵$S_w$；\n",
    "\n",
    "* 求解向量$w$，其中$w = S_w^{-1}(\\mu_1 - \\mu_2)$；\n",
    "\n",
    "* 将原始样本集投影到以$w$为基向量生成的低维空间中(1维)，投影后的样本集就是我们需要的样本集$D′$(1维特征)。\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 1.0\n"
     ]
    }
   ],
   "source": [
    "from sklearn import datasets\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "\n",
    "def shuffle_data(X, y, seed=None):\n",
    "    if seed:\n",
    "        np.random.seed(seed)\n",
    "\n",
    "    idx = np.arange(X.shape[0])\n",
    "    np.random.shuffle(idx)\n",
    "\n",
    "    return X[idx], y[idx]\n",
    "\n",
    "\n",
    "\n",
    "# 正规化数据集 X\n",
    "def normalize(X, axis=-1, p=2):\n",
    "    lp_norm = np.atleast_1d(np.linalg.norm(X, p, axis))\n",
    "    lp_norm[lp_norm == 0] = 1\n",
    "    return X / np.expand_dims(lp_norm, axis)\n",
    "\n",
    "\n",
    "# 标准化数据集 X\n",
    "def standardize(X):\n",
    "    X_std = np.zeros(X.shape)\n",
    "    mean = X.mean(axis=0)\n",
    "    std = X.std(axis=0)\n",
    "    \n",
    "    # 做除法运算时请永远记住分母不能等于0的情形\n",
    "    # X_std = (X - X.mean(axis=0)) / X.std(axis=0) \n",
    "    for col in range(np.shape(X)[1]):\n",
    "        if std[col]:\n",
    "            X_std[:, col] = (X_std[:, col] - mean[col]) / std[col]\n",
    "    \n",
    "    return X_std\n",
    "\n",
    "\n",
    "# 划分数据集为训练集和测试集\n",
    "def train_test_split(X, y, test_size=0.2, shuffle=True, seed=None):\n",
    "    if shuffle:\n",
    "        X, y = shuffle_data(X, y, seed)\n",
    "        \n",
    "    n_train_samples = int(X.shape[0] * (1-test_size))\n",
    "    x_train, x_test = X[:n_train_samples], X[n_train_samples:]\n",
    "    y_train, y_test = y[:n_train_samples], y[n_train_samples:]\n",
    "\n",
    "    return x_train, x_test, y_train, y_test\n",
    "\n",
    "\n",
    "def accuracy(y, y_pred):\n",
    "    y = y.reshape(y.shape[0], -1)\n",
    "    y_pred = y_pred.reshape(y_pred.shape[0], -1)\n",
    "    return np.sum(y == y_pred)/len(y)\n",
    "\n",
    "\n",
    "# 计算矩阵X的协方差矩阵\n",
    "def calculate_covariance_matrix(X, Y=np.empty((0,0))):\n",
    "    if not Y.any():\n",
    "        Y = X\n",
    "    n_samples = np.shape(X)[0]\n",
    "    covariance_matrix = (1 / (n_samples-1)) * (X - X.mean(axis=0)).T.dot(Y - Y.mean(axis=0))\n",
    "\n",
    "    return np.array(covariance_matrix, dtype=float)\n",
    "\n",
    "\n",
    "\n",
    "class BiClassLDA():\n",
    "    \"\"\"\n",
    "    线性判别分析分类算法(Linear Discriminant Analysis classifier). 既可以用来分类也可以用来降维.\n",
    "    此处实现二类情形(二类情形分类)\n",
    "    \"\"\"\n",
    "    def __init__(self):\n",
    "        self.w = None\n",
    "        \n",
    "        \n",
    "    def transform(self, X, y):\n",
    "        self.fit(X, y)\n",
    "        # Project data onto vector\n",
    "        X_transform = X.dot(self.w)\n",
    "        return X_transform\n",
    "\n",
    "    \n",
    "    def fit(self, X, y):\n",
    "        # Separate data by class\n",
    "        X = X.reshape(X.shape[0], -1)\n",
    "        \n",
    "        X1 = X[y == 0]\n",
    "        X2 = X[y == 1]\n",
    "        y = y.reshape(y.shape[0], -1)\n",
    "        \n",
    "        # 计算两个子集的协方差矩阵\n",
    "        S1 = calculate_covariance_matrix(X1)\n",
    "        S2 = calculate_covariance_matrix(X2)\n",
    "        Sw = S1 + S2\n",
    "        \n",
    "        # 计算两个子集的均值\n",
    "        mu1 = X1.mean(axis=0)\n",
    "        mu2 = X2.mean(axis=0)\n",
    "        mean_diff = np.atleast_1d(mu1 - mu2)\n",
    "        mean_diff = mean_diff.reshape(X.shape[1], -1)\n",
    "        \n",
    "        # 计算w. 其中w = Sw^(-1)(mu1 - mu2), 这里我求解的是Sw的伪逆, 因为Sw可能是奇异的\n",
    "        self.w = np.linalg.pinv(Sw).dot(mean_diff)\n",
    "        \n",
    "\n",
    "    def predict(self, X):\n",
    "        y_pred = []\n",
    "        for sample in X:\n",
    "            sample = sample.reshape(1, sample.shape[0])\n",
    "            h = sample.dot(self.w)\n",
    "            y = 1 * (h[0][0] < 0)\n",
    "            y_pred.append(y)\n",
    "        return y_pred\n",
    "\n",
    "\n",
    "def main():\n",
    "    # 加载数据\n",
    "    data = datasets.load_iris()\n",
    "    X = data.data\n",
    "    y = data.target\n",
    "\n",
    "    # 只取label=0和1的数据，因为是二分类问题\n",
    "    X = X[y != 2]\n",
    "    y = y[y != 2]\n",
    "    \n",
    "    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33)\n",
    "    \n",
    "    # 训练模型\n",
    "    lda = BiClassLDA()\n",
    "    lda.fit(X_train, y_train)\n",
    "    lda.transform(X_train, y_train)\n",
    "    \n",
    "    # 在测试集上预测\n",
    "    y_pred = lda.predict(X_test)\n",
    "    y_pred = np.array(y_pred)\n",
    "    accu = accuracy(y_test, y_pred)\n",
    "    print (\"Accuracy:\", accu)\n",
    "    \n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    main()\n"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}