846 Fundamental

Monte Carlo Pattern

Functional Programming

Tutorial Video

Text description (accessibility)

This video demonstrates the "Monte Carlo Pattern" functional Rust example. Difficulty level: Fundamental. Key concepts covered: Functional Programming. Some problems are too complex for exact algorithms but can be approximated by random sampling. Key difference from OCaml: | Aspect | Rust | OCaml |

Tutorial

The Problem

Some problems are too complex for exact algorithms but can be approximated by random sampling. Monte Carlo methods estimate quantities by sampling: estimate pi by generating random points in a unit square and counting those inside the unit circle; estimate integrals by averaging function values at random points; approximate Nash equilibria in game theory; simulate financial option pricing. Error decreases as O(1/sqrt(n)) — independent of dimensionality, which is the key advantage over deterministic quadrature for high-dimensional integration. Monte Carlo underpins financial risk models, particle physics simulation, and probabilistic algorithms.

🎯 Learning Outcomes

• Implement pi estimation: sample (x,y) uniformly in [-1,1]^2, count points where x^2+y^2 < 1

• Understand O(1/sqrt(n)) convergence: 4x more samples gives 2x more accuracy

• Implement Monte Carlo integration: estimate integral(f, a, b) by averaging f at random points

• Apply to high-dimensional integration where deterministic methods are exponentially expensive

• Recognize the importance of a good random number generator (LCG vs PCG vs Mersenne Twister)

Code Example

#![allow(clippy::all)]
//! Placeholder

(* Monte Carlo Methods in OCaml *)

(* Simple LCG PRNG returning float in [0, 1) *)
let state = ref 12345678

let rand_float () =
  state := (!state * 1664525 + 1013904223) land 0x7fffffff;
  float_of_int !state /. float_of_int 0x7fffffff

(* Pi estimation: sample n points in unit square, count those in unit circle *)
let estimate_pi (n : int) : float =
  let inside = ref 0 in
  for _ = 1 to n do
    let x = rand_float () and y = rand_float () in
    if x *. x +. y *. y <= 1.0 then incr inside
  done;
  4.0 *. float_of_int !inside /. float_of_int n

(* Monte Carlo integration: estimate ∫_a^b f(x) dx using n samples *)
let mc_integrate (f : float -> float) (a b : float) (n : int) : float =
  let sum = ref 0.0 in
  for _ = 1 to n do
    let x = a +. (b -. a) *. rand_float () in
    sum := !sum +. f x
  done;
  (b -. a) *. !sum /. float_of_int n

(* Monte Carlo: estimate E[X²] where X ~ Uniform[0, 1] *)
(* Exact answer: ∫₀¹ x² dx = 1/3 ≈ 0.333 *)

(* Monte Carlo: accept-reject sampling from a non-trivial distribution *)
(* Sample from f(x) ∝ sin(x) on [0, π] using uniform proposal *)
let sample_sin_distribution (n : int) : float list =
  let results = ref [] in
  let attempts = ref 0 in
  while List.length !results < n && !attempts < 100 * n do
    incr attempts;
    let x = rand_float () *. Float.pi in
    let u = rand_float () in
    if u < sin x then results := x :: !results
  done;
  !results

let () =
  Printf.printf "Pi estimation:\n";
  List.iter (fun n ->
    let pi_est = estimate_pi n in
    Printf.printf "  n=%7d: π ≈ %.6f  error = %.6f\n" n pi_est (abs_float (pi_est -. Float.pi))
  ) [100; 1000; 10000; 100000];

  Printf.printf "\nMC integration:\n";
  Printf.printf "  ∫₀¹ x² dx ≈ %.6f  (exact = 0.333333)\n"
    (mc_integrate (fun x -> x *. x) 0.0 1.0 100000);
  Printf.printf "  ∫₀^π sin(x) dx ≈ %.6f  (exact = 2.0)\n"
    (mc_integrate sin 0.0 Float.pi 100000)

Key Differences

Aspect	Rust	OCaml
RNG	`rand::RngCore` trait	`Random` module (LCG-based)
Testability	Seeded `SmallRng::seed_from_u64`	`Random.init seed`
Parallelism	`rayon::par_iter()`	`Domain.spawn` (OCaml 5.0)
Sample counting	`filter().count()`	Mutable counter in loop
Convergence	O(1/sqrt n)	Same
Higher dimensions	Same complexity advantage	Same

OCaml Approach

OCaml's Monte Carlo uses Random.float 2.0 -. 1.0 for uniform samples. The pi estimation is: let count = ref 0 in for _ = 1 to n do let x = ... and y = ... in if x*.x +. y*.y <= 1.0 then incr count done; 4.0 *. float !count /. float n. OCaml's Random.self_init () seeds from system entropy. The Seq module creates an infinite sample sequence for lazy Monte Carlo. OCaml's Float.sqrt and .*. operators handle the computation.

Full Source

#![allow(clippy::all)]
//! Placeholder

(* Monte Carlo Methods in OCaml *)

(* Simple LCG PRNG returning float in [0, 1) *)
let state = ref 12345678

let rand_float () =
  state := (!state * 1664525 + 1013904223) land 0x7fffffff;
  float_of_int !state /. float_of_int 0x7fffffff

(* Pi estimation: sample n points in unit square, count those in unit circle *)
let estimate_pi (n : int) : float =
  let inside = ref 0 in
  for _ = 1 to n do
    let x = rand_float () and y = rand_float () in
    if x *. x +. y *. y <= 1.0 then incr inside
  done;
  4.0 *. float_of_int !inside /. float_of_int n

(* Monte Carlo integration: estimate ∫_a^b f(x) dx using n samples *)
let mc_integrate (f : float -> float) (a b : float) (n : int) : float =
  let sum = ref 0.0 in
  for _ = 1 to n do
    let x = a +. (b -. a) *. rand_float () in
    sum := !sum +. f x
  done;
  (b -. a) *. !sum /. float_of_int n

(* Monte Carlo: estimate E[X²] where X ~ Uniform[0, 1] *)
(* Exact answer: ∫₀¹ x² dx = 1/3 ≈ 0.333 *)

(* Monte Carlo: accept-reject sampling from a non-trivial distribution *)
(* Sample from f(x) ∝ sin(x) on [0, π] using uniform proposal *)
let sample_sin_distribution (n : int) : float list =
  let results = ref [] in
  let attempts = ref 0 in
  while List.length !results < n && !attempts < 100 * n do
    incr attempts;
    let x = rand_float () *. Float.pi in
    let u = rand_float () in
    if u < sin x then results := x :: !results
  done;
  !results

let () =
  Printf.printf "Pi estimation:\n";
  List.iter (fun n ->
    let pi_est = estimate_pi n in
    Printf.printf "  n=%7d: π ≈ %.6f  error = %.6f\n" n pi_est (abs_float (pi_est -. Float.pi))
  ) [100; 1000; 10000; 100000];

  Printf.printf "\nMC integration:\n";
  Printf.printf "  ∫₀¹ x² dx ≈ %.6f  (exact = 0.333333)\n"
    (mc_integrate (fun x -> x *. x) 0.0 1.0 100000);
  Printf.printf "  ∫₀^π sin(x) dx ≈ %.6f  (exact = 2.0)\n"
    (mc_integrate sin 0.0 Float.pi 100000)

Deep Comparison

OCaml vs Rust: Monte Carlo Pattern

Overview

See the example.rs and example.ml files for detailed implementations.

Key Differences

Aspect	OCaml	Rust
Type system	Hindley-Milner	Ownership + traits
Memory	GC	Zero-cost abstractions
Mutability	Explicit ref	mut keyword
Error handling	Option/Result	Result<T, E>

See README.md for detailed comparison.

Exercises

Implement Monte Carlo integration for a 1D function and compare with the exact analytical result.

Extend to 10D integration of x1^2 + ... + x10^2 <= 1 to demonstrate the dimensionality advantage.

Use stratified sampling (divide the domain into strata, sample uniformly within each) and compare convergence with naive Monte Carlo.

Implement parallel Monte Carlo using Rayon with independent per-thread RNGs and measure linear speedup.

Estimate the variance of your pi estimator as a function of n and verify it matches the theoretical O(1/n) variance.

Open Source Repos

functional-rust

View the source for this example on GitHub — OCaml and Rust side by side in the repo.

Rust