gibbs.fs

let gaussLik (a: float array) (b: float array) =
    Math.Exp -0.5* (dist a b)

let sampleDiscrete (p: float array) =
    let d = new Distributions.ArbitraryDistribution()
    d.SetDistributionParameters(0, p)
    d.NextInt32()


// Samples the z variable for a single theta
let sampleZ alpha (C: float array) (theta: float array) (mu: float array array) (oldz: int) =
    C.[oldz] <- C.[oldz] - 1.0
    let p = Array.init C.Length (fun i -> C.[i] + alpha)
    let mutable total = 0.0
    for i in 0..C.Length-1 do
        let f = p.[i] * (gaussLik theta mu.[i])
        p.[i] <- f
        total <- total + f
    for i in 0..C.Length-1 do
        p.[i] <- p.[i] / total
    let z  = sampleDiscrete p
    C.[z] <- C.[z] + 1.0
    z
    
    
// Samples all the mus at once
let sampleMu (thetas: float array array) (z: int array) =
    let K = z.Length
    let D = thetas.[0].Length
    let T = thetas.Length
    let mus = Array.init K (fun i -> Array.init D (fun j -> 0.0))
    for k in 0..z.Length-1 do
        let mu = mus.[k]
        let mutable sumtheta = Array.init D (fun i -> 0.0)
        let mutable n = 0.0
        for i in 0..T-1 do
            if z.[i] = k then
                n <- n +  1.0
                sumtheta <- arraySum sumtheta thetas.[i]
        for i in 0..D-1 do
            let N = new Distributions.NormalDistribution(sumtheta.[i]/(n+1.0), 1.0/(n+1.0))
            mu.[i] <- N.NextDouble()