第四章：分 类

第四章：分 类

---

4.1 分类问题叙述

• 在线性回归模型中，我们假设相应变量Y是定量的，如：Y是从2001年到2019年中国总人口数量等。但是在很多情况下，响应变量Y 却是 定性的(qualitative)，例如：人的眼睛颜色是定性变量，取值为蓝色，棕色，绿色。定性变量也称为 分类(categorical) 变量，两者的统计意义是一样的。
• 现实生活中分类问题是比较常见的，下面举几个例子：
  (1). 某人来到急诊室，医生通过病人的一系列症状将病人归类到三类可能的病症中的一类。
  (2). 网上银行服务需要基于用户的IP地址，交易记录等信息辨别一个在线网上交易是否存在诈骗。
  (3). 一位生物学家对一定数量的患先天疾病的患者与没有患有先天疾病的人进行DNA排序信息分析，发现哪些的DNA突变是致病的，哪些是不致病的。
• 下面我们将用一个模拟的数据集Default(违约)数据集来阐述分类模型的概念。
  在该问题中，感兴趣的是一个人的年收入和月信用卡余额来预测其违约的状态数据集如图4-1所示：（橙色为违约状态，蓝色为未构成违约状态）
  
  

4.2 为什么线性回归不可用

在之前，我们已经用线性回归解决了很多问题，那我们自然就能想到线性回归能不能用于解决分类问题呢？那我们来探讨这个极具价值的问题吧：

• 我们先用线性回归来进行建模：$default={\beta }_0+{\beta }_{1}Balance+{\beta }_{2}Income$，只要输入Balance 和 Income 以及default的数据就能用最小二乘法估计出${\beta }_0,{\beta }_1$,设定预测的default>0.5就是违约反之不违约，感觉很完美的样子，但事实真的是这样吗？

• 我们仔细来看这个线性回归模型：
  (1).我们假设有一个穷人Lisa,他的Balance和Income都很小，那么有可能会导致default的值为负数，那么这个负数代表什么意义呢？
  
  （图4-2表示线性回归模型，在这个模型上有些default竟然是负的！）

  (2).当我们的分类变量是多类的时候，以0.5为界限划分分类就不可用了，那么我们应该怎么找到一个界限衡量多分类呢？

  基于以上问题，现在大家是否还觉得线性回归模型作为一个分类模型是否足够优秀呢？其实，为了解决以上的问题（1）我们来想想能不能将线性回归的结果default转化为区间[0:1]上，让default转变成一个违约的概率呢？下面我们来解决这个问题吧。

4.3 Logistic Regression（逻辑回归）

• 在推导逻辑回归之前，我们先来认识下一组函数，这组函数具有神奇的作用，可以将是实数轴上的数转换为[0:1]区间上的概率。
  首先，我们假设我们的线性回归模型为 $Y={\beta }_0+{\beta }_{1}X$，那么这个函数是如何将线性回归的结果转化为概率呢？这个函数就是logistic 函数，具体的形式为 $p(X)=\frac{e^{{\beta }_0+{\beta }_{1}X}}{1+e^{{\beta }_0+{\beta }_{1}X}}$，他的函数图像如图4-2：（左边是线性回归，右边是逻辑函数）

• 那么问题又来了，我们怎么去估计$p(X)=\frac{e^{{\beta }_0+{\beta }_{1}X}}{1+e^{{\beta }_0+{\beta }_{1}X}}$ 中的${\beta }_0,{\beta }_1$呢？显然，极大似然估计是个好办法，我们先写出$p(X)=\frac{e^{{\beta }_0+{\beta }_{1}X}}{1+e^{{\beta }_0+{\beta }_{1}X}}$的似然函数，如果不懂似然函数的可以用查看概率统计，这里直接用结论：我们写出他的似然函数$l({\beta }_0,{\beta }_1)=\prod _{i:y_i=1}p(x_i)\prod _{i:y_{i^{‘}}=0}(1-p(x_i))$，对这个似然函数求对数，对${\beta }_0,{\beta }_1$求偏导，求极值就可以得出系数啦，不过现实中还是使用相关的库求吧，手算太麻烦了，结果如下图所示：（这里没有将income也写进去）
  

• 到现在为止，我们已经推导出这个逻辑回归模型了，并且已经把所有的系数都估计出来啦，离问题解决只剩下最后一步了，那就是预测。其实，预测非常easy，只要将数据扔进去我们刚刚求出来的模型就好啦，下面我们来看看吧：
  加入Lisa的Balance=1000美元时，我们刚刚建立的模型是：$\hat{p}(X)=\frac{e^{\widehat{{\beta }_0}+\widehat{{\beta }_1}X}}{1+e^{\widehat{{\beta }_0}+\widehat{{\beta }_1}X}}$，我们把balance=1000代入X中，得出$\hat{p}(X)=0.00576$，违约的概率小于1%，可见她不会违约咯。

• 基于以上问题我们已经成功地将线性回归模型过渡到分类模型，让模型更加适合分类，我们成功地解决了问题(1)，但是还有点小瑕疵，那就是多个自变量的时候我们没有说明白，其实非常简单，就是将一元线性回归变成多元线性回归输入我们的逻辑函数，$p(X)=\frac{e^{{\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+....}}{1+e^{{\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+....}}$，就是这么简单咯。

• 值得注意的一点是：建立完模型记得看P-value的大小，P-value < 0.05才可以使用这个模型，如果不满足只能选择其他变量建立模型了。
  

• 对于问题(2),我们值得注意的是，逻辑回归在实际中不太用于多分类问题，因为实际效果不是很好，所以我们必须借助其他模型来解决这个问题，那让我们来解决这个遗留下来的问题吧。

4.4 线性判别分析（LDA）

• 在讨论如何解决多分类问题之前，我们先来说说贝叶斯的那些事吧。在概率统计的领域里有一条神奇的公式叫贝叶斯定理，具体的形式是：$P(Y=k|X=x)=\frac{{\pi }_k{f}_k(x)}{\sum _{l=1}^K{\pi }_l{f}_l(x)}$ ，我们先不要被公式的符号吓到，我们先来看看符号具体代表什么意思。我们假设观测有$K$类，${\pi }_k$为随机选择的观测来自第$k$类的 先验概率，也就是样本里面第$k$类的样本个数除以总样本的个数：${\pi }_k=\frac{n_k}{n}$。再来 $f_k(x)=P(X=x|Y=k)$，表示第$k$类观测的X的密度函数，说的直白一点就是在$Y=k$的样本里$X=x$的样本个数，即$f_k(x)=P(X=x|Y=k)=\frac{n_{(X=x,Y=k)}}{n_{(Y=k)}}$，最后，$\sum _{l=1}^K{\pi }_l{f}_l(x)=P(X=x)=\frac{n_{(X=x)}}{n}$，也就是样本中$X=x$的概率。

• 在讨论贝叶斯定理后，我们回到分类问题，这个定理跟我们的分类问题有什么关联呢？没错，这个公式$P(Y=k|X=x)=\frac{{\pi }_k{f}_k(x)}{\sum _{l=1}^K{\pi }_l{f}_l(x)}$给出了给定样本条件下，$Y=k$这个类别下的概率，这给分类问题提供了一条思路，那就是计算这个$P(Y=k|X=x)$，而且我们的逻辑回归就是这么干的，但是在$P(Y=k|X=x)=\frac{{\pi }_k{f}_k(x)}{\sum _{l=1}^K{\pi }_l{f}_l(x)}$这个公式中，分母$\sum _{l=1}^K{\pi }_l{f}_l(x)=P(X=x)$当样本给定的时候是一个与分类$k$无关的常数,所以我们的问题可以简化为只需要计算分子${\pi }_k{f}_k(x)$,进而比较哪个类别的概率最大就知道属于哪个类别了，因此我们的分类思路就出来啦，这个思路不同于逻辑回归，逻辑回归需要计算具体的$P(Y=k|X=x)$概率值，而我们现在的思路是通过贝叶斯定理计算贝叶斯定理的分子，比较分子最大的那个类别为最终类别。

• 在我们推导复杂算法之前，我们先推导下简单的当自变量个数只有一个的模型，即$p=1$的简单模型。我们记$P(Y=k|X=x)=\frac{{\pi }_k{f}_k(x)}{\sum _{l=1}^K{\pi }_l{f}_l(x)}$ 的分子为$g_k(x)={\pi }_k{f}_k(x)$。在这里，我们做个模型假设：假设$f_k(x)$服从正态分布，即$f_k(x)\sim N(\mu ,{\sigma }_k^2)$，而且每个${\sigma }_k^2={\sigma }^2$，同方差假设。因此$f_k(x)=\frac{1}{\sqrt{2\pi }{\sigma }_k}{e}^{-\frac{1}{2{\sigma }^2}(x-{\mu }_k{)}^2}$，最终我们的$g_k(x)={\pi }_k\frac{1}{\sqrt{2\pi }{\sigma }_k}{e}^{-\frac{1}{2{\sigma }^2}(x-{\mu }_k{)}^2}$,终于算出来啦。这个式子不是很好计算，我们对$g_k(x)$取个对数，令${\delta }_k(x)=ln(g_k(x))=ln{\pi }_k+\frac{\mu }{{\sigma }^2}x-\frac{{\mu }^2}{2{\sigma }^2}$，到这里我们的模型建立模型，我们只需要把位置的${\mu }_k$与${\sigma }^2$估计出来就好了。${\hat{\mu }}_k=\frac{1}{n_k}\sum _{i:y_i=k}{x}_i$，也就是当$y=k$这一类中$x$的平均值；${\hat{\sigma }}^2=\frac{1}{n-K}\sum _{k=1}^K\sum _{i:y_i=k}(x_i-{\hat{\mu }}_k{)}^2$，说白了就是计算每一类的方差，再求平均值。总结下上面的公式就是：
  $\{\begin{array}{l}{\delta }_k(x)=ln(g_k(x))=ln{\pi }_k+\frac{\mu }{{\sigma }^2}x-\frac{{\mu }^2}{2{\sigma }^2}\\ {\hat{\mu }}_k=\frac{1}{n_k}\sum _{i:y_i=k}{x}_i\\ {\hat{\sigma }}^2=\frac{1}{n-K}\sum _{k=1}^K\sum _{i:y_i=k}(x_i-{\hat{\mu }}_k{)}^2\end{array}$
  至此，我们的模型就建立完成了，我们只需要代入数据求出${\delta }_k(x)$，哪个$k$对应的${\delta }_k(x)$大，就是哪一类。
  （下图虚线是线性判别分析的决策边界，正态曲线哪边高样本就是哪一类）
  

• 就像逻辑回归一样，我们推到出了一个自变量的简单模型，就要泛化为多个自变量的线性判别分析了，即$p>1$。其实原理一样的，只是将一元正态分布扩展为多元正态分布：
  $f_k(x)=\frac{1}{2(\pi {)}^{\frac{p}{2}}|\Sigma {|}^{\frac{1}{2}}}{e}^{[-\frac{1}{2}(x-{\mu }_k{)}^T\Sigma (x-{\mu }_k)]}$
  $\widehat{{\mu }_k}=({\mu }_{k1},{\mu }_{k2},......,{\mu }_{kp}),\hat{\Sigma }=\frac{1}{p-1}\sum _{j=1}^p(x_j-\overline{x})(x_j-\overline{x}{)}^T$
  ${\delta }_k(x)=ln({\pi }_k{f}_k(x))=ln({\pi }_k)-(\frac{p}{2}ln(2\pi )+\frac{1}{2}ln(|\Sigma |))-\frac{1}{2}(x-{\mu }_k{)}^T{\Sigma }^{-}1(x-{\mu }_k)=x^T\hat{\Sigma }{\hat{\mu }}_k-\frac{1}{2}{\hat{\mu }}_k^T{\hat{\Sigma }}^{-1}{\hat{\mu }}_k+ln{\hat{\pi }}_k$

• 从上面的分析，我们至少能完整解决多分类问题了，那我们又开始想一个问题了，能不能优化这个算法呢？线性判别分析究竟存在什么问题呢？

4.5 二次判别分析（QDA）

• 为了优化线性判别分析，我们从线性判别分析的模型假设入手，我们假设$f_k(x)$是同方差的正态分布，然而现实中真的完全都是同方差吗？例如同是蓝色眼睛的人身高和体重和黑眼睛的身高体重的方差一定一样吗？为了打破这个假设，我们必须假设$f_k(x)$的方差各自不同，即我们假设${\Sigma }_1\ne {\Sigma }_2\ne ......\ne {\Sigma }_K$,那判别函数就变成${\delta }_k(x)=-\frac{1}{2}{x}^T{\Sigma }_k^{-1}x+x^T{\Sigma }_k^{-1}{\mu }_k-\frac{1}{2}{\mu }_k^T{\Sigma }_k^{-1}{\mu }_k+ln{\pi }_k$，我们仔细观察下这个表达式，就会知道 当不同 的$Cov(x_i,x_j)=0$时，刚好QDA模型就是朴素贝叶斯模型，所以朴素贝叶斯是简化版的QDA。
  （下图是QDA的绿色决策边界，可见QDA模型是非线性的二次函数决策边界）
  
  

4.6 分类方法的比较

• 在前面的探讨中，我们已经完美地解决了某些分类问题，并且我们从两个大的角度：计算$P(Y=k|X=x)$还是比较${\delta }_k(x)$的大小来决定某个样本的类别，也建立了三个分类模型：逻辑回归，线性判别分析LDA和二次判别分析QDA来解决分类问题，同时之前我们也探讨过了K近邻法（KNN）也可以解决某些分类问题，那么在这些方法之间有什么优劣之处呢？

• 逻辑回归 VS 线性判别分析LDA
  在回答这个问题之前，我们来重新对这两个模型的公式做些比较。首先是LDA，我们假设类别$K=2$，并且假设每种类别分别为$K=1，K=2$，设$P_1(X),P_2(X)=1-P_1(X)$分别为$K=1，K=2$的概率。对LDA模型下取对数发生比$log(\frac{P_1(X)}{1-P_1(X)})=log(\frac{P_1(X)}{P_2(X)})=c_0+c_{1}X$，其中$c_0,c_1$是${\mu }_0,{\mu }_1,{\sigma }^2$的函数。
  再来看逻辑回归，$log(\frac{P_1}{P_2})={\beta }_0+{\beta }_{1}X$，LDA与逻辑回归的对数发生比都是关于$X$的线性函数，因此他们的决策边界都是线性的，但由于逻辑回归使用极大似然估计计算出最优的${\beta }_0,{\beta }_1$，而LDA是用正态分布的均值和方差计算出来的，他们的拟合过程有差异，但是他们的结果应该是相近的，但并非是必然的。LDA假设观测服从每一类协方差矩阵都相同的正态分布（均值不一定相同），所以当样本对这个假设近似满足的时候，LDA应该会比逻辑回归的效果会更好，如果假设不满足，则逻辑回归的效果比LDA更好。

• KNN VS 逻辑回归与LDA
  当我们谈论LDA时，我们应该清楚KNN的分类原理跟我们刚刚谈论的完全不同，对观测$X=x$进行预测时，先找出该点最近的$K$个点，找出K个点中属于类别数最多的点作为该点的分类，说白了就是“少数服从多数”，因此KNN是一个彻底的非参数方法，也就是他没有假设决策边界的形状。因此，当决策边界高度非线性时，KNN会比逻辑回归和LDA效果更好，反之。

• 二次判别分析QDA VS KNN,逻辑回归，LDA
  QDA是非参数方法KNN与线性方法LDA，逻辑回归的一种折中方法，因为QDA假设得到的是一个二次决策边界，他比线性决策边界的模型适用性更广，虽然不如KNN方法的光滑程度（模型复杂度）高，但是QDA在固定训练数据量的问题上一般比KNN有更好的效果，因为QDA对决策边界的形状做了一定的假设。
  （下图可看到在不通类型的数据下，各类模型的错误率不一样。）
  

• 最后做个总结，当真实决策边界是线性的LDA和逻辑回归会比较好，当真实决策边界是一般非线性的QDA会比较好，当真实决策边界是更加复杂的，KNN会比较合适。

总而言之，选择压倒一切，没有万能的模型，只有适合的模型。

4.7 具体代码实现（R语言版本）

4.7.1 R语言版本

为什么使用R而不用python呢？因为这几个我们的模型都是传统的统计模型，统计模型在R这种专业的统计语言里比较全面，所以在面对统计模型时R是个十分不错的选择。

• 1.数据的初步探索：
  我们将会用到两个数据: Smarket数据和Purchase数据。这两个数据都在R包ISLR中，因此需要我们先用函数require()或者library()加载ISLR(如果电脑中没有安装ISLR，请先安装ISLR，然后再加载。首先，请先使用函数class(), names(), dim(), head(), summary()等初步探索该数据，可以对数据有一个直观的认识。例如：

require(ISLR)   # 加载ISIR包
names(Smarket)  # 查看Smarket的变量名
dim(Smarket)    # 查看Smarket数据集的维度
cor(Smarket[, -9])   # 查看各变量的相关性

如果大家对股票投资有一些了解，可以看出上面的结果和一般的预期是一致的，大部分变量之间的相关性都是非常小的。如果发现某些变量和当天的收益率相关性很强，投资者就可以利用这种相关性投资股票，获利丰厚了。现实是，在股票上赚钱是一件很困难的事情，因此也不太可能这么容易就可以找到和未来收益相关的变量。仔细看上面的结果，也可以发现“Year”和“Volume”之间的相关关系是“0.539”, 一个比较大的相关关系。这是为什么呢？

require("RColorBrewer")    # 设置颜色
miscolores = brewer.pal(8,"Set2")[1:8]  # 设置颜色
par(mar=c(5,5,0.1,0.1))   # 设置颜色
plot(Smarket[, "Volume"], col = miscolores[1], xlab = "Index", ylab = "Volume")  # 作出Volume散点图

可以发现，“Volume”总的趋势是在上升，说明该股票的交易越来越活跃。因此，“Year”和“Volume”可以计算出一个比较大的相关关系。

直观来看，Sarmket数据中有两个变量可以作为因变量，“Today”和“Direction”。“Today”表示股票当天的收益率，是定量变量。“Direction”表示股票当天是涨还是跌，是定性变量。因此，当“Today”是非负的，“Direction”等于“Up”；反之，当“Today”是负的，“Direction”等于“Down”。

attach(Smarket) # 引用Smarket数据，以便在后面不需要再同一个代码块内重复输入数据名字
# 把Today>0的变量direction设为up，反之为down
my_Direction <- ifelse(Today >= 0, "Up", "Down")
sum(my_Direction == Direction)
as.factor(my_Direction[1:5])
Direction[1:5]
detach() # 解除数据约定

当“Today”作为因变量时，我们需要建立的是回归模型；当“Direction”作为因变量时，我们需要建立分类模型。在这次的实验中，我们要学习的是分类方法的建模，因此把“Direction”作为因变量。

• 2.Logistic回归模型
  2.1 建立逻辑回归模型
  在R中，可以使用函数glm()建立logistic回归模型，函数glm()和函数lm()的使用方法有很多一样的地方。特别需要注意的是，函数glm()需要使用参数family = binomial来表示我们要建立的模型是logristic回归模型。因为, “glm”表示generalized linear model(广义线性模型), 函数glm()可以拟合这个系列的模型。同学们可以在R中输入?glm查看其它的参数family的选择。

glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Smarket, family = binomial) # 建立逻辑回归模型
glm.fit

我们建立的logistic模型，以Direction为因变量，以Lag1, Lag2, Lag3, Lag4, Lag5, Volume为自变量。输入glm.fit可以看到logistic模型的拟合系数。和函数lm()的结果一样， 我们也可以用names()可以看到glm.fit包含哪些信息：

names(glm.fit) # 我们可以看到很多结果，包括RSS，weight等

glm.fit$coefficients # 查看逻辑回归拟合出来的系数

我们也可以用函数summary()看到更多的信息：

summary(glm.fit)

函数summary()的结果主要包含了系数，系数的标准差，统计量，p值。从上面的结果，可以看到所有的自变量的p-value都比较大，说明这些自变量对Direction没有显著的作用。在这些自变量中，Lag1的p-value最小，虽然它的p-value还是比较大，等于0.15。Lag1的系数的值是负的，说明了如果前一天股票是下跌的，那么该股票今天上涨的可能性比较大；如果前一天股票是上涨的，那么该股票今天下跌的可能性比较大。

这模型glm.fit中，我们知道R在计算的过程中把Direction的Up或者Down变成了“1”。从上述的模型中我们无法得到这个信息。这时我们可以使用函数contrasts()查看R是怎么把定性变量变成傀儡变量的：

contrasts(Smarket$Direction)  # 查看R语言怎么设置定性变量的

在Smarkdet数据中，R把Up记为1， Down记为0。 函数predict()可以得到得到概率P(Y=1|X)。在这个例子中，函数predict()得到的是股票上涨的概率。

glm.probs <- predict(glm.fit, type = "response")
length(glm.probs)  # 查看结果的长度

glm.probs[1:6] # 查看结果的前6个

如果大家在R console中输入?predict.glm查看函数predict()的帮助文档，可以看到第一个参数是object, 该object是函数glm()的返回值，在我们这里也就是“glm.fit”； 第二个参数是newdata, 指的是我们要预测的自变量的值组成的数据框，该参数的缺省值是训练数据，即函数predict()的返回值是拟合值；第三个参数是type, 有三个选择，“link”, “response”, “terms”, 一般我们选择“response”, 表示我们要得到的结果是概率P(Y=1|X)。在这里，我们没有改变第二个参数newdata的值，因此，函数predict()默认使用训练数据作为参数newdata的输入，因此我们可以看到训练数据中1250个观测值的拟合概率P(Y=1|X)。结果“glm.probs”的第一个数字表示第一个观测值的预测概率为0.5070841。

如果我们已经知道昨天的成交量“Volume = 1.5”, 并且“Lag1 = 1, Lag2 = 1, Lag3 = 1, Lag4 = 1, Lag5 = 1”, 我们想预测今天的股票是涨还是跌，可以使用如下的方式：

glm.today <- predict(glm.fit, newdata = data.frame(Lag1 = 1, Lag2 = 1, Lag3 = 1, Lag4 = 1, Lag5 = 1, Volume = 1.5), type = "response")
glm.today

通过使用函数predict()，我们已经得到了概率值P(Y=1|X)。如果我们使用如下规则判断股票涨跌：如果概率值P(Y=1|X)大于等于0.5，那么我们判断今天的股票是涨,“Up”;如果如果概率值P(Y=1|X)小于0.5，那么我们判断今天的股票是跌,“Down”。可以使用下面两行代码实现上面的规则:

glm.pred <- rep("Down", 1250)
glm.pred[glm.probs >= 0.5] <- "Up"
glm.pred[1:6]

第一行代码表示产生一个向量，该向量的元素全部是“Down”, 向量的长度是1250。第二行代码glm.probs >= 0.5判断哪些观测点的“Up”的概率大于等于0.5。然后把glm.pred中对应的元素赋值为“Up”。然后我们可以使用函数table()得到一个confusion matrix:

attach(Smarket)
table(glm.pred, Direction)
detach()
(507 + 145) / 1250
mean(glm.pred == Smarket$Direction)

函数table()输出的结果就是confusion matrix了。在上面的代码中，函数table()的第一个参数是预测的股票涨跌，第二个参数是实际的股票涨跌。输出的结果confusion matrix，列表示预测的结果，行表示实际的结果。因此，confusion matrix中，145表示实际股票是“Down”,预测也为“Down”的观测点个数是145; 141表示实际股票是“Up”, 但是预测为“Down”的观测点个数是141；457表示实际股票是“Down”, 但是预测为“Up”的观测点个数是457；507表示实际股票是“Up”,预测也为“Up”的观测点个数507。因此正确的预测比例是0.5216。

在股票的预测中，仅仅从几天内的收益率与成交量的信息，用logistic 回归模型得到52%的成绩应该已经算不错的成绩了。但是，值得注意的是现在得到的这个正确率是训练数据的拟合正确率。该模型的训练数据的错误率是100%−52%=48%。我们说过，训练误差(training error)一般情况下是过于乐观的，也就是说，训练误差通常是小于真正的预测误差的。因此，现在得到的这个训练正确率有可能是大于真实的正确率的。
2.2Logistic回归模型的预测误差
为了得到更加准确的预测正确率的估计，一个自然的想法是，只用一部分数据建立Logistic回归模型，然后用剩下的数据计算模型的预测误差。在这个数据中，我们可以使用2001-2004的数据建立一个Logistic回归模型，然后再用建立的模型预测2005的股票涨跌，把预测的股票涨跌与真实的股票涨跌相比较，我们就可以得到模型的预测正确率或者错误率。

train <- (Smarket$Year < 2005)
train[1:10]
(!train)[1:10]
summary(train)

Smarket_2001TO2004 <- Smarket[train, ]
head(Smarket_2001TO2004)
tail(Smarket_2001TO2004)

Smarket_2005 <- Smarket[!train, ]
head(Smarket_2005)

train是一个布尔向量(Boolean vector), train的长度是1250，只包含FALSE或者“TRUE”。在这里，如果数据的Year小于2005，那么对应位置的train的值赋值为“TRUE”, 如果数据的Year大于2005，那么对应位置的train的值赋值为“FALSE”。在R中， !表示相反的，即!TRUE表示FALSE, !FALSE表示TRUE。因此, !train表示与train长度相同的向量，只是train向量中的TRUE，在!train中是FALSE; train向量中的FALSE，在!train中是TRUE。 表达式 Smarket[train, ]表示取数据框 Smarket中train为TRUE所对应的行。因此，Smarket_2001TO2004表示2001-2004的股票数据， Smarket_2005表示2005的股票数据。

这时我们可以使用2005前的股票数据建立Logistic回归模型, 然后再预测2005的股票涨跌。只需要在函数glm()中加入参数subset:

glm.fit <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume, data = Smarket, family = binomial, subset = train)
summary(glm.fit)

然后预测2005的股票涨跌:

glm.probs <- predict(glm.fit, newdata = Smarket_2005, type = "response")
length(glm.probs)

glm.preds <- rep("Down", length(glm.probs))
glm.preds[glm.probs >= 0.5] <- "Up"
glm.probs[1:6]
glm.preds[1:6]

现在，我们可以用函数table()得到confusion matrix:

table(glm.preds, Smarket_2005$Direction)
mean(glm.preds == Smarket_2005$Direction)
mean(glm.preds != Smarket_2005$Direction)

在R中，符号== 表示判断 ==两边的向量或者矩阵对应的元素是否相等；符号!=表示判断!=两边的向量或者矩阵对应的元素是否不相等。因此, mean(glm.preds == Smarket_2005$Direction)表示预测正确率；mean(glm.preds!=Smarke{t}_{2}005$Direction)表示预测错误率。

现在，预测正确率只有48%, 低于扔硬币的准确率！说明我们建立的Logistic回归模型不仅不能帮助我们预测股票走势，正确率反而低于随机的猜测。

有什么简单的办法可以帮助我们提高预测的准确率呢？我们认真看看上面summary(glm.fit)的结果，可以看到，所有的自变量的系数都不是显著的不等于0的，即，所有自变量的p值都是大于0.05的。只是Lag1和Lag2的p值相对其他的小一些。如果我们只选取Lag1和Lag2作为自变量建立Logistic回归模型效果会不会更好一些呢？

glm.sub <- glm(Direction ~ Lag1 + Lag2, data = Smarket_2001TO2004, family = binomial)
glm.sub.probs <- predict(glm.sub, newdata = Smarket_2005, type = "response")
glm.sub.preds <- rep("Down", length(glm.sub.probs))
glm.sub.preds[glm.sub.probs >= 0.5] <- "Up"
table(glm.sub.preds, Smarket_2005$Direction)
mean(glm.sub.preds == Smarket_2005$Direction)

这时我们的模型表现的比之前的好很多了，正确率可以达到56%。从confusion matrix, 我们也可以看到，当模型预测股票跌时，模型的正确率只有35/(35+35)=50%; 当模型预测股票涨时，模型的正确率有106/(106+76)=58%。所以，我们可以采用如下的策略投资股票：模型预测为跌时，我们不做任何的投资；模型预测为涨时，我们可以买进股票。这个策略的正确率可以达到58%!

• 3. LDA与QDA
  3.1 LDA线性判别分析
  在这一小节中，我们将学习使用函数lda()在R中实现LDA。函数lda()在R的包MASS中，并且函数lda()的使用方法和函数lm()或者glm()是非常相似的。

require(MASS)
lda.fit <- lda(Direction ~ Lag1 + Lag2, data = Smarket_2001TO2004)
lda.fit

lda.fit包含了很多有用的信息，比如，每一类的先验概率。上面结果中，0.491984表示在数据Smarket_2001TO2004中，有49.1984%天该股票是下跌的。

与函数lm()的结果类似，我们也可以使用函数predict()得到预测值:

lda.pred <- predict(lda.fit, Smarket_2005)
names(lda.pred)

结果lda.pred包含了最终的预测值，class和posterior。posterior表示预测为“Down”或者“Up”的概率。class是使用阈值0.5对posterior分类的结果，即，当posterior中，表示“Up”的概率大于0.5时，class = “Up”， 反之，class = “Down”。

lda.pred$class[1:10]
head(lda.pred$posterior)
post.pred <- rep("Down", nrow(Smarket_2005))
post.pred[lda.pred$posterior[, "Up"] > 0.5] <- "Up"
sum(post.pred == lda.pred$class)
nrow(Smarket_2005)

这里，同样的我们可以用函数table()计算confusion matrix:

table(lda.pred$class, Smarket_2005$Direction)
mean(lda.pred$class == Smarket_2005$Direction)

从上面结果可以看到，在数据Smarket中, 使用logistic回归模型和lda的效果是非常接近的。(印证了我们刚刚的理论推导)
3.2 QDA二次判别分析
在R中，QDA可以使用函数qda()实现，函数qda()同样包含在R包MASS中。而且使用的方法和函数lda()基本一样。

qda.fit <- qda(Direction ~ Lag1 + Lag2, data = Smarket_2001TO2004)
qda.fit

和函数lda()的结果相似，函数qda()的结果也包含了不同类的先验概率和每一个类下，自变量的平均值。可以使用函数predict()得到预测值:

qda.pred <- predict(qda.fit, Smarket_2005)
names(qda.pred)
table(qda.pred$class, Smarket_2005$Direction)
mean(qda.pred$class == Smarket_2005$Direction)

在这个数据中，QDA预测的正确率可以达到60%。这在股票的预测中是一个非常不错的结果。QDA的表现比LDA的表现好很多，说明了LDA假设不同类中的协方差矩阵相等在这个例子中可能是不成立的。

• 4. KNN
  在R中，可以使用包class中的函数knn()实现。函数knn()的语法和我们学过的lm(), lda(), qda()等有所不同。lm(), lda(), qda()等函数做预测都有两步：第一步用这些函数建立模型，第二步用函数predict()得到预测。而函数knn()需要输入4个参数:
  (1).一个矩阵或者数据框，包含了训练数据的自变量的信息
  (2).一个矩阵或者数据框，包含了测试数据的自变量的信息
  (3).一个向量，包含了训练数据的因变量
  (4).一个常数K, 最近的观测点的个数

首先对数据做一些处理:

train.x <- Smarket_2001TO2004[, c("Lag1", "Lag2")]
test.x <- Smarket_2005[, c("Lag1", "Lag2")]
train.y <- Smarket_2001TO2004[, "Direction"]

然后再把这些数据带入到函数knn()中:

require(class)  # 如果提示没有包“class”， 要先安装包“class”哦
knn.pred <- knn(train.x, test.x, train.y, k = 1)
table(knn.pred, Smarket_2005$Direction)
mean(knn.pred == Smarket_2005$Direction)

上面的KNN模型，我们使用了参数k=1, 得到的预测正确率是50%。说明模型的效果等同于随机的猜测。我们在课堂上结果，如果k的值很小，模型容易过拟合，最终的预测效果会变差。因此，我们可以尝试更大的k。

knn.pred <- knn(train.x, test.x, train.y, k = 3)
table(knn.pred, Smarket_2005$Direction)
mean(knn.pred == Smarket_2005$Direction)
knn.pred <- knn(train.x, test.x, train.y, k = 5)
table(knn.pred, Smarket_2005$Direction)
mean(knn.pred == Smarket_2005$Direction)

我们可以看到，当k=3时，模型的预测正确率有所提高，但是当k=5时，模型的预测正确率反而降低了。

4.8 结语

到现在为止，我们已经掌握了分类的基本模型，这些模型在我们面对一般的分类问题的时候已经绰绰有余了，下一步我们将要一起探索几个问题：
（1）既然我们的模型那么多，我们应该怎么在面对数据的时候选出最优秀的模型呢？是逻辑回归？LDA?QDA?KNN?
（2）当我们选定模型后怎么选择合适的变量使我们的模型最优呢？

参考文献

[1]加雷斯.詹姆斯，丹妮拉.威腾 .统计学习导论[M].北京:机械工业出版社,2015.6.
[2]李航.统计学习方法[M].北京:清华大学出版社,2012.3.
[3]杰克.万托布拉斯.Python数据科学手册[M].北京:人民邮电出版社,2018.2.
[4]富朗索瓦.肖莱.Python深度学习[M].北京:人民邮电出版社,2018.8