6 files changed, 69 insertions, 2 deletions

diff --git a/CHANGES b/CHANGES
index 64450e4..c6297eb 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,6 +1,8 @@
 6.0.0.0
 =======
 
+- added randomVector, gaussianSample, meanCov
+
 - added rankSVD, nullspaceSVD
 
 - rank, nullspacePrec, and economy svd defined in terms of ranksvd.
diff --git a/hmatrix.cabal b/hmatrix.cabal
index 4cd4402..bed36cf 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -129,6 +129,7 @@ library
                         Data.Packed.Convert,
                         Data.Packed.ST,
                         Data.Packed.Development
+                        Data.Packed.Random
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Signatures,
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 ()