MTH3007b Lecture 7

Me, in the lecture

zzzzz…

This session covers two distinct but related threads: Monte Carlo integration as a numerical quadrature technique, and stochastic processes starting with the random walk and building up to the Ornstein-Uhlenbeck process. Both rely heavily on the machinery of pseudo-random number generation.

Monte Carlo Integration

Monte Carlo integration estimates a definite integral by sampling the integrand at random points rather than at a regular grid.

Pseudo-Random Numbers

Pseudo-random numbers are generated deterministically by an algorithm but pass statistical tests for randomness. In Python, numpy.random provides a uniform distribution on $[a,b]$ via np.random.uniform(a, b, N).

A seed fixes the starting point of the random number generator, giving reproducible results:

Python

import numpy as np np.random.seed(0)

Output

Gaussian (normal) random numbers with mean $\mu$ and standard deviation $\sigma$ are generated by np.random.normal(mu, sigma, N).

1D Integration

To estimate $F={\int }_a^{b}f(x)dx$, draw $N$ uniform samples $x_i\in [a,b]$ and compute:

$$\boxed{F\approx (b-a)⟨f{⟩}_N}$$

where $⟨f{⟩}_N=\frac{1}{N}{\sum }_{i=1}^{N}f(x_i)$ is the sample mean. The error on this estimate is:

$$\boxed{{\sigma }_F=(b-a)\sqrt{\frac{⟨f^2{⟩}_N-⟨f{⟩}_N^2}{N}}}$$

The error scales as $O(N^{-1/2})$, so quadrupling the number of samples halves the error.

Higher-Dimensional Integration

For an integral over a $d$-dimensional hypervolume $V$:

$$F\approx V⟨f{⟩}_N$$

The error is still $O(N^{-1/2})$, regardless of the dimension $d$. This is the key advantage of Monte Carlo over grid-based quadrature, which degrades as dimension increases.

Example: Integrating $f(x)=x^2$ on $[0,1]$

Python

import numpy as np import matplotlib.pyplot as plt number_of_samples = 1000 np.random.seed(0) x_samples = np.random.uniform(0, 1, number_of_samples) f_values = x_samples ** 2 integral_estimate = (1 - 0) * np.mean(f_values) print("F =", integral_estimate) error_estimate = (1 - 0) * np.sqrt( (np.mean(f_values ** 2) - np.mean(f_values) ** 2) / number_of_samples ) print("error =", error_estimate)

Output

The exact answer is ${\int }_0^1{x}^{2}dx=1/3\approx 0.333$.

Stochastic Processes

Discrete Random Walk

The random walk is the simplest discrete stochastic process. At each step, a particle moves $+1$ or $-1$ with equal probability:

$$X_{n+1}=X_n\pm 1$$

The Wiener Process

The Wiener process $W(t)$ is the continuous-time limit of the random walk. It satisfies:

• $W(0)=0$
• $W(t)-W(s)\sim \mathcal{N}(0,t-s)$ for $t>s$
• Increments over non-overlapping intervals are independent

The numerical Euler-Maruyama scheme for the Wiener process is:

$$\boxed{W(t+\Delta t)=W(t)+\sqrt{\Delta t}Z}$$

where $Z\sim \mathcal{N}(0,1)$.

Python

import numpy as np import matplotlib.pyplot as plt time_step = 0.01 time_end = 10.0 number_of_steps = int(round(time_end / time_step)) W_path = np.zeros(number_of_steps + 1) t_values = np.zeros(number_of_steps + 1) for step_index in range(number_of_steps): t_values[step_index + 1] = t_values[step_index] + time_step W_path[step_index + 1] = W_path[step_index] + np.sqrt(time_step) * np.random.normal() plt.plot(t_values, W_path) plt.xlabel('t') plt.ylabel('W(t)') plt.tight_layout() plt.show()

Output

The Ornstein-Uhlenbeck Process

The Ornstein-Uhlenbeck process (OU process) adds a linear restoring force. The continuous equation is:

$$dV=-kVdt+dW$$

The numerical update rule is:

$$\boxed{V(t+\Delta t)=(1-k\Delta t)V(t)+\sqrt{\Delta t}Z}$$

Python

import numpy as np mean_reversion_rate = 1.0 time_step = 0.01 time_end = 10.0 number_of_steps = int(round(time_end / time_step)) velocity_path = np.zeros(number_of_steps + 1) t_values = np.zeros(number_of_steps + 1) velocity_path[0] = 0.0 for step_index in range(number_of_steps): t_values[step_index + 1] = t_values[step_index] + time_step velocity_path[step_index + 1] = ( (1 - mean_reversion_rate * time_step) * velocity_path[step_index] + np.sqrt(time_step) * np.random.normal() )

Output

The Langevin Equation

The Langevin equation is the physical framing of the OU process:

$$\frac{dV}{dt}=-kV+\xi (t)$$

where $\xi (t)$ is white noise with the correlation:

$$⟨\xi (t)\xi (t^{\prime })⟩=2\delta (t-t^{\prime })$$

Numerically, the white noise term is approximated as $\xi \to Z/\sqrt{\Delta t}$, which recovers the OU update rule above.

---

Pre-Lecture Notes from University Notes

• Monte Carlo integration: $F\approx (b-a)⟨f{⟩}_N$ with error $(b-a)\sqrt{(⟨f^2⟩-⟨f{⟩}^2)/N}$.
• Error is $O(N^{-1/2})$ regardless of dimension - decisive advantage over grid-based methods in high dimensions.
• np.random.seed(0) for reproducibility; np.random.uniform(a,b,N) for uniform samples.
• Wiener process: $W(0)=0$, increments $W(t)-W(s)\sim \mathcal{N}(0,t-s)$, independent increments.
• Euler-Maruyama update: $W(t+\Delta t)=W(t)+\sqrt{\Delta t}Z$, $Z\sim \mathcal{N}(0,1)$.
• OU process: $dV=-kVdt+dW$; update $(1-k\Delta t)V+\sqrt{\Delta t}Z$.
• Langevin: $\dot{V}=-kV+\xi (t)$ with $⟨\xi (t)\xi (t^{\prime })⟩=2\delta (t-t^{\prime })$; numerically $\xi \to Z/\sqrt{\Delta t}$.
• Next session: Dirac delta initial conditions for the diffusion equation, and first-passage time problems.