From dc39c11d2d2e75428e6394cf8542c8c3ff2cd887 Mon Sep 17 00:00:00 2001
From: Alberto Ruiz <aruiz@um.es>
Date: Sat, 14 Nov 2009 12:57:52 +0000
Subject: added module Random.hs

---
 lib/Data/Packed/Random.hs          | 52 ++++++++++++++++++++++++++++++++++++++
 lib/Data/Packed/Vector.hs          |  3 +--
 lib/Numeric/LinearAlgebra.hs       |  2 ++
 lib/Numeric/LinearAlgebra/Tests.hs | 11 ++++++++
 4 files changed, 66 insertions(+), 2 deletions(-)
 create mode 100644 lib/Data/Packed/Random.hs

(limited to 'lib')

diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs
new file mode 100644
index 0000000..56d038e
--- /dev/null
+++ b/lib/Data/Packed/Random.hs
@@ -0,0 +1,52 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Data.Packed.Vector
+-- Copyright   :  (c) Alberto Ruiz 2009
+-- License     :  GPL
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+--
+-- Random vectors and matrices.
+--
+-----------------------------------------------------------------------------
+
+module Data.Packed.Random (
+    RandDist(..),
+    randomVector,
+    gaussianSample,
+    meanCov,
+) where
+
+import Numeric.GSL.Vector
+import Data.Packed
+import Numeric.LinearAlgebra.Algorithms
+import Numeric.LinearAlgebra.Instances()
+import Numeric.LinearAlgebra.Interface
+
+-- | Obtains a matrix whose rows are pseudorandom samples from a multivariante
+-- Gaussian distribution.
+gaussianSample :: Int -- ^ seed
+               -> Int -- ^ number of rows
+               -> Vector Double -- ^ mean vector
+               -> Matrix Double -- ^ covariance matrix
+               -> Matrix Double -- ^ result
+gaussianSample seed n med cov = m where
+    (l,v) = eigSH' cov
+    c = dim l
+    meds = constant 1 n `outer` med
+    rs = reshape c $ randomVector seed Gaussian (c * n)
+    ds = sqrt (abs l)
+    m = rs <> (diag ds <> trans v) + meds
+
+------------ utilities -------------------------------
+
+-- | Compute mean vector and covariance matrix of the rows of a matrix.
+meanCov :: Matrix Double -> (Vector Double, Matrix Double)
+meanCov x = (med,cov) where
+    r    = rows x
+    k    = 1 / fromIntegral r
+    med  = constant k r <> x
+    meds = constant 1 r `outer` med
+    xc   = x - meds
+    cov  = (trans xc <> xc) / fromIntegral (r-1)
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 836f11a..21f51e5 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -21,8 +21,7 @@ module Data.Packed.Vector (
     vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex,
     mapVector, zipVector,
     fscanfVector, fprintfVector, freadVector, fwriteVector,
-    foldLoop, foldVector, foldVectorG,
-    RandDist(..), randomVector
+    foldLoop, foldVector, foldVectorG
 ) where
 
 import Data.Packed.Internal
diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs
index f92c40d..337c007 100644
--- a/lib/Numeric/LinearAlgebra.hs
+++ b/lib/Numeric/LinearAlgebra.hs
@@ -16,12 +16,14 @@ This module reexports the most comon functions (including "Numeric.LinearAlgebra
 -----------------------------------------------------------------------------
 module Numeric.LinearAlgebra (
     module Data.Packed,
+    module Data.Packed.Random,
     module Numeric.LinearAlgebra.Linear,
     module Numeric.LinearAlgebra.Algorithms,
     module Numeric.LinearAlgebra.Interface
 ) where
 
 import Data.Packed
+import Data.Packed.Random
 import Numeric.LinearAlgebra.Linear
 import Numeric.LinearAlgebra.Algorithms
 import Numeric.LinearAlgebra.Instances()
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index d7c3f99..097756e 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -131,6 +131,16 @@ rootFindingTest = TestList [ utest "root Hybrids" (fst sol1 ~~ [1,1])
 
 ---------------------------------------------------------------------
 
+randomTest = c :~1~: snd (meanCov dat) where
+    a = (3><3) [1,2,3,
+                2,4,0,
+               -2,2,1]
+    m = 3 |> [1,2,3]
+    c = a <> trans a
+    dat = gaussianSample 7 (10^6) m c
+
+---------------------------------------------------------------------
+
 rot :: Double -> Matrix Double
 rot a = (3><3) [ c,0,s
                , 0,1,0
@@ -227,6 +237,7 @@ runTests n = do
         , utest "polySolve" (polySolveProp [1,2,3,4])
         , minimizationTest
         , rootFindingTest
+        , utest "random" randomTest
         ]
     return ()
 
-- 
cgit v1.2.3