R cor有时会返回NaN

我一直在研究一些数据,这里有: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
> 

我想知道我在这里做错了什么,或者它是否是一个[n] [un]已知错误.如果是这样的话,我很感激,如果有人比我更有帮助我帮助我向CRAN报告.

我读了一篇非常古老的(2001)文章,其中cor.test表现出同样的行为(参见cor.test有时会产生NaN).

我很感激你的解释,因为我是R的nOOb谢谢!

Per 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  

查找结果("cor"):

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

---------- ### Second Edit ### --------

我重新启动了会话(我没有找到如何传递--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

我在新的会话中再次运行命令,仍然得到总和(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)
        }
    }
}

再次感谢.