7 判别分析

判别分析(discriminant analysis)属于“有监督（指导）学习”方法，从所谓“训练样本”经过分析计算得到一个判别规则，对新的数据可以利用判别规则判断新数据观测的类属。训练样本中既有用来分类的解释变量（自变量），又有真实的类属（标签，因变量）。

常用判别分析方法有距离判别、Fisher判别和Bayes判别，逻辑斯谛回归也经常用在判别问题中（尤其是两类的判别），分类树也是用于判别的方法。

当标签（因变量）只有两个类时，判别分析问题与假设检验问题有相似之处。假设检验问题更强调拒绝零假设的结论一定要可靠，所以更重视对立假设；两类的判别问题并不强调某个类别，或者按照先验概率、损失函数对不同类别施加不同的影响。

7.1 距离判别

7.1.1 距离

距离判别的思想是计算每个样品与各类中心的距离，把样品分到最近的一类。

常用距离为欧式距离与马氏(Mahalonobis)距离。

设\(\boldsymbol X=(X_1, X_2, \dots, X_p)^T\)与\(\boldsymbol Y = (Y_1, Y_2, \dots, Y_p)^T\) 是两个随机向量，有相同的协方差矩阵\(\Sigma\)，定义\(\boldsymbol X\)的观测值\(\boldsymbol x=(x_1, x_2, \dots, x_p)^T\) 与\(\boldsymbol Y\)的观测值\(\boldsymbol y=(y_1, y_2, \dots, y_p)^T\) 的马氏距离为 \[\begin{align} d(\boldsymbol x, \boldsymbol y) = \sqrt{ (\boldsymbol x - \boldsymbol y)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol y) } . \end{align}\] 当\(\Sigma\)为单位阵\(I\)时马氏距离与欧式距离相同。

7.1.1.1 马氏距离的解释

马氏距离排除了量纲的影响，并考虑到了不同变量之间的相关性。

设\(E \boldsymbol X = \boldsymbol\mu_{\boldsymbol X}\), \(E \boldsymbol Y = \boldsymbol\mu_{\boldsymbol Y}\), \(\text{Var}(\boldsymbol X) = \text{Var}(\boldsymbol Y)=\Sigma\)。令\(\boldsymbol U = \Sigma^{-1/2}(\boldsymbol X - \boldsymbol\mu_{\boldsymbol X})\), \(\boldsymbol V = \Sigma^{-1/2}(\boldsymbol Y - \boldsymbol\mu_{\boldsymbol Y})\), 则\(E\boldsymbol U = E\boldsymbol V = \boldsymbol 0\), \(\text{Var}(\boldsymbol U) = \text{Var}(\boldsymbol V) = I\)。

距离 \[\begin{align*} d^2(\boldsymbol X, \boldsymbol Y) =& (\boldsymbol X - \boldsymbol Y)^T \Sigma^{-1} (\boldsymbol X - \boldsymbol Y) \\ =& \left( \Sigma^{-1/2} \boldsymbol X - \Sigma^{-1/2} \boldsymbol Y \right)^T \left( \Sigma^{-1/2} \boldsymbol X - \Sigma^{-1/2} \boldsymbol Y \right) \\ =& \left\| \Sigma^{-1/2} \boldsymbol X - \Sigma^{-1/2} \boldsymbol Y \right\|^2 = \left\| \boldsymbol U - \boldsymbol V + \Sigma^{-1/2}(\boldsymbol\mu_{\boldsymbol X} - \boldsymbol\mu_{\boldsymbol Y}) \right\|^2 \end{align*}\] 所有的随机部分都已变成两个零均值、单位方差阵的随机变量的差。

考虑\(\Sigma=\text{diag}(\sigma_1^2, \sigma_2^2, \dots, \sigma_p^2)\)的情形。这时 \[\begin{align*} d^2(\boldsymbol x, \boldsymbol y) =& \sum_{j=1}^p \left(\frac{x_i - y_i}{\sigma_i} \right)^2, \end{align*}\] 不受量纲影响。

对一般正定协方差阵\(\Sigma\), 令\(\Lambda = \text{diag}(\sqrt{\sigma_{11}}, \sqrt{\sigma_{22}}, \dots, \sqrt{\sigma_{pp}})\), \(R = \Lambda^{-1} \Sigma \Lambda^{-1}\), 则\(R\)为\(\Sigma\)对应的相关系数阵。\(\Sigma = \Lambda R \Lambda\)。这时 \[\begin{align*} d^2(\boldsymbol x, \boldsymbol y) =& \left[ \Lambda^{-1}(\boldsymbol x - \boldsymbol y) \right]^T R^{-1} \left[ \Lambda^{-1}(\boldsymbol x - \boldsymbol y) \right] \end{align*}\] 可以看出\(\Lambda^{-1}(\boldsymbol x - \boldsymbol y)\)不受量纲影响。

7.1.2 两个总体的距离判别

设有两个总体\(G_1, G_2\)，均值分别为\(\boldsymbol\mu_1, \boldsymbol\mu_2\)， \(\boldsymbol\mu_1 \neq \boldsymbol\mu_2\)。有共同的协方差阵\(\Sigma\)。在\(\boldsymbol\mu_1, \boldsymbol\mu_2, \Sigma\)已知时，为了判断一个样品\(\boldsymbol x\)属于哪一个总体（类），可以用如下判别规则： \[\begin{aligned} \begin{cases} \text{判} \boldsymbol x \in G_1, & d(\boldsymbol x, \boldsymbol\mu_1) \leq d(\boldsymbol x, \boldsymbol\mu_2), \\ \text{判} \boldsymbol x \in G_2, & d(\boldsymbol x, \boldsymbol\mu_1) > d(\boldsymbol x, \boldsymbol\mu_2) . \end{cases} \end{aligned}\]

为判断上述条件，只要计算\(d^2(\boldsymbol x, \boldsymbol\mu_2) - d^2(\boldsymbol x, \boldsymbol\mu_1)\)的正负号，非负时判入\(G_1\)，负号时判入\(G_2\)。

7.1.2.1 判别函数

距离平方的差为 \[\begin{aligned} & d^2(\boldsymbol x, \boldsymbol\mu_2) - d^2(\boldsymbol x, \boldsymbol\mu_1) \\ =& \boldsymbol x^T \Sigma^{-1} \boldsymbol x - 2 \boldsymbol\mu_2^T \Sigma^{-1} \boldsymbol x + \boldsymbol\mu_2^T \Sigma^{-1} \boldsymbol\mu_2 \\ & - \boldsymbol x^T \Sigma^{-1} \boldsymbol x + 2 \boldsymbol\mu_1^T \Sigma^{-1} \boldsymbol x - \boldsymbol\mu_1^T \Sigma^{-1} \boldsymbol\mu_1 \\ =& 2(\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x + \boldsymbol\mu_2^T \Sigma^{-1} \boldsymbol\mu_2 - \boldsymbol\mu_1^T \Sigma^{-1} \boldsymbol\mu_1 \\ =& 2(\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x \\ & + \boldsymbol\mu_2^T \Sigma^{-1} \boldsymbol\mu_2 - (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol\mu_1 - \boldsymbol\mu_2^T \Sigma^{-1} \boldsymbol\mu_1 \\ =& 2(\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x \\ & - \boldsymbol\mu_2^T \Sigma^{-1} (\boldsymbol\mu_1 - \boldsymbol\mu_2) - (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol\mu_1 \\ =& 2(\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x - (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} (\boldsymbol\mu_1 + \boldsymbol\mu_2) \\ =& 2(\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \left(\boldsymbol x - \frac{\boldsymbol\mu_1 + \boldsymbol\mu_2}{2} \right) \end{aligned}\]

令\(\boldsymbol a = \Sigma^{-1}(\boldsymbol\mu_1 - \boldsymbol\mu_2)\), 则 \[\begin{align} W(\boldsymbol x) \stackrel{\triangle}{=} \frac12 (d^2(\boldsymbol x, \boldsymbol\mu_2) - d^2(\boldsymbol x, \boldsymbol\mu_1)) = \boldsymbol a^T \left(\boldsymbol x - \frac{\boldsymbol\mu_1 + \boldsymbol\mu_2}{2} \right) . \end{align}\] 这是\(\boldsymbol x\)的线性函数，称为线性判别函数。称\(\boldsymbol a\)为判别系数。对应的判别规则为: \[\begin{align} \begin{cases} \text{判} \boldsymbol x \in G_1, & \text{当} W(\boldsymbol x) \geq 0, \\ \text{判} \boldsymbol x \in G_2, & \text{当} W(\boldsymbol x) < 0 . \end{cases} \end{align}\]

方程 \(\left\{\boldsymbol x \in \mathbb R^p:\; \boldsymbol a^T (\boldsymbol x - \frac{\boldsymbol\mu_1 + \boldsymbol\mu_2}{2}) = 0 \right\}\) 是\(\mathbb R^p\)中超平面，穿过连接两个类中心的线段中点， \(\boldsymbol a\)是法线方向。超平面把空间\(\mathbb R^p\)分成判入\(G_1\)的半空间和判入\(G_2\)的半空间，法线正向的半空间判入\(G_1\)。

当\(p=1\)时，设\(G_1\)为\(\text{N}(\mu_1, \sigma^2)\), \(G_2\)为\(\text{N}(\mu_2, \sigma^2)\), 则：判别系数\(a = \sigma^{-2}(\mu_1 - \mu_2)\)；线性判别函数 \[\begin{align} W(x) = a \left(x - \frac{\mu_1 + \mu_2}{2} \right) = \frac{(\mu_1 - \mu_2) \left(x - \frac{\mu_1 + \mu_2}{2} \right)}{\sigma^2} . \end{align}\] 设\(\mu_1 < \mu_2\), 则\(a < 0\), 判别准则为 \[\begin{aligned} \begin{cases} \text{判} \boldsymbol x \in G_1, & \text{当} x \leq \frac12(\mu_1 + \mu_2), \\ \text{判} \boldsymbol x \in G_2, & \text{当} x > \frac12(\mu_1 + \mu_2) . \end{cases} \end{aligned}\]

要注意的是，以上的简单判别规则是在假定两类的协方差阵相等时得到的。如果协方差阵不相等，判别函数不再是线性的，分界点也不恰好在中点。总共只有两类时规则比较简单，可以得到一个判别函数。如果有多个类则需要计算多个判别函数。

7.1.2.2 距离判别的规则训练

对于实际的判别问题，真实的\(\boldsymbol\mu_1, \boldsymbol\mu_2, \Sigma\)都是未知的，可以从训练样本估计。

设\(\boldsymbol x^{(1)}_1, \boldsymbol x^{(1)}_2, \dots, \boldsymbol x^{(1)}_{n_1}\)是来自总体\(G_1\)的训练样本的值， \(\boldsymbol x^{(2)}_1, \boldsymbol x^{(2)}_2, \dots, \boldsymbol x^{(2)}_{n_2}\)是来自总体\(G_2\)的训练样本的值，则\(\boldsymbol\mu_1, \boldsymbol\mu_2\)分别有无偏估计 \[\begin{align*} \hat{\boldsymbol\mu}_1 = \frac{1}{n_1} \sum_{i=1}^{n_1} \boldsymbol x^{(1)}_i, \quad \hat{\boldsymbol\mu}_2 = \frac{1}{n_2} \sum_{i=1}^{n_2} \boldsymbol x^{(2)}_i, \end{align*}\] 共同的协方差阵\(\Sigma\)有如下无偏估计: \[\begin{align*} \hat\Sigma =& \frac{1}{n_1 + n_2 - 2} \left[ \sum_{i=1}^{n_1} \left(\boldsymbol x^{(1)}_i - \hat{\boldsymbol\mu}_1 \right) \left(\boldsymbol x^{(1)}_i - \hat{\boldsymbol\mu}_1 \right)^T \right. \\ & \left. + \sum_{i=1}^{n_2} \left(\boldsymbol x^{(2)}_i - \hat{\boldsymbol\mu}_2 \right) \left(\boldsymbol x^{(2)}_i - \hat{\boldsymbol\mu}_2 \right)^T \right] . \end{align*}\]

这时判别函数中的\(\boldsymbol\mu_1, \boldsymbol\mu_2, \Sigma\)分别用估计值 \(\hat{\boldsymbol\mu}_1, \hat{\boldsymbol\mu}_2, \hat\Sigma\)代替，变成 \[\begin{align*} \hat{\boldsymbol a} =& \hat\Sigma^{-1}(\hat{\boldsymbol\mu}_1 - \hat{\boldsymbol\mu}_2), \\ \quad W(\boldsymbol x) =& \hat{\boldsymbol a}^{T} \left( \boldsymbol x - \frac{\hat{\boldsymbol\mu}_1 + \hat{\boldsymbol\mu}_2}{2} \right) . \end{align*}\] 当\(W(\boldsymbol x) \geq 0\)时判\(\boldsymbol x \in G_1\), 否则判入\(G_2\)。

7.1.2.3 二次判别

当两个总体\(G_1\)和\(G_2\)的协方差阵不相同，分别为\(\Sigma_1\)和\(\Sigma_2\)时，判别函数为 \[\begin{align*} W^*(\boldsymbol x) =& d^2(\boldsymbol x, \boldsymbol\mu_2) - d^2(\boldsymbol x, \boldsymbol\mu_1) \\ =& (\boldsymbol x - \boldsymbol\mu_2)^T \Sigma_2^{-1} (\boldsymbol x - \boldsymbol\mu_2) - (\boldsymbol x - \boldsymbol\mu_1)^T \Sigma_1^{-1} (\boldsymbol x - \boldsymbol\mu_1) . \end{align*}\] 这是\(\boldsymbol x\)的多元二次多项式函数，称为二次判别函数。

相应判别规则仍为 \[\begin{align} \begin{cases} \text{判} \boldsymbol x \in G_1, & \text{当} W^*(\boldsymbol x) \geq 0, \\ \text{判} \boldsymbol x \in G_2, & \text{当} W^*(\boldsymbol x) < 0 . \end{cases} \end{align}\] 从样本中估计\(\Sigma_1, \Sigma_2\)时分别估计即可。

7.1.2.4 例5.1（昆虫雌雄判别）

P.54例5.1。某种昆虫的体长与翅长用来判别性别。雌虫(总体\(G_1\))平均值\(\boldsymbol\mu_1 = (6, 5)^T\), 雄虫（总体\(G_2\)）平均值\(\boldsymbol\mu_2 = (8, 6)^T\)，雄虫体型较大。两个总体共同的协方差阵为 \(\Sigma=\left(\begin{array}{cc} 9 & 2 \\ 2 & 4 \end{array}\right)\)。

现有某虫测量值为\(\boldsymbol x = (7.2, 5.6)^T\)，用距离判别法做判别，判断其性别。用R程序实现。

假定两个总体有相同协方差阵时判别函数为线性函数，否则判别函数是二次函数。

7.1.2.4.1 线性判别

线性判别函数，简化公式：

W1 <- function(x, mu1, mu2, Sigma){
    a <- solve(Sigma, mu1 - mu2)
    mu3 <- (mu1 + mu2)/2

    sum(a * (x - mu3))
}

线性判别函数，直接用马氏距离的差计算，应该与简化公式得到完全相同的判别函数：

W2 <- function(x, mu1, mu2, Sigma){
    0.5*(mahalanobis(x, mu2, Sigma)
         - mahalanobis(x, mu1, Sigma))
}

取判别函数为W1:

W <- W1

判别函数大于等于零时判入第一类，否则判入第二类。

用线性判别法建立雌虫与雄虫的判别公式，只需要输入两个中心，与共同的协方差阵：

mu1 <- c(6, 5) # 雌虫平均值
mu2 <- c(8, 6) # 雄虫平均值
Sigma <- rbind(c(9,2), c(2,4))

当测量值为\((7.2, 5.6)\)做判别：

x1 <- c(7.2, 5.6)
y1 <- W(x1, mu1, mu2, Sigma)
cat('当测量值为：', x1, '时，线性判别函数值=', y1, '\n')

## 当测量值为： 7.2 5.6 时，线性判别函数值= -0.053125

判别函数值-0.053125判入第二类，雄虫。

当测量值为\((6.3, 4.9)\)时做判别：

x2 <- c(6.3, 4.9)
y2 <- W(x2, mu1, mu2, Sigma)
cat('当测量值为：', x2, '时，线性判别函数值=', y2, '\n')

## 当测量值为： 6.3 4.9 时，线性判别函数值= 0.225

判别函数值为0.225, 判入第一类, 雌虫。

7.1.2.4.2 二次判别

假设方差阵不同，这时计算马氏距离的差，判别函数为二次函数，称为二次判别。

二次判别函数。

Wsqr <- function(x, mu1, mu2, Sigma1, Sigma2){
    0.5*(mahalanobis(x, mu2, Sigma2)
         - mahalanobis(x, mu1, Sigma1))
}

输入两个总体的均值与方差阵：

mu1 <- c(6, 5) # 雌虫
mu2 <- c(8, 6) # 雄虫
Sigma1 <- rbind(c(9,2), c(2,4))
Sigma2 <- rbind(c(6,2), c(2,3))

当测量值为\((7.2, 5.6)\)的二次判别计算：

x1 <- c(7.2, 5.6)
y1 <- Wsqr(x1, mu1, mu2, Sigma1, Sigma2)
cat('当测量值为：', x1, '时，二次判别函数值=', y1, '\n')

## 当测量值为： 7.2 5.6 时，二次判别函数值= -0.03848214

判别函数值-0.03848为负数，判入第二类, 雄虫。

当测量值为\((6.3, 4.9)\)时的二次判别计算：

x2 <- c(6.3, 4.9)
y2 <- Wsqr(x2, mu1, mu2, Sigma1, Sigma2)
cat('当测量值为：', x2, '时，二次判别函数值=', y2, '\n')

## 当测量值为： 6.3 4.9 时，二次判别函数值= 0.2928795

判别函数值0.2929为非负数, 判入第一类, 雌虫。

7.1.3 多个总体的距离判别

设有\(k\)个总体\(G_1, G_2, \dots, G_k\), 分别有平均值\(\boldsymbol\mu_1, \boldsymbol\mu_2, \dots, \boldsymbol\mu_k\)，有共同的协方差阵\(\Sigma\)。对于新的测量值\(\boldsymbol x \in \mathbb R^p\), 计算它到每个总体均值的马氏距离 \[\begin{align*} d^2(\boldsymbol x, \boldsymbol\mu_j) =& (\boldsymbol x - \boldsymbol\mu_j)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol\mu_j) \\ =& \boldsymbol x^T \Sigma^{-1} \boldsymbol x - 2 \boldsymbol\mu_j^T \Sigma^{-1} \boldsymbol x + \boldsymbol\mu_j^T \Sigma^{-1} \boldsymbol\mu_j . \end{align*}\]

对不同的总体，距离中第一项是相同的。记\(\boldsymbol a_j = \Sigma^{-1} \boldsymbol\mu_j\), \(c_j = -\frac12 \boldsymbol\mu_j^T \Sigma^{-1} \boldsymbol\mu_j\), 则 \[\begin{align*} d^2(\boldsymbol x, \boldsymbol\mu_j) = \boldsymbol x^T \Sigma^{-1} \boldsymbol x - 2 (\boldsymbol a_j^T \boldsymbol x + c_j) . \end{align*}\]

令 \[\begin{align*} W_j(\boldsymbol x) = \boldsymbol a_j^T \boldsymbol x + c_j, \ j=1,2,\dots,k \end{align*}\] 把\(\boldsymbol x\)判入使得\(W_j(\boldsymbol x)\)最大的那一类。

当参数\(\boldsymbol\mu_1, \boldsymbol\mu_2, \dots, \boldsymbol\mu_k\)和\(\Sigma\)未知时，可以从训练样本中估计。

当各总体的协方差阵不相同时，直接使用\(\boldsymbol x\)到各总体的马氏距离，把\(\boldsymbol x\)判入马氏距离最小的那一类。这时的判别函数是\(\boldsymbol x\)的多元二次多项式函数，称为二次判别函数。

7.2 Fisher判别法

Fisher判别法的思想是把多维数据投影到一维直线上，使得同类数据尽量接近，异类数据尽量分开。从方差分析角度看，是组内变差要小，组间变差要大。

也要利用距离判别。属于确定性判别法，有线性判别、非线性判别和典型判别等多种常用方法。这里仅介绍线性判别法。

7.2.1 两总体Fisher判别

设有两个总体\(G_1\)和\(G_2\)。求一个向量\(\boldsymbol a \in \mathbb R^p\), 使得两个总体投影到这个向量上后距离最远。设\(G_1, G_2\)的均值分别为\(\boldsymbol\mu_1, \boldsymbol\mu_2\)，有共同的协方差阵\(\Sigma\)。

令\(\boldsymbol X^{(1)}\)来自\(G_1\), \(\boldsymbol X^{(2)}\)来自\(G_2\), \(Y_1 = \boldsymbol a^T \boldsymbol X^{(1)}\), \(Y_2 = \boldsymbol a^T \boldsymbol X^{(2)}\), 则 \[\begin{align*} E Y_1 =& \boldsymbol a^T \boldsymbol\mu_1, \quad E Y_2 = \boldsymbol a^T \boldsymbol\mu_2, \\ \text{Var}(Y_1) =& \text{Var}(Y_2) = \boldsymbol a^T \Sigma \boldsymbol a . \end{align*}\]

对投影后的两个均值计算一维的马氏距离，参数\(\boldsymbol a\)的取法应令如下一维的马氏距离最大： \[\begin{align} \max_{\boldsymbol a \in \mathbb R^p} \frac{ \left[\boldsymbol a^T(\boldsymbol\mu_1 - \boldsymbol\mu_2) \right]^2 } { \boldsymbol a^T \Sigma \boldsymbol a } . \tag{7.1} \end{align}\]

定理7.1 当\(\boldsymbol a = c \Sigma^{-1}(\boldsymbol\mu_1 - \boldsymbol\mu_2)\)时（\(c\neq 0\)为常数），式(7.1)达到最大。

特别地，当\(c=1\)时，线性函数 \[\begin{align} y = \boldsymbol a^T \boldsymbol x = (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x \end{align}\] 称为Fisher线性判别函数。

证明:

设\(\boldsymbol d \in \mathbb R^p\), 求 \[\begin{align} \max_{\boldsymbol a \in \mathbb R^p} \frac{ \left[\boldsymbol a^T \boldsymbol d \right]^2 } { \boldsymbol a^T \Sigma \boldsymbol a } . \tag{7.2} \end{align}\] 令\(\boldsymbol b = \Sigma^{\frac{1}{2}} \boldsymbol a\), 则\(\boldsymbol a = \Sigma^{-\frac{1}{2}} \boldsymbol b\)， \[\begin{aligned} \frac{ \left[\boldsymbol a^T \boldsymbol d \right]^2 } { \boldsymbol a^T \Sigma \boldsymbol a } =& \frac{[\boldsymbol b^T \Sigma^{-\frac{1}{2}} \boldsymbol d]^2} {\boldsymbol b^T \boldsymbol b} \\ =& \left[ \frac{\boldsymbol b}{\| \boldsymbol b \|} \Sigma^{-\frac{1}{2}} \boldsymbol d \right]^2 \\ \leq& \left[ \left\| \frac{\boldsymbol b}{\| \boldsymbol b \|} \right\| \left\| \Sigma^{-\frac{1}{2}} \boldsymbol d \right\| \right]^2, \end{aligned}\] 等号成立当且仅当存在常数\(c \neq 0\)使得 \[ \boldsymbol b = c \Sigma^{-\frac{1}{2}} \boldsymbol d, \] 这时(7.2)达到最大，由\(\boldsymbol b = \Sigma^{\frac{1}{2}} \boldsymbol b\)知这时 \[ \boldsymbol a = \Sigma^{-\frac{1}{2}} \boldsymbol b = c \Sigma^{-1} \boldsymbol d, \] 取\(\boldsymbol d = (\boldsymbol \mu_1 - \boldsymbol \mu_2)\)就证明了定理。

○○○○○○

取 \[\begin{align*} \mu_y =& \frac12(E Y_1 + E Y_2) = \frac12 \boldsymbol a^T (\boldsymbol\mu_1 + \boldsymbol\mu_2) \\ =& \frac12 (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} (\boldsymbol\mu_1 + \boldsymbol\mu_2), \end{align*}\] 则 \[\begin{align*} EY_1 - \mu_y =& \boldsymbol a^T \left( \boldsymbol\mu_1 - \frac{\boldsymbol\mu_1 - \boldsymbol\mu_2}{2} \right) \\ =& \frac12 (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} (\boldsymbol\mu_1 - \boldsymbol\mu_2) > 0, \\ EY_2 - \mu_y =& \boldsymbol a^T \left( \boldsymbol\mu_2 - \frac{\boldsymbol\mu_1 - \boldsymbol\mu_2}{2} \right) \\ =& -\frac12 (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} (\boldsymbol\mu_1 - \boldsymbol\mu_2) < 0 . \end{align*}\] 于是得Fisher线性判别规则为 \[\begin{align*} \begin{cases} \text{判} \boldsymbol x \in G_1, & \text{如果} y = (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x \geq \mu_y, \\ \text{判} \boldsymbol x \in G_2, & \text{如果} y = (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x < \mu_y . \end{cases} \end{align*}\]

如果记线性函数 \[\begin{align*} W(\boldsymbol x) =& (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \boldsymbol x - \mu_y \\ =& (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \left( \boldsymbol x - \frac{\boldsymbol\mu_1 + \boldsymbol\mu_2}{2} \right), \end{align*}\] 则判别规则为 \[\begin{align*} \begin{cases} \text{判} \boldsymbol x \in G_1, & \text{如果} W(\boldsymbol x) \geq 0, \\ \text{判} \boldsymbol x \in G_2, & \text{如果} W(\boldsymbol x) < 0 . \end{cases} \end{align*}\]

这与§1的在协方差相等条件下马氏距离判别法得到的线性判别函数相同。

与§1类似，如果参数\(\boldsymbol\mu_1, \boldsymbol\mu_2, \Sigma\)未知，可以从训练样本估计，估计公式与§1相同。

7.2.2 多总体Fisher判别

有多个总体时需要找多个投影方向，建立多个判别函数。设有\(k\)个总体\(G_1, G_2, \dots, G_k\)，均值为\(\boldsymbol\mu_1, \boldsymbol\mu_2, \dots, \boldsymbol\mu_k\), 有共同的协方差阵\(\Sigma\)。设随机向量\(\boldsymbol X^{(j)}\)来自总体\(G_j\), \(\boldsymbol a \in \mathbb R^p\), 随机变量\(Y_j = \boldsymbol a^T \boldsymbol X^{(j)}\)。记 \[\begin{align*} \bar{\boldsymbol\mu} = \frac{1}{k} \sum_{j=1}^k \boldsymbol\mu_j, \quad G = \sum_{j=1}^k (\boldsymbol\mu_j - \bar{\boldsymbol\mu}) (\boldsymbol\mu_j - \bar{\boldsymbol\mu})^T . \end{align*}\] 令\(\bar Y = \frac{1}{k}(Y_1 + Y_2 + \dots + Y_k)\), 则 \[\begin{align*} \mu_y = E \bar Y = \boldsymbol a^T \bar{\boldsymbol\mu} . \end{align*}\]

考虑\(EY_j\)与\(\mu_y\)的一维马氏距离平方，即 \[\begin{align*} \frac{(EY_j - \mu_y)^2}{\text{Var}(Y_j)} = \frac{[\boldsymbol a^T (\boldsymbol\mu_j - \bar{\boldsymbol\mu})]^2} { \boldsymbol a^T \Sigma \boldsymbol a }, \ j=1,2,\dots, k, \end{align*}\] 令这\(k\)个平方距离的和最大： \[\begin{align} \max_{\boldsymbol a \in \mathbb R^p} \frac{\sum_{j=1}^k \boldsymbol a^T (\boldsymbol\mu_j - \bar{\boldsymbol\mu}) (\boldsymbol\mu_j - \bar{\boldsymbol\mu})^T \boldsymbol a} { \boldsymbol a^T \Sigma \boldsymbol a } = \frac{\boldsymbol a^T G \boldsymbol a} { \boldsymbol a^T \Sigma \boldsymbol a } . \tag{7.3} \end{align}\] 则\(k\)个类投影到\(\boldsymbol a\)方向上时类间距离最远。

设\(\Sigma^{-1} G\)的非零特征值为 \(\lambda_1 \geq \lambda_2 \geq \dots \lambda_s\) (\(s \leq \min(k-1, p)\)), \(\boldsymbol e_1, \boldsymbol e_2, \dots, \boldsymbol e_s\)为相应的特征向量，且满足\(\boldsymbol e^T \Sigma \boldsymbol e = 1\)。则当\(\boldsymbol a = \boldsymbol e_1\)时，类间平方距离(7.3)达到最大，称\(\boldsymbol e_1^T \boldsymbol x\)为第一判别函数。 \(\boldsymbol a = \boldsymbol e_2\)是在约束\(\text{Cov}(\boldsymbol e_1^T \boldsymbol X, \boldsymbol a^T \boldsymbol X)=0\) 下使得(7.3)达到最大的解，称\(\boldsymbol e_2^T \boldsymbol x\)为第二判别函数。类似地，\(\boldsymbol a = \boldsymbol e_j\)是在约束 \(\text{Cov}(\boldsymbol e_i^T \boldsymbol X, \boldsymbol a^T \boldsymbol X)=0\)(\(i=1,2,\dots,j-1\)) 下使得(7.3)达到最大的解，称\(\boldsymbol e_j^T \boldsymbol x\)为第\(j\)判别函数。

(证明略去)

当参数\(\boldsymbol\mu_1, \dots, \boldsymbol\mu_k, \Sigma\)未知时，可以从训练样本中估计。多总体的Fisher判别规则比较复杂。一种做法是，把各线性判别函数值看作是数据降维后的值，对降维后的数据使用距离判别法。

7.2.2.1 例5.2（鸢尾花数据线性判别）

R中iris数据，三个类（Species），每类50个样品。每个观测有花萼长、宽与花瓣长、宽计四个测量值。以此作为训练集。

Fisher判别法将数据投影到分辨力最高的方向上，进行判别。在R中用MASS::lda()函数，也称线性判别。用predict()进行判别。用table()列表训练准确程度。注意：需要使用prior=指定一个各类相等的先验概率，如这里需要prior=rep(1/3, 3)。

数据的情况：

str(iris)

## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

knitr::kable(head(iris))

Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa

用150个样本点估计判别函数：

library(MASS)
lda1 <- lda(
  Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
  prior=rep(1/3,3), data=iris)

分类结果显示：

print(lda1)

## Call:
## lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width, 
##     data = iris, prior = rep(1/3, 3))
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa            5.006       3.428        1.462       0.246
## versicolor        5.936       2.770        4.260       1.326
## virginica         6.588       2.974        5.552       2.026
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8293776  0.02410215
## Sepal.Width   1.5344731  2.16452123
## Petal.Length -2.2012117 -0.93192121
## Petal.Width  -2.8104603  2.83918785
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9912 0.0088

结果包括公式、先验概率、各组均值向量、第一和第二线性判别系数、两个判别函数对区分各总体的贡献比例等。

可以对参与估计的150个样本点每个进行判别，用predict()函数：

res1 <- predict(lda1)
print(res1)

## $class
##   [1] setosa     setosa     setosa     setosa     setosa     setosa    
##   [7] setosa     setosa     setosa     setosa     setosa     setosa    
##  [13] setosa     setosa     setosa     setosa     setosa     setosa    
##  [19] setosa     setosa     setosa     setosa     setosa     setosa    
##  [25] setosa     setosa     setosa     setosa     setosa     setosa    
##  [31] setosa     setosa     setosa     setosa     setosa     setosa    
##  [37] setosa     setosa     setosa     setosa     setosa     setosa    
##  [43] setosa     setosa     setosa     setosa     setosa     setosa    
##  [49] setosa     setosa     versicolor versicolor versicolor versicolor
##  [55] versicolor versicolor versicolor versicolor versicolor versicolor
##  [61] versicolor versicolor versicolor versicolor versicolor versicolor
##  [67] versicolor versicolor versicolor versicolor virginica  versicolor
##  [73] versicolor versicolor versicolor versicolor versicolor versicolor
##  [79] versicolor versicolor versicolor versicolor versicolor virginica 
##  [85] versicolor versicolor versicolor versicolor versicolor versicolor
##  [91] versicolor versicolor versicolor versicolor versicolor versicolor
##  [97] versicolor versicolor versicolor versicolor virginica  virginica 
## [103] virginica  virginica  virginica  virginica  virginica  virginica 
## [109] virginica  virginica  virginica  virginica  virginica  virginica 
## [115] virginica  virginica  virginica  virginica  virginica  virginica 
## [121] virginica  virginica  virginica  virginica  virginica  virginica 
## [127] virginica  virginica  virginica  virginica  virginica  virginica 
## [133] virginica  versicolor virginica  virginica  virginica  virginica 
## [139] virginica  virginica  virginica  virginica  virginica  virginica 
## [145] virginica  virginica  virginica  virginica  virginica  virginica 
## Levels: setosa versicolor virginica
## 
## $posterior
##           setosa   versicolor    virginica
## 1   1.000000e+00 3.896358e-22 2.611168e-42
## 2   1.000000e+00 7.217970e-18 5.042143e-37
## 3   1.000000e+00 1.463849e-19 4.675932e-39
## 4   1.000000e+00 1.268536e-16 3.566610e-35
## 5   1.000000e+00 1.637387e-22 1.082605e-42
## 6   1.000000e+00 3.883282e-21 4.566540e-40
## 7   1.000000e+00 1.113469e-18 2.302608e-37
## 8   1.000000e+00 3.877586e-20 1.074496e-39
## 9   1.000000e+00 1.902813e-15 9.482936e-34
## 10  1.000000e+00 1.111803e-18 2.724060e-38
## 11  1.000000e+00 1.185277e-23 3.237084e-44
## 12  1.000000e+00 1.621649e-18 1.833201e-37
## 13  1.000000e+00 1.459225e-18 3.262506e-38
## 14  1.000000e+00 1.117219e-19 1.316642e-39
## 15  1.000000e+00 5.487399e-30 1.531265e-52
## 16  1.000000e+00 1.261505e-27 2.268705e-48
## 17  1.000000e+00 6.754338e-25 3.868271e-45
## 18  1.000000e+00 4.223741e-21 1.224313e-40
## 19  1.000000e+00 1.774911e-22 2.552153e-42
## 20  1.000000e+00 2.593237e-22 5.792079e-42
## 21  1.000000e+00 1.274639e-19 4.357774e-39
## 22  1.000000e+00 1.465999e-20 1.987241e-39
## 23  1.000000e+00 6.569280e-25 7.769177e-46
## 24  1.000000e+00 8.912348e-15 9.178624e-32
## 25  1.000000e+00 1.070702e-15 1.167516e-33
## 26  1.000000e+00 2.497339e-16 5.710269e-35
## 27  1.000000e+00 3.967732e-17 4.378624e-35
## 28  1.000000e+00 1.548165e-21 1.595360e-41
## 29  1.000000e+00 9.271847e-22 6.297955e-42
## 30  1.000000e+00 9.665144e-17 2.977974e-35
## 31  1.000000e+00 2.299936e-16 7.182666e-35
## 32  1.000000e+00 1.975404e-19 2.788334e-38
## 33  1.000000e+00 7.100041e-27 2.216408e-48
## 34  1.000000e+00 1.610295e-28 2.743783e-50
## 35  1.000000e+00 1.205219e-17 1.277245e-36
## 36  1.000000e+00 1.597186e-21 9.033772e-42
## 37  1.000000e+00 1.939869e-24 1.662808e-45
## 38  1.000000e+00 3.310234e-23 7.004971e-44
## 39  1.000000e+00 4.190242e-17 6.991441e-36
## 40  1.000000e+00 1.769359e-20 3.541694e-40
## 41  1.000000e+00 1.063014e-21 2.003866e-41
## 42  1.000000e+00 2.174217e-11 1.213781e-28
## 43  1.000000e+00 1.540753e-18 1.305719e-37
## 44  1.000000e+00 8.940589e-16 1.315511e-32
## 45  1.000000e+00 1.616206e-17 3.205992e-35
## 46  1.000000e+00 1.714743e-16 7.172435e-35
## 47  1.000000e+00 2.083089e-22 2.289783e-42
## 48  1.000000e+00 2.793482e-18 2.629539e-37
## 49  1.000000e+00 2.597560e-23 9.820820e-44
## 50  1.000000e+00 2.322258e-20 4.241757e-40
## 51  1.969732e-18 9.998894e-01 1.105878e-04
## 52  1.242878e-19 9.992575e-01 7.425297e-04
## 53  2.088263e-22 9.958069e-01 4.193053e-03
## 54  2.198898e-22 9.996423e-01 3.576502e-04
## 55  4.213678e-23 9.955903e-01 4.409655e-03
## 56  8.127287e-23 9.985020e-01 1.497982e-03
## 57  3.549900e-22 9.858346e-01 1.416542e-02
## 58  5.007065e-14 9.999999e-01 1.119811e-07
## 59  5.683334e-20 9.998781e-01 1.218649e-04
## 60  1.241039e-20 9.995027e-01 4.973085e-04
## 61  1.956628e-18 9.999986e-01 1.420841e-06
## 62  5.968900e-20 9.992294e-01 7.705716e-04
## 63  2.716128e-18 9.999988e-01 1.220169e-06
## 64  1.184445e-23 9.943267e-01 5.673286e-03
## 65  5.574931e-14 9.999984e-01 1.649215e-06
## 66  2.369511e-17 9.999573e-01 4.268212e-05
## 67  8.429328e-24 9.806471e-01 1.935289e-02
## 68  2.505072e-16 9.999991e-01 9.151716e-07
## 69  1.670352e-27 9.595735e-01 4.042653e-02
## 70  1.341503e-17 9.999967e-01 3.296105e-06
## 71  7.408118e-28 2.532282e-01 7.467718e-01
## 72  9.399292e-17 9.999907e-01 9.345291e-06
## 73  7.674672e-29 8.155328e-01 1.844672e-01
## 74  2.683018e-22 9.995723e-01 4.277469e-04
## 75  7.813875e-18 9.999758e-01 2.421458e-05
## 76  2.073207e-18 9.999171e-01 8.290530e-05
## 77  6.357538e-23 9.982541e-01 1.745936e-03
## 78  5.639473e-27 6.892131e-01 3.107869e-01
## 79  3.773528e-23 9.925169e-01 7.483138e-03
## 80  9.555338e-12 1.000000e+00 1.910624e-08
## 81  1.022109e-17 9.999970e-01 3.007748e-06
## 82  9.648075e-16 9.999997e-01 3.266704e-07
## 83  1.616405e-16 9.999962e-01 3.778441e-06
## 84  4.241952e-32 1.433919e-01 8.566081e-01
## 85  1.724514e-24 9.635576e-01 3.644242e-02
## 86  1.344746e-20 9.940401e-01 5.959931e-03
## 87  3.304868e-21 9.982223e-01 1.777672e-03
## 88  2.034571e-23 9.994557e-01 5.443096e-04
## 89  5.806986e-18 9.999486e-01 5.137101e-05
## 90  5.981190e-21 9.998183e-01 1.816870e-04
## 91  5.878614e-23 9.993856e-01 6.144200e-04
## 92  5.399006e-22 9.980934e-01 1.906591e-03
## 93  3.559507e-18 9.999887e-01 1.128570e-05
## 94  2.104146e-14 9.999999e-01 1.135016e-07
## 95  4.700877e-21 9.996980e-01 3.020226e-04
## 96  1.584328e-17 9.999817e-01 1.826327e-05
## 97  2.802293e-19 9.998892e-01 1.108315e-04
## 98  1.626918e-18 9.999536e-01 4.640488e-05
## 99  7.638378e-11 1.000000e+00 1.867332e-08
## 100 4.679301e-19 9.999269e-01 7.305863e-05
## 101 7.503075e-52 7.127303e-09 1.000000e+00
## 102 5.213802e-38 1.078251e-03 9.989217e-01
## 103 1.231264e-42 2.592826e-05 9.999741e-01
## 104 1.537499e-38 1.068139e-03 9.989319e-01
## 105 6.242501e-46 1.812963e-06 9.999982e-01
## 106 4.209281e-49 6.656263e-07 9.999993e-01
## 107 3.797837e-33 4.862025e-02 9.513797e-01
## 108 1.352176e-42 1.395463e-04 9.998605e-01
## 109 1.323390e-42 2.235313e-04 9.997765e-01
## 110 3.453358e-46 1.727277e-07 9.999998e-01
## 111 5.452660e-32 1.305353e-02 9.869465e-01
## 112 1.182560e-37 1.673875e-03 9.983261e-01
## 113 5.204321e-39 2.006352e-04 9.997994e-01
## 114 1.269953e-40 1.948672e-04 9.998051e-01
## 115 1.685361e-45 1.000455e-06 9.999990e-01
## 116 5.141640e-40 2.605493e-05 9.999739e-01
## 117 1.909820e-35 6.083553e-03 9.939164e-01
## 118 1.207799e-44 1.503799e-06 9.999985e-01
## 119 3.181265e-59 1.317279e-09 1.000000e+00
## 120 1.598511e-33 2.207990e-01 7.792010e-01
## 121 1.119461e-42 6.451865e-06 9.999935e-01
## 122 3.038170e-37 8.272676e-04 9.991727e-01
## 123 6.032879e-50 9.509838e-07 9.999990e-01
## 124 1.951261e-31 9.711942e-02 9.028806e-01
## 125 1.956408e-39 8.836845e-05 9.999116e-01
## 126 1.109337e-36 2.679310e-03 9.973207e-01
## 127 7.841997e-30 1.883675e-01 8.116325e-01
## 128 7.964690e-30 1.342431e-01 8.657569e-01
## 129 6.190641e-44 1.303681e-05 9.999870e-01
## 130 1.406448e-32 1.036823e-01 8.963177e-01
## 131 4.108129e-42 1.442338e-04 9.998558e-01
## 132 1.555697e-36 5.198047e-04 9.994802e-01
## 133 1.320330e-45 3.014091e-06 9.999970e-01
## 134 1.283891e-28 7.293881e-01 2.706119e-01
## 135 1.926560e-35 6.602253e-02 9.339775e-01
## 136 1.271083e-45 2.152818e-06 9.999978e-01
## 137 3.038963e-44 8.881859e-07 9.999991e-01
## 138 4.605973e-35 6.165648e-03 9.938344e-01
## 139 4.538634e-29 1.925262e-01 8.074738e-01
## 140 2.140232e-36 8.290895e-04 9.991709e-01
## 141 6.570902e-45 1.180810e-06 9.999988e-01
## 142 6.202588e-36 4.276398e-04 9.995724e-01
## 143 5.213802e-38 1.078251e-03 9.989217e-01
## 144 1.073945e-45 1.028519e-06 9.999990e-01
## 145 4.048249e-46 2.524984e-07 9.999997e-01
## 146 4.970070e-39 7.473361e-05 9.999253e-01
## 147 4.616611e-36 5.898784e-03 9.941012e-01
## 148 5.548962e-35 3.145874e-03 9.968541e-01
## 149 1.613687e-40 1.257468e-05 9.999874e-01
## 150 2.858012e-33 1.754229e-02 9.824577e-01
## 
## $x
##            LD1          LD2
## 1    8.0617998  0.300420621
## 2    7.1286877 -0.786660426
## 3    7.4898280 -0.265384488
## 4    6.8132006 -0.670631068
## 5    8.1323093  0.514462530
## 6    7.7019467  1.461720967
## 7    7.2126176  0.355836209
## 8    7.6052935 -0.011633838
## 9    6.5605516 -1.015163624
## 10   7.3430599 -0.947319209
## 11   8.3973865  0.647363392
## 12   7.2192969 -0.109646389
## 13   7.3267960 -1.072989426
## 14   7.5724707 -0.805464137
## 15   9.8498430  1.585936985
## 16   9.1582389  2.737596471
## 17   8.5824314  1.834489452
## 18   7.7807538  0.584339407
## 19   8.0783588  0.968580703
## 20   8.0209745  1.140503656
## 21   7.4968023 -0.188377220
## 22   7.5864812  1.207970318
## 23   8.6810429  0.877590154
## 24   6.2514036  0.439696367
## 25   6.5589334 -0.389222752
## 26   6.7713832 -0.970634453
## 27   6.8230803  0.463011612
## 28   7.9246164  0.209638715
## 29   7.9912902  0.086378713
## 30   6.8294645 -0.544960851
## 31   6.7589549 -0.759002759
## 32   7.3749525  0.565844592
## 33   9.1263463  1.224432671
## 34   9.4676820  1.825226345
## 35   7.0620139 -0.663400423
## 36   7.9587624 -0.164961722
## 37   8.6136720  0.403253602
## 38   8.3304176  0.228133530
## 39   6.9341201 -0.705519379
## 40   7.6882313 -0.009223623
## 41   7.9179372  0.675121313
## 42   5.6618807 -1.934355243
## 43   7.2410147 -0.272615132
## 44   6.4144356  1.247301306
## 45   6.8594438  1.051653957
## 46   6.7647039 -0.505151855
## 47   8.0818994  0.763392750
## 48   7.1867690 -0.360986823
## 49   8.3144488  0.644953177
## 50   7.6719674 -0.134893840
## 51  -1.4592755  0.028543764
## 52  -1.7977057  0.484385502
## 53  -2.4169489 -0.092784031
## 54  -2.2624735 -1.587252508
## 55  -2.5486784 -0.472204898
## 56  -2.4299673 -0.966132066
## 57  -2.4484846  0.795961954
## 58  -0.2226665 -1.584673183
## 59  -1.7502012 -0.821180130
## 60  -1.9584224 -0.351563753
## 61  -1.1937603 -2.634455704
## 62  -1.8589257  0.319006544
## 63  -1.1580939 -2.643409913
## 64  -2.6660572 -0.642504540
## 65  -0.3783672  0.086638931
## 66  -1.2011726  0.084437359
## 67  -2.7681025  0.032199536
## 68  -0.7768540 -1.659161847
## 69  -3.4980543 -1.684956162
## 70  -1.0904279 -1.626583496
## 71  -3.7158961  1.044514421
## 72  -0.9976104 -0.490530602
## 73  -3.8352593 -1.405958061
## 74  -2.2574125 -1.426794234
## 75  -1.2557133 -0.546424197
## 76  -1.4375576 -0.134424979
## 77  -2.4590614 -0.935277280
## 78  -3.5184849  0.160588866
## 79  -2.5897987 -0.174611728
## 80   0.3074879 -1.318871459
## 81  -1.1066918 -1.752253714
## 82  -0.6055246 -1.942980378
## 83  -0.8987038 -0.904940034
## 84  -4.4984664 -0.882749915
## 85  -2.9339780  0.027379106
## 86  -2.1036082  1.191567675
## 87  -2.1425821  0.088779781
## 88  -2.4794560 -1.940739273
## 89  -1.3255257 -0.162869550
## 90  -1.9555789 -1.154348262
## 91  -2.4015702 -1.594583407
## 92  -2.2924888 -0.332860296
## 93  -1.2722722 -1.214584279
## 94  -0.2931761 -1.798715092
## 95  -2.0059888 -0.905418042
## 96  -1.1816631 -0.537570242
## 97  -1.6161564 -0.470103580
## 98  -1.4215888 -0.551244626
## 99   0.4759738 -0.799905482
## 100 -1.5494826 -0.593363582
## 101 -7.8394740  2.139733449
## 102 -5.5074800 -0.035813989
## 103 -6.2920085  0.467175777
## 104 -5.6054563 -0.340738058
## 105 -6.8505600  0.829825394
## 106 -7.4181678 -0.173117995
## 107 -4.6779954 -0.499095015
## 108 -6.3169269 -0.968980756
## 109 -6.3277368 -1.383289934
## 110 -6.8528134  2.717589632
## 111 -4.4407251  1.347236918
## 112 -5.4500957 -0.207736942
## 113 -5.6603371  0.832713617
## 114 -5.9582372 -0.094017545
## 115 -6.7592628  1.600232061
## 116 -5.8070433  2.010198817
## 117 -5.0660123 -0.026273384
## 118 -6.6088188  1.751635872
## 119 -9.1714749 -0.748255067
## 120 -4.7645357 -2.155737197
## 121 -6.2728391  1.649481407
## 122 -5.3607119  0.646120732
## 123 -7.5811998 -0.980722934
## 124 -4.3715028 -0.121297458
## 125 -5.7231753  1.293275530
## 126 -5.2791592 -0.042458238
## 127 -4.0808721  0.185936572
## 128 -4.0770364  0.523238483
## 129 -6.5191040  0.296976389
## 130 -4.5837194 -0.856815813
## 131 -6.2282401 -0.712719638
## 132 -5.2204877  1.468195094
## 133 -6.8001500  0.580895175
## 134 -3.8151597 -0.942985932
## 135 -5.1074897 -2.130589999
## 136 -6.7967163  0.863090395
## 137 -6.5244960  2.445035271
## 138 -4.9955028  0.187768525
## 139 -3.9398530  0.614020389
## 140 -5.2038309  1.144768076
## 141 -6.6530868  1.805319760
## 142 -5.1055595  1.992182010
## 143 -5.5074800 -0.035813989
## 144 -6.7960192  1.460686950
## 145 -6.8473594  2.428950671
## 146 -5.6450035  1.677717335
## 147 -5.1795646 -0.363475041
## 148 -4.9677409  0.821140550
## 149 -5.8861454  2.345090513
## 150 -4.6831543  0.332033811

下面对比用判别函数判别的结果与真实的类属的关系：

newG <- res1$class
tab <- table(iris$Species, newG)
knitr::kable(tab)

	setosa	versicolor	virginica
setosa	50	0	0
versicolor	0	48	2
virginica	0	1	49

仅分错了3例，是很难得的。把分错的3例找出来:

knitr::kable(cbind(iris, newG)[iris$Species != newG,])

	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species	newG
71	5.9	3.2	4.8	1.8	versicolor	virginica
84	6.0	2.7	5.1	1.6	versicolor	virginica
134	6.3	2.8	5.1	1.5	virginica	versicolor

7.3 贝叶斯判别

距离判别与Fisher判别都不考虑某一总体实际出现的概率，也不考虑错判后的损失差别。贝叶斯判别可以解决这两个问题，而且方法容易解释。贝叶斯判别属于概率判别法，通过计算属于某类的后验概率，并把待判样品判入后验概率最大的类来进行分类; 也可以指定一个损失函数，把待判样品判入使得总平均损失最小的类。

7.3.1 两总体的贝叶斯判别

设有两个总体\(G_1\)和\(G_2\), 分别有概率密度函数\(f_1(\boldsymbol x)\)和\(f_2(\boldsymbol x)\) （概率质量函数也可以）， \(\boldsymbol x \in \mathcal X \subset \mathbb R^p\)。设\(R_1 \cup R_2 = \mathcal X\), \(R_1 \cap R_2 = \emptyset\)，则\(\boldsymbol x \in R_1\)和\(\boldsymbol x \in R_2\)可以分别作为把\(\boldsymbol x\) 判入\(G_1\)和\(G_2\)的判别准则；反之，判别准则有对应的\(R_1\)和\(R_2\)集合。

7.3.1.1 用平均误判损失最小判别

设随机向量\(\boldsymbol X\)来自\(G_1\)，则\(\boldsymbol X\)被错判入\(G_2\)的概率为 \[\begin{align*} P(2|1) = P(X \in R_2) = \int_{R_2} f_1(\boldsymbol x) \,d\boldsymbol x . \end{align*}\] 类似地，把属于\(G_2\)的观测错判入\(G_1\)的概率为 \[\begin{align*} P(1|2) = \int_{R_1} f_2(\boldsymbol x) \,d\boldsymbol x . \end{align*}\]

而正确判别的概率为 \[\begin{align*} P(1|1) =& \int_{R_1} f_1(\boldsymbol x) \,d\boldsymbol x, \\ P(2|2) =& \int_{R_2} f_2(\boldsymbol x) \,d\boldsymbol x . \end{align*}\]

设一个观测来自\(G_1\)的先验概率为\(p_1\)，来自\(G_2\)的先验概率为\(p_2\)(\(p_2 = 1 - p_1\))，令随机变量\(I\)分布为\(P(I=1)=p_1=1-P(I=2)\)，用\(I=1\)表示观测来自\(G_1\), \(I=2\)表示观测来自\(G_2\)。则 \[\begin{align*} P(\boldsymbol X \text{来自}G_1\text{且被正确判入}G_1) =& P(I=1, \boldsymbol X \in R_1) \\ =& P(I=1) P(\boldsymbol X \in R_1 | I=1) \\ =& p_1 P(1|1) \end{align*}\]

类似推理可得 \[\begin{align*} P(\boldsymbol X \text{来自}G_1\text{且被正确判入}G_1) =& p_1 P(1|1) \\ P(\boldsymbol X \text{来自}G_2\text{且被正确判入}G_2) =& p_2 P(2|2) \\ P(\boldsymbol X \text{来自}G_2\text{但被错误判入}G_1) =& p_2 P(1|2) \\ P(\boldsymbol X \text{来自}G_1\text{但被错误判入}G_2) =& p_1 P(2|1) \end{align*}\] 设\(L(j|i)\)是把一个来自第\(i\)类的样品判入第\(j\)类的损失，则\(L(1|1)=L(2|2)=0\), \(L(2|1)>0, L(1|2)>0\)。

定义平均误判损失(ECM, Expected cost of misclassification)为 \[\begin{align*} \text{ECM}(R_1, R_2) = L(2|1) P(2|1) p_1 + L(1|2) P(1|2) p_2 . \end{align*}\] 判别规则的一种合理选择是选\(R_1\)和\(R_2\)使得平均误判损失最小。这样得到的空间划分为 \[\begin{align*} R_1 =& \left\{ \boldsymbol x \,\left|\, \frac{f_1(\boldsymbol x)}{f_2(\boldsymbol x)} \geq \frac{L(1|2)}{L(2|1)} \cdot \frac{p_2}{p_1} \right.\right\} \\ R_2 =& \left\{ \boldsymbol x \,\left|\, \frac{f_1(\boldsymbol x)}{f_2(\boldsymbol x)} < \frac{L(1|2)}{L(2|1)} \cdot \frac{p_2}{p_1} \right.\right\} \end{align*}\]

上面的判别法则适用于有密度的情形。如果是离散分布，在规则取等号的时候应该单独考虑。

平均误判损失可以表示为 \[\begin{align*} \text{ECM}(R_1, R_2) = \ & \int_{R_2} L(2|1) p_1 f_1(\boldsymbol x) \,d\boldsymbol x \\ + & \int_{R_1} L(1|2) p_2 f_2(\boldsymbol x) \,d\boldsymbol x , \end{align*}\]

为了让两个积分的和最小，应该取积分区域\(R_2\)为第一个被积函数比第二个被积函数小的区域，取积分区域\(R_1\)为第二个被积函数比第一个被积函数小的区域。

7.3.1.2 用后验概率最大判别

在给定\(\boldsymbol X = \boldsymbol x\)后，反推\(\boldsymbol x\)来自哪一类，即求\(P(I=1|\boldsymbol X=\boldsymbol x)\)和\(P(I=2|\boldsymbol X=\boldsymbol x)\)。类似于逆概率公式有如下后验概率 \[\begin{align*} P(I=1|\boldsymbol X=\boldsymbol x) =& \frac{p_1 f_1(\boldsymbol x)}{p_1 f_1(\boldsymbol x) + p_2 f_2(\boldsymbol x)}, \\ P(I=2|\boldsymbol X=\boldsymbol x) =& \frac{p_2 f_2(\boldsymbol x)}{p_1 f_1(\boldsymbol x) + p_2 f_2(\boldsymbol x)} . \end{align*}\] 可以把\(\boldsymbol x\)判入后验概率最大的那一类，这时\(R_1 = \{ \boldsymbol x | p_1 f_1(\boldsymbol x) \geq p_2 f_2(\boldsymbol x) \}\), \(R_2 = \{ \boldsymbol x | p_1 f_1(\boldsymbol x) < p_2 f_2(\boldsymbol x) \}\)。当损失函数\(L(2|1)=L(1|2)\)时，平均错判损失最小法则与后验概率最大法则是一致的。

7.3.1.3 正态总体的贝叶斯判别

贝叶斯判别需要知道\(G_1\)和\(G_2\)的分布。考虑两个总体为多元正态分布的情形。设\(G_1\)为\(\text{N}(\boldsymbol\mu_1, \Sigma_1)\), 设\(G_2\)为\(\text{N}(\boldsymbol\mu_2, \Sigma_2)\)。如果\(\Sigma_1 = \Sigma_2 = \Sigma\)(有共同的协方差阵)，则两个密度为 \[\begin{align*} f_j(\boldsymbol x) = (2\pi)^{-p/2} |\Sigma|^{-1/2} \exp\left\{ -\frac12 (\boldsymbol x - \boldsymbol\mu_j)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol\mu_j) \right\} . \end{align*}\]

最小平均误判损失准则中的\(\frac{f_1(\boldsymbol x)}{f_2(\boldsymbol x)} \geq c\)要求，很容易变成 \[\begin{align*} \frac12(\boldsymbol x - \boldsymbol\mu_2)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol\mu_2) - \frac12(\boldsymbol x - \boldsymbol\mu_1)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol\mu_1) \geq \beta \end{align*}\] 这又等价于（参见§1） \[\begin{align*} (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \left( \boldsymbol x - \frac{\boldsymbol\mu_1 + \boldsymbol\mu_2}{2} \right) \geq \beta, \end{align*}\] 其中 \[\begin{align*} \beta = \log\left( \frac{L(1|2)}{L(2|1)} \cdot \frac{p_2}{p_1} \right) . \end{align*}\]

定义 \[\begin{align*} W(x) = (\boldsymbol\mu_1 - \boldsymbol\mu_2)^T \Sigma^{-1} \left( \boldsymbol x - \frac{\boldsymbol\mu_1 + \boldsymbol\mu_2}{2} \right), \end{align*}\] 令\(R_1 = \{\boldsymbol x | W(\boldsymbol x) \geq \beta \}\), \(R_2 = \{\boldsymbol x | W(\boldsymbol x) < \beta \}\)。

\(W(\boldsymbol x)\)与距离判别和Fisher判别的判别函数形式相同，并且当\(p_1=p_2\), \(L(2|1)=L(1|2)\)时\(\beta=0\), 判别规则也相同。

在\(\Sigma_1 \neq \Sigma_2\)时，可以类似得到判别函数\(W(\boldsymbol x)\)，这时\(W(\boldsymbol x)\)是二次函数。

在未知总体参数时可以从训练样本中估计参数然后代入到判别函数中。

这样的方法是参数方法，也可以用非参数方法估计两个总体的密度。

7.3.2 多总体的贝叶斯判别

设有\(G_1, G_2, \dots, G_k\)共\(k\)个总体， \(G_j\)有密度\(f_j(\boldsymbol x)\)。利用后验概率最大准则，取 \[\begin{align*} R_j = \left\{ \boldsymbol x \left| p_j f_j(\boldsymbol x) = \max_{1 \leq i \leq k} p_i f_i(\boldsymbol x) \right.\right\} . \end{align*}\] 与最小平均错判损失准则相比，后验概率最大准则相当于假定所有错判损失都相等。

7.3.2.1 多元正态多总体的贝叶斯后验概率判别

设\(G_j\)为\(\text{N}(\boldsymbol\mu_j, \Sigma_j)\), \(j=1,2,\dots,k\)。有共同协方差阵时(\(\Sigma_1 = \dots = \Sigma_k = \Sigma\)), 总体\(G_j\)的密度函数为 \[\begin{align*} f_j(\boldsymbol x) = (2\pi)^{-p/2} |\Sigma|^{-1/2} \exp\left\{ -\frac12 (\boldsymbol x - \boldsymbol\mu_j)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol\mu_j) \right\} . \end{align*}\] 令 \[\begin{align*} d_j^2(\boldsymbol x) = (\boldsymbol x - \boldsymbol\mu_j)^T \Sigma^{-1} (\boldsymbol x - \boldsymbol\mu_j) - 2\ln(p_j), \end{align*}\] 称\(d_j^2(\cdot)\)为\(\boldsymbol x\)到第\(j\)类的广义平方距离。则后验概率最大准则的空间划分为 \[\begin{align*} R_j = \left\{ \boldsymbol x | d_j^2(\boldsymbol x) = \min_{i=1,2,\dots,k} d_i^2(\boldsymbol x) \right\}. \end{align*}\]

当各总体协方差阵不相同时，只要修改广义平方距离\(d_j^2(\boldsymbol x)\)定义为 \[\begin{align*} d_j^2(\boldsymbol x) = (\boldsymbol x - \boldsymbol\mu_j)^T \Sigma_j^{-1} (\boldsymbol x - \boldsymbol\mu_j) + \ln|\Sigma_j| - 2\ln(p_j) . \end{align*}\]

总体参数未知时，可以从训练样本中估计。

7.3.2.2 例5.3（有无春旱的判别分析）

P.63例5.3。 14年的有无春旱的数据，用两个综合性气象指标x1, x2作为自变量。 G为组别（1春旱，2无春旱），x1, x2为两个综合预报指标。在R中用lda()可以进行贝叶斯线性判别（假定正态总体且协方差阵相同），用qda()进行贝叶斯二次判别（假定正态总体且协方差阵不同），用prior=可以指定先验概率，不指定先验概率时自动从训练样本中估计先验概率。用后验概率最大准则进行判别。均值、协方差阵从训练样本中估计。

贝叶斯线性判别的程序：

library(MASS)
load('springDrought.RData')
d <- springDrought
lda1 <- lda(
  G ~ x1 + x2,
  prior=c(6,8)/14,
  data=d)

在上面lda()调用中先验分布用数据中旱与不旱的比例来估计。

判别分析结果概况:

print(lda1)

## Call:
## lda(G ~ x1 + x2, data = d, prior = c(6, 8)/14)
## 
## Prior probabilities of groups:
##         1         2 
## 0.4285714 0.5714286 
## 
## Group means:
##         x1        x2
## 1 25.31667 -2.416667
## 2 22.02500 -1.187500
## 
## Coefficients of linear discriminants:
##           LD1
## x1 -0.6312826
## x2  1.0020661

对训练数据集的判别效果:

res1 <- predict(lda1)
newG <- res1$class
tab <- table(d$G, newG); print(tab)

##    newG
##     1 2
##   1 5 1
##   2 0 8

只有一例错判。

print(prop.table(tab))

##    newG
##              1          2
##   1 0.35714286 0.07142857
##   2 0.00000000 0.57142857

cat('正确判别率=', sum(diag(prop.table(tab))), '\n')

## 正确判别率= 0.9285714

显示错判的训练数据:

print(cbind(d, '判别结果'=newG)[d$G != newG, ])

##   G   x1   x2 判别结果
## 4 1 23.5 -1.9        2

每个训练样例属于每组的后验概率:

print(cbind('实际组别'=d$G, res1$post))

##    实际组别            1            2
## 1         1 0.9386546174 6.134538e-02
## 2         1 0.9303445828 6.965542e-02
## 3         1 0.9999448424 5.515761e-05
## 4         1 0.4207076326 5.792924e-01
## 5         1 0.9892508267 1.074917e-02
## 6         1 0.9999925582 7.441831e-06
## 7         2 0.0007277911 9.992722e-01
## 8         2 0.0026045742 9.973954e-01
## 9         2 0.0008227369 9.991773e-01
## 10        2 0.0585597189 9.414403e-01
## 11        2 0.0349605147 9.650395e-01
## 12        2 0.0005620155 9.994380e-01
## 13        2 0.0038092358 9.961908e-01
## 14        2 0.0012325974 9.987674e-01

7.4 类别不平衡问题

在两类的判别中，设两类用\(1\)和\(-1\)标记，称为正例和反例。如果两类的实际比例相差太大，比如正例、反例比例为1:99，则总是预测反例也有99%的精度，但这就失去了建立判别法的意义。利用贝叶斯判别增加先验概率可以部分地缓解这个问题，但是某一类过少还是会使得判别规则难以准确反映实际情况。多类的判别也有类似问题。

考虑正例比例很小的情形。一种解决方法是欠采样法，即删除部分反例，使得两类基本平衡。欠采样法如果随机删除反例，可能会丢失重要信息，典型的EasyEnsemble算法则利用了集成学习机制充分利用这些信息。

比例不平衡的另一种解决方法是过采样法，设法生成更多的正例。但是，不能简单地重复正例，这会造成严重的过度拟合，因为重复的正例会低估正例预测的误差。典型的过采样算法SMOTE利用了插值的方法生成不重复的正例。

7.5 判别效果评价指标

如何判断一种判别方法的效果？以二分类的判别为例，并用医学检验问题作说明。在医学检验中，将有某种病的人正确检出的概率叫做敏感度(sensitivity)，将无病的人正确地判为无病的概率叫做特异度(specificity)。将这样的术语对应到统计假设检验问题： \[ H_0: \text{无病} \longleftrightarrow H_a: \text{有病} \] 第一类错误是\(H_0\)成立但拒绝\(H_0\)的错误，概率为\(\alpha\)；第二类错误是\(H_a\)成立但承认\(H_0\)的错误，概率为\(\beta\)。则敏感度为\(1-\beta\)（检验功效），特异度为\(1-\alpha\)。

称\(\alpha\)为假阳性率(FPR)，等于1减去特异度；称敏感度为真阳性率(TPR)，也称为召回率(recall)。

设图7.1中TN表示真阴性例数，即无病判为无病个数； FN表示假阴性例数，即有病判为无病例数； FP表示假阳性例数，即无病判为有病例数； TP表示真阳性例数，即有病判为有病例数。这时，可以估计 \[\begin{aligned} \text{TPR(真阳性率)} =& 1 - \beta = \text{Sensitivity(敏感度)} = \text{Recall(查全率)} = \frac{TP}{TP + FN} ; \\ \text{FPR(假阳性率)} =& \alpha = \frac{FP}{FP + TN} ; \\ \text{TNR(真阴性率)} =& 1-\alpha = \text{Specificity(特异度)} = \frac{TN}{TN + FP}; \\ \text{FNR(假阴性率)} =& \beta = \frac{FN}{FN + TP} . \end{aligned}\]

另外， \[\begin{aligned} \text{Precision(查准率)} =& \text{Positive predictive value} = 1 - \text{FDR} = \frac{TP}{TP + FP}; \\ \text{FDR(错误发现率)} =& \frac{FP}{FP + TP} = 1 - \text{Precision}. \end{aligned}\]

F1指标是查准率与查全率的如下“调和平均”： \[ \frac{1}{\text{F1}} = \frac{1}{2} \left( \frac{1}{\text{Precision}} + \frac{1}{\text{Recall}} \right) . \]

图7.1: 医学诊断四种情况表格

在判别分析中，以\(-1\)表示无病，称为反例或者阴性， \(1\)表示有病，称为正例或者阳性。设判别函数为\(f(\boldsymbol x)\)，任意取分界值\(t\)，都可以按\(f(\boldsymbol x) > t\)和\(f(\boldsymbol x) < t\)将因变量值（标签）判为1或者-1。对每个\(t\)，都可以在训练集上计算对应的假阳性率（健康人判入有病的比例）和真阳性率（病人判为有病的概率）。 \(t\)值取得越大，相当于检验水平\(\alpha\)取值越小，越不容易判为有病，这时假阳性率低，真阳性率也低； \(t\)值取得越小，越容易判为有病，相当于检验水平\(\alpha\)取值很大，这时假阳性率很高，真阳性率也很高。最好是对比较大的分界点\(t\)，当假阳性率很低时，真阳性率就很高。以给定\(t\)后的假阳性率为横坐标，真阳性率为纵坐标作曲线，称为ROC(Receiver Operating Curve，受试者工作特征)曲线。曲线越靠近左上角越好。可以用ROC曲线下方的面积大小度量判别方法的优劣，面积越大越好。记这个面积为AUC(Area Under the Curve)。

R的扩展包ROCR可以用来进行ROC分析。其中的prediction(predicitons, labels, label.ordering)三个参数在二值判别中：

predictions为训练集每个样本点的判别函数值（不能是类别）；
labels是训练集每个样本点的真实类别，是因子类型；
label.ordering将labels的两个水平指定次序，分别对应于代表\(-1\)和代表\(+1\)的两个因变量水平。

用performance()函数生成ROC曲线所需的坐标，用plot()函数作图。

例如，对iris数据集，以versicolor为1， setosa和virginica为-1，用Fisher线性判别法，做ROC曲线：

library(MASS)
library (ROCR)
iris2 <- iris
iris2[["Species"]] <- ifelse(
  iris[["Species"]] == "versicolor", 1, -1)
lda1 <- lda(
  Species ~ Sepal.Length + Sepal.Width + 
    Petal.Length + Petal.Width, data=iris2)
pred1 <- predict(
  lda1, newdata=iris2,
  decision.values = TRUE)
decf1 <- pred1$posterior[,"1"]
predobj <- prediction(
  decf1, iris2[["Species"]],
  label.ordering=c("-1", "1"))
perf1 <- performance(predobj, "tpr", "fpr")
auc1 <- performance(
       predobj, "auc")@y.values[[1]]
plot(perf1, main="iris中versicolor用Fisher线性判别的ROC",
     sub=paste("AUC=", auc1))

plotROC包可以利用ggplot2制作ROC曲线。

7.6 案例分析与R实现：各省经济增长差异判别分析

P.65， 1994年30个省、直辖市、自治区影响经济增长差异的四个制度变量数据，用来判定各省的经济增长组别，有两组。 G为组别（1较好，2较差），x1, x2, x3, x4为自变量。 x1是经济增长率（%）， x2是非国有化水平（%）， x3是开放度（%）， x4是市场化程度（%）。分别用两总体距离判别法、Fisher判别法和贝叶斯判别法进行判别分析，对留出的三个省（28–30号观测）用建立的判别规则进行检验判别。（原始数据中这三个观测也录入了类别，实际应该是缺失值）。

载入训练数据和待判样本：

load('provincesDevelope.RData')
d <- provincesDevelope[1:27,] # 训练样本
d.test <- provincesDevelope[28:30,] # 待判样本

7.6.1 (1) 距离判别法

分别计算两个组的协方差阵:

S1 <- var(d[d$G==1,c('x1', 'x2', 'x3', 'x4')])
S2 <- var(d[d$G==2,c('x1', 'x2', 'x3', 'x4')])

组1的协方差阵:

print(S1)

##          x1        x2        x3       x4
## x1 11.10055  26.82577  29.48864  8.30330
## x2 26.82577 115.07050  62.42816 47.93081
## x3 29.48864  62.42816 452.02607 39.19863
## x4  8.30330  47.93081  39.19863 51.32260

组2的协方差阵:

print(S2)

##           x1        x2         x3        x4
## x1  9.023833  22.09105   4.138192   2.67120
## x2 22.091050 276.71312 -38.842827  80.50661
## x3  4.138192 -38.84283  35.373310 -12.10012
## x4  2.671200  80.50661 -12.100123  72.85081

协方差阵差别较大，用不同协方差阵的二次判别函数：

mu1 <- colMeans(d[d$G==1,c('x1', 'x2', 'x3', 'x4')])
mu2 <- colMeans(d[d$G==2,c('x1', 'x2', 'x3', 'x4')])
W1 <- function(x, mu1, mu2, Sigma1, Sigma2){
    0.5*(mahalanobis(x, mu2, Sigma2)
         - mahalanobis(x, mu1, Sigma1))
}

用以上二次判别函数对训练样本判别：

w.train <- W1(as.matrix(d[,c('x1', 'x2', 'x3', 'x4')]),
               mu1, mu2, S1, S2)
newG.train <- ifelse(w.train >= 0, 1, 2)

二次距离判别的训练结果:

tab <- table(d$G, newG.train); print(tab)

##    newG.train
##      1  2
##   1 10  1
##   2  0 16

错判的训练样例

print(cbind(d, '错判入'=newG.train)[d$G != newG.train,])

##    province G x1    x2    x3    x4 错判入
## 10   广  西 1 16 57.11 12.57 60.91      2

训练样本只错判了广西。

对三个待判省的应用：

w.test <- W1(as.matrix(d.test[,c('x1', 'x2', 'x3', 'x4')]),
               mu1, mu2, S1, S2)
newG.test <- ifelse(w.test >= 0, 1, 2)
print(cbind(d.test[,'province'], '判别类'=newG.test))

##             判别类
## 28 "江  苏" "1"   
## 29 "安  徽" "2"   
## 30 "陕  西" "2"

江苏1，安徽、陕西2。

7.6.1.1 距离判别法：用qda作距离判别

在贝叶斯判别法选择相等先验和相等误判损失时，就是距离判别法。

res1 <- qda(
  G ~ x1 + x2 + x3 + x4,
  prior=c(0.5, 0.5), # 相等先验
  data=d)
print(res1)

## Call:
## qda(G ~ x1 + x2 + x3 + x4, data = d, prior = c(0.5, 0.5))
## 
## Prior probabilities of groups:
##   1   2 
## 0.5 0.5 
## 
## Group means:
##         x1       x2        x3     x4
## 1 15.73636 65.02818 25.149091 74.350
## 2 11.56250 40.10625  9.228125 58.105

用qda()作二次距离判别对训练数据的判别效果：

newG.train <- predict(res1)$class
tab <- table(d$G, newG.train); print(tab)

##    newG.train
##      1  2
##   1 10  1
##   2  0 16

错判的训练样例：

print(cbind(d, '错判入'=newG.train)[d$G != newG.train,])

##    province G x1    x2    x3    x4 错判入
## 10   广  西 1 16 57.11 12.57 60.91      2

只错判了广西。

7.6.2 (2) Fisher线性判别法

用lda()作Fisher线性判别。不用先验。

lda1 <- lda(
  G ~ x1 + x2 + x3 + x4,
  prior=rep(1/2, 2), # 不用先验
  data=d)
print(lda1)

## Call:
## lda(G ~ x1 + x2 + x3 + x4, data = d, prior = rep(1/2, 2))
## 
## Prior probabilities of groups:
##   1   2 
## 0.5 0.5 
## 
## Group means:
##         x1       x2        x3     x4
## 1 15.73636 65.02818 25.149091 74.350
## 2 11.56250 40.10625  9.228125 58.105
## 
## Coefficients of linear discriminants:
##            LD1
## x1 -0.06034498
## x2 -0.01661878
## x3 -0.02532111
## x4 -0.08078449

用lda()作Fisher判别的训练结果:

newG.train <- predict(lda1)$class
tab <- table(d$G, newG.train); print(tab)

##    newG.train
##      1  2
##   1 10  1
##   2  0 16

错判的训练样例：

print(cbind(d, '错判入'=newG.train)[d$G != newG.train,])

##    province G x1    x2    x3    x4 错判入
## 10   广  西 1 16 57.11 12.57 60.91      2

只错判了广西。

Fisher线性判别对三个待判省：

newG.test <- predict(lda1, newdata=d.test)$class
print(cbind(d.test[,'province'], '判别类'=newG.test))

##               判别类
## [1,] "江  苏" "1"   
## [2,] "安  徽" "1"   
## [3,] "陕  西" "2"

江苏、安徽1，陕西2，线性判别法不如二次判别法。

7.6.3 (3) 贝叶斯判别法（线性）

用lda()作贝叶斯线性判别，不输入先验，不指定先验时自动从训练样本中估计。

lda2 <- lda(
  G ~ x1 + x2 + x3 + x4,
  data=d)
print(lda2)

## Call:
## lda(G ~ x1 + x2 + x3 + x4, data = d)
## 
## Prior probabilities of groups:
##         1         2 
## 0.4074074 0.5925926 
## 
## Group means:
##         x1       x2        x3     x4
## 1 15.73636 65.02818 25.149091 74.350
## 2 11.56250 40.10625  9.228125 58.105
## 
## Coefficients of linear discriminants:
##            LD1
## x1 -0.06034498
## x2 -0.01661878
## x3 -0.02532111
## x4 -0.08078449

训练结果:

newG.train <- predict(lda2)$class
tab <- table(d$G, newG.train); print(tab)

##    newG.train
##      1  2
##   1 10  1
##   2  0 16

错判的训练样例：

print(cbind(d, '错判入'=newG.train)[d$G != newG.train,])

##    province G x1    x2    x3    x4 错判入
## 10   广  西 1 16 57.11 12.57 60.91      2

人为输入估计的先验。应该与lda2相同:

lda3 <- lda(
  G ~ x1 + x2 + x3 + x4,
  prior=c(11,16)/27,
  data=d)
print(lda3)

## Call:
## lda(G ~ x1 + x2 + x3 + x4, data = d, prior = c(11, 16)/27)
## 
## Prior probabilities of groups:
##         1         2 
## 0.4074074 0.5925926 
## 
## Group means:
##         x1       x2        x3     x4
## 1 15.73636 65.02818 25.149091 74.350
## 2 11.56250 40.10625  9.228125 58.105
## 
## Coefficients of linear discriminants:
##            LD1
## x1 -0.06034498
## x2 -0.01661878
## x3 -0.02532111
## x4 -0.08078449

newG.train <- predict(lda3)$class
tab <- table(d$G, newG.train); print(tab)

##    newG.train
##      1  2
##   1 10  1
##   2  0 16

print(cbind(d, '错判入'=newG.train)[d$G != newG.train,])

##    province G x1    x2    x3    x4 错判入
## 10   广  西 1 16 57.11 12.57 60.91      2