r - 在 R 中不使用循环平滑序列

Question

我正在 R 中实现一篇学术论文(请参阅引文末尾)的统计方法。我认为有一种方法可以在不使用循环的情况下执行其中一个步骤，但我无法决定如何攻击它。

此方法对具有三个变量的数据框进行操作:x、n 和 p。仅当所有 i 的 p[i] <= p[i+1] 时，它才能运行。如果一对点违反了这一点，则通过将 p[i] 和 p[i+1] 设置为等于其加权平均值来平滑它们 (n[i]*p[i]+n[i+1]*p[i+1])/(n[i]+n[i+1]) 迭代此平滑，直到 p_i 成为非递减序列。

这种平滑的问题在于:a) R 中的循环形式不好，b) 如果连续有多个点，使得 p_i > p_(i+1) >= p_(i+2)，则方法可能无法终止或需要很长时间才能收敛。例如，如果发生这样的序列:

x  n  p
2  10 0.6
5  10 0.5
10 10 0.5

平滑会将 p 的前两个值设置为 0.55，然后将后两个值设置为 0.525，然后将前两个值设置为 0.5325，依此类推并永远循环(或者如果我幸运的话达到显着性极限)无数次迭代)。应该有一种数学上等效但更有效的方法来做到这一点，通过识别相邻的递减数据点并将它们作为一个组进行平均，但我不确定如何在 R 中实现这一点。

如果您需要更多背景知识，相关论文是 Martin A. Hamilton, Rosemarie C. Russo, Robert V. Thurston. "Trimmed Spearman-Karber method for estimating median lethal concentrations in toxicity bioassays." Environ. Sci. Technol., 1977, 11 (7), pp 714–719 。我指的是第 716 页的“第一步”部分。

Answer 1

据我了解该算法，您需要找到 p 减少的位置，并从每个位置开始，找出(累积)加权平均值减少的时间，以便 p 可以逐 block 更新。我不明白如何在没有某种循环的情况下完成此操作。某些解决方案可能会将循环隐藏在 lapply 或等效项下，但恕我直言，这是足够复杂的算法之一，我更喜欢一个好的旧循环。您可能会损失一点效率，但代码可读性很好。我的尝试，使用 while 循环:

smooth.p <- function(df) {

   while (any(diff(df$p) < 0)) {

      # where does it start decreasing
      idx <- which(diff(df$p) < 0)[1]

      # from there, compute the cumulative weighted average
      sub <- df[idx:nrow(df), ]
      cuml.wavg <- cumsum(sub$n * sub$p) / cumsum(sub$n)

      # and see for how long it is decreasing
      bad.streak.len <- rle(diff(cuml.wavg) <= 0)$lengths[1]

      # these are the indices for the block to average
      block.idx <- seq(from = idx, length = bad.streak.len + 1)

      # compute and apply the average p
      df$p[block.idx] <- sum(df$p[block.idx] * df$n[block.idx]) /
                     sum(df$n[block.idx])
   }
   return(df)
}

这里是一些数据，包括您建议的粗略补丁:

df <- data.frame(x = 1:9,
                 n = rep(1, 9),
                 p = c(0.1, 0.3, 0.2, 0.6, 0.5, 0.5, 0.8, 1.0, 0.9))
df
#   x n   p
# 1 1 1 0.1
# 2 2 1 0.3
# 3 3 1 0.2
# 4 4 1 0.6
# 5 5 1 0.5
# 6 6 1 0.5
# 7 7 1 0.8
# 8 8 1 1.0
# 9 9 1 0.9

输出:

smooth.p(df)
#   x n         p
# 1 1 1 0.1000000
# 2 2 1 0.2500000
# 3 3 1 0.2500000
# 4 4 1 0.5333333
# 5 5 1 0.5333333
# 6 6 1 0.5333333
# 7 7 1 0.8000000
# 8 8 1 0.9500000
# 9 9 1 0.9500000