pub fn prs<F>(f: F, n: usize, (a, b): (f64, f64), m: usize, eps: f64) -> anyhow::Result<Vec<f64>>
where
    F: Fn(f64) -> f64 + Copy,
{
    let mut rng = rand::rng();

    let mut result = vec![0f64; n];

    let w = WeightedUniform::from_max_pool_1d(f, (a, b), m, eps)?;

    let mut initial_x = w.sample(n);
    let mut left_num = n;

    while left_num > 0 {
        for &x in initial_x.iter() {
            let weight = w.weight_at(x);
            if weight <= 0f64 {
                continue;
            } else {
                let y = rng.random_range(0f64..=weight);

                if y <= f(x) {
                    result[n - left_num] = x;
                    left_num -= 1;
                    if left_num == 0 {
                        return Ok(result);
                    }
                }
            }
        }
        initial_x = w.sample(left_num);
    }
    panic!("Error: failed to generate {} samples", n);
}
Piecewise Rejection Sampling

# Arguments
* `f` - Function to sample (unnormalized function is allowed)
* `n` - Number of samples
* `(a, b)` - Range of sampling
* `m` - Number of pieces
* `eps` - Epsilon for max pooling

# Examples
```
use peroxide::fuga::*;

fn main() -> Result<(), Box<dyn Error>> {
let f = |x: f64| {
if (0f64..=2f64).contains(&x) {
-(x - 1f64).powi(2) + 1f64
} else {
0f64
}
};

let samples = prs(f, 1000, (-1f64, 3f64), 200, 1e-4)?;
samples.mean().print(); // near 1

Ok(())
}
src/statistics/rand.rs