1. 马氏距离定义

2. 源代码

输入mahanobis 即可得到代码：

> mahalanobis
function (x, center, cov, inverted = FALSE, ...) 
{
    x <- if (is.vector(x)) 
        matrix(x, ncol = length(x))
    else as.matrix(x)
    if (!isFALSE(center)) 
        x <- sweep(x, 2L, center)
    if (!inverted) 
        cov <- solve(cov, ...)
    setNames(rowSums(x %*% cov * x), rownames(x))
}
<bytecode: 0x0000000003e30778>
<environment: namespace:stats>

3. 案例操作

数据：7个变量（x1,x2,……x7）；12个样本
目标：计算每个样本到总体的马氏距离

数据录入：

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)
)

1. 直接带入函数计算马氏距离：

classX1 <- as.matrix(classX1) # 源代码中有把数据转化为矩阵的步骤，这里重新定义classX1,只是用来计算均值和方差.
n1 <- nrow(classX1)
S  <- (n1-1)*var(classX1)
mu1<-colMeans(classX1)
mahalanobis(classX1,mu1,S)
 [1] 0.3018766 0.1541574 0.2973262 0.3084365 0.7051001 0.7761428 0.7316278 0.6880709 0.6616244 0.9029164 0.7775262 0.6951945

2. 源代码 主要步骤学习：

function (x, center, cov, inverted = FALSE, ...)  
{
    # cov <- make.positive.definite (cov) 若矩阵不可逆
    x <- if (is.vector(x))    
        matrix(x, ncol = length(x))
    else as.matrix(x)  # 将X转换为矩阵形式
    
    if (!isFALSE(center))  # 如果样本不是均值(中心)，则用sweep函数计算（x-mu）
        x <- sweep(x, 2L, center) #2L:按照列的方向计算 12列
    if (!inverted)  # 默认值为 inverted = FALSE
        cov <- solve(cov, ...)   
    setNames(rowSums(x %*% cov * x), rownames(x))
}

【注释】 sweep () 函数 ：

 > classX1
       x1  x2  x3  x4 x5   x6  x7
 [1,] 6.6  39 1.0 6.0  6 0.12  20
 [2,] 6.6  39 1.0 6.0 12 0.12  20
 [3,] 6.1  47 1.0 6.0  6 0.08  12
 [4,] 6.1  47 1.0 6.0 12 0.08  12
 [5,] 8.4  32 2.0 7.5 19 0.35  75
 [6,] 7.2   6 1.0 7.0 28 0.30  30
 [7,] 8.4 113 3.5 6.0 18 0.15  75
 [8,] 7.5  52 1.0 6.0 12 0.16  40
 [9,] 7.5  52 3.5 7.5  6 0.16  40
[10,] 8.3 113 0.0 7.5 35 0.12 180
[11,] 7.8 172 1.0 3.5 14 0.21  45
[12,] 7.8 172 1.5 3.0 15 0.21  45
> mu1 = colMeans(classX1)
 x1         x2         x3         x4         x5         x6          
 7.3583333 73.6666667  1.4583333  6.0000000 15.2500000  0.1716667 
x7
49.5000000 
> x<-sweep(TrnX1,2L,mu1)  #用classX1的每一个值 减去mu1  
> x
              x1        x2          x3   x4    x5          x6    x7
 [1,] -0.7583333 -34.66667 -0.45833333  0.0 -9.25 -0.05166667 -29.5
 [2,] -0.7583333 -34.66667 -0.45833333  0.0 -3.25 -0.05166667 -29.5
 [3,] -1.2583333 -26.66667 -0.45833333  0.0 -9.25 -0.09166667 -37.5
 [4,] -1.2583333 -26.66667 -0.45833333  0.0 -3.25 -0.09166667 -37.5
 [5,]  1.0416667 -41.66667  0.54166667  1.5  3.75  0.17833333  25.5
 [6,] -0.1583333 -67.66667 -0.45833333  1.0 12.75  0.12833333 -19.5
 [7,]  1.0416667  39.33333  2.04166667  0.0  2.75 -0.02166667  25.5
 [8,]  0.1416667 -21.66667 -0.45833333  0.0 -3.25 -0.01166667  -9.5
 [9,]  0.1416667 -21.66667  2.04166667  1.5 -9.25 -0.01166667  -9.5
[10,]  0.9416667  39.33333 -1.45833333  1.5 19.75 -0.05166667 130.5
[11,]  0.4416667  98.33333 -0.45833333 -2.5 -1.25  0.03833333  -4.5
[12,]  0.4416667  98.33333  0.04166667 -3.0 -0.25  0.03833333  -4.5

【注释】 对于不可逆矩阵，使其正定：

方法1（稍稍改变矩阵的eigen value 使其正定）：

# Method by Higham 1988  
make.positive.definite = function(m, tol)
{
  if (!is.matrix(m))
    m = as.matrix(m)
  
  d = dim(m)[1]
  if (dim(m)[2] != d)
    stop("Input matrix is not square!") #判断是否是方阵
  
  es = eigen(m, symmetric = TRUE)
  esv = es$values
  
  if (missing(tol))
    tol = d * max(abs(esv)) * .Machine$double.eps
  delta =  2 * tol
  
  tau = pmax(0, delta - esv)
  dm = es$vectors %*% diag(tau, d) %*% t(es$vectors)
  
  return(m +  dm)

【注释】计算马氏距离

  setNames(rowSums(x %*% cov * x), rownames(x))

这里我们是计算，x中每一个样本到总体的距离，所以相当于要把每一个一位向量抽离出来，与方差矩阵相乘。因此后面一个乘号不是矩阵相乘，最后在用加法按照行的方向进行求和，得出12个（样本）距离总体的距离。

原理解析：

上面的公式计算的是一个样本x到总体X的距离，我们假设：

可以得到原方程（一个距离值） =

当我们把两个的样本放在同一个矩阵中时，我们假设：

我们会发现，如果完全按照矩阵相乘，我们最后会得到一个2*2的矩阵，然而我们的计算步骤应该是这样的：

由此我们会发现，

  setNames(rowSums(x %*% cov * x), rownames(x))

对代码的改动以及步骤是：（此处x=t(x)）

1. 不要再对第二个x进行转置
2. 不要对第二个乘号*进行举证运算
3. 按照行的方向，对列进行求和