使用 RStorm 计算多个 data.table 列的 Welford 方差答案

【问题标题】：Computing Welford's variance for multiple data.table columns using RStorm使用 RStorm 计算多个 data.table 列的 Welford 方差
【发布时间】：2015-12-31 20:10:54
【问题描述】：

鉴于以下data.tabledt：

我想使用 Welford's Method 和 RStorm 包工具计算列 a 和 b 按列 i 分组的运行方差。我遵循page 4 of RStorm's vignette 上的示例并通读了introductory paper on RStorm，但我无法弄清楚如何使其工作。这是我的代码：

library(RStorm)
dt = data.table(i=c(1,2,2,3,3,3,4,4,4,4), a=c(1,1:2,1:3,1:4), b=c(NA,NA,2,1,2,NA,NA,2,3,NA)
in_cols = c('a','b')
out_cols <- paste0(in_cols, '.var.Welford')
## Calculaing variance using Welford's method
## See: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance
## See: "RStorm: Developing and Testing Streaming Algorithms in R", R Journal Vol 6/1
var.Welford <- function(x, ...) {
    x <- as.numeric(x[1])
    params <- GetHash("params2")
    if (!is.data.frame(params)) {
        params <- list()
        params$M <- params$S <- params$n <- 0
    }
    x <- ifelse(is.na(x), params$M, x)
    n <- params$n + 1
    delta <- (x - params$M)
    M <- params$M + ( delta / (n + 1) )
    S <- params$S + delta*(x - M)
    SetHash("params2", data.frame(n=n,M=M,S=S))
    var <- ifelse(n > 1, S / (n-1), 0)
    TrackRow("var.Welford", data.frame(var = var))
}
computeVarWelford <- function(x) {
    topology <- Topology(as.data.frame(x=as.data.frame(x)))
    topology <- AddBolt(topology, Bolt(var.Welford, listen = 0))
    result <- RStorm(topology)
    # GetTrack('var.Welford', result)
    result$track$var.Welford
}

## Execute:
dt[, eval(out_cols) := lapply(.SD, function(x) {return(as.list(computeVarWelford(x))[1])})
, by=i, .SDcols = in_cols]

执行上面的代码会将dt 转换为：

    i a  b                       a.var.Welford                       b.var.Welford
 1: 1 1 NA                                   0                                   0
 2: 2 1 NA                                 0,2                   0.000000,2.666667
 3: 2 2  2                                 0,2                   0.000000,2.666667
 4: 3 1  1                         0.0,2.0,2.5                               0,2,1
 5: 3 2  2                         0.0,2.0,2.5                               0,2,1
 6: 3 3 NA                         0.0,2.0,2.5                               0,2,1
 7: 4 1 NA 0.000000,2.000000,2.500000,3.333333 0.000000,2.666667,3.375000,2.250000
 8: 4 2  2 0.000000,2.000000,2.500000,3.333333 0.000000,2.666667,3.375000,2.250000
 9: 4 3  3 0.000000,2.000000,2.500000,3.333333 0.000000,2.666667,3.375000,2.250000
10: 4 4 NA 0.000000,2.000000,2.500000,3.333333 0.000000,2.666667,3.375000,2.250000

从结果中可以清楚地看出，每个 (column,group) 对的整个方差列表都被复制到该 (column,group) 对的 每个元素 中，而不是映射到该（列、组）对的所有元素。这是我真正想要的：

    i a  b     a.var.Welford        b.var.Welford
 1: 1 1 NA     0                    0
 2: 2 1 NA     0                    0
 3: 2 2  2     2                    2.666667
 4: 3 1  1     0.0                  0
 5: 3 2  2     2.0                  2
 6: 3 3 NA     2.5                  1
 7: 4 1 NA     0.000000             0.000000
 8: 4 2  2     2.000000             2.666667
 9: 4 3  3     2.500000             3.375000
10: 4 4 NA     3.333333             2.250000

我真的希望有一个简单的解决方法，但我无法在我的一生中解决这个问题。每次我尝试我认为应该有效的方法时，我最终都会收到来自data.table 的错误说

j=list(...) 中的所有项目都应该是原子向量或列表。如果你是尝试类似 j=list(.SD,newcol=mean(colA)) 然后使用 := by 改为分组（更快），或者之后 cbind 或合并。

我理解这意味着我在lapply(.SD, FUN) 代码中尝试的任何FUN 的返回值的维度与data.table 期望的该组列的维度不对应。

非常感谢任何和所有的帮助。

编辑：好的，解决方案非常简单。我觉得我好笨。但这是以后可能需要的人的答案

## Make sure to use [[]] at the end. My problem came entirely down to using [].
dt[, eval(out_cols) := lapply(.SD, function(x) {return(as.list(computeVarWelford(x))[[1]])})
   , by=i, .SDcols = in_cols]

这就像一个魅力。我得到了我需要的东西：

    i a  b a.var.Welford b.var.Welford
 1: 1 1 NA      0.000000      0.000000
 2: 2 1 NA      0.000000      0.000000
 3: 2 2  2      2.000000      2.666667
 4: 3 1  1      0.000000      0.000000
 5: 3 2  2      2.000000      2.000000
 6: 3 3 NA      2.500000      1.000000
 7: 4 1 NA      0.000000      0.000000
 8: 4 2  2      2.000000      2.666667
 9: 4 3  3      2.500000      3.375000
10: 4 4 NA      3.333333      2.250000

【问题讨论】：

@VeerendraGadekar ：这将给我 k^2 行，每个组最初有 k 行。 cSplit 将分割每一行。我只需要从每个组中拆分第一行，然后将拆分值放入该组的剩余行中。
或者你可以像这样使用rbindlist rbindlist(lapply(split(data, data$i), function(x){cSplit(x[1,], c('a.var.Welford', 'b.var.Welford'), ',', 'long')}))
@VeerendraGadekar 是的，它是0。我修正了错字。您上面的两个代码语句都可以工作....几乎。我在上面的问题中添加了结果作为编辑。如何摆脱奇怪的解析/评估错误？
我没有得到你在那里显示的内容，但你可以使用listCol_l(output, c('a.var.Welford', 'b.var.Welford')) 摆脱它，其中输出将是上述命令检查this 的结果以获取更多信息
@VeerendraGadekar 感谢您的努力。我能够在上游解决我的问题，因此我不再需要每行解析那些奇怪的列表结果。

标签： r statistics data.table grouping lapply

【解决方案1】：

编辑：我不再需要下面的 hack 解决方案。这是解决此问题的代码（注意 [[]] 而不是 [] 修复）：

dt[, eval(out_cols) := lapply(.SD, function(x) {return(as.list(computeVarWelford(x))[[1]])})
   , by=i, .SDcols = in_cols]

OLD：好的，所以我终于找到了让它工作的方法。但我觉得这条路很丑。我现在将接受这个作为我的答案，但是如果有人有更好的解决方案，我很乐意听到它，如果它比我的更好，我会接受它作为这个问题的答案。

解决方案：

out_cols_fixed <- paste0(out_cols, '.fixed')
dt[,eval(out_cols_fixed) := lapply(.SD, function(x) { return(x[1][[1]]) }), by=i, .SDcols = out_cols]
dt[,eval(out_cols) := NULL]
setnames(dt, old = out_cols_fixed, new = out_cols)

dt 所需的结果：

    i a  b a.var.Welford b.var.Welford
 1: 1 1 NA      0.000000      0.000000
 2: 2 1 NA      0.000000      0.000000
 3: 2 2  2      2.000000      2.666667
 4: 3 1  1      0.000000      0.000000
 5: 3 2  2      2.000000      2.000000
 6: 3 3 NA      2.500000      1.000000
 7: 4 1 NA      0.000000      0.000000
 8: 4 2  2      2.000000      2.666667
 9: 4 3  3      2.500000      3.375000
10: 4 4 NA      3.333333      2.250000

我先尝试了以下方法，但没有成功。谁能解释一下为什么？

dt[,eval(out_cols) := lapply(.SD, function(x) { return(x[1][[1]]) }), by=i, .SDcols = out_cols]

运行上面的行时出现以下错误：

[.data.table(dt, , :=(eval(out_cols), lapply(.SD, function(x) { : RHS ('double') 的类型必须与 LHS ('list') 匹配。到 check 和 coerce 对最快的性能影响太大案例。要么更改目标列的类型，要么强制 RHS of := 你自己（例如，使用 1L 而不是 1）

【讨论】：