biobabble作者群有一个小伙伴提出来的问题。

没错，凡是给biobabble投稿的小伙伴都会被拉入作者群。

话说有一个矩阵，我们让它某个值缺失，算出来的距离和原来的竟然是一样的：

> dat = matrix(1:10, 2)
> dist(dat)
         1
2 2.236068
> dat[1,1] = NA
> dist(dat)
         1
2 2.236068

这有点反直觉，实际上你再搞个缺失值，它算出来还是一样的：

> dat[1,3] = NA
> dist(dat)
         1
2 2.236068

但以我的经验，我通常不敢轻易说人家有bug，我们得确认一下。

于是我读了R的源代码中欧式距离的代码：

static double R_euclidean(double *x, int nr, int nc, int i1, int i2)
{
    double dev, dist;
    int count, j;

    count= 0;
    dist = 0;
    for(j = 0 ; j < nc ; j++) {
	if(both_non_NA(x[i1], x[i2])) {
	    dev = (x[i1] - x[i2]);
	    if(!ISNAN(dev)) {
		dist += dev * dev;
		count++;
	    }
	}
	i1 += nr;
	i2 += nr;
    }
    if(count == 0) return NA_REAL;
    if(count != nc) dist /= ((double)count/nc);
    return sqrt(dist);
}

把NA干掉？

首先从代码上看，缺失值会被去掉，这就是count计数看有多少有效值的作用。如果我们简单粗爆地把NA值去掉，我们会算出来什么值呢？

  1. 完整数据的情况下，距离是2.236068
  2. dat[1,1] = NA的情况下，距离是2
  3. dat[1,c(1,3)] = NA的情况下，距离是1.732051

从这里我们可以看到，缺失值越多，距离是越小的，这个从欧式距离的公式上也反映出来，所有的量都是正值，少了就小了。所以如果我们简单地去NA值，那么对距离的估算会低估。

怎么补救？

这就是倒数第二句干的事情：

    if(count != nc) dist /= ((double)count/nc);

如果有缺少，用已知的点的均值来估计。也就是用mean((qi-pi)^2)去代替缺失点。

不难验证一下，前面第二种情况时距离是2，那么被一下缺失值：

> sqrt(2^2 / (4/5))
[1] 2.236068

出来的就是第一种情况时的值了。

如果你看不明白的话，我再拆一下给你看：

  • 2^2 是总和
  • 2^2/4 是均值
  • 2^2 + 2^2/4 就是补了缺失值的总和
  • 2^2 / (4/5)等同于上面的值
  • sqrt(2^2 / (4/5))算出来就是距离，请对照公式来看一下。

为什么能够完美估计原来的值？

这也是一开头让大家疑惑的地方，为什么做为估计出来的值能和真实的值一样？

这其实是个巧合，真的是巧合。你如果是一般的随机数拿来试试，肯定无法一模一样的值。

当然对于一个determinant的计算，巧合也是算出来的，必定是巧合，原因必须存在于数值之中。

如果你对dat算相关系数的话，就明白了。我们对原始值算一下：

> cor(dat[1,], dat[2,])
[1] 1

能够完美估算出距离，其实等同于能够完美复原缺失点，而这种情况只在于数据的相关系数等于一的情况下，数据的相关性越好，也就能够复原得越好，这很好理解。

欧式距离对于outlier是很敏感的，你可以试试不同的值去替换掉原始值，再看看缺失值的情况下，计算出来的距离和不缺失的情况下，差距有多少。

写到这里，我想到了欧式距离可以拿来检测outlier啊。比如我们每次让一个点缺失，比如说我们有100个点，就可以计算出100个距离，这些距离各不相等（正常情况下数据相关系统不会刚刚好是1），它们和真实值的差值也可以算出来，这100个差值里面，差值较大的，相对应的缺失点，就是outlier。