判别分析（Discriminant Analysis, DA） - ricket-sjtu/bi028 GitHub Wiki

1. 判别分析

判别分析是判别个体所属群体的一种统计学方法，常用的判别分析方法主要包括：

基于距离的判别分析方法
基于贝叶斯的判别分析方法
Fisher判别分析方法（Fisher's discriminant analysis, FDA）

1.1 贝叶斯判别分析

先验概率：对分析对象的先验知识：$\pi_i$
通过样本分析得到后验概率分布：$p(\pi_i | x) = \frac{\pi_i f(x|\pi_i)}{\sum_{j=1}^g \pi_j f(x | \pi_j)}$
通过后验概率分布进行统计推断：$\operatorname{argmax}_i p(\pi_i | x)$ (后验概率最大化) 实际上就是使得平均误判损失（误判概率与误判损失的加权和）ECM达到最小的过程。

1.2 判别分析的过程

收集数据
设置先验概率（prior）:

等配先验（equal prior）：$\hat{\pi}_i = 1/g$；
任意先验（arbitrary prior）：$\sum_i \hat{\pi}_i = 1$
经验先验（estimated prior）：$\hat{\pi}_i = n_i / N$

使用Bartlett检验确定方差-协方差矩阵（variance-covariance matrix）是否同质（homogeneous）：

如果同质，则采用线性判别分析（linear discriminant analysis, LDA）；
如果异质，则采用二次判别分析（quadratic discriminant analysis, QDA）；

确定条件概率密度函数的参数，必须满足以下假设：

第$i$组数据具有相同的平均向量$\mu_i$；
第$i$组数据具有相同的方差-协方差矩阵$\Sigma$；
个体的独立性（independence）；
数据呈多元正态分布（multivariate normal distribution）

计算判别函数（discriminant function）
使用**交叉验证（cross validation）**估计误分类概率（misclassification probabilities）
为未知分组样本确定其分组

1.2.1 线性判别分析（Linear discriminant analysis, LDA）

假定在总体$i$中，$X \sim N(\mu_i, \Sigma_p)$多元正态分布： $$ f(x | \pi_i) = \frac{1}{(2\pi)^{p/2} |\Sigma|^{1/2}} \exp \left[ (x-\mu_i)^T \Sigma^{-1} (x-\mu_i)\right] $$ 因此分组就是最大化$\pi_i (f(x | \pi_i))$，也可以转化为$\log [\pi_i (f(x | \pi_i))] $。

线性打分函数（Linear Score Function）

$$ s^L_i(\mathbf{X}) = -\frac{1}{2}\mathbf{\mu_i^T \Sigma^{-1}\mu_i + \mu_i^T \Sigma^{-1} x}+ \log \pi_i = d_{i0}+\sum_{j=1}^{p}d_{ij}x_j + \log \pi_i $$ 其中

$d_{i0} = -\frac{1}{2}\mathbf{\mu_i^T \Sigma^{-1}\mu_i}$；
$d_{ij}$是$\mu_i^T \Sigma^{-1}$的第$j$各元素

LDA步骤

设置先验概率：$p_i = Pr(\pi_i); i=1, 2, \ldots, g$；
计算每个组的均值：$\mu_i = E(X | \pi_i) = \bar{x}_i, i=1,2, \ldots, g$；
计算方差-协方差矩阵：$\Sigma = var(X | \pi_i) = \mathbf{S}p = \frac{\sum{i=1}^{g}(n_i-1)\mathbf{S}i}{\sum{i=1}^{g}(n_i-1)}$；
估计线性打分函数：$\hat{s}^L_i(\mathbf{x}) = -\frac{1}{2}\mathbf{\bar{x}'_i S^{-1}_p \bar{x}_i +\bar{x}'i S^{-1}p x } + \log{\hat{p}i} = \hat{d}{i0} + \sum{j=1}^{p}\hat{d}{ij}x_j + \log{p}_i$，其中

$\hat{d}_{i0} = -\frac{1}{2}\mathbf{\bar{x}'_i S^{-1}_p \bar{x}_i}$；
$\hat{d}_{ij}$是$\mathbf{\bar{x}'_iS^{-1}_p}$的第$j$个元素；

判别$\mathbf{x}$为令$\hat{s}^L_i(\mathbf{x})$最大的$i$。

1.2.2 二次判别分析（Quadratic discriminant analysis, QDA）

QDA与LDA相似，也要定义一个二次打分函数（Quadratic score function）: $$ s_i^Q (\mathbf{x}) = -\frac{1}{2}\log{|\mathbf{\Sigma_i}|}-\frac{1}{2}{\mathbf{(x-\mu_i)^T \Sigma_i^{-1}(x - \mu_i)}}+\log{p_i} $$ 不同的是，这里$\Sigma_i$的采用估计值$\mathbf{S}_i$替代。

1.2.3 两个总体的贝叶斯判别分析方法

discriminant.bayes <- function(TrnX1, TrnX2, rate=1, TstX=NULL, var.equal=FALSE) {
  #@TrnX1, TrnX2: data.frame, training set for class 1 and class 2, respectively
  #@rate = p1/p2
  #@TstX: data.frame, test set
  #@var.equal: TRUE/FALSE, the two class has equal variance
  if (is.null(TstX)) TstX <- rbind(TrnX1, TrnX2)
  if (is.vector(TstX)) TstX <- t(as.matrix(TstX))
  if (!is.matrix(TstX)) TstX <- as.matrix(TstX)
  if (!is.matrix(TrnX1)) TrnX1 <- as.matrix(TrnX1)
  if (!is.matrix(TrnX2)) TrnX2 <- as.matrix(TrnX2)
  nx <- nrow(TstX)
  blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE, 
                  dimnames=list("blong", 1:nx))
  mu1 <- colMeans(TrnX1)
  mu2 <- colMeans(TrnX2)
  if (var.equal) {
    S <- var(rbind(TrnX1, TrnX2))
    beta <- 2*log(rate)
    w <- mahalanobis(TstX, mu2, S) - mahalanobis(TstX, mu1, S)
  } else {
    S1 <- var(TrnX1)
    S2 <- var(TrnX2)
    beta <- 2*log(rate) + log(det(S1)/det(S2))
    w <- mahalanobis(TstX, mu2, S2) - mahalanobis(TstX, mu1, S1)
  }
  blong <- (w <= beta) + 1
  blong
}

TrnX1<-matrix(c(24.8, 24.1, 26.6, 23.5, 25.5, 27.4, 
                -2.0, -2.4, -3.0, -1.9, -2.1, -3.1),
              ncol=2)
TrnX2<-matrix(c(22.1, 21.6, 22.0, 22.8, 22.7, 21.5, 22.1, 21.4,
               -0.7, -1.4, -0.8, -1.6, -1.5, -1.0, -1.2, -1.3),
             ncol=2)
p <- nrow(TrnX1)
q <- nrow(TrnX2)
discriminant.bayes(TrnX1, TrnX2, rate=q/p,var.equal=TRUE)
discriminant.bayes(TrnX1, TrnX2, rate=q/p)

1.2.3 多总体贝叶斯判别分析

multiclass.discriminant.bayes <- function(TrnX, TrnG, p=rep(1, length(levels(TrnG))), TstX=NULL, var.equal=FALSE) {
  if (!is.factor(TrnG)) {
    mx <- nrow(TrnX); mg <- nrow(TrnG)
    TrnX <- rbind(TrnX, TrnG)
    TrnG <- factor(rep(1:2, c(mx, mg)))
  }
  if (is.null(TstX)) TstX <- TrnX
  if (is.vector(TstX)) TstX <- as.matrix(TstX)
  else if (!is.matrix(TstX)) TstX <- as.matrix(TstX)
  if (!is.matrix(TrnX)) TrnX <- as.matrix(TrnX)

  nx <- nrow(TstX)
  blong <- matrix(rep(0, nx), nrow=1, dimnames=list("blong", 1:nx))

  g <- length(levels(TrnG))
  mu <- matrix(0, nrow=g, ncol=ncol(TrnX))
  for (i in 1:g) {
    mu[i, ] <- colMeans(TrnX[TrnG==i, ])
  }
  D <- matrix(0, nrow=g, ncol=nx)
  if (var.equal) {
    for (i in 1:g) {
      d2 <- mahalanobis(TstX, mu[i, ], var(TrnX))
      D[i,] <- d2 - 2*log(p[i])
    }
  } else {
    for (i in 1:g) {
      S <- var(TrnX[TrnG==i,])
      d2 <- mahalanobis(TstX, mu[i, ], S)
      D[i,] <- d2 - 2*log(p[i]) - log(det(S)) 
    }
  }
  for (j in 1:nx) blong[j] <- which.min(D[,j])

  blong
}

data(iris)
X<-iris[,1:4]
G<-gl(3,50)
multiclass.discriminant.bayes(X, G, var.equal=T)
multiclass.discriminant.bayes(X, G)