MonteCarloZoo.jl

MonteCarloZoo

A broad collection of Monte Carlo algorithms implemented in Julia.

BoxMullerTransformSampler()

A Box-Muller transform sampler.

The Box-Muller transform is a technique for generating independent, standard, normally distributed samples. It involves sampling from a polar coordinate system in which the angle of the sample is uniformly distrubition and the radius is distributed as $\text{exp}(\tfrac{1}{2})$ (equivalent to $\chi^2_1$). It can then be shown that the cartesian coordinates of these samples have independent standard normal distributions.

The full algorithm is given here:

• Sample $U_1, U_2 \stackrel{\text{i.i.d.}}{\sim} \text{Unif}(0, 1)$
• Compute $Z_1 = \sqrt{-2 \log U_1}\cos(2\pi U_2)$, $Z_2 = \sqrt{-2 \log U_1}\sin(2\pi U_2)$

Arguments

None

Examples

# Generating 10 standard normal samples
s = BoxMullerTransformSampler()
samples = sample(s, 10)

InverseTransformSampler(F_inv)

An inverse transform sampler.

Inverse transform sampling is based on the result that when $U \sim \text{Unif}(0, 1)$, $X = F^{-1}(U)$ will be distributed according to the CDF $F$. For non-decreasing CDFs we can replace $F^{-1}$ with $F$'s generalised inverse $F^{-}(u) = \inf{x : F(x) > u}$ and use the same method.

Arguments

• F_inv::Function: the inverse CDF of the target distrubition.

Notes

• Because of the symmetry of a uniform random sample, one can replace $1-u$ in the inverse CDF with $u$.

Examples

# Generating 10 standard exponential random variables
F_inv(u) = -log(u)  # using symmetry of U (see notes)
s = InverseTransformSampler(F_inv)
samples = sample(s, 10)

LinearCongruentialGenerator(a, c, m, seed)

A linear congruential generator.

A linear congruential generator (LCG) is a primitive psueudo-random number generator (PRNG), producing samples that approximate a $\text{Uniform}(0, 1)$ distribution. Starting from a seed $X_0$, a sequence is generated according to the recurence relation,

\[X_{n+1} = (aX_n + c) \mod m\]

which are then turned into samples $Z_n = \frac{X_n}{m} \in [0, 1)$.

Although the generator is favoured for esoteric purposes due to its simplicity and speed, in practice it is highly flawed and can result in samples that have minimal random structure when viewed through an appropriate projection (as in Entacher). The generator is also highly sensitive to the choice of $a$, $c$, and $m$.

Arguments

• a::Integer: the "multiplier"; must be between 0 and m, exclusive.
• c::Integer: the "increment"; must be between 0 and m-1, inclusive.
• m::Integer: the "modulus"; must be non-negative proposal distribution.
• seed::Integer: the seed value; must be between 0 and m-1 inclusive

Examples

# Generating 10 uniform random variables
# Parameters taken from "Numerical Recipes"
a = 1664525
c = 1013904223
m = 2^32
seed = 1729
s = LinearCongruentialGenerator(a, c, m, seed)
samples = sample(s, 10)

# Batch sampling
s = LinearCongruentialGenerator(a, c, m, seed)
samples1 = sample(s, 10)
samples2 = sample(s, 10)
# samples1 ≠ samples2 in general

NormalNormalTransformer(normal_sampler, new_params, old_params=(0.0, 1.0))

Generate normal random variables from another normal sampler.

Transform samples from a normal distribution sampler into samples from a normal distribution with different parameters. This uses the result that if $Z$ is a standard normal random variable ($Z \sim \text{Normal}(0, 1)$) then $X = \sigma Z + \mu$ has distribution $\text{Normal}(\mu, \sigma^2)$.

Arguments

• normal_sampler::sampler: a sampler producing normally distributed samples.
• new_params::Tuple{Real, Real}: a tuple (mu, sigma2) giving the mean and variance of the target distribution
• old_params::Tuple{Real, Real}=(0.0, 1.0): a tuple (mu, sigma2) giving the mean and variance of the normal sampler. Defaults to a standard normal distribution.

Examples

# Transforming ten standard normal random variables into N(1, 2)
s = BoxMullerTransformSampler()
t = NormalNormalTransformer(s, (1, 2))
transformed_samples = sample(t, 10)

RejectionSampler(f, g, proposal_sampler, M)

A rejection sampler.

Rejection sampling involves sampling from a target distribution $f$ using a proposal distribution $g$ which we are able to directly sample from. In order for the resulting samples to be distributed according to $f$ we require that $\frac{f(x)}{g(x)}$ is bounded from above by some constant $M$ (M) for all $x$ in the support of $f$.

Under these conditions, the rejection sampling algorithm given below will generate independent samples from $f$:

• Sample $X \sim g$
• Accept $X$ with probability $\frac{f(X)}{Mg(X)}$

Arguments

• f::Function: density function for the target distribution.
• g::Function: density function for the proposal distribution.
• proposal_sampler::Function: a function that generates samples from the proposal distribution.
• M::Real: a scaling constant that bounds the ratio of the target and proposal density.

Notes

• No checks are made to ensure that the scale constant is valid or that the proposal sampler and density match. If these conditions are not met, the resulting samples will not be distributed according to the target density.
• The expected acceptance rate for the sampler is the reciprocal of M.

Examples

# Generating 10 standard normal samples using a Cauchy proposal
f(x) = exp(-x^2 / 2) / sqrt(2π)
g(x) = 1 / (π * (1 + x^2))
proposal_sampler = InverseTransformSampler(u -> tan(π * (u - 0.5)))
M = sqrt(2π / ℯ)
s = RejectionSampler(f, g, proposal_sampler, M)
samples = sample(s, 10)

sample(s::BoxMullerTransformSampler, N)

Generate samples using a Box-Muller transform sampler.

Arguments

• s::BoxMullerTransformSampler: a Box-Muller transform sampler.
• N::Integer: the number of samples to generate.

Returns

• samples::Array{Float64}: a 2-D array in which each column is a sample.

sample(s::InverseTransformSampler, N)

Generate samples using an inverse transform sampler.

Arguments

• s::InverseTransformSampler: an inverse transform sampler.
• N::Integer: the number of samples to generate.

Returns

• samples::Array{Float64}: a 2-D array in which each column is a sample.

sample(s::LinearCongruentialGenerator, N)

Generate samples using a linear congruential generator.

Arguments

• s::LinearCongruentialGenerator: a linear congruential generator.
• N::Integer: the number of samples to generate.

Returns

• samples::Array{Float64}: a 2-D array in which each column is a sample.

sample(s::NormalNormalTransformer, N)

Transform normal samples into normal samples with new parameters.

Arguments

• s::NormalNormalTransformer: a normal-normal transformer.
• N::Integer: the number of samples to generate.

Returns

• samples::Array{Float64}: a 2-D array in which each column is a sample.

sample(s::RejectionSampler, N)

Generate samples using a rejection sampler.

Arguments

• s::RejectionSampler: a rejection sampler.
• N::Integer: the number of samples to generate.

Returns

• samples::Array{Float64}: a 2-D array in which each column is a sample.