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

use crate::matrix::Matrix;
use crate::vector::Vector;
use crate::{Kernel, Mean};
use rand::distributions::Distribution;
use rand::Rng;
use rand_distr;

pub use self::mvn::MultivariateNormal;
pub use self::posterior::GaussianProcessPosterior;
pub use self::prior::GaussianProcessPrior;

mod mvn;
mod posterior;
mod prior;

#[derive(Debug)]
pub struct Posterior {
    inner: rand_distr::Normal<f64>,
}
impl Posterior {
    pub fn new<X, F, M, K>(x: &X, xs: &[X], f: F, mean: M, kernel: K) -> Option<Self>
    where
        X: Clone,
        F: Fn(&X) -> f64,
        M: Mean<X>,
        K: Kernel<X>,
    {
        let cov0 = Matrix::covariance(&[x.clone()], xs, &kernel);
        let cov1_i = Matrix::covariance(xs, xs, &kernel).inverse()?;

        let m0 = mean.mean(x);
        let means0 = xs
            .iter()
            .map(|x| f(x) - mean.mean(x) + m0)
            .collect::<Vector>();
        let means =
            cov0.clone().into_inner() * (cov1_i.clone().into_inner() * means0.clone().into_inner());
        let mean = means.get((0, 0))?;

        let cov2 = kernel.kernel(x, x);
        let cov3 = Matrix::covariance(xs, &[x.clone()], &kernel);
        let cov4 = cov0.into_inner() * (cov1_i.into_inner() * cov3.into_inner());
        let variance = cov2 - cov4.get((0, 0))?;

        Some(Self {
            inner: rand_distr::Normal::new(*mean, variance.sqrt()).ok()?,
        })
    }
}
impl Distribution<f64> for Posterior {
    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
        self.inner.sample(rng)
    }
}