Yesterday I described my experiences implementing multivariate Gaussian distributions in Julia. I’ve since repeated the exercise in Clojure, for any time you might need to do such a thing on the JVM. I’m using core.matrix for array operations.

μ, again is just (mean x).

To calculate the covariance matrix:

(defn covariance
  "Calculate covariance matrix"
  [x]
  (let [mu (mean x)]
    (div
     (inner-product (sub x mu) (transpose (sub x mu)))
     (row-count x))))

This is just Σ = (X - μ)’(X - μ) / n (where n is the number of examples).

And now that we have μ and Σ, we can use this function for probability:

(defn probability
  "Calculate probability given mean, covariance matrix, and X"
  [mu sigma x]
  (let [x-minus-mu (sub x mu)
        n (row-count mu)]
    (*  (exp (* -0.5 (esum
                      (emul (transpose x-minus-mu)
                            (inverse sigma)
                            x-minus-mu))))
        (Math/pow (* 2 Math/PI) (* n -0.5))
        (Math/pow (det sigma) -0.5))))

I’ve wrapped the code for this up in a library, available on Github.