simsam 0.1.0

Sample from custom discrete and continuous distributions (SciPy-like API)
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83

use crate::error::{BuildError, SampleError};
use crate::multivar::support::HyperRect;
use crate::multivar::traits::{HasSupportNd, PdfNd};
use rand::{Rng, RngExt};

#[derive(Debug, Clone, Copy)]
pub struct RejectionOptions {
    /// Maximum number of proposal draws per sample.
    pub max_trials: usize,
}

impl Default for RejectionOptions {
    fn default() -> Self {
        Self { max_trials: 100_000 }
    }
}

/// Generic N-D rejection sampler on a bounded hyper-rectangle.
///
/// Requires an upper bound `pdf_max` such that `pdf(x) <= pdf_max` on the whole support.
pub struct RejectionSamplerNd<P> {
    pdf: P,
    support: HyperRect,
    pdf_max: f64,
    opts: RejectionOptions,
}

impl<P> RejectionSamplerNd<P>
where
    P: PdfNd + HasSupportNd,
{
    pub fn new(pdf: P, pdf_max: f64, opts: RejectionOptions) -> Result<Self, BuildError> {
        if !pdf_max.is_finite() || pdf_max <= 0.0 {
            return Err(BuildError::InvalidPdfMax);
        }
        let support = pdf.support().clone();
        support.validate()?;
        Ok(Self {
            pdf,
            support,
            pdf_max,
            opts,
        })
    }

    pub fn dim(&self) -> usize {
        self.support.dim()
    }

    pub fn support(&self) -> &HyperRect {
        &self.support
    }

    pub fn pdf_max(&self) -> f64 {
        self.pdf_max
    }

    pub fn sample(&self) -> Result<Vec<f64>, SampleError> {
        self.sample_with_rng(&mut rand::rng())
    }

    pub fn sample_with_rng<R: Rng + ?Sized>(&self, rng: &mut R) -> Result<Vec<f64>, SampleError> {
        let dim = self.dim();
        let mut x = vec![0.0; dim];

        for _ in 0..self.opts.max_trials {
            for i in 0..dim {
                let u: f64 = rng.random();
                x[i] = self.support.lo[i] + u * (self.support.hi[i] - self.support.lo[i]);
            }
            let fx = self.pdf.pdf(&x);
            if !fx.is_finite() || fx < 0.0 {
                continue;
            }
            let u: f64 = rng.random();
            if u * self.pdf_max <= fx {
                return Ok(x);
            }
        }
        Err(SampleError::RejectionSamplingFailed)
    }
}