From b1e29801c38258bee71e13cbbc3815af6f25842f Mon Sep 17 00:00:00 2001
From: Simon Byrne <simonbyrne@gmail.com>
Date: Thu, 21 Jan 2016 12:43:47 +0000
Subject: [PATCH] add covariance for multivariate mixtures

---
 src/mixtures/mixturemodel.jl | 28 ++++++++++++++++++++++++++++
 test/mixture.jl              |  1 +
 2 files changed, 29 insertions(+)

diff --git a/src/mixtures/mixturemodel.jl b/src/mixtures/mixturemodel.jl
index 1944e3288a..e46b261daa 100644
--- a/src/mixtures/mixturemodel.jl
+++ b/src/mixtures/mixturemodel.jl
@@ -122,6 +122,34 @@ function var(d::UnivariateMixture)
     return v
 end
 
+function cov(d::MultivariateMixture)
+    K = ncomponents(d)
+    p = probs(d)
+    m = zeros(length(d))
+    md = zeros(length(d))
+    V = zeros(length(d),length(d))
+
+    for i = 1:K
+        pi = p[i]
+        if pi > 0.0
+            c = component(d, i)
+            BLAS.axpy!(pi, mean(c), m)
+            BLAS.axpy!(pi, cov(c), V)
+        end
+    end
+    for i = 1:K
+        pi = p[i]
+        if pi > 0.0
+            c = component(d, i)
+            # todo: use more in-place operations
+            md = mean(c) - m
+            BLAS.axpy!(pi, md*md', V)
+        end
+    end
+    return V
+end
+
+
 
 #### show
 
diff --git a/test/mixture.jl b/test/mixture.jl
index 6a91935f0f..198cfbeed3 100644
--- a/test/mixture.jl
+++ b/test/mixture.jl
@@ -118,6 +118,7 @@ function test_mixture(g::MultivariateMixture, n::Int, ns::Int)
     @test isa(Xs, Matrix{Float64})
     @test size(Xs) == (length(g), ns)
     @test_approx_eq_eps vec(mean(Xs, 2)) mean(g) 0.01
+    @test_approx_eq_eps cov(Xs, vardim=2) cov(g) 0.01
 end