多元统计之判别分析

发表于： 更新于： 分类于：
字数： 1377 阅读：≈ 7分钟 浏览： 评论：

概述

  • 判别分析是判别样品所属类型的一种分析方法，是在分类确定的条件下，根据某一研究对象的各种特征值判别其类型归属问题的一种多变量统计分析方法。
  • 判别分析于聚类分析的功能差不多，区别在于，聚类分析之前，没有人知道具体的是怎么分的类，分了哪几大类。而判别分析是已经把类别给分好，要做的是把没有分好类的数据观测，按照之前分好的类再进行分类。这里不同于生活中常见的分类先有具体的分类逻辑（这里叫做判别函数）。所以判别分的难点在于先由分好类的数据观测找到一个或者多个判别函数，然后对未进行分类的观测按照该判别公式进行分类。

分类

  • 按判别的总体数来区分: 两个总体判别分析和多总体判别分析
  • 按区分不同总体所用的数学模型来分: 线性判别和非线性判别
  • 常用的几种判别分析方法: 距离判别法、Fisher判别法、Bayes判别法和逐步判别法
  • 按判别时所处理的变量方法不同: 逐步判别和序贯判别等

距离判别法

马氏距离

  • 定义点到总体G的马氏距离为 $$D^2(X, G) = (X - \mu)’\Sigma^{-1}(X - \mu)$$
    • $当\Sigma = I(单位矩阵)时, 即为欧氏距离的情形$

距离判别的思想和方法

两个总体的距离判别

  • 问题:
    • 设有协方差矩阵相等的两个总体$G_1$和$G_2$, 其均值分别是$\mu_1$和$\mu_2$, 对于一个新的样品X, 判断它来自哪个总体
  • 思路:
    • 求新样品到两个总体的马氏距离的差值$W(X)$, 若$W(X) \le 0, 则X \in G_1; 反之X \in G_2$, 其中$\bar{\mu} = \frac{\mu_1 + \mu_2}{2}$是两个总体均值的平均值, $\alpha = \Sigma^{-1}(\mu_1 - \mu_2)$
    $$
\begin{aligned}
  W(X) & = D^2(X, G_1) - D^2(X, G_2) \\
       & = -2(X - \frac{\mu_1 + \mu_2}{2})'\Sigma^{-1}(\mu_1 - \mu_2) \\
       & = -2(X - \bar{\mu})'\alpha  \\
       & = -2\alpha'(X - \bar{\mu})
\end{aligned} \tag{2.2.1-1}
$$
    • 记$W(X) = \alpha’(X - \bar{\mu})$, 则判别规则
    $$
\begin{cases}
    X \in G_1, \quad 如果W(X) \ge 0 \\
    X \in G_2, \quad 如果W(X) < 0 \\
\end{cases} \tag{2.2.1-2}
$$
    这里称W(X)为两总体距离判别的判别函数, 由于它是X的线性函数, 故又称为线性判别函数, $\alpha$称为判别系数
实际应用
  • 在实际应用中, 总体的均值和协方差矩阵一般是未知的, 可由样本均值和样本协方差矩阵分别进行估计。
  • $设X_1^{(1)},\cdots,X_{n_1}^{(1)}$来自总体$G_1$的样本, $设X_2^{(2)},\cdots,X_{n_2}^{(2)}$是来自总体$G_2$的样本, $\mu_1$和$\mu_2$的一个无偏估计分别为 $$\bar{X}^{(1)} = \frac{1}{n_1}\sum_{i=1}^{n_1} X_i^{(1)} 和 \bar{X}^{(1)} = \frac{1}{n_2}\sum_{i=1}^{n_2} X_i^{(2)}$$ $\Sigma$的一个联合无偏估计为 $$\hat{\Sigma} = \frac{1}{n_1 + n_2 -2}(S_1 + S_2)$$ 其中$$S_\alpha = \sum_{1}^{n_\alpha}(X_i^{(\alpha)} - \bar{X}^{(\alpha)})(X_i^{(\alpha)} - \bar{X}^{(\alpha)})', \alpha=1,2\\$$ 此时, 两总体距离判别的判别函数为$$\hat{W}(X) = \hat{\alpha}'(X - \bar{X});$$其中$$\bar{X} = \frac{\bar{X}^{(1)} + \bar{X}^{(2)}}{2}, \hat{\alpha} = \hat{\Sigma}^{-1}(\bar{X}^{(1)} - \bar{X}^{(2)});$$判别规则为
    $$
  \begin{cases}
    X \in G_1, \quad 如果\hat{W}(X) \ge 0 \\
    X \in G_2, \quad 如果\hat{W}(X) < 0 \\
  \end{cases}
$$
  • 示例
## state_train.csv ##
# state,life expectancy,literacy ,class
# USA,76,99,1
# Japan,79.5,99,1
# Switzerland,78,99,1
# Argentina,72.1,95.9,1
# UAE,73.8,77.7,1
# Bulgaria,71.2,93,0
# Cuba,75.3,94.9,0
# Paraguay,70,91.2,0
# Georgia,72.8,99,0
# South Africa,62.9,80.6,0

## state_test.csv ##
# state,life expectancy,literacy 
# China,68.5,79.3
# Romania,69.9,96.9
# Greece,77.6,93.8
# Columbia,69.3,90.3

train = read.csv("state_train.csv", sep = ",", header = T)
test = read.csv("state_test.csv", sep = ",", header = T)
t_ed = subset(train, train$class==1)[,-4]
t_ing = subset(train, train$class==0)[,-4]
# class=1
u1 = apply(t_ed[,-1], 2, mean)
s1 = cov(t_ed[,-1])
# calss=0
u2 = apply(t_ing[,-1], 2, mean)
s2 = cov(t_ing[,-1])
# 联合无偏估计
s = (s1*4 + s2*4) / (5+5-2)
ubar = (u1 + u2)/2
# 判别系数
alpha = solve(s)%*%(u1 - u2)

# 判别函数
discriminate = function(test) {
  class = NULL
  w.value = NULL
  n = dim(test)[1]
  state = test$state
  data = as.matrix(test[,-1])
  
  for (i in 1:n) {
    # 读取每一行的数据
    x = data[i,]
    w = t(alpha)%*%(x - ubar)
    w.value[i] = w
    if (w < 0) {
      # 属于第二类
      class[i] = 2
    } else {
      class[i] = 1
    }
  }
  result = data.frame(state, w.value, class)
  return(result)
}

discriminate(test)
discriminate(t_ed)
discriminate(t_ing)

多个总体的距离判别

  • 问题:

    • 设有k个总体$G_1,G_2,\cdots,G_k$, 其均值和协方差矩阵分别是$\mu_1,\mu_2,\cdots,\mu_k$和$\Sigma_1,\Sigma_2,\cdots,\Sigma_k$, 而且$\Sigma_1 = \Sigma_2 = \cdots = \Sigma_k = \Sigma$。对于一个新的样品X, 要判断它来自哪个总体

  • 思路:

    • 与两个总体的距离判别问题解决思路一样。计算新样品X到每一个总体的距离, 即
    $$
  \begin{aligned}
    D^2(X, G_\alpha) & = (X - \mu_\alpha)'\Sigma^{-1}(X - \mu_\alpha) \\
                     & = X'\Sigma^{-1}X - \mu_\alpha'(X\Sigma^{-1})' -\mu_\alpha'\Sigma^{-1}X + \mu'_\alpha\Sigma^{-1}\mu_\alpha\\
                     & = X'\Sigma^{-1}X - 2\mu'_\alpha\Sigma^{-1}X + \mu'_\alpha\Sigma^{-1}\mu_\alpha \\
                     & = X'\Sigma^{-1}X - 2(I'_\alpha X + C_\alpha)
  \end{aligned}
$$

    其中$$I_\alpha = \Sigma^{-1}\mu_\alpha, C_\alpha = -\frac{1}{2}\mu'_\alpha\Sigma^{-1}\mu_\alpha, \alpha=1,2,\cdots,k。$$

    • 线性判别函数为:
    $$W_\alpha(X) = I'_\alpha X + C_\alpha, \quad \alpha=1,2,\cdots,k$$

    相应的判别规则为

    $$X \in G_i \quad 如果W_i(X) = \max\limits_{1 \le \alpha le k}(I'_\alpha X + C_a)$$

实际应用

  • 针对实际问题, 当$\mu_1,\mu_2,\cdots,\mu_k$和$\Sigma$均未知时, 可以通过相应的样本值来替代
  • 设$X_1^{(\alpha)},\cdots,X_{n_\alpha}^{\alpha}$是来自总体$G_\alpha$中的样本$(\alpha=1,2,\cdots,k)$,则$\mu_\alpha(1,2,\cdots,k)和\Sigma$可估计为
$$
\bar{X}^{(\alpha)} = \frac{1}{n_\alpha}\sum_{i=1}^{n_\alpha}X_i^{(\alpha)}, \alpha=1,2,\cdots,k \\
\hat{\Sigma} = \frac{1}{n - k}\sum_{\alpha=1}^kS_\alpha,n=n_1 + n_2 + \cdots + n_k, \\
S_\alpha =\sum_{i=1}^{n_\alpha}(X_i^{(\alpha)} - \bar{X}^{(\alpha)})(X_i^{(\alpha)} - \bar{X}^{(\alpha)})',\alpha = 1,2,\cdots,k
$$
  • 如果总体$G_1,G_2,\cdots,G_k$的协方差矩阵分别是$\Sigma_1,\Sigma_2,\cdots,\Sigma_k$, 而且它们不全相等, 则计算X到各总体的马氏距离,即 $$D^2(D,G_\alpha) = (X - \mu_\alpha)’\Sigma_\alpha^{-1}(X - \mu_\alpha) \quad \alpha=1,2,\cdots,k$$ 判别规则为
$$X \in G_i \quad 如果D^2(X, G_i) = \min\limits_{1 \le \alpha \le k}D^2(X, G_\alpha)$$
当$\mu_1,\mu_2,\cdots,\mu_k和\Sigma_1,\Sigma_2,\cdots,\Sigma_k$均未知时, $\mu_\alpha(\alpha=1,2,\cdots,k)$的估计同前, $\Sigma_\alpha(\alpha=1,2,\cdots,k)$的估计为
$$\hat{\Sigma}_\alpha = \frac{1}{n_\alpha - 1}S_\alpha, \alpha=1,2,\cdots,k$$
  • 示例
x = data.frame(iris)

# 计算参数
sps = as.vector(unique(x$Species))
for (i in 1:3) {
  # 分类
  assign(paste("g", i, sep = ""), subset(x, x$Species==sps[i])[,-5])
  group = get(paste("g", i, sep = ""))
  # 计算均值
  assign(paste("u", i, sep = ""), apply(group, 2, mean))
  # 协方差矩阵
  assign(paste("s", i, sep = ""), cov(group))
}
# 联合无偏估计
s = (s1 + s2 + s3)*49/(150 - 3)
for (i in 1:3) {
  ubar = get(paste("u", i, sep = ""))
  assign(paste("I", i, sep = ""), solve(s)%*%ubar)
  assign(paste("C", i, sep = ""), t(ubar)%*%solve(s)%*%ubar*(-1/2))
}

# 判别函数
discriminant = function(test) {
  class = NULL
  n = dim(test)[1]
  data = as.matrix(test[,-5])
  # 待分类的n个样品
  for (i in 1:n) {
    df = data[i,]
    w = NULL
    # 计算每个样品到每一个总体的距离
    for (j in 1:3) {
      iv = get(paste("I", j, sep = ""))
      cv = get(paste("C", j, sep = ""))
      w[j] = t(iv)%*%df + cv
    }
    class[i] = sps[which.max(w)]
  }
  result = data.frame(test$Species, class)
  return(result)
}

result = discriminant(x)
table(result)

贝叶斯判别法

  • 在判别分析时, 考虑先验概率, 并利用Bayes公式导出后验概率(posterior probability), 即各个样品属于每一类的概率。判别的准则是按后验概率大小归类。

基本思想

  • 问题:
    • 问题：设有k个总体$G_1,G_2,\cdots,G_k$，其各自的分布密度函数$f_1(x),f_2(x),\cdots,G_k$互不相同的，假设个总体各自出现的概率分别为$q_1,q_2,\cdots,q_k$（先验概率），$q > 0$，$\sum_{i=1}^kq_i=1$。假设已知若将本来属于$G_i$总体的样品错判到总体$G_j$时造成的损失为$C(j|i)，i,j=1,2,\cdots,k$。在这样的情形下，对于新的样品判断其来自哪个总体。
  • 思路
    • 设k个总体$G_1,G_2,\cdots,G_k$相应的p维样本空间为$R_1,R_2,\cdots.R_k$, 即为一个划分, 故我们可以简记一个判别规则为$R = R(R_1,R_2,\cdots,R_k)$. 在规则R下, 将属于$G_i$的样品错判到$G_j$的概率为 $$P(j|i,R) = \int_{R_j}f_i(x)dx \quad i,j=1,2,\cdots,k \quad i \neq j$$ 此概率值等于总体$G_i$的分布密度在区域$R_j$上的面积.
    • 判别规则R对总体$G_i$而言样品错判后所造成的平均损失为 $$r(i|R) = \sum_{j=1}^k[C(j|i)P(j|i,R)] \quad i=1,2,\cdots,k$$ 其中$C(i|i) = 0$
    • 由于$k$个总体$G_1,G_2,\cdots,G_k$出现的先验概率分别为$q_1,q_2,\cdots,q_k$, 则用规则R来进行判别所造成的总平均损失为 $$\begin{aligned} g(R) & = \sum_{i=1}^kq_ir(i,R) \ & = \sum_{i=1}^kq_i\sum_{j=1}^kC(j|i)P(j|i,R) \end{aligned} $$
    • Bayes判别法则, 就是要选择$R_1,R_2,\cdots,R_k$, 使得g(R)达到极小.

基本方法

  • 已知样品X来自总体$G_i$的先验概率为$q_i,i=1,2,\cdots,k$, 则在规则R下, 误判的总平均损失为
    $$
\begin{aligned}
  g(R) & = \sum_{i=1}^kq_i\sum_{j=1}^kC(j|i)P(j|i,R) \\
       & = \sum_{i=1}^kq_i\sum_{j=1}^kC(j|i)\int_{R_j}f_i(x)dx \\
       & = \sum_{j=1}^k \int_{R_j}(\sum_{i=1}^kq_iC(j|i)f_i(x))dx \\
       & = \sum_{j=1}^k\int_{R_j}h_j(x)dx \\
\end{aligned}
$$
    $$其中, \quad h_j(x) = \sum_{i=1}^kq_iC(j|i)f_i(x) \tag{3.2.1}$$
  • 如果空间$R^p$有另一种划分$$R^* = (R_1^*, R_2^*, \cdots, R_k^*)$$, 则它的总平均损失为
    $$g(R^*) = \sum_{j=1}^k\int_{R_j^*}h_j(x)dx$$
    则两种划分下的总平均损失之差为
    $$g(R) - g(R^*) = \sum_{i=1}^k\sum_{j=1}^k\int_{R_i \cap R_i^*}[h_i(x) - h_j(x)]dx \tag{3.2.2}$$
  • 由$R_i$的定义, 在$R_i$上$h_i(x) \le h_j(x)$对于一切j成立, 故3.2.2式小于或等于零, 说明$R_1,R_2,\cdots,R_k$确能使总平均损失达到极小, 它是Bayes判别的解.
  • 这样, 我们以Bayes判别的思想得到的划分$R= (R_1,R_2,\cdots,R_k)$为 $$R_i={x|h_i(x) = \min \limits_{1 \le j \le k}h_j(x)} \quad i=1,2,\cdots,k$$

性质

  • 当k=2时, 由3.2.1可得
    $$
h_1(x) = q_2C(1|2)f_2(x) \quad h_2(x) = q_1C(2|1)f_1(x) \\
令 V(x) = \frac{f_1(x)}{f_2(x)}, \quad d = \frac{q_2C(1|2)}{q_1C(2|1)}
$$
    判别规则表示为
    $$
  \begin{cases}
    x \in G_1, \quad 当V(x) \ge d \\
    x \in G_2 , \quad 当V(x) < d
  \end{cases} \tag{3.3.1}
$$
    若$f_1(x)$与$f_2(x)$分别为$N(\mu_1, \Sigma)$和$N(\mu_2, \Sigma)$, 那么
    $$
  \begin{aligned}
    V(x) & = \frac{f_1(x)}{f_2(x)}  \\
         & = exp\{-\frac{1}{2}(x - \mu_1)'\Sigma^{-1}(x - \mu_1) + \frac{1}{2}(x - \mu_2)'\Sigma^{-1}(x - \mu_2)\}  \\
         & = exp\{[x - (\mu_1 + \mu_2)/2]'\Sigma^{-1}(\mu_1 - \mu_2)\} \\
         & = exp W(x)
  \end{aligned}
$$
    判定样品X来自该总体时, 判别规则为
    $$
  \begin{cases}
    x \in G_1, \quad 当W(x) \ge lnd \\
    x \in G_2 , \quad 当W(x) < lnd
  \end{cases} \tag{3.3.2}
$$
    其中$W(x)$由2.2.1-1定义

注意

  • 在实际应用中要确定先验概率是很困难的事
  • 可将样本的各类样本含量之构成比作为先验概率的估计
  • 在无法确定时, 可取各类的先验概率相等, 此时Bayes判别失去其优越性

示例 

# bayes判别
# prior 指定先验概率
ld = lda(iris[,1:4], iris[,5], prior = c(0.2, 0.3, 0.5))
ld
Z = predict(ld)
newG = Z$class
tab = table(newG, Species)
sum(diag(prop.table(tab)))

费歇(Fisher)判别法

基本思想

  • 通过将多维数据投影到某个方向上，投影的原则是将总体与总体之间尽可能的放开，然后再选择合适的判别规则，将新的样品进行分类判别.
  • 从$k$个总体中抽取具有$p$个指标的样品观测数据，借助方差分析的思想构造一个线性判别函数 $$U(X) = u_1X_1 + u_2X_2 + \cdots + u_pX_p \tag{4.1-1}$$ 其中系数$u = (u_1,u_2,\cdots,u_p)’$确定的原则是使得总体之间区别最大，而使每个总体内部的离差最小。

判别函数的构造

针对两个总体的情况

  • 假设有两个总体$G_1,G_2$，其均值分别为$\mu_1$和$\mu_2$，协方差矩阵为$\Sigma_1$和$\Sigma_2$。当$X \in G_i$时，我们可以求出$u’X$的均值和方差，即
    $$
E(u'X) = E(u'X | G_i) = u'E(X|G_i) = u'\mu_i \triangleq \bar{u}_i, \quad i=1,2 \\
D(u'X) = D(u'X | G_i) = u'D(X|G_i)u = u'\Sigma_iu \triangleq \sigma_i^2, \quad i=1,2 \\
$$
  • 在求线性判别函数时，尽量使得总体之间差异大，也就是要求$u’\mu_1 - u’\mu_2$尽可能的大，即$\bar{\mu}_1 - \bar{\mu}_2$变大；同时要求每一个总体内的离差平方和最小，即$\sigma_1^2 + \sigma_2^2$，则我们可以建立一个目标函数 $$\varPhi(u) = \frac{\bar{u}_1 - \bar{u}_2}{\sigma_1^2 + \sigma_2^2} \tag{4.2.1-2}$$ 这样，我们就将问题转化为，寻找$u$使得目标函数$\varPhi$达到最大。从而可以构造出所要求的线性判别函数。

针对多个总体的情况

  • 假设有k个总体$G_1,G_2,\cdots,G_k$，其均值和协方差矩阵分别为$\mu_i$和$\Sigma_i(>0)$。当$X \in G_i$时，我们可以求出$u’X$的均值和方差，即
    $$
E(u'X) = E(u'X | G_i) = u'E(X|G_i) = u'\mu_i, \quad i=1,2,\cdots,k \\
D(u'X) = D(u'X | G_i) = u'D(X|G_i)u = u'\Sigma_iu, \quad i=1,2,\cdots,k \\
$$
    令
    $$
b = \sum_{i=1}^k(u'\mu_i - u'\bar{\mu})^2 \\
e = \sum_{i=1}^ku'\Sigma_iu = u'(\sum_{i=1}^k\Sigma_i)u = u'Eu \\
其中\bar{u} = \frac{1}{k}\sum_{i=1}^k\mu_i, \quad E = \sum_{i=1}^k\Sigma_i
$$
    b相当于一元方差分析中的组间差, e相当于组内差, 应用方差分析的思想, 选择u使得目标函数 $$\varPhi(u) = \frac{b}{e} \tag{4.2.2-1}$$ 达到极大.
    • 线性判别函数$u’X$, 对于一个新的样品$X$可以构造一个判别规则, 如果 $$|u’X - u’\mu_j| = \min \limits_{1 \le i \le k}|u’X - u’\mu_i| \tag{4.2.2-2}$$ 则判定$X$来自总体$G_j$

线性判别函数的求法

  • 针对多个总体的情形, 设$X$为$p$维空间的样品, 那么 $$\bar{\mu} = \frac{1}{k}\sum_{i=1}^k\mu_i = \frac{1}{k}M'1$$ 其中
    $$ M = 
\begin{pmatrix}
  \mu_{11} & \mu_{21} & \cdots & \mu_{p1} \\
  \mu_{12} & \mu_{22} & \cdots & \mu_{p2} \\
  \vdots   & \vdots   & \ddots & \vdots   \\
  \mu_{1k} & \mu_{2k} & \cdots & \mu_{pk}
\end{pmatrix} = 
\begin{pmatrix}
  \mu'_1 \\
  \mu'_2 \\
  \vdots \\
  \mu'_k
\end{pmatrix} \quad 1 =
\begin{pmatrix}
  1 \\
  1 \\
  \vdots \\
  1
\end{pmatrix}
$$

注意到

$$M'M = 
  \begin{pmatrix}
    \mu_1 & \mu_2 & \cdots & \mu_k
  \end{pmatrix}
  \begin{pmatrix}
    \mu'_1 \\
    \mu'_2 \\
    \vdots \\
    \mu'_k
  \end{pmatrix} = \sum_{i=1}^k\mu_i\mu'_i
$$

从而

$$
  \begin{aligned}
    b & = \sum_{i=1}^k(u'\mu_i - u'\bar{\mu})^2 \\
      & = u'\sum_{i=1}^k(\mu_i - \bar{\mu})(\mu_i - \bar{\mu})u \\
      & = u'[\sum_{i=1}^k\mu_i\mu'_i - k\bar{\mu}\bar{\mu}']u \\
      & = u'(M'M - \frac{1}{k}M'11'M)u  \\
      & = u'M'(I - \frac{1}{k}J)Mu \\
      & = u'Bu
  \end{aligned}
$$

这里,

$$
B=M'(I - \frac{1}{k}J)M, \quad I_{p \times p}为p \times p的单位阵, \quad J = \begin{pmatrix}
  1 & \cdots & 1  \\
    &  \ddots   \\
  1 & \cdots & 1
\end{pmatrix}
$$

即有 $$\varPhi(u) = \frac{u’Bu}{u’Eu} \tag{4.3-1}$$ 求得使$(4.3-1)$式达到极大的$u$.

  • 为确保解的唯一性, 不妨设$u’Eu = 1$, 转化为在该条件下, 求u使得$u’Bu$式达到极大。考虑目标函数 $$ \varphi(u) = u’Bu - \lambda(u’Eu - 1) \tag{4.3-2} $$ 求导有
$$
  \left\{
  \begin{aligned}
    \frac{\partial \varphi}{\partial u} & = 2(B - \lambda E)u = 0  \qquad (a)\\
    \frac{\partial \varphi}{\partial u} & = u'Eu - 1 = 0 \qquad (b)
  \end{aligned} \right. \tag{4.3-3}
$$

对$(4.3-3)$的(a)式两边同乘$u’$, 有$u’Bu = \lambda u’Eu = \lambda$
从而, $u’Bu$的极大值为$\lambda$. 再用$E^{-1}$左乘$(4.3-3)$的b式, 有 $$(E^{-1}B - \lambda I)u = 0 \tag{4.3-4}$$ 由$(4.3-4)$式说明$\lambda$为$E^{-1}B$特征值, u为$E^{-1}B$的特征向量。在此最大特征值所对应的特征向量$u = (u_1, u_2, \cdots, u_p)’$即为所求。

示例 

library(MASS)

data(iris)
attach(iris)

# fisher判别
# 构建公式
#ld = lda(Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width)
ld = lda(iris[,1:4], iris[,5])
ld
# 预测
Z = predict(ld)
newG = Z$class

result = cbind(Species, newG, Z$x)
tab = table(newG, Species)
sum(diag(prop.table(tab)))
plot(ld)