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

extern crate nalgebra;
extern crate rand;

use self::nalgebra::{DMatrix, DVector};
use self::rand::Rng;
use dist::{InvWishart, MvGaussian};
use std::io;
use traits::*;

/// Common conjugate prior on the μ and Σ parameters in the Multivariate
/// Gaussian, Ν(μ, Σ)
///
/// Ν(μ, Σ) ~ NIW(μ<sub>0</sub>, κ<sub>0</sub>, ν, Ψ) implies
/// μ ~ N(μ<sub>0</sub>, Σ/k<sub>0</sub>) and
/// Σ ~ W<sup>-1</sup>(Ψ, ν)
///
/// # Example
///
/// Draw a Multivariate Gaussian from GIW
///
/// ```
/// # extern crate rv;
/// extern crate rand;
/// extern crate nalgebra;
///
/// use nalgebra::{DMatrix, DVector};
/// use rv::prelude::*;
///
/// let mu = DVector::zeros(3);
/// let k = 1.0;
/// let df = 3;
/// let scale = DMatrix::identity(3, 3);
///
/// let niw = NormalInvWishart::new(mu, k, df, scale).unwrap();
///
/// let mut rng = rand::thread_rng();
///
/// let mvg: MvGaussian = niw.draw(&mut rng);
/// ```
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde_support", derive(Serialize, Deserialize))]
pub struct NormalInvWishart {
    /// The mean of μ, μ<sub>0</sub>
    pub mu: DVector<f64>,
    /// A scale factor on Σ, κ<sub>0</sub>
    pub k: f64,
    /// The degrees of freedom, ν > |μ| - 1
    pub df: usize,
    /// The positive-definite scale matrix, Ψ
    pub scale: DMatrix<f64>,
}

impl NormalInvWishart {
    /// Create a new `NormalInvWishart` distribution
    pub fn new(
        mu: DVector<f64>,
        k: f64,
        df: usize,
        scale: DMatrix<f64>,
    ) -> io::Result<Self> {
        let dims = mu.len();
        let err = if k <= 0.0 {
            Some("k must be > 0.0")
        } else if df < dims {
            Some("df must be >= the number of dimensions")
        } else if !scale.is_square() {
            Some("scale must be square")
        } else if dims != scale.nrows() {
            Some("dimensions in mu and scale must match")
        } else if scale.clone().cholesky().is_none() {
            Some("scale is not positive definite")
        } else {
            None
        };

        match err {
            Some(msg) => Err(io::Error::new(io::ErrorKind::InvalidInput, msg)),
            None => Ok(NormalInvWishart { mu, k, df, scale }),
        }
    }
}

// TODO: We might be able to make things faster by storing the InvWishart
// because each time we create it, it clones and validates the parameters.
impl Rv<MvGaussian> for NormalInvWishart {
    fn ln_f(&self, x: &MvGaussian) -> f64 {
        let m = self.mu.clone();
        let sigma = x.cov.clone() / self.k;
        let mvg = MvGaussian::new(m, sigma).unwrap();
        let iw = InvWishart::new(self.scale.clone(), self.df).unwrap();
        mvg.ln_f(&x.mu) + iw.ln_f(&x.cov)
    }

    fn draw<R: Rng>(&self, mut rng: &mut R) -> MvGaussian {
        let iw = InvWishart::new(self.scale.clone(), self.df).unwrap();
        let sigma = iw.draw(&mut rng);

        let mvg =
            MvGaussian::new(self.mu.clone(), sigma.clone() / self.k).unwrap();
        let mu = mvg.draw(&mut rng);

        MvGaussian::new(mu, sigma).unwrap()
    }
}

impl Support<MvGaussian> for NormalInvWishart {
    fn supports(&self, x: &MvGaussian) -> bool {
        let p = self.mu.len();
        x.mu.len() == p && x.cov.clone().cholesky().is_some()
    }
}

impl ContinuousDistr<MvGaussian> for NormalInvWishart {}