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

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

根据 this book 我需要计算“统计效率低下”，或者粗略地计算系列中数据的相关时间。为此，我必须将系列划分为不同长度的块，并且对于每个块长度 (t_b)，块平均值的方差 (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。当然，我需要分析的数据是第二列。

最佳答案

假设 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。

关于r - 统计效率低下(块平均)，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/12694295/