gen-rs 0.2.3

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

use rand::rngs::ThreadRng;
use super::{Distribution,normal};
use std::f64::consts::PI;
use nalgebra::{DVector,DMatrix};


/// Multivariate Gaussian distribution type
pub struct MvNormal { }

/// Instantiation of the Multivariate Gaussian distribution
pub const mvnormal: MvNormal = MvNormal { };

impl Distribution<DVector<f64>,(DVector<f64>,DMatrix<f64>)> for MvNormal {
    fn logpdf(&self, x: &DVector<f64>, params: (DVector<f64>,DMatrix<f64>)) -> f64 {
        let (mu, cov) = params;
        let k = mu.len() as f64;
        let cov_det = cov.determinant();
        let cov_inv = cov.clone().try_inverse().unwrap();
        let centered_x = x - mu;
        let mahalanobis_squared = (centered_x.transpose() * cov_inv * centered_x).trace();
        -(k*(2.*PI).ln() + cov_det.ln() + mahalanobis_squared)/2.
    }

    fn random(&self, rng: &mut ThreadRng, params: (DVector<f64>,DMatrix<f64>)) -> DVector<f64> {
        let (mu, cov) = params;
        let transform: DMatrix<f64>;
        match cov.clone().cholesky() {
            Some(c) => {
                transform = c.l();
            },
            None => {
                let decomp = cov.clone().symmetric_eigen();
                transform = decomp.eigenvectors * DMatrix::from_diagonal(&decomp.eigenvalues.map(|v| v.sqrt()));
            }
        }
        let samples = transform * &mu.map(|_| normal.random(rng, (0.,1.))) + mu;
        samples
    }
}