r - 统计效率低下( block 平均)

Question

我有一系列数据，这些数据是通过分子动力学模拟获得的，因此在时间上是连续的，并且在某种程度上是相关的。我可以将平均值计算为数据的平均值，我想估计与以这种方式计算的平均值相关的误差。

根据 this book我需要计算“统计效率低下”，或者粗略地计算系列数据的相关时间。为此，我必须将系列划分为不同长度的 block ，并且对于每个 block 长度 (t_b)， block 平均值的方差 (v_b)。那么，如果整个级数的方差为v_a(即t_b=1时的v_b)，我就得求出(t_b*v_b/v_a)t_b趋于无穷大的极限，也就是无效率s .

则均值误差为sqrt(v_a*s/N)，其中N为总点数。所以，这意味着每 s 个点只有一个不相关。

我想这可以用 R 来完成，也许已经有一些包可以做到这一点，但我是 R 的新手。谁能告诉我该怎么做？我已经找到了如何读取数据序列并计算均值和方差。

数据样本，根据要求:

# t(ps) dH/dl(kJ/mol)
0.0000 582.228
0.0100 564.735
0.0200 569.055
0.0300 549.917
0.0400 546.697
0.0500 548.909
0.0600 567.297
0.0700 638.917
0.0800 707.283
0.0900 703.356
0.1000 685.474
0.1100 678.07
0.1200 687.718
0.1300 656.729
0.1400 628.763
0.1500 660.771
0.1600 663.446
0.1700 637.967
0.1800 615.503
0.1900 605.887
0.2000 618.627
0.2100 587.309
0.2200 458.355
0.2300 459.002
0.2400 577.784
0.2500 545.657
0.2600 478.857
0.2700 533.303
0.2800 576.064
0.2900 558.402
0.3000 548.072

... 这一直持续到 500 ps。当然，我需要分析的数据是第二列。

Answer 1

假设 x 保存数据序列(例如，来自第二列的数据)。

v = var(x)
m = mean(x)
n = length(x)

si = c()
for (t in seq(2, 1000)) {
    nblocks = floor(n/t)
    xg = split(x[1:(nblocks*t)], factor(rep(1:nblocks, rep(t, nblocks))))
    v2 = sum((sapply(xg, mean) - m)**2)/nblocks
    #v rather than v1
    si = c(si, t*v2/v)
}
plot(si)

下图是我从我的一些时间序列数据中得到的。当 si 的曲线变得近似平坦(斜率 = 0)时，您就有了 t_b 的下限。参见 http://dx.doi.org/10.1063/1.1638996