3

私は一連のデータを持っています。これらは分子動力学シミュレーションによって得られたものであり、したがって時間的に連続しており、ある程度相関しています。データの平均として平均を計算できます。この方法で計算された平均に関連するエラーを推定したいと思います。

この本によると、「統計的非効率性」、またはシリーズ内のデータの大まかな相関時間を計算する必要があります。このために、シリーズをさまざまな長さのブロックに分割し、ブロックの長さ (t_b) ごとに、ブロック平均の分散 (v_b) を分割する必要があります。次に、シリーズ全体の分散が v_a (つまり、t_b=1 の場合の v_b) である場合、t_b は無限大になる傾向があるため、(t_b*v_b/v_a) の制限を取得する必要があり、それが非効率 s です。 .

次に、平均の誤差は sqrt(v_a*s/N) です。ここで、N は点の総数です。したがって、これは、s ポイントごとに 1 つだけが無相関であることを意味します。

これは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 まで続きます。もちろん、分析する必要があるデータは 2 番目の列です。

4

3 に答える 3

2

一連xのデータ (たとえば、2 番目の列のデータ) を保持しているとします。

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 に答える