我一直在处理一些数据,这些数据可以在这里找到: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>
再次感谢您。
sum(is.nan(x))
等于零。(1) 尝试在干净的R会话中开始(如果可能,请使用--vanilla
); (2)sessionInfo()
的结果?(3)find("cor")
的结果? - Ben BolkersessionInfo()
输出中出现非推荐包的外观表明您还没有这样做。 - Ben Bolker