3

我有一个 OpenBUGS 模型，它使用随时间 (x.values) 观察到的数据 (y.values) 来模拟多次运行 (~100000)，并为每次运行提供新的 y 值 (y.est) 估计值。观察到的数据显示出从最大值明显下降。

我想跟踪每次运行从最大丰度 (T.max) 下降到最大丰度的 10% (T.10%) 所需的时间长度。因为最大丰度值随运行而变化，该最大值的 10% 也将随运行而变化，因此 T.10% 将随运行而变化。

设置一个参数来存储 T.max 很容易，它不会因运行而异，因为最大值远远大于任何其他值。

我想不通的是如何存储 y-est 值和 T.10% 的交集。

我的第一次尝试是使用以下函数确定每个 y-est 值是高于还是低于 T.10% step()：

above.below[i] <- step(T.10% - y.est[i])

这会为每个 y.est 值生成一串一和零（例如，0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1 , 1, 1, 1, 1, 1 等）如果每次运行只是从最大值连续下降到最小值，我可以使用该rank()函数来确定有多少above.below[i]值出现在 T.10% 以上：

decline.length <- rank(above.below[1:N], 0)

在此示例中，decline.length将等于上述字符串中“0”的数量，即 9。不幸的是，y-est 值偶尔会在下降到 T.10% 以下后显示增长期。因此，above.below值向量可以如下所示： 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 , 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 1 等。因此，decline.length给定向量中的后续 0，将等于 14 而不是 9。

我想要做的是弄清楚如何在above.below第一个“1”之前只存储“0”的数量；above.below[1:10]而不是above.below[1:N]. above.below不幸的是，第一个“1”并不总是出现在第 10 个时间步长，所以我需要在模拟过程中使每次运行的变化范围达到最大。

我正在努力在 OpenBUGS 中实现这一点，因为它是一种非过程语言，但我认为它可以做到，我只是不知道这样做的诀窍。我希望更熟悉step()和rank()功能的人可以提供一些专家建议。

非常感谢任何指导！

4

1 回答 1

2

提供给我的两个解决方案：

1）计算每个时间步的累积和：

for (i in 1:N){
    above.below[i] <- step(T.10% - y.est[i])
    cum.above.below[i] <- sum(above.below[1:i])
}

decline.length <- rank(cum.above.below[1:N], 0)

2）直接计算每一年是在阈值之上还是之下，没有1和0：

for(i in 1:N){
    above.below[i] <- step(T.10% - y.est[i])
    dummy[i] <- above.below[i] * i + (1 - above.below[i]) * (N+1)
}

decline.length <- ranked(dummy[], 1)

当为 1时dummy也是如此，当为 0 时为 N+1 。iabove.belowdummyabove.below

于 2014-10-17T14:23:12.373 回答