一、概述

距离判别法是最简单、最直观的一种判别方法，该方法适用于连续型随机变量的判别，对变量的概率分布没有限制。
原理：计算待测点与各类的距离，取最短者为其所属分类。值得注意的是，距离的衡量有很多种方式，这里采用的是马氏距离。

二、马氏距离

1.欧式距离与马氏距离

通常，我们所定义的距离是欧式距离。若x,y是n维空间中的两个点，则x与y的距离为：

python马氏距离判别分析 spss马氏距离判别_协方差


但在统计分析与计算中，欧式距离就不适用了。从以下例子可以看出：

python马氏距离判别分析 spss马氏距离判别_方差_02


图中A点，若用欧式距离判断，离μ1更近；但是从概率角度来判断，情况并非如此。经计算，A点与μ1的距离为1.66σ1，μ2的距离为1.17σ2，离μ2更近，统计学上的意义则为A点有更大的概率属于右边的分布。

直观上，可以理解为欧式距离是绝对距离，而马氏距离是考虑了随机变量的方差的一种相对距离，体现了概率的思想。

2.马氏距离的定义

设x,y是服从均值为μ，协方差阵为Σ的总体X中抽取的样本，则x,y两点的马氏距离为：

python马氏距离判别分析 spss马氏距离判别_二分类_03


样本点x与总体X的马氏距离为：

python马氏距离判别分析 spss马氏距离判别_python马氏距离判别分析_04


公式中的x，y，μ都是向量。

三、判别准则与判别函数

在这里，讨论两个总体的距离判断，分别讨论两个总体协方差阵相同和协方差阵不同的情况。

设总体X1和X2的均值向量分别为μ1和μ2，协方差真分别为Σ1和Σ2，给定一个样本点x，要判断x来自于哪一个总体。

其思路为，分别计算样本点离两个样本的中心点的距离，然后比较两个距离的大小，从而判断其分类。具体来看：

（1）两个总体协方差阵相同

判别函数为：

python马氏距离判别分析 spss马氏距离判别_算法_05


判别准则为：

python马氏距离判别分析 spss马氏距离判别_算法_06


（2）两个总体协方差阵不同

判别函数为：

python马氏距离判别分析 spss马氏距离判别_算法_07


判别准则为：

python马氏距离判别分析 spss马氏距离判别_算法_06

（与上述相同）

注：这里以R语言实现为重点，故具体的数学推导过程不作详解。

四、R程序

由于R中没有专门的函数用于距离判别法的计算（可能存在，只是我不知道~），故需要编写自定义函数，将其写成R程序。
具体如下：

discriminiant.distance <-function
  (TrnX1,TrnX2,TstX=NULL,var.equal=F)
{
  #双分类距离判别方法，TrnX1、TrnX2代表训练集 TstX代表测试集。var.equal如为T表示训练样本的协方差
  #相同,默认值为不相同
  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=T,dimnames=list("blong",1:nx))
  mu1 <- colMeans(TrnX1)
  mu2 <- colMeans(TrnX2)
  if (var.equal == TRUE)
  {
    S <- var(rbind(TrnX1, TrnX2)) 
    w <- mahalanobis(TstX,mu2,S) - mahalanobis(TstX,mu1,S)
  }
  else
  {
    S1 <- var(TrnX1); S2 <- var(TrnX2)
    w <- mahalanobis(TstX,mu2, S2) - mahalanobis(TstX,mu1,S1)
  }
  for (i in 1:nx)
  {
    if (w[i] > 0) blong[i] <- 1 else blong[i] <-2
  }  
  blong
}

五、实例

在研究砂基液化问题中，选了7个因子，抽取35个样本，其中12个已液化，23个未液化。

#读取数据
> classX1<-data.frame(
+   x1=c(6.60, 6.60, 6.10, 6.10, 8.40, 7.2, 8.40, 7.50,
+        7.50, 8.30, 7.80, 7.80),
+   x2=c(39.00,39.00, 47.00, 47.00, 32.00, 6.0, 113.00, 52.00,
+        52.00,113.00,172.00,172.00),
+   x3=c(1.00, 1.00, 1.00, 1.00, 2.00, 1.0, 3.50, 1.00,
+        3.50, 0.00, 1.00, 1.50),
+   x4=c(6.00, 6.00, 6.00, 6.00, 7.50, 7.0, 6.00, 6.00,
+        7.50, 7.50, 3.50, 3.00),
+   x5=c(6.00, 12.00, 6.00, 12.00, 19.00, 28.0, 18.00, 12.00,
+        6.00, 35.00, 14.00, 15.00),
+   x6=c(0.12, 0.12, 0.08, 0.08, 0.35, 0.3, 0.15, 0.16,
+        0.16, 0.12, 0.21, 0.21),
+   x7=c(20.00,20.00, 12.00, 12.00, 75.00, 30.0, 75.00, 40.00,
+        40.00,180.00, 45.00, 45.00))
> classX2<-data.frame(
+     x1=c(8.40, 8.40, 8.40, 6.3, 7.00, 7.00, 7.00, 8.30,
+          8.30, 7.2, 7.2, 7.2, 5.50, 8.40, 8.40, 7.50,
+          7.50, 8.30, 8.30, 8.30, 8.30, 7.80, 7.80),
+     x2=c(32.0 ,32.00, 32.00, 11.0, 8.00, 8.00, 8.00,161.00,
+          161.0, 6.0, 6.0, 6.0, 6.00,113.00,113.00, 52.00,
+          52.00, 97.00, 97.00,89.00,56.00,172.00,283.00),
+     x3=c(1.00, 2.00, 2.50, 4.5, 4.50, 6.00, 1.50, 1.50,
+          0.50, 3.5, 1.0, 1.0, 2.50, 3.50, 3.50, 1.00,
+          1.00, 0.00, 2.50, 0.00, 1.50, 1.00, 1.00),
+     x4=c(5.00, 9.00, 4.00, 7.5, 4.50, 7.50, 6.00, 4.00,
+          2.50, 4.0, 3.0, 6.0, 3.00, 4.50, 4.50, 6.00,
+          7.50, 6.00, 6.00, 6.00, 6.00, 3.50, 4.50),
+     x5=c(4.00, 10.00, 10.00, 3.0, 9.00, 4.00, 1.00, 4.00,
+          1.00, 12.0, 3.0, 5.0, 7.00, 6.00, 8.00, 6.00,
+          8.00, 5.00, 5.00,10.00,13.00, 6.00, 6.00),
+     x6=c(0.35, 0.35, 0.35, 0.2, 0.25, 0.25, 0.25, 0.08,
+          0.08, 0.30, 0.3, 0.3, 0.18, 0.15, 0.15, 0.16,
+          0.16, 0.15, 0.15, 0.16, 0.25, 0.21, 0.18),
+     x7=c(75.00,75.00, 75.00, 15.0,30.00, 30.00, 30.00, 70.00,
+          70.00, 30.0, 30.0, 30.0,18.00, 75.00, 75.00, 40.00,
+          40.00,180.00,180.00,180.00,180.00,45.00,45.00))
> source("D:\\program files\\RStudio\\sources\\discriminiant.distance.R")
> discriminiant.distance(classX1, classX2, var.equal=TRUE)
      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
blong 1 1 1 1 1 1 1 1 2  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  1  1  2  2  2  2  2  2
> discriminiant.distance(classX1, classX2)
      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
blong 1 1 1 1 1 1 1 1 2  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2

在认为两总体协方差矩阵相同的情况下，将训练样本回代进行判别，有三个点判错；在认为两总体协方差矩阵不相同的情况下，将训练样本回代进行判别，有一个点判错。

六、多分类问题的距离判别

1.概述

由于距离判别的本质就是计算待分类点与各个总体的马氏距离，因此很容易由二分类推广到多分类问题。

假设样本共有k类，分别为X1，X2……Xk。

若这k类总体方差相等，则可以用全部样本的方差作为总体方差的估计值Σ；若这k类总体方差不相等，则用各自的样本计算样本方差作为各自的总体方差Σj的估计值。

相应的判别准则为：

python马氏距离判别分析 spss马氏距离判别_方差_09

2.R程序

用同样的办法编写R程序，具体如下：

distinguish.distance1 <-function
(TrnX, TrnG, TstX=NULL, var.equal= F)  #将所有的训练样本都保存在TranX中，将训练样本的分类保存在TrnG中
{
  #以下这段程序确保如果按二分类问题中分别按种类一和种类二输入训练集，该R程序也能通过转换后兼容
  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 <- t(as.matrix(TstX))

  if (!is.matrix(TstX)) TstX <- as.matrix(TstX)
  if (!is.matrix(TrnX)) TrnX <- as.matrix(TrnX)
  #设置保存结果集的矩阵
  nx <- nrow(TstX)
  belong <- matrix(rep(0,nx),nrow=1,dimnames=list("belong",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)
      D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX))
  }
  else
  {
    for (i in 1:g)
      D[i,] <- mahalanobis(TstX, mu[i,],var(TrnX[TrnG ==i,]))
  }
  #将待测点划分到与之距离最小的类中
  for ( j in 1:nx)
  {
    dmin <- Inf
    for (i in 1:g)
      if (D[i,j] < dmin)
      {
        dmin <- D[i,j]; belong[j] <- i
      }
  }
  #返回结果
  belong
}

3.实例

这里采用自带的Iris数据集进行实例演示。

> X <- iris[,1:4]
> G <- gl(3,50)
> source("D:\\program files\\RStudio\\sources\\discriminiant.distance1.R")
> distinguish.distance1(X,G)
       1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
belong 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
       38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
belong  1  1  1  1  1  1  1  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  3
       72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
belong  2  3  2  2  2  2  2  2  2  2  2  2  3  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2   2   3   3   3   3
       105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
belong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
       131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
belong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3

计算结果显示，有三个样本判错，判别正确率为147/150 = 98%。

七、过程中遇到的问题

这里特别强调一点，就是在用source函数时，如果不指定路径，会特别容易出错，如下：

> source("discriminiant.distance.R")
Error in file(filename, "r", encoding = encoding) : 
  cannot open the connection
In addition: Warning message:
In file(filename, "r", encoding = encoding) :
  cannot open file 'discriminiant.distance.R': No such file or directory

在网上找了很久发现很多人都遇到了此问题，但是没有解决办法，后来意识到是路径问题，也就是在用function创建自定义函数时，不知道系统把.R函数放在哪个路径中，我甚至尝试在电脑上搜索也没能找到，不知是何原因。后来我的解决方法是先用cat()函数创建一个空的.R程序，该程序会存在于默认文件夹里（可以用getwd()查询），然后将编好的函数粘贴在.R中，再次调用就成功了。具体如下：

> cat("",file = "discriminiant.distance.R")
#期间将编好的函数粘贴在.R中
> source("discriminiant.distance.R") #调用成功