平滑序列而不使用 R 中的循环

Question

我正在用 R 中的一篇学术论文（参见末尾的引文）实现一种统计方法。我认为有一种方法可以在不使用循环的情况下执行其中一个步骤，但我在决定如何攻击时遇到了麻烦它。

此方法对具有三个变量的数据框进行操作：x、n 和 p。只有当 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

smooth 会将 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 块按块。如果没有某种循环，我看不出如何做到这一点。某些解决方案可能会将循环隐藏在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

Answer 2

按照上面的 Glen_b，Hamilton 论文中描述的内容相当于 CRAN 包中的 gpava isotone。