3

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

根据这本书，我需要计算“统计效率低下”，或者大致是系列中数据的相关时间。为此，我必须将系列划分为不同长度的块，并且对于每个块长度 (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。当然，我需要分析的数据是第二列。

4

3 回答 3

2

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

时间序列的统计效率低下

于 2012-11-29T08:21:25.180 回答
1

有几种不同的方法可以计算统计效率低下或积分自相关时间。在 R 中，最简单的方法是使用 CODA 包。它们有一个函数 ，effectiveSize它为您提供有效样本量，即样本总数除以统计无效率。均值标准差的渐近估计量是sd(x)/sqrt(effectiveSize(x))。

require('coda')
n_eff = effectiveSize(x)
于 2014-12-23T09:03:21.397 回答
1

好吧，提出问题永远不会太晚，不是吗？当我自己做一些分子模拟时，我确实解决了这个问题，但还没有看到这个线程。我发现 Allen & Tildesley 实际提出的方法与现代错误分析方法相比似乎有些过时。本书的其余部分足够好，值得一看。

虽然 Sunhwan Jo 对块平均方法的回答是正确的，但关于错误分析，您可以在此处找到其他方法，例如 jacknife 和 bootstrap 方法（彼此密切相关）：http ://www.helsinki.fi/~rummukai/lectures/montecarlo_oulu /讲座/mc_notes5.pdf

简而言之，使用 bootstrap 方法，您可以从数据中制作一系列随机人工样本，并在新样本上计算您想要的值。我写了一小段 Python 代码来处理一些数据（不要忘记导入 numpy 或我使用的函数）：

def Bootstrap(data):
    B  = 100 # arbitraty number of artificial samplings
    es = 0.
    means = numpy.zeros(B)
    sizeB = data.shape[0]/4 # (assuming you pass a numpy array)
                            # arbitrary bin-size proportional to the one of your
                            # sampling.
    for n in range(B):
        for i in range(sizeB):
        # if data is multi-column array you may have to add the one you use
        # specifically in randint, else it will give you a one dimension array.
        # Check the doc.
            means[n] = means[n] + data[numpy.random.randint(0,high=data.shape[0])] # Assuming your desired value is the mean of the values
                                             # Any calculation is ok.
            means[n] = means[n]/sizeB
    es = numpy.std(means,ddof = 1)
    return es

我知道它可以升级，但这是第一次尝试。使用您的数据，我得到以下信息：

Mean              =  594.84368
Std               =   66.48475
Statistical error =    9.99105

我希望这可以帮助任何在数据统计分析中遇到这个问题的人。如果我错了或其他任何事情（第一篇文章，我不是数学家），欢迎任何更正。

于 2016-03-29T09:09:27.163 回答