Monte Carlo Simulations in Physics

Three short Monte Carlo simulations, each checked against an exact result and analysed with least-squares fits and error estimates. Written in Python with NumPy, SciPy and Matplotlib.

Script	Physics	Key result
random_walk.py	Diffusion	Slope of ⟨x²⟩ vs N = 1.002 ± 0.011 (theory: 1)
radioactive_decay.py	Exponential decay	Half-life recovered to within 0.1%; scatter falls as 1/√N₀
bump_hunt.py	Resonance search	125 GeV peak fitted at 124.85 ± 0.22 GeV

Quick start

pip install -r requirements.txt
python3 random_walk.py
python3 radioactive_decay.py
python3 bump_hunt.py

Each script prints its results and saves a plot in the repository's plots/ folder, whichever directory you run it from. Random seeds are fixed, so you get the same numbers shown here.

To run the tests (theory checks, output shapes and values, reproducibility, parameter checks, and the documented numbers below):

python3 -m pytest

---

1. Random walk and diffusion

10,000 walkers each take 1,000 steps of +1 or −1. After $N$ steps the position is $x = s_1 + s_2 + \dots + s_N$ with each $s_i = \pm 1$. Squaring and averaging, the cross terms $\langle s_i s_j \rangle$ vanish and each $s_i^2 = 1$, so the mean squared distance grows exactly linearly with the number of steps:

$$\langle x^2 \rangle = N, \qquad x_{\text{rms}} = \sqrt{N}$$

By the central limit theorem, the final positions are approximately Gaussian with variance $N$.

• Fitted slope: 1.002 ± 0.011 (theory: 1)
• Variance of final positions: 1008 (theory: 1000)

Note on the error bar. Every point on the $\langle x^2 \rangle$ curve uses the same walkers, so the points are strongly correlated and the standard least-squares error is far too small. Instead, the walkers are split into 10 independent groups, each group is fitted separately, and the uncertainty comes from the scatter between groups.

2. Radioactive decay

Start with $N_0$ nuclei. In each time step, every surviving nucleus decays with probability $p = 0.01$, so the number decaying in one step is a binomial random number, $\text{Binomial}(N_{\text{survivors}}, p)$. On average this gives

$$N(t) = N_0\, e^{-\lambda t}, \qquad \lambda = -\ln(1 - p), \qquad t_{1/2} = \frac{\ln 2}{\lambda}$$

so the true half-life is 68.97 steps. The half-life is extracted with a simple unweighted least-squares fit of that curve (see Limitations).

One experiment with N₀ = 1,000 gives a fitted half-life of 68.91 steps (true: 68.97).

Repeating each experiment 200 times shows how the result scatters:

N₀	Fitted half-life (mean ± spread)	Relative spread
100	69.23 ± 7.66	11.1%
1,000	69.04 ± 2.52	3.7%
10,000	68.93 ± 0.72	1.0%

Each tenfold increase in $N_0$ shrinks the spread by about a factor of $\sqrt{10} \approx 3$, as expected for counting statistics, where the relative spread scales as $1/\sqrt{N_0}$.

3. Toy resonance search ("bump hunt")

A simplified version of how a new particle shows up in collider data. 10,000 background events follow a falling exponential mass spectrum between 100 and 160 GeV. Signal events form a Gaussian peak at 125 GeV with a 2 GeV width. The masses are histogrammed in 1 GeV bins and fitted with background + peak, using $\sqrt{n_i}$ as the error on bin $i$. The fit model, in events per GeV, is

$$f(m) = b_0\, e^{-(m - 100)/\tau} + \frac{s}{\sqrt{2\pi}\,\sigma}\, e^{-(m - \mu)^2 / 2\sigma^2}$$

with background normalisation $b_0$ and slope $\tau$, signal yield $s$, peak mass $\mu$ and width $\sigma$ all free.

Signal events	Fitted mass (GeV)	Fitted signal	Rough S/√B
500	124.85 ± 0.22	574 ± 55	14.9
200	125.33 ± 0.53	136 ± 41	4.0
50	124.92 ± 1.48	56 ± 43	1.5

With 500 signal events the peak is obvious. With 200 it is still visible, but the mass uncertainty more than doubles. With 50 it disappears into the background: the fitted signal is consistent with zero, and the fit lands near 125 GeV mainly because that is where it was told to start looking.

$S/\sqrt{B}$ is a rough toy estimate: fitted signal events $S$ divided by the square root of fitted background events $B$ within $\mu \pm 2\sigma$ of the peak. It is not a formal discovery significance or a p-value (see Limitations).

Limitations

These are teaching-scale simulations. The main statistical simplifications:

• Random-walk Gaussian is an approximation. The exact distribution of final positions is binomial and only takes even values after an even number of steps. The Gaussian curve is its central-limit-theorem approximation, which is close for 1,000 steps but not exact, especially far out in the tails.
• Decay fitting is illustrative. The fit treats every point of the survivor curve as independent with equal weight. In reality the points are correlated (each count is the previous count minus new decays, so a fluctuation carries forward to every later point) and their scatter changes over time, so the fit's own error estimate would not be reliable and is not used. The spread quoted above comes from repeating the experiment 200 times instead. A rigorous analysis would use a maximum-likelihood fit of the decay times.
• $S/\sqrt{B}$ is a rough toy estimate. It ignores the uncertainty on the background shape, uses $\sqrt{n_i}$ errors that are only approximate in low-count bins, and does not account for the look-elsewhere effect (searching for a peak at any mass, not only 125 GeV). It is not a formal discovery significance and not a p-value. Real searches use likelihood-ratio tests.

Methods used

Monte Carlo sampling (uniform, binomial, exponential and Gaussian random numbers), histogramming, linear and non-linear least-squares fitting (numpy.polyfit, scipy.optimize.curve_fit), and statistical error estimation.

License

MIT