Step Optimised Rejection Sampling

Fast and efficient sampling from general univariate probability density functions. Implements a rejection sampling approach designed to take advantage of modern CPU caches and minimise evaluation of the target density for most samples. Many standard densities are internally implemented in 'C' for high performance, with general user defined densities also supported. A paper describing the methodology will be released soon.

📦 {stors} R package

StORS (Step Optimised Rejection Sampling) is an optimised rejection sampling method designed for efficient sampling from univariate probability density functions. It constructs a proposal distribution using computationally efficient steps that adapt around the target density.

Built-in distributions can have any type of unbounded tails, as Inverse Transform sampling is used when necessary. However, for user-provided distributions, if unbounded, the tails must be log-concave since Adaptive Rejection Sampling (ARS) is applied in this case.

Installation

The stable binary release can be install from CRAN:

install.packages("stors")

The latest development binary can be installed from R-universe:

install.packages('stors', repos = c('https://ahmad-alqabandi.r-universe.dev', 'https://cloud.r-project.org'))

Or from GitHub:

# install.packages("pak")
pak::pak("ahmad-alqabandi/stors")

Example

Sampling from a Built-in Normal Distribution

The first example demonstrates sampling from a built-in Normal distribution using a pre-optimised proposal that is constructed at package installation time for efficiency.

library("stors")

# Generate 1000 samples from the standard normal distribution
x <- srnorm(1000)

# Visualize the density
hist(x)

Re-optimizing the Built-in Normal Proposal

If needed, the user can re-optimize the proposal for the normal distribution. This can be useful for tasks such as sampling from a truncated normal distribution.

# Re-optimize the proposal for a truncated normal distribution (e.g., between -2 and 2)
srnorm_optimize(xl = -2, xr = 2)

# Generate new samples using the optimised proposal
x <- srnorm(1000)

# Visualize the truncated normal density
hist(x)

Building a Proposal for a Custom Distribution

Users can build a custom proposal by passing a probability density function (PDF) to build_proposal(). The following example demonstrates building a proposal for a Pareto distribution.

# Define the Pareto PDF
dpareto <- function(x, alpha) {
  alpha / x^(alpha + 1)
}

# Build a proposal distribution for Pareto
g <- build_proposal(f = dpareto,
                    modes = 1,
                    lower = 1,
                    upper = Inf,
                    steps = 256,
                    alpha = 2)

# Print and plot the proposal
print(g)
plot(g)

# Create a sampler from the proposal
rpareto_sampler <- build_sampler(g)

# Generate 1000 samples from the Pareto distribution
x <- rpareto_sampler(1000)

# Visualize the sampled distribution
hist(x)

For more details, please see the Getting Started vignette