Julia package for fast flexible adaptive rejection sampling
Switch branches/tags
Nothing to show
Clone or download
mauriciogtec
Latest commit a314e5c Feb 20, 2018

README.md

Build Status AdaptiveRejectionSampling AdaptiveRejectionSampling

Coverage Status codecov

MIT Licence

AdaptiveRejectionSampling

using AdaptiveRejectionSampling
using Plots

Sampling from a shifted normal distribution

# Define function to be sampled
μ, σ = 1.0, 2.0
f(x) = exp(-0.5(x - μ)^2 / σ^2) / sqrt(2pi * σ^2) 
support = (-Inf, Inf)

# Build the sampler and simulate 10,000 samples
sampler = RejectionSampler(f, support, max_segments = 5)
@time sim = run_sampler!(sampler, 10000);
  0.010434 seconds (192.15 k allocations: 3.173 MiB)

Let's verify the result

x = -linspace(-10.0, 10.0, 100)
envelop = eval_envelop.(sampler.envelop, x)
target = f.(x)

histogram(sim, normalize = true, label = "Histogram")
plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])

Let's try a Gamma

α, β = 5.0, 2.0
f(x) = β^α * x^-1) * exp(-β*x) / gamma(α)
support = (0.0, Inf)

sampler = RejectionSampler(f, support)
@time sim = run_sampler!(sampler, 10000) 

x = linspace(0.0, 5.0, 100)
envelop = eval_envelop.(sampler.envelop, x)
target = f.(x)

histogram(sim, normalize = true, label = "Histogram")
plot!(x, [target envelop], width = 2, label = ["Gamma(α, β)" "Envelop"])
  0.007299 seconds (182.00 k allocations: 3.027 MiB)

Truncated distributions and unknown normalisation constant

We don't to provide an exact density--it will sample up to proportionality--and we can do truncated distributions

α, β = 5.0, 2.0
f(x) = β^α * x^-1) * exp(-β*x) / gamma(α)
support = (1.0, 3.5)

sampler = RejectionSampler(f, support)
@time sim = run_sampler!(sampler, 10000) 

x = linspace(0.01, 8.0, 100)
envelop = eval_envelop.(sampler.envelop, x)
target = f.(x)

histogram(sim, normalize = true, label = "histogram")
plot!(x, [target envelop], width = 2, label = ["target density" "envelop"])
  0.007766 seconds (181.82 k allocations: 3.024 MiB)