1

我rdist用来计算矩阵与其自身之间的欧几里得距离：

> m = matrix(c(1,1,1,2,2,2,3,4,3),nrow=3, ncol=3)
> m
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    4
[3,]    1    2    3
library(fields)
> rdist(m)
      [,1]  [,2]  [,3]
[1,] 1e-10 1e+00 1e-10
[2,] 1e+00 1e-10 1e+00
[3,] 1e-10 1e+00 1e-10

让我感到困惑的是，我认为它应该在对角线上有 0（肯定一个向量到自身的距离是 0？），出于同样的原因，它应该在比较第一行和第三行的地方有 0。相反，我看到的值（1e-10）看起来很大，是数字噪声。怎么了？

编辑：rdist来自包fields。

4

1 回答 1

4

首先1e-10是简单1*10^-10的0.0000000001，所以数值上非常接近0（因为它是平方根的结果，所以计算中的实际误差是数量级的1e-20）。是不是“太大”了？好吧，库是用 fortran 编写的，并且专注于速度，所以它是完全可以接受的。如果你分析确切的代码，你会发现它是如何计算的：

# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
"rdist" <- function(x1, x2) {
    if (!is.matrix(x1)) 
        x1 <- as.matrix(x1)
    if (missing(x2)) 
        x2 <- x1
    if (!is.matrix(x2)) 
        x2 <- as.matrix(x2)
    d <- ncol(x1)
    n1 <- nrow(x1)
    n2 <- nrow(x2)
    par <- c(1/2, 0)
    temp <- .Fortran("radbas", nd = as.integer(d), x1 = as.double(x1), 
        n1 = as.integer(n1), x2 = as.double(x2), n2 = as.integer(n2), 
        par = as.double(par), k = as.double(rep(0, n1 * n2)))$k
    return(matrix(temp, ncol = n2, nrow = n1))
}

确切的答案隐藏在 fortran 文件中（radfun.f调用 from radbas.f），您可以在其中找到该行

if( dtemp.lt.1e-20) dtemp =1e-20

它将小（甚至 0）值视为1e-20，在取平方根后得到1e-10。似乎动机是通过使用值的对数来加速该过程（因此，平方根只是简单地除以 2），这当然没有为0.

于 2013-09-04T18:05:27.403 回答