距离判别法

距离判别法：距离判别方是通过计算待测点到各个分类的距离，在根据计算出距离的大小，进行判别该待测点属于那个分类。但是距离的计算是通过马氏距离进行计算的，而不是我们平常几何中用的欧式距离。

欧式距离的定义：

       欧氏距离是最易于理解的一种距离计算方法，源自欧氏空间中两点间的距离公式。

(1)二维平面上两点a(x1,y1)与b(x2,y2)间的欧氏距离：

 

(2)三维空间两点a(x1,y1,z1)与b(x2,y2,z2)间的欧氏距离：

 

(3)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的欧氏距离：

 

　　也可以用表示成向量运算的形式：

马氏距离的定义：

       有M个样本向量X1~Xm，协方差矩阵记为S，均值记为向量μ，则其中样本向量X到u的马氏距离表示为：

 

       而其中向量Xi与Xj之间的马氏距离定义为：

       若协方差矩阵是单位矩阵（各个样本向量之间独立同分布）,则公式就成了：

为什么使用马氏距离不使用欧式距离：

从上图两个分类服从正态分布u1和u2，  u1 均值为0 方差为1   、 u2 均值为4 方差为2。A点在x轴的值为1.66

按照欧式距离进行计算A 点跟接近u1分类，他离u1 的平均值跟接近

按照马氏距离进行计算，我们对两个分类进行密度函数计算。

正态分布概率密度计算公式为：

计算A属于u1分类概率密度=0.10

计算A 属于u2分类概率密度=0.33

上图也很好的解释为什么要考虑使用马氏距离，上图黄色椭圆部分代表整个分类在直角坐标系的排布情况，u代表分类的中心。图中有圆上有两个点分别是A和B ，使用欧式距离计算到u的距离是一样。但是从图中可以看出A点已经快出整个分类的范围，而B却在中心很密的位置。

距离判别法公式：

这里要讨论两种情况

情况1：分类来自一个总体，这样分类的协方差是一样的。

情况2：分类是来不同的总体，这样他们的协方差就不一样。

所以这里要考虑协方差相同的情况和不相同的情况

公式推理：

设总体x1 和x2,他们的均值分别是u1和u2,协方差分别是 Σ1和Σ2.

x 为待测样本

情况1，协方差相同

μ1 ≠ μ2, Σ1 = Σ2 = Σ

均值和协方差计算

化简得

情况2，协方差不相同

μ1  ≠μ2, Σ1  ≠ Σ2 

R 实现代码

#TrnX1类别1  TrnX2类别2 TstX测试集 var.equal方差是否一样

discriminiant.distance <- function(TrnX1, TrnX2, TstX = NULL, var.equal = FALSE){

  if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)

  if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))

  else if (is.matrix(TstX) != TRUE)

    TstX <- as.matrix(TstX)

  if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)

  if (is.matrix(TrnX2) != TRUE) 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 == TRUE || var.equal == T){

    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

}

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)

)

discriminiant.distance(classX1, classX2, var.equal=TRUE)

discriminiant.distance(classX1, classX2)

上述函数方法是解决两个分类的情况。如果出现多分类情况，可以通过分类的距离进行排序，取距离最近一个分类，即为当前待测变量的分类。

实现代码如下：

distinguish.distance <- function(TrnX, TrnG, TstX = NULL, var.equal = FALSE){

  if ( is.factor(TrnG) == FALSE){

    mx <- nrow(TrnX); mg <- nrow(TrnG)

    TrnX <- rbind(TrnX, TrnG)

    TrnG <- factor(rep(1:2, c(mx, mg)))

  }

  if (is.null(TstX) == TRUE) TstX <- TrnX

  if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))

  else if (is.matrix(TstX) != TRUE)

    TstX <- as.matrix(TstX)

  if (is.matrix(TrnX) != TRUE) 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 == TRUE || var.equal == T){

    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]; blong[j] <- i

      }

  }

  blong

}