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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

use nalgebra::{dmatrix, dvector, DVector};
use rand::rngs::ThreadRng;
use std::any::Any;
use std::fs::{create_dir_all, write};
use std::sync::Arc;

use modppl::prelude::*;

mod pointed_model;
use pointed_model::types_2d::Bounds;
use pointed_model::{DriftProposal, PointedModel};

mod dyngenfns;
use dyngenfns::{
    add_or_remove_param_proposal, hierarchical_drift_proposal, hierarchical_model, read_coeffs,
};
use dyngenfns::{pointed_2d_drift_proposal, pointed_2d_model};

#[test]
fn test_metropolis_hastings_handcoded() -> std::io::Result<()> {
    create_dir_all("../data")?;

    const NUM_ITERS: u32 = 25000;

    let model = PointedModel {
        obs_cov: dmatrix![1., -3./5.; -3./5., 2.],
    };
    let proposal = DriftProposal {
        drift_cov: dmatrix![0.25, 0.; 0., 0.25],
    };
    let bounds = Bounds {
        xmin: -5.,
        xmax: 5.,
        ymin: -5.,
        ymax: 5.,
    };
    let obs = dvector![0., 0.];

    let constraints = (None, Some(obs));

    let (mut trace, _) = model.generate(bounds, constraints);
    for iter in 0..NUM_ITERS {
        dbg!(iter);
        let (new_trace, accepted) = modppl::mh(&model, trace, &proposal, ());
        dbg!(accepted);
        trace = new_trace;
        let data = trace.data.0.clone().unwrap();
        let json = format!("[{},{}]", data[0], data[1]);
        write(format!("../data/mh_trace_{}.json", iter), json)?;
    }

    Ok(())
}

#[test]
pub fn test_metropolis_hastings_dyngenfn() -> std::io::Result<()> {
    create_dir_all("../data")?;

    const NUM_ITERS: u32 = 25000;

    let bounds = Bounds {
        xmin: -5.,
        xmax: 5.,
        ymin: -5.,
        ymax: 5.,
    };
    let obs = dvector![0., 0.];

    let mut observations = Trie::new();
    observations.observe("obs", Arc::new(obs) as Arc<dyn Any + Send + Sync>);

    let mut trace = pointed_2d_model
        .generate((bounds, dmatrix![1., -3./5.; -3./5., 2.]), observations)
        .0;
    for iter in 0..NUM_ITERS {
        dbg!(iter);
        let (new_trace, accepted) = mh(
            &pointed_2d_model,
            trace,
            &pointed_2d_drift_proposal,
            dmatrix![0.25, 0.; 0., 0.25],
        );
        dbg!(accepted);
        trace = new_trace;
        let data = trace.data.read::<DVector<Real>>("latent");
        let json = format!("[{},{}]", data[0], data[1]);
        write(format!("../data/mh_trace_{}.json", iter), json)?;
    }

    Ok(())
}

#[test]
pub fn test_metropolis_hastings_hierarchical() -> std::io::Result<()> {
    create_dir_all("../data")?;

    let mut rng = ThreadRng::default();

    let xs = vec![-5., -4., -3., -2., -1., 0., 1., 2., 3., 4., 5.];

    let mut observations = DynTrie::new();
    let (a, b, c) = (0.3, 0.4, 0.5);
    let ys = xs
        .iter()
        .map(|x| a + b * x + c * x * x + normal.random(&mut rng, (0., 0.1)))
        .collect::<Vec<Real>>();
    write(
        "../data/hierarchical_data.json",
        format!("[{:?}, {:?}]", xs, ys),
    )?;
    ys.into_iter().enumerate().for_each(|(i, y)| {
        observations.observe(&format!("(y, {})", i), Arc::new(y));
    });

    let mut trace = hierarchical_model.generate(xs, observations).0;
    let mut all_coeffs = vec![];
    for _ in 0..100 {
        let (new_trace, _) = mh(
            &hierarchical_model,
            trace,
            &add_or_remove_param_proposal,
            (),
        );
        trace = new_trace;
        all_coeffs.push(read_coeffs(&trace));
        for _ in 0..3 {
            let (new_trace, _) = mh(
                &hierarchical_model,
                trace,
                &hierarchical_drift_proposal,
                0.1,
            );
            trace = new_trace;
            all_coeffs.push(read_coeffs(&trace));
        }
        for _ in 0..10 {
            let (new_trace, _) = mh(
                &hierarchical_model,
                trace,
                &hierarchical_drift_proposal,
                0.01,
            );
            trace = new_trace;
            all_coeffs.push(read_coeffs(&trace));
        }
        write(
            "../data/hierarchical_model.json",
            format!("{:?}", all_coeffs),
        )?;
    }
    Ok(())
}