R cor有时会返回NaN

9

我一直在处理一些数据,这些数据可以在这里找到:Dropbox 的 csv 文件(请使用它来复制错误时要友善)。

当我运行以下代码时:

t<-read.csv("120.csv")
x<-NULL
for (i in 1:100){
  x<-c(x,cor(t$nitrate,t$sulfate,use="na.or.complete"))
}
sum(is.nan(x))

我得到了最后一个表达式的随机值,通常在55到60左右。我期望cor能给出可重复的结果,因此我期望x是一个长度为100的向量,由相同的值组成。例如,请参见两次独立运行的输出:
> x<-NULL; for (i in 1:100){x<-c(x,cor(t$nitrate,t$sulfate,use="na.or.complete"))}
> sum(is.nan(x))
[1] 62
> head(x,10)
 [1]       NaN       NaN 0.2967441       NaN 0.2967441       NaN       NaN       NaN
 [9] 0.2967441       NaN
> x<-NULL; for (i in 1:100){x<-c(x,cor(t$nitrate,t$sulfate,use="na.or.complete"))}
> sum(is.nan(x))
[1] 52
> head(x,10)
 [1] 0.2967441       NaN       NaN       NaN       NaN 0.2967441 0.2967441       NaN
 [9] 0.2967441 0.2967441
> 

我不确定我是否做错了什么,或者这是一个未知的错误。如果是这样的话,我希望比我更熟练的人能帮助我向CRAN报告。
我读到了一个非常古老的(2001年)帖子,其中cor.test表现出相同的行为(请参见cor.test produces NaN sometimes)。
由于我是R的新手,所以非常感谢您的解释。谢谢!
根据Ben的建议:
> sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252    LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
[5] LC_TIME=Spanish_Colombia.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] stringr_0.6.2     digest_0.6.4      RCurl_1.95-4.3    bitops_1.0-6      qpcR_1.4-0        Matrix_1.1-4      robustbase_0.91-1 rgl_0.95.1157    
 [9] minpack.lm_1.1-8  MASS_7.3-35       plyr_1.8.1        swirl_2.2.16      ggplot2_1.0.0     lattice_0.20-29  

loaded via a namespace (and not attached):
 [1] colorspace_1.2-4 DEoptimR_1.0-2   grid_3.1.1       gtable_0.1.2     httr_0.5         labeling_0.3     munsell_0.4.2    proto_0.3-10     Rcpp_0.11.3     
[10] reshape2_1.4     scales_0.2.4     testthat_0.9.1   tools_3.1.1      yaml_2.1.13  

find("cor")的结果:

> find("cor")
[1] "package:stats"

---------- ### 第二次编辑 ###--------

我重新启动了会话(我没有找到如何传递--vanilla参数。我正在使用Rstudio),这是新的sessionInfo:

> sessionInfo()
R version 3.1.1 (2014-07-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252    LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
[5] LC_TIME=Spanish_Colombia.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] tools_3.1.1

我在新会话中再次运行命令,但仍然得到sum(is.nan(x))=52 :(

以防有用:

> cor
function (x, y = NULL, use = "everything", method = c("pearson", 
    "kendall", "spearman")) 
{
    na.method <- pmatch(use, c("all.obs", "complete.obs", "pairwise.complete.obs", 
        "everything", "na.or.complete"))
    if (is.na(na.method)) 
        stop("invalid 'use' argument")
    method <- match.arg(method)
    if (is.data.frame(y)) 
        y <- as.matrix(y)
    if (is.data.frame(x)) 
        x <- as.matrix(x)
    if (!is.matrix(x) && is.null(y)) 
        stop("supply both 'x' and 'y' or a matrix-like 'x'")
    if (!(is.numeric(x) || is.logical(x))) 
        stop("'x' must be numeric")
    stopifnot(is.atomic(x))
    if (!is.null(y)) {
        if (!(is.numeric(y) || is.logical(y))) 
            stop("'y' must be numeric")
        stopifnot(is.atomic(y))
    }
    Rank <- function(u) {
        if (length(u) == 0L) 
            u
        else if (is.matrix(u)) {
            if (nrow(u) > 1L) 
                apply(u, 2L, rank, na.last = "keep")
            else row(u)
        }
        else rank(u, na.last = "keep")
    }
    if (method == "pearson") 
        .Call(C_cor, x, y, na.method, FALSE)
    else if (na.method %in% c(2L, 5L)) {
        if (is.null(y)) {
            .Call(C_cor, Rank(na.omit(x)), NULL, na.method, method == 
                "kendall")
        }
        else {
            nas <- attr(na.omit(cbind(x, y)), "na.action")
            dropNA <- function(x, nas) {
                if (length(nas)) {
                  if (is.matrix(x)) 
                    x[-nas, , drop = FALSE]
                  else x[-nas]
                }
                else x
            }
            .Call(C_cor, Rank(dropNA(x, nas)), Rank(dropNA(y, 
                nas)), na.method, method == "kendall")
        }
    }
    else if (na.method != 3L) {
        x <- Rank(x)
        if (!is.null(y)) 
            y <- Rank(y)
        .Call(C_cor, x, y, na.method, method == "kendall")
    }
    else {
        if (is.null(y)) {
            ncy <- ncx <- ncol(x)
            if (ncx == 0) 
                stop("'x' is empty")
            r <- matrix(0, nrow = ncx, ncol = ncy)
            for (i in seq_len(ncx)) {
                for (j in seq_len(i)) {
                  x2 <- x[, i]
                  y2 <- x[, j]
                  ok <- complete.cases(x2, y2)
                  x2 <- rank(x2[ok])
                  y2 <- rank(y2[ok])
                  r[i, j] <- if (any(ok)) 
                    .Call(C_cor, x2, y2, 1L, method == "kendall")
                  else NA
                }
            }
            r <- r + t(r) - diag(diag(r))
            rownames(r) <- colnames(x)
            colnames(r) <- colnames(x)
            r
        }
        else {
            if (length(x) == 0L || length(y) == 0L) 
                stop("both 'x' and 'y' must be non-empty")
            matrix_result <- is.matrix(x) || is.matrix(y)
            if (!is.matrix(x)) 
                x <- matrix(x, ncol = 1L)
            if (!is.matrix(y)) 
                y <- matrix(y, ncol = 1L)
            ncx <- ncol(x)
            ncy <- ncol(y)
            r <- matrix(0, nrow = ncx, ncol = ncy)
            for (i in seq_len(ncx)) {
                for (j in seq_len(ncy)) {
                  x2 <- x[, i]
                  y2 <- y[, j]
                  ok <- complete.cases(x2, y2)
                  x2 <- rank(x2[ok])
                  y2 <- rank(y2[ok])
                  r[i, j] <- if (any(ok)) 
                    .Call(C_cor, x2, y2, 1L, method == "kendall")
                  else NA
                }
            }
            rownames(r) <- colnames(x)
            colnames(r) <- colnames(y)
            if (matrix_result) 
                r
            else drop(r)
        }
    }
}
<bytecode: 0x0000000008ce0158>
<environment: namespace:stats>

再次感谢您。

2
如果有人能添加“cor”标签,我将不胜感激。我的声望(仍然少于1500)不允许我添加新标签,而我认为这个标签对于遇到同样问题的人非常关键。谢谢! - PavoDive
1
我无法复制;我总是得到sum(is.nan(x))等于零。(1) 尝试在干净的R会话中开始(如果可能,请使用--vanilla); (2)sessionInfo()的结果?(3)find("cor")的结果? - Ben Bolker
我怀疑这是你方面的问题,而不是一般性的错误,但这是一个非常好的问题。 - Ben Bolker
越来越好奇了。我唯一剩下的建议是重复上面的第1点(确保您在干净的R会话中开始) - sessionInfo()输出中出现非推荐包的外观表明您还没有这样做。 - Ben Bolker
1
希望有人能在运行于64位Windows上的3.1.1上测试这个奇怪的问题。虽然出现bug的可能性非常非常小,但我已经没有其他解释了。 - Ben Bolker
显示剩余3条评论
3个回答

6

以下是几点注释:

  • 没有人能够重现您的问题
  • 这不可能是120.csv文件出了问题,它是完好无损的。
  • 实际上,使用另一个 use=".." 选项只是一种权宜之计
  • R源代码中的底层C代码在检测值是否为NA或NaN时使用ISNAN(.),而这个函数将转到(系统内部的)C库的isnan(.)函数。
  • 你(而且只有你)有时会得到NaN,因为在某些情况下ISNAN(.)并没有返回“true”,而浮点数算法却在与NA一起计算,并正确地返回NaN。

作为一个“老”R核心成员,我可以向您保证,在R的核心计算中,ISNAN(.)在很多基本位置都被使用,而对于您来说,它有时似乎不能检测到NA/NaN,以至于它们传播到结果中非常棘手。 正如Duncan Murdoch在回答您的R错误报告https://bugs.r-project.org/bugzilla/show_bug.cgi?id=16058时所说,这必须是您特定的“系统”出了问题。 我假设您只是从CRAN下载了R,即使对于R 3.1.2,您仍然遇到问题,我倾向于认为您的系统软件(Windows)或 - 不太可能 - 您的硬件必须轻微地损坏/损坏。


1
感谢您的回复。我原本计划在错误报告中链接到讨论,但您已经提前做了!Duncan建议可能存在msvcrt.dll库的问题或浮点处理器出现故障。我读到了一些类似的Octave构建投诉,这让我再次考虑迁移到Linux或OSX。如果这确实是一个库的问题,我相信有些人并不知道他们的系统中存在这个问题。我希望有人能解决这个问题。同时,如您所建议的,我会将自己的回复标记为解决方法。 - PavoDive

1

Giovanni,对我来说它运行得很好。也许你应该尝试更改cor的参数为=“complete.obs”,看看是否有所帮助。此外,你应该检查你的CSV文件是否损坏。

希望这能有所帮助。


1
发布您的R版本和平台对于这种类型的答案会很有帮助。 - Ben Bolker
@Danish:感谢您的建议。它仍然产生不一致的结果。另一方面,如果csv文件已损坏,则会产生一致的意外(或错误)结果。真正让我困惑的是不一致性。 - PavoDive

0
虽然我感到困惑,但我开始尝试使用cor函数的不同选项。我发现如果我使用use="pairwise.complete.obs",就可以得到一致的结果:cor(t$nitrate,t$sulfate,use="pairwise.complete.obs")
> x<-NULL; for (i in 1:100){x<-c(x,cor(t$nitrate,t$sulfate,use="pairwise.complete.obs"))};x
  [1] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [12] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [23] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [34] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [45] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [56] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [67] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [78] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
 [89] 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441 0.2967441
[100] 0.2967441

我仍然不明白其他用户传递的use选项为什么没有引起奇怪的行为。


正如@MartinMächler所说,这更像是一个解决方法而不是一个解决方案,但可能会对某些人有所帮助,因此我将其标记为答案。 - PavoDive

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接