R中基于不等式的条件累积和答案

【问题标题】：Inequality based conditional cumulative sum in RR中基于不等式的条件累积和
【发布时间】：2013-11-14 22:50:59
【问题描述】：

我一直在尝试编写这样的代码：对于每个

到目前为止，我想出的最好方法是使用循环。这是一个示例

y=rnorm(10)
x=c(1,1,1,2,2,2,3,3,3,4)
z=c(5,5,6,6,7,7,8,8,9,9)
data=data.frame(y,x,z)
n=10

s=rep(NA,length(unique(x))*length(unique(z)))
dim(s)=c(length(unique(x)),length(unique(z)))
for (i in 1:length(unique(x))){
  for (j in 1:length(unique(z))){
       s[i,j]=sum(y*as.numeric((x<=unique(x)[i]))*
                    as.numeric((z<=unique(z)[j])))
}
}

这样的输出还可以，但是当我的尺寸增加时，这变得低效。因为对于给定的 z，这看起来像是一个条件累积和，所以我 100% 确信有一种更有效的方法可以做到这一点，而无需循环。

你们有什么建议吗？如果我没有 z，我知道我可以使用 data.table：

 s=data[order(x)][,lapply(.SD, sum),by=c("x"), .SDcols=c("y")]
  s=s[,lapply(.SD, cumsum), .SDcols=c("y")]

但有多个索引（x 和 z，而不仅仅是 x）我无法制定程序。

【问题讨论】：

如果您显示一小部分输入输出数据和所需的输出，这将更容易理解。
同意。我现在正在这样做

标签： r data.table plyr

【解决方案1】：

我认为您不需要data.table，因为您为每个组使用了整个“y”。这会更容易通过一些线性代数来完成：

t(y*outer(x, unique(x), '<=')) %*% outer(z, unique(z), '<=')
          [,1]       [,2]       [,3]       [,4]       [,5]
[1,] 0.3538152  0.1762013  0.1762013  0.1762013  0.1762013
[2,] 0.3538152 -0.7308157 -1.2421102 -1.2421102 -1.2421102
[3,] 0.3538152 -0.7308157 -1.2421102 -1.1770919 -1.8315592
[4,] 0.3538152 -0.7308157 -1.2421102 -1.1770919 -4.1171477

这是您的 3 维代码版本：

set.seed(1)
y=rnorm(10)
x=c(1,1,1,2,2,2,3,3,3,4)
z=c(5,5,6,6,7,7,8,8,9,9)
w=c(7,7,8,8,9,9,10,10,11,11)
n=10

s=rep(NA,length(unique(w))*length(unique(z))*length(unique(x)))
dim(s)=c(length(unique(w)),length(unique(z)), length(unique(x)))
for (i in 1:length(unique(w))) {
  for (j in 1:length(unique(z))) {
    for (k in 1:length(unique(x))) {
       s[i,j, k]=sum(y*as.numeric((w<=unique(w)[i]))*
                    as.numeric((z<=unique(z)[j]))*
                    as.numeric((x<=unique(x)[k])))
    }
  }
}

以下是您可以使用与我之前的答案相同的想法来完成此操作的方法：

t1 <- outer(x, unique(x), '<=')
t2 <- outer(z, unique(z), '<=')
t3 <- outer(w, unique(w), '<=')
lapply(seq_along(unique(x)), function(idx) t(y*t1[,idx]*t2) %*% t3)

这里输出的是一个列表（而不是数组），但是输出是一样的，你可以用“s”来比较结果。你应该可以从这里拿走它。

【讨论】：

如果我添加其他索引，比如说4个呢？这就是为什么我倾向于使用 data.table，不是吗？
就像我在答案中所说的那样，您对每个 i,j 都使用了整个“y”。但是 data.table 的 .SD 只会为“by”中的组提供“y”的值。我不知道如何做到这一点。使用我认为的当前方法将其扩展到 4 个以上的索引很简单。
您介意举个例子，您将如何将其扩展到 3 维而不是 2 维？因为这需要使用数组（3 维），而不是标准矩阵（2 维）。非常感谢
@user17645，希望此编辑对您有所帮助。您可以将其扩展到您的数据维度。
@user17645 鉴于这似乎可以有效地解决您的问题，可以扩展到 n 维（这不是 OP 的一部分 - 直到您编辑它 - 您很幸运 Arun 有无限的耐心！）我希望您通过按下绿色对勾来接受答案，以表达您对为此付出的所有努力的感谢？

【解决方案2】：

遵循@Arun 参数，我设法嵌套了两个 lapply 函数，以将解决方案推广到更高维度。

lapply(seq_along(unique(x)), function(idx){lapply(seq_along(unique(r)),
                                              function(idr) t(y*t1[,idx]*t2)%*%(t3
                                              *t4[,idr]))})

为了添加其他维度，我将继续嵌套 lapply 函数。有没有更清洁的方法来做到这一点？

【讨论】：