Sample from the multivariate normal distribution with a given vector-valued $\mu$ and covariance matrix $\Sigma$. For example, the chart below shows samples from the bivariate normal distribution. The dimension of the mean $n$ is statically checked to be compatible with the dimension of the covariance matrix $n \times n$.

Example code to generate the chart:

{-# LANGUAGE DataKinds #-}

import qualified Graphics.Rendering.Chart as C
import Graphics.Rendering.Chart.Backend.Diagrams

import Data.Random.Distribution.Static.MultivariateNormal

import qualified Data.Random as R
import Data.Random.Source.PureMT
import Control.Monad.State
import Numeric.LinearAlgebra.Static

nSamples :: Int
nSamples = 10000

sigma1, sigma2, rho :: Double
sigma1 = 3.0
sigma2 = 1.0
rho = 0.5

singleSample :: R.RVarT (State PureMT) (R 2)
singleSample = R.sample $ Normal (vector [0.0, 0.0])
               (sym $ matrix [ sigma1, rho * sigma1 * sigma2
                             , rho * sigma1 * sigma2, sigma2])

multiSamples :: [R 2]
multiSamples = evalState (replicateM nSamples $ R.sample singleSample) (pureMT 3)

pts = map f multiSamples
  where
    f z = (x, y)
      where
        (x, t) = headTail z
        (y, _) = headTail t

chartPoint pointVals n = C.toRenderable layout
  where

    fitted = C.plot_points_values .~ pointVals
              $ C.plot_points_style  . C.point_color .~ opaque red
              $ C.plot_points_title .~ "Sample"
              $ def

    layout = C.layout_title .~ "Sampling Bivariate Normal (" ++ (show n) ++ " samples)"
           $ C.layout_y_axis . C.laxis_generate .~ C.scaledAxis def (-3,3)
           $ C.layout_x_axis . C.laxis_generate .~ C.scaledAxis def (-3,3)

           $ C.layout_plots .~ [C.toPlot fitted]
           $ def

diagMS = do
  denv <- defaultEnv C.vectorAlignmentFns 600 500
  return $ fst $ runBackend denv (C.render (chartPoint pts nSamples) (500, 500))