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 +--
 2 files changed, 53 insertions(+), 2 deletions(-)
 create mode 100644 lib/Data/Packed/Random.hs

(limited to 'lib/Data/Packed')

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
-- 
cgit v1.2.3