rejection-sampler 0.1.4

A statistical computing toolkit for validating proposal distributions and computing optimal rejection-sampling constants.

rejection-sampler

A small Python package for validating rejection sampling setups and computing the optimal rejection constant (M). The package supports both callable Python functions and symbolic SymPy expressions for target and proposal probability density functions (PDF).

In rejection sampling, optimal M is the smallest constant such that, for all x in the target support:

target_pdf(x) <= M * proposal_pdf(x)

Installation

pip install rejection-sampler

or

uv add rejection-sampler

Usage

Import the main function:

from rejection_sampler import find_optimal_M

The function supports two ways of PDF definition: callable input and symbolic input (with SymPy). In general, symbolic input has a much higher chance of success. When calling the function, 4 arguments must be provided: target_pdf, target_support, proposal_pdf, proposal_support.

Example 1: Callable input

def target_pdf(x):
    return 2 * x if 0 <= x <= 1 else 0.0

def proposal_pdf(x):
    return 1.0 if 0 <= x <= 1 else 0.0

M = find_optimal_M(
    target_pdf=target_pdf,
    target_support=(0.0, 1.0),
    proposal_pdf=proposal_pdf,
    proposal_support=(0.0, 1.0),
)

print(M)

Example 2: SymPy input

import sympy as sp

x = sp.Symbol("x", real=True)

target_pdf = 2 * x
proposal_pdf = sp.Integer(1)

M = find_optimal_M(
    target_pdf=target_pdf,
    target_support=(0, 1),
    proposal_pdf=proposal_pdf,
    proposal_support=(0, 1),
)

print(M)

Infinite support

For numerical inputs with infinite support, the argument bounds must be provided (see Note for more info on bounds):

import numpy as np
from rejection_sampler import find_optimal_M

def target_pdf(x):
    return np.exp(-0.5 * x * x) / np.sqrt(2 * np.pi)

def proposal_pdf(x):
    return 1.0 / (np.pi * (1 + x * x))

M = find_optimal_M(
    target_pdf=target_pdf,
    target_support=(-sp.oo, sp.oo),
    # or use (-np.inf, np.inf)
    # or use (-float("inf"), float("inf")) for infinite support
    proposal_pdf=proposal_pdf,
    proposal_support=(-np.inf, np.inf),
    bounds=(-10.0, 10.0),
)

print(M)

Parameters

• target_pdf: target probability density function, either callable or SymPy expression
• target_support: support of the target PDF
• proposal_pdf: proposal probability density function, either callable or SymPy expression
• proposal_support: support of the proposal PDF
• error: numerical tolerance for validation
• bounds: search interval for numerical optimization for pdfs with infinite support

Note

• When writing mathematical expressions (eg. exp, log, sqrt, inf), use SymPy or NumPy instead of the built-in math module.
• When checking complicated PDFs, SymPy input has a much higher rate of success. Always use SymPy if possible.
• For infinite-support callable inputs, bounds defines the finite interval used for numerical optimization. It should contain the maximum of target_pdf / proposal_pdf. The returned M is optimized only over this interval and does not prove global domination outside it, so selecting appropriate bounds may require domain knowledge or experimentation.
• Known issue: For certain PDFs, there might be overflowing/underflowing issues. For example, when sampling the standard Gumbel distribution for maxima with Cauchy distribution, the calculation overflows when bounds are too large, and an appropriate bounds is (-10, 10).

Source code

• GitHub: https://github.com/HankTaiwan869/rejection-sampler
• Issues are welcome.

License

This project is licensed under the MIT License.