logo标签列表

在矩阵中带有重置的计数列数据

r
17
我正在收集有关我的猫便便数量的数据，并将其制成矩阵：
m <- cbind(fluffy=c(1.1,1.2,1.3,1.4),misterCuddles=c(0.9,NA,1.1,1.0))
row.names(m) <- c("2013-01-01", "2013-01-02", "2013-01-03","2013-01-04")

这给了我这个：

           fluffy misterCuddles
2013-01-01    1.1           0.9
2013-01-02    1.2            NA
2013-01-03    1.3           1.1
2013-01-04    1.4           1.0

每天，我想知道每只猫已经连续便便了多少天。因此，生成的矩阵应该如下所示：

           fluffy misterCuddles
2013-01-01      1             1
2013-01-02      2             0
2013-01-03      3             1
2013-01-04      4             2
有没有一种高效的方法可以做到这一点？cumsum 函数执行的是类似的操作，但它是一个基本函数，我无法修改它来适应我的需求。
我可以使用 for 循环并存储计数，如下所示：
m.output <- matrix(nrow=nrow(m),ncol=ncol(m))
for (column in 1:ncol(m)) {
  sum <- 0
  for (row in 1:nrow(m)) {
    if (is.na(m[row,column])) sum <- 0
    else sum <- sum + 1

    m.output[row,column] <- sum
  }
}

这是最有效的方法吗？我有很多只猫，并且记录了多年的粪便数据。我能通过某种方式按列并行化处理它吗？

数字是如何工作的？为什么Fluffy（例如）会拉1.1次大便？ - gung - Reinstate Monica
2
@gung感谢您花时间充分理解问题。这些是十分之一便士的权重。因此，1.1意味着它的重量相当于11个便士。 - dvmlls
哦，鉴于 Fluffy 在1月4日拉了1.4克的屎，你怎么知道它（？）拉了4次呢？Mister Cuddles在3号拉了1.1克一次，但在4号拉了两次，每次1.0克；我没看出来规律。 - gung - Reinstate Monica
@gung 在1月4日，Fluffy已经连续4天拉屎了。因为Mister Cuddles在1月2日没有拉屎，所以1月3日是她第一次连续拉屎的日子。 - dvmlls
我明白了，所以问题是，“如何计算每组NA值之间的非NA值数量”，对吗？ - gung - Reinstate Monica
1
@gung 更像是“你如何找到自上一个NA值以来的天数”？ - dvmlls
7个回答
11

这里所有的答案都太复杂了（包括我自己之前的回答，如下所述）。Reduce系列的答案只是在单个函数调用中掩盖了一个for循环。我喜欢Roland和Ananda的答案，但我认为两者都有点过于复杂。

因此，这里有一个简单的向量化解决方案：

reset <- function(x) {
    s <- seq_along(x)
    s[!is.na(x)] <- 0
    seq_along(x) - cummax(s)
}

> apply(m, 2, reset)
     fluffy misterCuddles
[1,]      1             1
[2,]      2             0
[3,]      3             1
[4,]      4             2

它也适用于Roland的示例：

m2 <- cbind(fluffy=c(NA,1.1,1.2,1.3,1.4,1.0,2),
           misterCuddles=c(NA,1.3,2,NA,NA,1.1,NA))

> apply(m2, 2, reset)
     fluffy misterCuddles
[1,]      0             0
[2,]      1             1
[3,]      2             2
[4,]      3             0
[5,]      4             0
[6,]      5             1
[7,]      6             0

之前的代码没有向量化，但仍然可行：

pooprun <- function(x){
    z <- numeric(length=length(x))
    count <- 0
    for(i in 1:length(x)){
        if(is.na(x[i]))
            count <- 0
        else
            count <- + count + 1
        z[i] <- count
    }
    return(z)
}
apply(m, 2, pooprun)

> apply(m, 2, pooprun)
     fluffy misterCuddles
[1,]      1             1
[2,]      2             0
[3,]      3             1
[4,]      4             2

基准测试

在这里，我只是把每个人的回答都包装在一个函数调用中（根据他们的名字）。

> library(microbenchmark)
> microbenchmark(alexis(), hadley(), thomas(), matthew(), thomasloop(), usobi(), ananda(), times=1000)
Unit: microseconds
         expr     min       lq   median       uq       max neval
     alexis()   1.540   4.6200   5.3890   6.1590   372.185  1000
     hadley()  87.755   92.758   94.298  96.6075  1767.012  1000
     thomas()  92.373  99.6860 102.7655 106.6140   315.223  1000
    matthew() 128.168 136.2505 139.7150 145.4880  5196.344  1000
 thomasloop() 133.556 141.6390 145.1030 150.4920 84131.427  1000
      usobi() 148.182 159.9210 164.7320 174.1620  5010.445  1000
     ananda() 720.507 742.4460 763.6140 801.3335  5858.733  1000

以下是Roland示例数据的结果：

> microbenchmark(alexis(), hadley(), thomas(), matthew(), thomasloop(), usobi(), ananda(), times=1000)
Unit: microseconds
         expr     min       lq   median       uq      max neval
     alexis()   2.310   5.3890   6.1590   6.9290   75.438  1000
     hadley()  75.053   78.902   80.058   83.136 1747.767  1000
     thomas()  90.834  97.3770 100.2640 104.3050  358.329  1000
    matthew() 139.715 149.7210 154.3405 161.2680 5084.728  1000
 thomasloop() 144.718 155.4950 159.7280 167.4260 5182.103  1000
      usobi() 177.048 188.5945 194.3680 210.9180 5360.306  1000
     ananda() 705.881 729.9370 753.4150 778.8175 8226.936  1000
注意：Alexis和Hadley的解决方案在我的机器上定义为函数需要一些时间，而其他解决方案可以直接使用，但除此之外，Alexis的解决方案是明显的赢家。
注意：Alexis和Hadley的解决方案在我的机器上定义为函数需要一些时间，而其他解决方案可以直接使用，但除此之外，Alexis的解决方案是明显的赢家。
1
你的向量化解决方案不起作用。请在一个包含多个NA值的向量上尝试它。 - Roland
1
在我的2500x2500矩阵上尝试微基准测试：在Windows上，您的解决方案比Ananda的快10倍，而Ananda的解决方案又比Usobi的快10倍。在Linux上使用64个核心运行mclapply：您的速度略高于Ananda（3.1秒对3.8秒），两者都比Usobi的速度快约两倍（7.7秒）。 - dvmlls
5

这应该可以工作。请注意，您的每只猫都是独立的个体，因此您可以将数据框转换为列表并使用mclapply，它使用并行处理方法。

count <- function(y,x){
  if(is.na(x)) return(0)
  return (y + 1)
}

oneCat = m[,1]

Reduce(count,oneCat,init=0,accumulate=TRUE)[-1]

编辑：以下是完整的答案。

count <- function(x,y){
 if(is.na(y)) return(0)
 return (x + 1)
}

mclapply(as.data.frame(m),Reduce,f=count,init=0,accumulate=TRUE)

EDIT2: 主要的问题是在开头会出现多余的0。
result = mclapply(as.data.frame(m),Reduce,f=count,init=0,accumulate=TRUE)
finalResult = do.call('cbind',result)[-1,]
rownames(finalResult) = rownames(m)

做得到。
这很好，你应该“应用”它，而且你似乎没有使用“mclapply”，尽管谈论了它...？ - Thomas
此外，在底层，这仍然只是一个“for”循环。 - Thomas
主要的点是使用Reduce函数并设置参数accumulate=TRUE。 - Usobi
finalSet = do.call('cbind',resultOfMcLapply) - Usobi
2
@Usobi 不，它与 cbind(a,b,c,d) 相同，因为它将列表内容传递给函数参数。 - Roland
显示剩余5条评论
4

另一种选择，类似于 @Usobi 使用 Reduce，但采用略有不同的方法：

apply(!is.na(m), 2, Reduce, f=function(x,y) if (y) x + y else y, accumulate=TRUE)
#      fluffy misterCuddles
# [1,]      1             1
# [2,]      2             0
# [3,]      3             1
# [4,]      4             2
哦，这样更好，我的方法确实在第一行得到了额外的0。 - Usobi
4

我曾经从 这里 保存了一小段代码，可以几乎完美地解决类似这样的问题：

countReset <- function(x) {
  x[!is.na(x)] <- 1
  y <- ave(x, rev(cumsum(rev(is.na(x)))), FUN=cumsum)
  y[is.na(y)] <- 0
  y
}
apply(m, 2, countReset)
#            fluffy misterCuddles
# 2013-01-01      1             1
# 2013-01-02      2             0
# 2013-01-03      3             1
# 2013-01-04      4             2
这段程序相关内容有些晦涩，所以我要确保自己理解了其中发生的事情。在NA值上进行cumsum操作会将非NA值分组。然后ave函数会对每个组运行cumsum操作。 - dvmlls
@davez0r，是的。您可以在此处阅读Bill Dunlap的解释（他还分享了一种更简单的方法来创建组）：http://r.789695.n4.nabble.com/partial-cumsum-tp899789p899795.html - A5C1D2H2I1M1N2O1R2T1
注意：折断了。在我的情况下，我需要调用 rev 吗？我不认为需要。测试中... - dvmlls
4

由于我正在适应使用.Call，这里有另一个看起来行得通且可能很快的想法（不过别相信我的话，我的技能不可靠！）：

library(inline)  #use "inline" package for convenience

f <- cfunction(sig = c(R_mat = "numeric", R_dims = "integer"), body = '
 R_len_t *dims = INTEGER(R_dims);
 R_len_t rows = dims[0], cols = dims[1];
 double *mat = REAL(R_mat);

 SEXP ans;
 PROTECT(ans = allocMatrix(INTSXP, rows, cols));
 R_len_t *pans = INTEGER(ans);

 for(int ic = 0; ic < cols; ic++)
  {
   pans[0 + ic*rows] = ISNA(mat[0 + ic*rows]) ? 0 : 1;

   for(int ir = 1; ir < rows; ir++)
    {
     if(ISNA(mat[ir + ic*rows]))
      {
       pans[ir + ic*rows] = 0;
      }else
      {
       if(!ISNA(mat[(ir - 1) + ic*rows]))
        {
         pans[ir + ic*rows] = pans[(ir - 1) + ic*rows] + 1;
        }else
        {
         pans[ir + ic*rows] = 1;
        }
      }
    }
  }

 UNPROTECT(1);

 return(ans);
')

f(m, dim(m))
#     [,1] [,2]
#[1,]    1    1
#[2,]    2    0
#[3,]    3    1
#[4,]    4    2
f(mm, dim(mm))   #I named Roland's matrix, mm ; I felt that I had to pass this test!
#     [,1] [,2]
#[1,]    0    0
#[2,]    1    1
#[3,]    2    2
#[4,]    3    0
#[5,]    4    0
#[6,]    5    1
#[7,]    6    0
1
@Thomas：太酷了！然而，在更大的矩阵上计时我的函数，速度的相对差异会降低。例如，看看 mat = matrix(as.vector(m), nrow = 1e3, ncol = 1e3)；microbenchmark(f(m, dim(m)), apply(m, 2, reset), times = 10)；microbenchmark(f(mat, dim(mat)), apply(mat, 2, reset), times = 10)。 - alexis_laz
1
还是非常快的。 - Thomas
3

所以解决这个问题有两个部分：

  1. 一个函数，接受每只猫的向量并返回一个向量，在每个日期告诉我自上次NA以来过了多少天
  2. 一个函数，接受一个NxM矩阵并返回一个NxM矩阵，对每列应用函数（1）

对于（2），我从 @Usobi 的答案进行了改编：

daysSinceLastNA <- function(matrix, vectorFunction, cores=1) {
  listResult <- mclapply(as.data.frame(matrix), vectorFunction, mc.cores=cores)
  result <- do.call('cbind', listResult)
  rownames(result) <- rownames(matrix)
  result
}

对于（1），我有两个解决方案：
@ananda-mahto的解决方案：
daysSinceLastNA_1 <- function(vector) {
  vector[!is.na(vector)] <- 1
  result <- ave(vector, rev(cumsum(rev(is.na(vector)))), FUN=cumsum)
  result[is.na(result)] <- 0
  result
}

@Usobi的解决方案：
daysSinceLastNA_2 <- function(vector) {
  reduction <- function(total, additional) ifelse(is.na(additional), 0, total + 1)
  Reduce(reduction, vector, init=0, accumulate=TRUE)[-1]
}

那么我像这样调用它们：

> system.time(result1 <- daysSinceLastNA (test, daysSinceLastNA_1 ))
   user  system elapsed 
   5.40    0.01    5.42 
> system.time(result2 <- daysSinceLastNA (test, daysSinceLastNA_2 ))
   user  system elapsed 
  58.02    0.00   58.03

在我的测试数据集上，大约是一个2500x2500矩阵，第一种方法比第二种方法快一个数量级。

如果我在具有64个核心的Linux上运行，解决方案（1）需要2秒钟，而解决方案（2）需要6秒钟。

非常有趣。那个向量化的算法真的值得记住。 - Usobi
1
使用 library(microbenchmark) 获取更可靠的基准测试结果... 因为 system.time 只会考虑一次执行，所以它的结果不够稳定。 - Thomas
3

对于这种可以使用for循环轻松解决的问题，我认为Rcpp是一个非常自然的答案。

library(Rcpp)

cppFunction("NumericVector cumsum2(NumericVector x) {
  int n = x.length();
  NumericVector out(x);

  for(int i = 0; i < n; ++i) {
    if (NumericVector::is_na(x[i]) || i == 0) {
      x[i] = 0;
    } else {
      x[i] = x[i - 1] + 1;
    }
  }

  return out;
}")

这段代码需要比等效的R代码更多的簿记工作，但函数的大部分是一个非常简单的 for 循环。

然后您可以像任何其他矢量化函数一样在 R 中应用：

m2 <- cbind(
  fluffy=c(NA,1.1,1.2,1.3,1.4,1.0,2),
  misterCuddles=c(NA,1.3,2,NA,NA,1.1,NA)
)

apply(m2, 2, cumsum2)

当然，你可以让C++代码迭代矩阵的列，但我认为既然这在R中已经很容易表达了，那么你最好使用内置工具。

1
应用于原始数据，这会产生错误的结果（您应该删除 || i==0 逻辑）。我还将其添加到我的答案中的基准测试中。 - Thomas
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接