• incanter

1.2.3-SNAPSHOT

sample-mvn

incanter.stats

• (sample-mvn size & options)
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

Source incanter/stats.clj:288 top

```(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)))))))```
Vars in incanter.stats/sample-mvn: > defn let map
Used in 0 other vars