(defn sample-mvn
" Returns a sample of the given size from a Multivariate Normal
  distribution. This is equivalent to R's mvtnorm::rmvnorm function.

  Arguments:
    size -- the size of the sample to return

  Options:
    :mean (default (repeat (ncol sigma) 0))
    :sigma (default (identity-matrix (count mean)))


  Examples:

    (use '(incanter core stats charts))
    (def mvn-samp (sample-mvn 1000 :mean [7 5] :sigma (matrix [[2 1.5] [1.5 3]])))
    (covariance mvn-samp)
    (def means (map mean (trans mvn-samp)))

    ;; plot scatter-plot of points
    (def mvn-plot (scatter-plot (sel mvn-samp :cols 0) (sel mvn-samp :cols 1)))
    (view mvn-plot)
    ;; add centroid to plot
    (add-points mvn-plot [(first means)] [(second means)])

    ;; add regression line to scatter plot
    (def x (sel mvn-samp :cols 0))
    (def y (sel mvn-samp :cols 1))
    (def lm (linear-model y x))
    (add-lines mvn-plot x (:fitted lm))


  References:
    http://en.wikipedia.org/wiki/Multivariate_normal

"
([^Integer size & options]
   (let [opts (when options (apply assoc {} options))
         mean (or (:mean opts)
                  (if (:sigma opts)
                    (repeat (ncol (:sigma opts)) 0)
                    [0]))
         sigma (or (:sigma opts)
                   (identity-matrix (count mean)))
         p (count mean)
         chol (decomp-cholesky sigma)
         norm-samp (mmult (matrix (sample-normal (* size p)) p) chol)
        ]
     (if (> (nrow norm-samp) 1)
       (matrix (map #(plus % (trans mean)) norm-samp))
       (matrix (plus norm-samp (trans mean)))))))