metromc 0.2.0

Markov chain Monte Carlo sampling using the Independence Metropolis-Hastings algorithm.
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

use std::ops::RangeInclusive;
use tinyrand::Rand;
use crate::Pdf;

pub struct Config<R: Rand, P: Pdf> {
    pub rand: R,
    pub dist: P,
    pub range: RangeInclusive<f64>,
}

pub struct Sampler<R: Rand, P: Pdf> {
    config: Config<R, P>,
    current: f64,
    prob_current: f64,
    span: f64,
}
impl<R: Rand, P: Pdf> Sampler<R, P> {
    pub fn new(config: Config<R, P>) -> Self {
        let span  = config.range.end() - config.range.start();
        let current = (config.range.start() + config.range.end()) / 2.0;
        let prob_current = config.dist.prob(current);
        Self {
            config,
            current,
            prob_current,
            span,
        }
    }

    #[inline(always)]
    pub fn next_sample(&mut self) -> f64 {
        let next = random_f64(&mut self.config.rand) * self.span + self.config.range.start();
        let prob_next = self.config.dist.prob(next);
        //println!("current: {:.3}, prob_current: {:.3}, next: {next:.3}, prob_next: {prob_next:.3}", self.current, self.prob_current);
        if prob_next >= self.prob_current {
            self.current = next;
            self.prob_current = prob_next;
        } else {
            let threshold = prob_next / self.prob_current;
            let rand = random_f64(&mut self.config.rand);
            if rand < threshold {
                self.current = next;
                self.prob_current = prob_next;
            }
        }
        self.current
    }
}

impl<R: Rand, P: Pdf> Iterator for Sampler<R, P> {
    type Item = f64;

    #[inline(always)]
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.next_sample())
    }
}

#[inline(always)]
fn random_f64(rand: &mut impl Rand) -> f64 {
    rand.next_u64() as f64 / u64::MAX as f64
}