friedrich 0.6.0

Gaussian Process Regression.
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

use crate::conversion::Input;
use nalgebra::{DMatrix, DVector};
use rand::{Rng, RngExt};
use rand_distr::StandardNormal;
use std::marker::PhantomData;

/// Multivariate Normal Distribution.
///
/// This class is meant to be produced by the `sample_at` method of the gaussian process
/// and can be used to sample the process at a given point / set of points :
///
/// ```rust
/// # use friedrich::gaussian_process::GaussianProcess;
/// // Trains a model.
/// let training_inputs = vec![vec![0.8], vec![1.2], vec![3.8], vec![4.2]];
/// let training_outputs = vec![3.0, 4.0, -2.0, -2.0];
/// let gp = GaussianProcess::default(training_inputs, training_outputs);
///
/// // Computes the distribution at some new coordinates.
/// let new_inputs = vec![vec![1.], vec![2.]];
/// let sampler = gp.sample_at(&new_inputs);
///
/// // Samples from the distribution.
/// let mut rng = rand::rng();
/// println!("samples a vector : {:?}", sampler.sample(&mut rng));
/// ```
///
/// Note that the output type is a function of the input of `sample_at`, the method can be used on a vector of vectors as well as a single row :
///
/// ```rust
/// # use friedrich::gaussian_process::GaussianProcess;
/// # // Trains a model.
/// # let training_inputs = vec![vec![0.8], vec![1.2], vec![3.8], vec![4.2]];
/// # let training_outputs = vec![3.0, 4.0, -2.0, -2.0];
/// # let gp = GaussianProcess::default(training_inputs, training_outputs);
/// // Computes the distribution at some new coordinate.
/// let new_input = vec![1.];
/// let sampler = gp.sample_at(&new_input);
///
/// // Samples from the distribution.
/// let mut rng = rand::rng();
/// println!("samples a value : {}", sampler.sample(&mut rng));
/// ```
pub struct MultivariateNormal<T: Input>
{
    mean: DVector<f64>,
    cholesky_covariance: DMatrix<f64>,
    input_type: PhantomData<T>
}

impl<T: Input> MultivariateNormal<T>
{
    /// Produces a new multivariate gaussian with the given parameters.
    pub fn new(mean: DVector<f64>, covariance: DMatrix<f64>) -> Self
    {
        let cholesky_covariance =
            covariance.cholesky().expect("MultivariateNormal: Cholesky decomposition failed!").unpack();
        MultivariateNormal { mean, cholesky_covariance, input_type: PhantomData }
    }

    /// Outputs the mean of the distribution.
    pub fn mean(&self) -> T::OutVector
    {
        T::from_dvector(&self.mean)
    }

    /// Takes a random number generator and uses it to sample from the distribution.
    pub fn sample<RNG: Rng>(&self, rng: &mut RNG) -> T::OutVector
    {
        let normal = DVector::from_fn(self.mean.nrows(), |_, _| rng.sample(StandardNormal));
        let sample = &self.mean + &self.cholesky_covariance * normal;
        T::from_dvector(&sample)
    }
}