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

use crate::MultivariateDistribution;
use opensrdk_linear_algebra::Matrix;
use rand::prelude::*;
use rand_distr::StandardNormal;

pub struct MultivariateNormal {
    mean: Vec<f64>,
    l_cov: Matrix,
}

impl MultivariateNormal {
    /// # Multivariate normal
    /// `L` is needed as second argument under decomposition `Sigma = L * L^T`
    pub fn new(mean: Vec<f64>, l_cov: Matrix) -> Self {
        Self { mean, l_cov }
    }

    pub fn from(mean: Vec<f64>, cov: Matrix) -> Result<Self, String> {
        if mean.len() != cov.rows() {
            return Err("dimension mismatch".to_owned());
        }

        let decomposed = {
            let (u, sigma, _) = cov.gesvd()?;

            u * sigma.dipowf(0.5)
        };

        Ok(Self::new(mean, decomposed))
    }

    pub fn mean(&self) -> &[f64] {
        &self.mean
    }

    pub fn l_cov(&self) -> &Matrix {
        &self.l_cov
    }

    pub fn cov(&self) -> Matrix {
        &self.l_cov * self.l_cov.t()
    }
}

impl MultivariateDistribution for MultivariateNormal {
    fn sample(&self, thread_rng: &mut ThreadRng) -> Result<Vec<f64>, String> {
        let z = (0..self.l_cov.rows())
            .into_iter()
            .map(|_| thread_rng.sample(StandardNormal))
            .collect::<Vec<_>>();

        let y = Matrix::col(self.mean.clone()).gemm(&self.l_cov, &Matrix::col(z), 1.0, 1.0)?;

        Ok(y.elems())
    }
}

#[cfg(test)]
mod tests {
    #[test]
    fn it_works() {
        assert_eq!(2 + 2, 4);
    }
}