(defn sample-mvn-params
" Returns samples of means (sampled from an mvn distribution) and vectorized covariance
  matrices (sampled from an inverse-wishart distribution) for the given mvn data.

  Arguments:
    size -- the number of samples to return
    y -- the data used to estimate the parameters


  Returns map with following fields:
    :means
    :sigmas


  Examples:

    (use '(incanter core stats bayes charts))
    (def y (sample-mvn 500 :mean [0 0] :sigma (identity-matrix 2)))
    (def samp (sample-mvn-params 1000 y))

    (map mean (trans (:means samp)))
    (symmetric-matrix (map mean (trans (:sigmas samp))) :lower false)

    (view (histogram (sel (:means samp) :cols 0) :x-label \"mean 1\"))
    (view (histogram (sel (:means samp) :cols 1) :x-label \"mean 2\"))
    (view (histogram (sel (:sigmas samp) :cols 1) :x-label \"covariance\"))
    (view (histogram (sel (:sigmas samp) :cols 0) :x-label \"variance 1\"))
    (view (histogram (sel (:sigmas samp) :cols 2) :x-label \"variance 2\"))

    (map #(quantile % :probs [0.025 0.0975]) (trans (:means samp)))
    (map #(quantile % :probs [0.025 0.0975]) (trans (:sigmas samp)))




    (use '(incanter core stats bayes charts))
    (def y (sample-mvn 500 :sigma (symmetric-matrix [10 5 10]) :mean [5 2]))
    (def samp (sample-mvn-params 1000 y))
    (symmetric-matrix (map mean (trans (:sigmas samp))) :lower false)
    (map mean (trans (:means samp)))


 

"
([^Integer size y & options]
    (let [opts (when options (apply assoc {} options))
          means (map mean (trans y))
          n (count y)
          S (reduce plus
                    (map #(mmult (minus (to-list %) means)
                                 (trans (minus (to-list %) means)))
                         y))
          sigma-samp (matrix (for [_ (range size)]
                               (half-vectorize (sample-inv-wishart :df (dec n) :scale (solve S)))))
          mu-samp (matrix (for [sigma sigma-samp]
                            (sample-mvn 1
                                        :mean means
                                        :sigma (div (symmetric-matrix sigma :lower false) n))))
          ]
  {:means mu-samp :sigmas sigma-samp})))