Metadata-Version: 2.4
Name: rejection-sampler
Version: 0.1.3
Summary: A statistical computing toolkit for validating proposal distributions and computing optimal rejection-sampling constants.
Keywords: rejection-sampling,sampling,computational-statistics,statistical-computing,scipy,sympy,probability,pdf,numerical-optimization
Author: Pin Hao Wang
Author-email: Pin Hao Wang <hank86030@gmail.com>
License-Expression: MIT
License-File: LICENSE
Requires-Dist: numpy>=2.4.6
Requires-Dist: scipy>=1.17.1
Requires-Dist: sympy>=1.14.0
Requires-Python: >=3.13
Project-URL: Issues, https://github.com/HankTaiwan869/rejection-sampler/issues
Project-URL: Repository, https://github.com/HankTaiwan869/rejection-sampler
Description-Content-Type: text/markdown

# 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
```bash
pip install rejection-sampler
```
or
```bash
uv add rejection-sampler
```
# Usage
Import the main function:
```python
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
```python
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
```python
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`):
```python
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 inputs, choose **appropriately large** `bounds`. `bounds` should contain where the maximum of target_pdf / proposal_pdf is. Choose a larger `bounds` if possible. For certain cases, the appropriate `bounds` might require some trial and errors to figure out.
- **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.