This iterative series converges and I don't understand why

Question

I am learning about the Metropolis Monte Carlo methods for applications in statistical thermodynamics and I wrote a small portion of code for the purposes of figuring something out. Turns out the dynamics of my loop were more interesting than expected. The loop is:

import numpy as np
tot = 0
k = 0.5
n = 100000
for i in range(n):
    rng = np.random.rand()
    if rng > k:
        tot+=rng
        k = rng
    else:
        tot-=rng
        k = rng
print(tot)

When you run this code for large n, it seems to converge to the figure 1/6*n (so for n = 1000, you would expect to get a value that falls around 166).

I understand why the output value is positive, since for a given random number that is greater than the preceding k value, you can't get a subsequent value for 'rng' that takes away more than what was added. This effectively introduces a positive pressure on the system.

Furthermore, if a large number (a number close to one) was just added, you'll likely end up with a number taken away (since k becomes rng, which was large (close to one) so the chances of the next random number being greater than it is low, so the chances of the subsequent contribution to the 'tot' value being negative is high). This constraint imposes a 'braking force' on the system.

Honestly when I first saw this I thought it would be trivial to decipher, but what appeared to be simple at first has left me scratching my head. If someone that's better than me at maths could find an analytical solution to this (in the limit that n is large), I'd be very happy to see it.

Lastly, perhaps this post is in the wrong thread. If so, please could you let me know of a more suitable location to post it!

Answer 1

Since k is updated regardless, it’s just uniformly distributed. The probability of moving down is of course $k$, with an expected change of $-k/2$; the up possibility is similar. As such,

$$\begin{align} \mathrm{E}(\Delta tot)&=\int_0^1 \left(\frac{-k}2k+\frac{1+k}2(1-k)\right)\,dk\\ &=\int_0^1 \left(\frac12-k^2\right)\,dk=\frac12-\frac13=\frac16 \end{align}$$

Intuitive understanding might be aided by recognizing that the expected change when $k=1/2$ (as initially) is $1/4$, since the two directions are equally likely but the upward movement is on average three times as large.