modppl 0.3.1

a experimental library for probabilistic programming in Rust.
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
84
85
86
87
88
89
90
91
92
93
94
95
96

use nalgebra::{dmatrix, dvector};
use rand::rngs::ThreadRng;

use modppl::{ParticleSystem, Real};

mod hmm;

#[test]
fn test_hmm_forward_algorithm() {
    let prior = dvector![0.4, 0.6];
    let emission_dists = dmatrix![0.1, 0.9; 0.7, 0.3].transpose();
    let transition_dists = dmatrix![0.5, 0.5; 0.2, 0.8].transpose();
    let obs = vec![1, 0];

    let mut true_marginal_likelihood = 0.0;
    // z = [0, 0]
    true_marginal_likelihood += prior[0]
        * emission_dists[(obs[0] as usize, 0)]
        * transition_dists[(0, 0)]
        * emission_dists[(obs[1] as usize, 0)];
    // z = [0, 1]
    true_marginal_likelihood += prior[0]
        * emission_dists[(obs[0] as usize, 0)]
        * transition_dists[(1, 0)]
        * emission_dists[(obs[1] as usize, 1)];
    // z = [1, 0]
    true_marginal_likelihood += prior[1]
        * emission_dists[(obs[0] as usize, 1)]
        * transition_dists[(0, 1)]
        * emission_dists[(obs[1] as usize, 0)];
    // z = [1, 1]
    true_marginal_likelihood += prior[1]
        * emission_dists[(obs[0] as usize, 1)]
        * transition_dists[(1, 1)]
        * emission_dists[(obs[1] as usize, 1)];

    let empirical_marginal_likelihood =
        hmm::hmm_forward_alg(prior, emission_dists, transition_dists, &obs);
    approx::assert_abs_diff_eq!(
        empirical_marginal_likelihood,
        true_marginal_likelihood,
        epsilon = Real::EPSILON.sqrt()
    );
}

#[test]
fn test_particle_filter() -> std::io::Result<()> {
    let rng = ThreadRng::default();
    const NUM_PARTICLES: usize = 10000;

    let prior = dvector![0.2, 0.3, 0.5];
    let emission_matrix = dmatrix![
        0.1, 0.2, 0.7;
        0.2, 0.7, 0.1;
        0.7, 0.2, 0.1
    ]
    .transpose();
    let transition_matrix = dmatrix![
        0.4, 0.4, 0.2;
        0.2, 0.3, 0.5;
        0.9, 0.05, 0.05
    ]
    .transpose();
    let params = hmm::HMMParams::new(
        prior.clone(),
        emission_matrix.clone(),
        transition_matrix.clone(),
    );

    let model = hmm::HMM::new(params);

    let data = vec![0, 0, 1, 2];
    let expected = hmm::hmm_forward_alg(prior, emission_matrix, transition_matrix, &data).ln();

    let mut filter = ParticleSystem::new(model, NUM_PARTICLES, rng);

    let store = hmm::ParamStore {};
    let mut data_it = data.into_iter();
    filter.init_step(store, (vec![None], vec![data_it.next()]));
    println!("T = {}", 1);

    for (t, obs) in data_it.enumerate() {
        println!("T = {}", t + 2); // time is 1-indexed and init_step used "1 => observation" (1 + 1 = 2)
        let constraints = (vec![None], vec![Some(obs)]);
        filter = filter.step(constraints);
        let ess = filter.effective_sample_size();
        dbg!(ess);
        let log_weight = filter.resample();
        dbg!(log_weight);
    }
    let lml_estimate = filter.log_marginal_likelihood_estimate();

    approx::assert_abs_diff_eq!(lml_estimate, expected, epsilon = 0.03);

    Ok(())
}