Ten million random walkers miss Pólya's constant by six error bars
A random walk on the 3D integer grid steps to one of its six neighbours, each with the same chance. Pólya proved in 1921 that such a walk doesn't always come back to where it started (in 1D and 2D it does). The chance that it ever returns is about 0.340537. I tried to get that number by brute force: 10 million walkers, each allowed 65,536 steps. The noise came out small, ±0.00015, but the answer was 0.33962, six error bars too low. Every walker that would have come home after step 65,536 counted as one that never did, and there are a lot of them.
The run
A short JavaScript loop: a 128-bit generator (sfc32), and each 32-bit draw below 612 becomes 12 steps. It records the step at which each walker first lands back on the origin. Eighty jobs of 125,000 walkers ran nine at a time, about 3.6 CPU-hours over 25 minutes.
To check the simulation I also computed the exact chance of returning by step n. The chance of being at the origin at step 2m splits into the steps that go along the x axis (a 1D walk) and the rest (a walk in the y–z plane, which is two independent 1D walks rotated 45°). First returns then come from the renewal equation u = 1 + f ∗ u. That's about 20 lines of numpy and takes 25 seconds up to step 65,536.
| steps n | simulated P(back by n) | exact | difference / error bar | simulated + 0.287/√n |
|---|---|---|---|---|
| 16 | 0.27053 | 0.27028 | +1.8 | 0.34227 |
| 64 | 0.30503 | 0.30485 | +1.2 | 0.34090 |
| 256 | 0.32284 | 0.32263 | +1.4 | 0.34077 |
| 1,024 | 0.33178 | 0.33157 | +1.4 | 0.34074 |
| 4,096 | 0.33626 | 0.33605 | +1.4 | 0.34074 |
| 16,384 | 0.33849 | 0.33830 | +1.3 | 0.34073 |
| 65,536 | 0.33962 | 0.33942 | +1.3 | 0.34074 |
The simulation matches the exact numbers at every length. Its excess of about 1.3 error bars is a single bit of luck carried along: the column is cumulative, and most of the excess is already there at step 2 (0.16681 against 1/6, +1.2).
The tail you can't wait out
The walkers still out at step n come home slowly. The chance of standing on the origin at step 2m falls like m−3/2, so the return chance still missing after n steps falls only like 1/√n. Working the constant out gives 0.287/√n: (1−p)² × 4√2 × (3/4π)3/2. At 65,536 steps that's 0.00112, seven times the noise. Each extra digit costs 100 times more steps per walker, so running longer doesn't really help.
Adding the formula back fixes it. The last column flattens out by about 1,000 steps and stays at 0.34074 ± 0.00015, which is within 1.3 error bars of the true 0.340537, the same lucky +1.3 as before. Applied to the exact column, the same correction gives 0.340537 to six places.
Try it
The same walk runs below, in your browser. Each walker either comes home or gets cut off at the step limit. The blue line is the share back by each step, with a band of two error bars. The dashed line is the true chance of ever coming back. In 3D, the orange line adds the 0.287/√n tail back on. Leave it running and watch the error bars shrink while the gap stays put. Then cut the step limit and watch the gap grow.
3D, limit + tail: –
1D and 2D come back for sure, but 2D takes its time: with 65,536 steps only about 77% of walkers are home. The missing share there falls like 1/log n, so no step limit you could afford gets close.
So what
In a Monte Carlo run, the error bar only covers the luck of the draw. It says nothing about stopping each walker early, and here that cut was the bigger error by a factor of seven. Adding more walkers would have shrunk the error bar around the wrong number. If you know how the tail falls off, add it back. If you don't, at least plot the estimate against the cutoff and see whether it has stopped moving. Here it hadn't, even at 65,536 steps.
Method: one machine, 10,000,000 walkers in 80 independently seeded streams, cutoff 216 steps. Error bars are binomial standard errors. The true value is Watson's 1939 closed form.