diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Author name here (c) 2016
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Author name here nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,4 @@
+# Multivariate Distributions in `random-fu`
+
+Presently, this package adds multivariate normal distributions,
+implemented using `hmatrix`.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/diagrams/src_Data_Random_Distribution_MultivariateNormal_diagM.svg b/diagrams/src_Data_Random_Distribution_MultivariateNormal_diagM.svg
new file mode 100644
# file too large to diff: diagrams/src_Data_Random_Distribution_MultivariateNormal_diagM.svg
diff --git a/diagrams/src_Data_Random_Distribution_Static_MultivariateNormal_diagMS.svg b/diagrams/src_Data_Random_Distribution_Static_MultivariateNormal_diagMS.svg
new file mode 100644
# file too large to diff: diagrams/src_Data_Random_Distribution_Static_MultivariateNormal_diagMS.svg
diff --git a/random-fu-multivariate.cabal b/random-fu-multivariate.cabal
new file mode 100644
--- /dev/null
+++ b/random-fu-multivariate.cabal
@@ -0,0 +1,40 @@
+name:                random-fu-multivariate
+version:             0.1.1.1
+synopsis:            Multivariate distributions for random-fu
+description:         Please see README.md
+homepage:            https://github.com/fpco/random-fu-multivariate
+license:             BSD3
+license-file:        LICENSE
+author:              Dominic Steinitz, Jacob West
+maintainer:          dominic@steinitz.org
+copyright:           (c) 2016 FP Complete Corporation
+category:            Math
+build-type:          Simple
+cabal-version:       >=1.10
+extra-source-files:  README.md, diagrams/*.svg
+extra-doc-files:     diagrams/*.svg
+
+source-repository head
+  type:     git
+  location: https://github.com/fpco/random-fu-multivariate
+
+library
+  default-language:  Haskell2010
+  hs-source-dirs:    src
+  exposed-modules:   Data.Random.Distribution.MultivariateNormal
+                   , Data.Random.Distribution.Static.MultivariateNormal
+                   , Data.Random.Distribution.MultiNormal
+  ghc-options:       -Wall
+  build-depends:     base >= 4.7 && < 5
+                   , random-fu
+                   , hmatrix
+                   , mtl
+
+test-suite random-fu-multivariate-test
+  default-language:  Haskell2010
+  hs-source-dirs:    test
+  main-is:           Spec.hs
+  type:              exitcode-stdio-1.0
+  ghc-options:       -Wall -threaded -rtsopts -with-rtsopts=-N
+  build-depends:     base >= 4.7 && < 5
+                   , random-fu-multivariate
diff --git a/src/Data/Random/Distribution/MultiNormal.hs b/src/Data/Random/Distribution/MultiNormal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/MultiNormal.hs
@@ -0,0 +1,109 @@
+--------------------------------------------------------------
+--- An implementation of multivariate normal distributions ---
+--------------------------------------------------------------
+{-
+Written by: Dominic Steinitz, Jacob West
+Last modified: 2016-07-27
+
+Summary: Multivariate normal distributions are necessary for Kalman
+filters and smoothers.  However, strictly speaking, the functionality
+provided here should exist elsewhere, perhaps in the package:
+random-fu.
+-}
+
+---------------------------
+--- File header pragmas ---
+---------------------------
+{-# LANGUAGE RecordWildCards #-}       -- Used by multiNormalRV, multiNormalConstant
+                                       -- and multiNormalQuadraticForm
+{-# LANGUAGE MultiParamTypeClasses #-} -- Necessary for Distribution instance
+{-# LANGUAGE FlexibleInstances #-}     -- Necessary for Show instance
+{-# LANGUAGE TypeFamilies #-}          -- Necessary for MultiNormal definition
+
+------------------------
+--- Module / Exports ---
+------------------------
+module Data.Random.Distribution.MultiNormal
+       (
+         MultiNormal(..)
+       , inv
+       )
+       where
+---------------
+--- Imports ---
+---------------
+import Control.Monad (replicateM, when)
+import Data.Maybe (fromMaybe, fromJust)
+import Data.Random
+import GHC.TypeLits
+import Numeric.LinearAlgebra.Static
+
+import qualified Numeric.LinearAlgebra as LA
+
+------------------------
+--- Helper Functions ---
+------------------------
+-- Matrix inverse: for some reason, this isn't built into the
+-- static interface; warning: no error handling
+
+-- WARNING: Needs better error handling
+inv :: KnownNat n => Sq n -> Sq n
+inv = fromMaybe (error "Failed attempting to invert non-invertible matrix.") .
+      flip linSolve eye
+
+----------------------------------------
+--- Multivariate Normal Distrubtions ---
+----------------------------------------
+-- This probably belongs elsewhere, maybe Data.Random, but that would
+-- create a dependence on Numric.LinearAlgebra which I believe is not
+-- there now and may be undesirable.
+
+data family MultiNormal k :: *
+data instance KnownNat n => MultiNormal (R n) =
+  MultiNormal { mu :: (R n), cov :: (Sym n) }
+
+--- Show Instance ---
+instance KnownNat n => Show (MultiNormal (R n)) where
+  show MultiNormal {..} = "Normal " ++ show mu ++ " " ++ show cov
+
+--- Distribution Instance ---
+instance KnownNat n => Distribution MultiNormal (R n) where
+  rvar = multiNormalRV
+
+-- WARNING: Needs better error handling
+multiNormalRV :: KnownNat n => MultiNormal (R n) -> RVarT m (R n)
+multiNormalRV MultiNormal {..} = do
+  let (vals, vecs) = eigensystem cov
+  when (any (<0) (LA.toList $ unwrap vals))
+    (error "Covariance matrix is not positive semi-definite.")
+
+  let lSqrt = diag (fromJust . create $ LA.cmap sqrt (extract vals))
+      bigA  = tr vecs <> lSqrt
+  
+  gnoise <- replicateM (size mu) (rvarT StdNormal)
+  return $ mu + bigA #> (vector gnoise)
+
+--- PDF Instance ---
+instance KnownNat n => PDF MultiNormal (R n) where
+  pdf    = multiNormalPDF
+  logPdf = multiNormalLogPDF
+
+multiNormalPDF :: KnownNat n => MultiNormal (R n) -> R n -> Double
+multiNormalPDF mn pt =
+  multiNormalConstant mn * exp (multiNormalQuadraticForm mn pt)
+
+multiNormalLogPDF :: KnownNat n => MultiNormal (R n) -> R n -> Double
+multiNormalLogPDF mn pt =
+  multiNormalConstant mn + multiNormalQuadraticForm mn pt
+
+multiNormalConstant :: KnownNat n => MultiNormal (R n) -> Double
+multiNormalConstant MultiNormal {..} = recip . sqrt $ (2*pi)^n * detCov
+  where
+    n = size mu
+    detCov = LA.det . extract . unSym $ cov
+
+multiNormalQuadraticForm :: KnownNat n => MultiNormal (R n) -> R n -> Double
+multiNormalQuadraticForm MultiNormal {..} pt = (diff LA.<.> invCov LA.#> diff) / (-2)
+  where
+    diff = extract (mu - pt)
+    invCov = extract . inv . unSym $ cov
diff --git a/src/Data/Random/Distribution/MultivariateNormal.hs b/src/Data/Random/Distribution/MultivariateNormal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/MultivariateNormal.hs
@@ -0,0 +1,124 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Random.Distribution.MultivariateNormal
+-- Copyright   :  (c) 2016 FP Complete Corporation
+-- License     :  MIT (see LICENSE)
+-- Maintainer  :  dominic@steinitz.org
+--
+-- 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.
+--
+-- <<diagrams/src_Data_Random_Distribution_MultivariateNormal_diagM.svg#diagram=diagM&height=600&width=500>>
+--
+-- Example code to generate the chart:
+--
+-- > import qualified Graphics.Rendering.Chart as C
+-- > import Graphics.Rendering.Chart.Backend.Diagrams
+-- >
+-- > import Data.Random.Distribution.MultivariateNormal
+-- >
+-- > import qualified Data.Random as R
+-- > import Data.Random.Source.PureMT
+-- > import Control.Monad.State
+-- > import qualified Numeric.LinearAlgebra.HMatrix as LA
+-- >
+-- > nSamples :: Int
+-- > nSamples = 10000
+-- >
+-- > sigma1, sigma2, rho :: Double
+-- > sigma1 = 3.0
+-- > sigma2 = 1.0
+-- > rho = 0.5
+-- >
+-- > singleSample :: R.RVarT (State PureMT) (LA.Vector Double)
+-- > singleSample = R.sample $ Normal (LA.fromList [0.0, 0.0])
+-- >                (LA.sym $ (2 LA.>< 2) [ sigma1, rho * sigma1 * sigma2
+-- >                                      , rho * sigma1 * sigma2, sigma2])
+-- >
+-- > multiSamples :: [LA.Vector Double]
+-- > multiSamples = evalState (replicateM nSamples $ R.sample singleSample) (pureMT 3)
+-- > pts = map (f . LA.toList) multiSamples
+-- >   where
+-- >     f [x, y] = (x, y)
+-- >     f _      = error "Only pairs for this chart"
+-- >
+-- >
+-- > 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
+-- >
+-- > diagM = do
+-- >   denv <- defaultEnv C.vectorAlignmentFns 600 500
+-- >   return $ fst $ runBackend denv (C.render (chartPoint pts nSamples) (500, 500))
+--
+-----------------------------------------------------------------------------
+
+{-# LANGUAGE TypeFamilies          #-}
+{-# LANGUAGE FlexibleInstances     #-}
+{-# LANGUAGE FlexibleContexts      #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+
+module Data.Random.Distribution.MultivariateNormal
+    ( Normal(..)
+    ) where
+
+import           Data.Random.Distribution
+import qualified Numeric.LinearAlgebra.HMatrix as H
+import           Control.Monad
+import qualified Data.Random as R
+import           Foreign.Storable ( Storable )
+import           Data.Maybe ( fromJust )
+
+normalMultivariate :: H.Vector Double -> H.Herm Double -> R.RVarT m (H.Vector Double)
+normalMultivariate mu bigSigma = do
+  z <- replicateM (H.size mu) (rvarT R.StdNormal)
+  return $ mu + bigA H.#> (H.fromList z)
+  where
+    (vals, bigU) = H.eigSH bigSigma
+    lSqrt = H.diag $ H.cmap sqrt vals
+    bigA = bigU H.<> lSqrt
+
+data family Normal k :: *
+
+data instance Normal (H.Vector Double) = Normal (H.Vector Double) (H.Herm Double)
+
+instance Distribution Normal (H.Vector Double) where
+  rvar (Normal m s) = normalMultivariate m s
+
+normalPdf :: (H.Numeric a, H.Field a, H.Indexable (H.Vector a) a, Num (H.Vector a)) =>
+             H.Vector a -> H.Herm a -> H.Vector a -> a
+normalPdf mu sigma x = exp $ normalLogPdf mu sigma x
+
+normalLogPdf :: (H.Numeric a, H.Field a, H.Indexable (H.Vector a) a, Num (H.Vector a)) =>
+                 H.Vector a -> H.Herm a -> H.Vector a -> a
+normalLogPdf mu bigSigma x = - H.sumElements (H.cmap log (diagonals dec))
+                              - 0.5 * (fromIntegral (H.size mu)) * log (2 * pi)
+                              - 0.5 * s
+  where
+    dec = fromJust $ H.mbChol bigSigma
+    t = fromJust $ H.linearSolve (H.tr dec) (H.asColumn $ x - mu)
+    u = H.cmap (\v -> v * v) t
+    s = H.sumElements u
+
+diagonals :: (Storable a, H.Element t, H.Indexable (H.Vector t) a) =>
+             H.Matrix t -> H.Vector a
+diagonals m = H.fromList (map (\i -> m H.! i H.! i) [0..n-1])
+  where
+    n = max (H.rows m) (H.cols m)
+
+instance PDF Normal (H.Vector Double) where
+  pdf (Normal m s) = normalPdf m s
+  logPdf (Normal m s) = normalLogPdf m s
diff --git a/src/Data/Random/Distribution/Static/MultivariateNormal.hs b/src/Data/Random/Distribution/Static/MultivariateNormal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Static/MultivariateNormal.hs
@@ -0,0 +1,146 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Random.Distribution.Static.MultivariateNormal
+-- Copyright   :  (c) 2016 FP Complete Corporation
+-- License     :  MIT (see LICENSE)
+-- Maintainer  :  dominic@steinitz.org
+--
+-- 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\).
+--
+-- <<diagrams/src_Data_Random_Distribution_Static_MultivariateNormal_diagMS.svg#diagram=diagMS&height=600&width=500>>
+--
+-- 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))
+--
+-----------------------------------------------------------------------------
+
+{-# OPTIONS_GHC -Wall                     #-}
+{-# OPTIONS_GHC -fno-warn-name-shadowing  #-}
+{-# OPTIONS_GHC -fno-warn-type-defaults   #-}
+{-# OPTIONS_GHC -fno-warn-unused-do-bind  #-}
+{-# OPTIONS_GHC -fno-warn-missing-methods #-}
+{-# OPTIONS_GHC -fno-warn-orphans         #-}
+
+{-# LANGUAGE MultiParamTypeClasses        #-}
+{-# LANGUAGE TypeFamilies                 #-}
+{-# LANGUAGE ScopedTypeVariables          #-}
+{-# LANGUAGE DataKinds                    #-}
+
+module Data.Random.Distribution.Static.MultivariateNormal
+    ( Normal(..)
+    ) where
+
+import           Data.Random hiding ( StdNormal, Normal )
+import qualified Data.Random as R
+import           Control.Monad.State ( replicateM )
+import qualified Numeric.LinearAlgebra.HMatrix as H
+import           Numeric.LinearAlgebra.Static
+                 ( R, vector, extract, Sq, Sym, col,
+                   tr, linSolve, uncol, chol, (<.>),
+                   ℝ, (<>), diag, (#>), eigensystem
+                 )
+import          GHC.TypeLits ( KnownNat, natVal )
+import          Data.Maybe ( fromJust )
+
+
+normalMultivariate :: KnownNat n =>
+                      R n -> Sym n -> RVarT m (R n)
+normalMultivariate mu bigSigma = do
+  z <- replicateM (fromIntegral $ natVal mu) (rvarT R.StdNormal)
+  return $ mu + bigA #> (vector z)
+  where
+    (vals, bigU) = eigensystem bigSigma
+    lSqrt = diag $ mapVector sqrt vals
+    bigA = bigU <> lSqrt
+
+mapVector :: KnownNat n => (ℝ -> ℝ) -> R n -> R n
+mapVector f = vector . H.toList . H.cmap f . extract
+
+sumVector :: KnownNat n => R n -> ℝ
+sumVector x = x <.> 1
+
+data family Normal k :: *
+
+data instance Normal (R n) = Normal (R n) (Sym n)
+
+instance KnownNat n => Distribution Normal (R n) where
+  rvar (Normal m s) = normalMultivariate m s
+
+normalLogPdf :: KnownNat n =>
+                R n -> Sym n -> R n -> Double
+normalLogPdf mu bigSigma x = - sumVector (mapVector log (diagonals dec))
+                             - 0.5 * (fromIntegral $ natVal mu) * log (2 * pi)
+                             - 0.5 * s
+  where
+    dec = chol bigSigma
+    t = uncol $ fromJust $ linSolve (tr dec) (col $ x - mu)
+    u = mapVector (\x -> x * x) t
+    s = sumVector u
+
+normalPdf :: KnownNat n =>
+             R n -> Sym n -> R n -> Double
+normalPdf mu sigma x = exp $ normalLogPdf mu sigma x
+
+diagonals :: KnownNat n => Sq n -> R n
+diagonals = vector . H.toList . H.takeDiag . extract
+
+instance KnownNat n => PDF Normal (R n) where
+  pdf (Normal m s) = normalPdf m s
+  logPdf (Normal m s) = normalLogPdf m s
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,2 @@
+main :: IO ()
+main = putStrLn "Test suite not yet implemented"