random-fu-multivariate (empty) → 0.1.1.1

raw patch · 10 files changed

+457/−0 lines, 10 filesdep +basedep +hmatrixdep +mtlsetup-changed

Dependencies added: base, hmatrix, mtl, random-fu, random-fu-multivariate

Files

LICENSE +30/−0
README.md +4/−0
Setup.hs +2/−0
diagrams/src_Data_Random_Distribution_MultivariateNormal_diagM.svg too large to diff
diagrams/src_Data_Random_Distribution_Static_MultivariateNormal_diagMS.svg too large to diff
random-fu-multivariate.cabal +40/−0
src/Data/Random/Distribution/MultiNormal.hs +109/−0
src/Data/Random/Distribution/MultivariateNormal.hs +124/−0
src/Data/Random/Distribution/Static/MultivariateNormal.hs +146/−0
test/Spec.hs +2/−0

+ LICENSE view

@@ -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.

+ README.md view

@@ -0,0 +1,4 @@+# Multivariate Distributions in `random-fu`++Presently, this package adds multivariate normal distributions,+implemented using `hmatrix`.

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

+ diagrams/src_Data_Random_Distribution_MultivariateNormal_diagM.svg view

file too large to diff

+ diagrams/src_Data_Random_Distribution_Static_MultivariateNormal_diagMS.svg view

file too large to diff

+ random-fu-multivariate.cabal view

@@ -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

+ src/Data/Random/Distribution/MultiNormal.hs view

@@ -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

+ src/Data/Random/Distribution/MultivariateNormal.hs view

@@ -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

+ src/Data/Random/Distribution/Static/MultivariateNormal.hs view

@@ -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

+ test/Spec.hs view

@@ -0,0 +1,2 @@+main :: IO ()+main = putStrLn "Test suite not yet implemented"