(defn sample-wishart
"
  Returns a p-by-p symmetric distribution drawn from a Wishart distribution

  Options:
    :p (default 2) -- number of dimensions of resulting matrix
    :df (default p) -- degree of freedoms (aka n), df <= p
    :scale (default (identity-matrix p)) -- positive definite matrix (aka V)

  Examples:
    (use 'incanter.stats)
    (sample-wishart :df 10  :p 4)

    ;; calculate the mean of 1000 wishart matrices, should equal (mult df scale)
    (div (reduce plus (for [_ (range 1000)] (sample-wishart :p 4))) 1000)


  References:
    http://en.wikipedia.org/wiki/Wishart_distribution#

"
  ([& options]
    (let [opts (when options (apply assoc {} options))
          scale (or (:scale opts) (when (:p opts) (identity-matrix (:p opts))))
          p (count scale)
          df (or (:df opts) p)
          diagonal (for [i (range 1 (inc p))]
                     (pow (sample-chisq 1 :df (inc (- df i))) 1/2))
          mat (diag diagonal)
          indices (for [i (range p) j (range p) :when (< j i)] [i j])
          _ (doseq [indx indices] (.set mat (first indx) (second indx) (sample-normal 1)))
          chol (decomp-cholesky scale)
          x (mmult chol mat (trans mat) (trans chol))]
        x)))