統計無效率（塊平均值）

我有一系列數據，這些數據是通過分子動力學模擬獲得的，因此在時間上是連續的並且在一定程度上相關。我可以將平均值計算爲數據的平均值，我想估計與這種方式計算的平均值相關的誤差。統計無效率（塊平均值）

根據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是點的總數。所以，這意味着每個點只有一個是不相關的。

我認爲這可以用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。當然，我需要分析的數據是第二欄。

來源

2012-10-02 Jellby

你可以發佈你創建的幾行示例數據嗎？這將幫助我們創造更好的解決方案。 – TARehman

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

statistical inefficiency of a time series

來源

2012-11-29 08:21:25

我不認爲這是OP所要求的。 –

@RomanLuštrik如果你認爲你知道OP要求什麼，只需發佈你的答案。 –

有幾種不同的方法來計算統計效率低下，或集成的自相關時間。 R中最簡單的就是CODA軟件包。它們有一個函數effectiveSize，它給出了有效的樣本大小，即樣本總數除以統計無效率。平均值的標準偏差的漸近估計量爲sd(x)/sqrt(effectiveSize(x))。

require('coda') 
n_eff = effectiveSize(x)

來源

2014-12-23 09:03:21

好吧，提出問題永遠不會太晚，不是嗎？由於我自己正在做一些分子模擬，所以我一步步解決了這個問題，但是並沒有看到這個線程。我發現艾倫& Tildesley提出的方法與現代錯誤分析方法相比似乎有些過時。儘管如此，本書的其餘部分還是很值得一看。

雖然Sunhwan喬的答案是正確的關於塊平均線法，關於錯誤的分析，你可以找到其他的方法，如jacknife和引導方法（密切相關，彼此）位置：http://www.helsinki.fi/~rummukai/lectures/montecarlo_oulu/lectures/mc_notes5.pdf

總之，自舉方法，您可以根據您的數據製作一系列隨機的人造樣本，並計算您想要的新樣本的值。我寫了一小段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-29 09:09:27

統計無效率（塊平均值）

回答

相關問題