Block average analysis for correlated data

In this example, a random walk is performed in one dimension with periodic boundaries. The step size is chosen randomly uniform in [-1, 1], and the repeating cell is 25. Use the blocking method to find the standard deviation of the mean from uncorrelated data, and to estimate the correlation time (i.e., the number of steps until the particle could be anywhere in the repeating cell).

For 1e6 steps, the error on the average position is about 0.08 with about 68% confidence, and a correlation time of roughly 2^10 steps.

import unittest
import feasst as fst

class TestBlockAv(unittest.TestCase):
    def test(self, pbc_length=25, step_size=1, num_steps=int(1e6)):
        av_position = fst.MakeAccumulator(fst.args({"max_block_operations": "20"}))
        av_position2 = fst.MakeAccumulator()
        rng = fst.MakeRandomMT19937()
        position = rng.uniform_real(-pbc_length/2, pbc_length/2)
        for step in range(num_steps):
            position += step_size*rng.uniform_real(-1, 1)
            if position < -pbc_length/2:
                position += pbc_length
            elif position > pbc_length/2:
                position -= pbc_length
            av_position.accumulate(position)
            av_position2.accumulate(position)
        for op in range(av_position.max_block_operations()):
            print(op, av_position.block_stdev(op))
        self.assertAlmostEqual(av_position2.block_stdev(), 0.08, delta=0.05)

Run the tests.

unittest.main(argv=[''], verbosity=2, exit=False)

test (__main__.TestBlockAv) ...

0 0.007238781090395475
1 0.009938366458832426
2 0.01369574098784891
3 0.018803556782505806
4 0.02556235262330736
5 0.034206826191159716
6 0.044640794340539275
7 0.05637749798772618
8 0.06655949625147244
9 0.07343506780491323
10 0.07698166109738189
11 0.07860233781069917
12 0.08191576015790028
13 0.08328572894610582
14 0.08253846353334518
15 0.07150823096402974
16 0.06180560019559975
17 0.054061765376740345
18 0.027105746508845847
19 0.0

ok

----------------------------------------------------------------------
Ran 1 test in 2.208s

OK

<unittest.main.TestProgram at 0x7f7d181d6520>

Did this tutorial work as expected? Did you find any inconsistencies or have any comments? Please contact us. Any feedback is appreciated!