summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
Diffstat (limited to 'packages')
-rw-r--r--packages/base/hmatrix-base.cabal2
-rw-r--r--packages/base/src/C/vector-aux.c49
-rw-r--r--packages/base/src/Data/Packed/IO.hs2
-rw-r--r--packages/base/src/Data/Packed/Numeric.hs8
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Base.hs20
-rw-r--r--packages/base/src/Numeric/LinearAlgebra/Random.hs80
-rw-r--r--packages/base/src/Numeric/Vectorized.hs30
-rw-r--r--packages/hmatrix/src/Numeric/Container.hs2
-rw-r--r--packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs61
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-aux.c4
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs4
11 files changed, 184 insertions, 78 deletions
diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal
index 30619e1..b884518 100644
--- a/packages/base/hmatrix-base.cabal
+++ b/packages/base/hmatrix-base.cabal
@@ -23,6 +23,7 @@ library
23 binary, 23 binary,
24 array, 24 array,
25 deepseq, 25 deepseq,
26 random,
26 storable-complex, 27 storable-complex,
27 vector >= 0.8 28 vector >= 0.8
28 29
@@ -56,6 +57,7 @@ library
56 Numeric.Matrix 57 Numeric.Matrix
57 Data.Packed.Internal.Numeric 58 Data.Packed.Internal.Numeric
58 Numeric.LinearAlgebra.Util.Convolution 59 Numeric.LinearAlgebra.Util.Convolution
60 Numeric.LinearAlgebra.Random
59 Numeric.Conversion 61 Numeric.Conversion
60 62
61 C-sources: src/C/lapack-aux.c 63 C-sources: src/C/lapack-aux.c
diff --git a/packages/base/src/C/vector-aux.c b/packages/base/src/C/vector-aux.c
index f1bb371..5b9c171 100644
--- a/packages/base/src/C/vector-aux.c
+++ b/packages/base/src/C/vector-aux.c
@@ -695,3 +695,52 @@ int saveMatrix(char * file, char * format, KDMAT(a)){
695 OK 695 OK
696} 696}
697 697
698////////////////////////////////////////////////////////////////////////////////
699
700// http://c-faq.com/lib/gaussian.html
701double gaussrand()
702{
703 static double V1, V2, S;
704 static int phase = 0;
705 double X;
706
707 if(phase == 0) {
708 do {
709 double U1 = (double)rand() / RAND_MAX;
710 double U2 = (double)rand() / RAND_MAX;
711
712 V1 = 2 * U1 - 1;
713 V2 = 2 * U2 - 1;
714 S = V1 * V1 + V2 * V2;
715 } while(S >= 1 || S == 0);
716
717 X = V1 * sqrt(-2 * log(S) / S);
718 } else
719 X = V2 * sqrt(-2 * log(S) / S);
720
721 phase = 1 - phase;
722
723 return X;
724}
725
726int random_vector(int seed, int code, DVEC(r)) {
727 srand(seed);
728 int k;
729 switch (code) {
730 case 0: { // uniform
731 for (k=0; k<rn; k++) {
732 rp[k] = (double)rand()/RAND_MAX;
733 }
734 OK
735 }
736 case 1: { // gaussian
737 for (k=0; k<rn; k++) {
738 rp[k] = gaussrand();
739 }
740 OK
741 }
742
743 default: ERROR(BAD_CODE);
744 }
745}
746
diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs
index db03d5f..94fb30a 100644
--- a/packages/base/src/Data/Packed/IO.hs
+++ b/packages/base/src/Data/Packed/IO.hs
@@ -149,6 +149,8 @@ apparentCols s = f . dropWhile null . map words . lines <$> readFile s
149 f [] = 0 149 f [] = 0
150 f (x:_) = length x 150 f (x:_) = length x
151 151
152
153-- | load a matrix from an ASCII file formatted as a 2D table.
152loadMatrix :: FilePath -> IO (Matrix Double) 154loadMatrix :: FilePath -> IO (Matrix Double)
153loadMatrix f = do 155loadMatrix f = do
154 v <- vectorScan f 156 v <- vectorScan f
diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs
index 6036e8c..d130ecd 100644
--- a/packages/base/src/Data/Packed/Numeric.hs
+++ b/packages/base/src/Data/Packed/Numeric.hs
@@ -84,7 +84,7 @@ linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n
84class Contraction a b c | a b -> c 84class Contraction a b c | a b -> c
85 where 85 where
86 infixl 7 <.> 86 infixl 7 <.>
87 {- | Matrix product, matrix vector product, and dot product 87 {- | Matrix product, matrix - vector product, and dot product
88 88
89Examples: 89Examples:
90 90
@@ -223,11 +223,7 @@ instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
223 223
224-------------------------------------------------------------------------------- 224--------------------------------------------------------------------------------
225 225
226{- | alternative operator for '(\<.\>)' 226-- | alternative unicode symbol (25c7) for the contraction operator '(\<.\>)'
227
228x25c7, white diamond
229
230-}
231(◇) :: Contraction a b c => a -> b -> c 227(◇) :: Contraction a b c => a -> b -> c
232infixl 7 ◇ 228infixl 7 ◇
233(◇) = (<.>) 229(◇) = (<.>)
diff --git a/packages/base/src/Numeric/LinearAlgebra/Base.hs b/packages/base/src/Numeric/LinearAlgebra/Base.hs
index 395c84a..8d44d26 100644
--- a/packages/base/src/Numeric/LinearAlgebra/Base.hs
+++ b/packages/base/src/Numeric/LinearAlgebra/Base.hs
@@ -35,10 +35,13 @@ module Numeric.LinearAlgebra.Base (
35 -- , 5.0, 5.0, 7.0 ] 35 -- , 5.0, 5.0, 7.0 ]
36 -- 36 --
37 37
38 -- * Products 38 -- * Matrix product
39 (<.>), 39 (<.>),
40 40
41 -- | The matrix product is also implemented in the "Data.Monoid" instance for Matrix, where 41 -- | This operator can also be written using the unicode symbol ◇ (25c7).
42 --
43
44 -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where
42 -- single-element matrices (created from numeric literals or using 'scalar') 45 -- single-element matrices (created from numeric literals or using 'scalar')
43 -- are used for scaling. 46 -- are used for scaling.
44 -- 47 --
@@ -48,9 +51,9 @@ module Numeric.LinearAlgebra.Base (
48 -- [ 1.0, 4.0, 0.0 51 -- [ 1.0, 4.0, 0.0
49 -- , 4.0, 10.0, 0.0 ] 52 -- , 4.0, 10.0, 0.0 ]
50 -- 53 --
51 -- mconcat uses 'optimiseMult' to get the optimal association order. 54 -- 'mconcat' uses 'optimiseMult' to get the optimal association order.
52 55
53 (◇), 56 -- * Other products
54 outer, kronecker, cross, 57 outer, kronecker, cross,
55 scale, 58 scale,
56 sumElements, prodElements, absSum, 59 sumElements, prodElements, absSum,
@@ -114,15 +117,15 @@ module Numeric.LinearAlgebra.Base (
114 -- * Norms 117 -- * Norms
115 norm1, norm2, normInf, pnorm, NormType(..), 118 norm1, norm2, normInf, pnorm, NormType(..),
116 119
117 -- * Correlation and Convolution 120 -- * Correlation and convolution
118 corr, conv, corrMin, corr2, conv2, 121 corr, conv, corrMin, corr2, conv2,
119 122
120 -- * Random arrays 123 -- * Random arrays
121 124
122 -- | rand, randn, RandDist(..), randomVector, gaussianSample, uniformSample 125 RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample,
123 126
124 -- * Misc 127 -- * Misc
125 meanCov, peps, relativeError, haussholder, optimiseMult, udot 128 meanCov, peps, relativeError, haussholder, optimiseMult, udot, Seed, (◇)
126) where 129) where
127 130
128import Numeric.LinearAlgebra.Data 131import Numeric.LinearAlgebra.Data
@@ -132,6 +135,7 @@ import Numeric.Vector()
132import Data.Packed.Numeric 135import Data.Packed.Numeric
133import Numeric.LinearAlgebra.Algorithms 136import Numeric.LinearAlgebra.Algorithms
134import Numeric.LinearAlgebra.Util 137import Numeric.LinearAlgebra.Util
138import Numeric.LinearAlgebra.Random
135 139
136 140
137 141
diff --git a/packages/base/src/Numeric/LinearAlgebra/Random.hs b/packages/base/src/Numeric/LinearAlgebra/Random.hs
new file mode 100644
index 0000000..b36c7a3
--- /dev/null
+++ b/packages/base/src/Numeric/LinearAlgebra/Random.hs
@@ -0,0 +1,80 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.LinearAlgebra.Random
4-- Copyright : (c) Alberto Ruiz 2009-14
5-- License : BSD3
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9-- Random vectors and matrices.
10--
11-----------------------------------------------------------------------------
12
13module Numeric.LinearAlgebra.Random (
14 Seed,
15 RandDist(..),
16 randomVector,
17 gaussianSample,
18 uniformSample,
19 rand, randn
20) where
21
22import Numeric.Vectorized
23import Data.Packed.Numeric
24import Numeric.LinearAlgebra.Algorithms
25import System.Random(randomIO)
26
27
28-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
29-- Gaussian distribution.
30gaussianSample :: Seed
31 -> Int -- ^ number of rows
32 -> Vector Double -- ^ mean vector
33 -> Matrix Double -- ^ covariance matrix
34 -> Matrix Double -- ^ result
35gaussianSample seed n med cov = m where
36 c = dim med
37 meds = konst' 1 n `outer` med
38 rs = reshape c $ randomVector seed Gaussian (c * n)
39 m = rs `mXm` cholSH cov `add` meds
40
41-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
42-- uniform distribution.
43uniformSample :: Seed
44 -> Int -- ^ number of rows
45 -> [(Double,Double)] -- ^ ranges for each column
46 -> Matrix Double -- ^ result
47uniformSample seed n rgs = m where
48 (as,bs) = unzip rgs
49 a = fromList as
50 cs = zipWith subtract as bs
51 d = dim a
52 dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
53 am = konst' 1 n `outer` a
54 m = fromColumns (zipWith scale cs dat) `add` am
55
56-- | pseudorandom matrix with uniform elements between 0 and 1
57randm :: RandDist
58 -> Int -- ^ rows
59 -> Int -- ^ columns
60 -> IO (Matrix Double)
61randm d r c = do
62 seed <- randomIO
63 return (reshape c $ randomVector seed d (r*c))
64
65-- | pseudorandom matrix with uniform elements between 0 and 1
66rand :: Int -> Int -> IO (Matrix Double)
67rand = randm Uniform
68
69{- | pseudorandom matrix with normal elements
70
71>>> disp 3 =<< randn 3 5
723x5
730.386 -1.141 0.491 -0.510 1.512
740.069 -0.919 1.022 -0.181 0.745
750.313 -0.670 -0.097 -1.575 -0.583
76
77-}
78randn :: Int -> Int -> IO (Matrix Double)
79randn = randm Gaussian
80
diff --git a/packages/base/src/Numeric/Vectorized.hs b/packages/base/src/Numeric/Vectorized.hs
index a2d7f70..5aebb14 100644
--- a/packages/base/src/Numeric/Vectorized.hs
+++ b/packages/base/src/Numeric/Vectorized.hs
@@ -17,7 +17,8 @@ module Numeric.Vectorized (
17 FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ, 17 FunCodeV(..), vectorMapR, vectorMapC, vectorMapF, vectorMapQ,
18 FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ, 18 FunCodeSV(..), vectorMapValR, vectorMapValC, vectorMapValF, vectorMapValQ,
19 FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ, 19 FunCodeVV(..), vectorZipR, vectorZipC, vectorZipF, vectorZipQ,
20 vectorScan, saveMatrix 20 vectorScan, saveMatrix,
21 Seed, RandDist(..), randomVector
21) where 22) where
22 23
23import Data.Packed.Internal.Common 24import Data.Packed.Internal.Common
@@ -308,7 +309,13 @@ vectorScan s = do
308foreign import ccall unsafe "saveMatrix" c_saveMatrix 309foreign import ccall unsafe "saveMatrix" c_saveMatrix
309 :: CString -> CString -> TM 310 :: CString -> CString -> TM
310 311
311saveMatrix :: FilePath -> String -> Matrix Double -> IO () 312{- | save a matrix as a 2D ASCII table
313-}
314saveMatrix
315 :: FilePath
316 -> String -- ^ \"printf\" format (e.g. \"%.2f\", \"%g\", etc.)
317 -> Matrix Double
318 -> IO ()
312saveMatrix name format m = do 319saveMatrix name format m = do
313 cname <- newCString name 320 cname <- newCString name
314 cformat <- newCString format 321 cformat <- newCString format
@@ -317,4 +324,23 @@ saveMatrix name format m = do
317 free cformat 324 free cformat
318 return () 325 return ()
319 326
327--------------------------------------------------------------------------------
328
329type Seed = Int
330
331data RandDist = Uniform -- ^ uniform distribution in [0,1)
332 | Gaussian -- ^ normal distribution with mean zero and standard deviation one
333 deriving Enum
334
335-- | Obtains a vector of pseudorandom elements (use randomIO to get a random seed).
336randomVector :: Seed
337 -> RandDist -- ^ distribution
338 -> Int -- ^ vector size
339 -> Vector Double
340randomVector seed dist n = unsafePerformIO $ do
341 r <- createVector n
342 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
343 return r
344
345foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
320 346
diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs
index b6e797b..2b61d90 100644
--- a/packages/hmatrix/src/Numeric/Container.hs
+++ b/packages/hmatrix/src/Numeric/Container.hs
@@ -16,7 +16,7 @@ module Numeric.Container (
16 meanCov 16 meanCov
17) where 17) where
18 18
19import Data.Packed.Numeric 19import Data.Packed.Numeric hiding (saveMatrix, loadMatrix)
20import Numeric.LinearAlgebra.IO 20import Numeric.LinearAlgebra.IO
21import Numeric.LinearAlgebra.Random hiding (Seed) 21import Numeric.LinearAlgebra.Random hiding (Seed)
22import Numeric.LinearAlgebra.Util(meanCov) 22import Numeric.LinearAlgebra.Util(meanCov)
diff --git a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
index 43e46c5..8bd0cc2 100644
--- a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
+++ b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
@@ -9,15 +9,15 @@
9----------------------------------------------------------------------------- 9-----------------------------------------------------------------------------
10 10
11module Numeric.GSL.LinearAlgebra ( 11module Numeric.GSL.LinearAlgebra (
12 RandDist(..), randomVector, 12 randomVector,
13 saveMatrix, 13 saveMatrix,
14 fwriteVector, freadVector, fprintfVector, fscanfVector 14 fwriteVector, freadVector, fprintfVector, fscanfVector
15) where 15) where
16 16
17import Data.Packed 17import Data.Packed
18import Numeric.LinearAlgebra.Base(RandDist(..))
18import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) 19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
19 20
20import Data.Complex
21import Foreign.Marshal.Alloc(free) 21import Foreign.Marshal.Alloc(free)
22import Foreign.Ptr(Ptr) 22import Foreign.Ptr(Ptr)
23import Foreign.C.Types 23import Foreign.C.Types
@@ -28,10 +28,6 @@ fromei x = fromIntegral (fromEnum x) :: CInt
28 28
29----------------------------------------------------------------------- 29-----------------------------------------------------------------------
30 30
31data RandDist = Uniform -- ^ uniform distribution in [0,1)
32 | Gaussian -- ^ normal distribution with mean zero and standard deviation one
33 deriving Enum
34
35-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed. 31-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
36randomVector :: Int -- ^ seed 32randomVector :: Int -- ^ seed
37 -> RandDist -- ^ distribution 33 -> RandDist -- ^ distribution
@@ -39,10 +35,10 @@ randomVector :: Int -- ^ seed
39 -> Vector Double 35 -> Vector Double
40randomVector seed dist n = unsafePerformIO $ do 36randomVector seed dist n = unsafePerformIO $ do
41 r <- createVector n 37 r <- createVector n
42 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector" 38 app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL"
43 return r 39 return r
44 40
45foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV 41foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV
46 42
47-------------------------------------------------------------------------------- 43--------------------------------------------------------------------------------
48 44
@@ -105,56 +101,7 @@ fwriteVector filename v = do
105 101
106foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV 102foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
107 103
108type PF = Ptr Float --
109type PD = Ptr Double -- 104type PD = Ptr Double --
110type PQ = Ptr (Complex Float) --
111type PC = Ptr (Complex Double) --
112type TF = CInt -> PF -> IO CInt --
113type TFF = CInt -> PF -> TF --
114type TFV = CInt -> PF -> TV --
115type TVF = CInt -> PD -> TF --
116type TFFF = CInt -> PF -> TFF --
117type TV = CInt -> PD -> IO CInt -- 105type TV = CInt -> PD -> IO CInt --
118type TVV = CInt -> PD -> TV --
119type TVVV = CInt -> PD -> TVV --
120type TFM = CInt -> CInt -> PF -> IO CInt --
121type TFMFM = CInt -> CInt -> PF -> TFM --
122type TFMFMFM = CInt -> CInt -> PF -> TFMFM --
123type TM = CInt -> CInt -> PD -> IO CInt -- 106type TM = CInt -> CInt -> PD -> IO CInt --
124type TMM = CInt -> CInt -> PD -> TM --
125type TVMM = CInt -> PD -> TMM --
126type TMVMM = CInt -> CInt -> PD -> TVMM --
127type TMMM = CInt -> CInt -> PD -> TMM --
128type TVM = CInt -> PD -> TM --
129type TVVM = CInt -> PD -> TVM --
130type TMV = CInt -> CInt -> PD -> TV --
131type TMMV = CInt -> CInt -> PD -> TMV --
132type TMVM = CInt -> CInt -> PD -> TVM --
133type TMMVM = CInt -> CInt -> PD -> TMVM --
134type TCM = CInt -> CInt -> PC -> IO CInt --
135type TCVCM = CInt -> PC -> TCM --
136type TCMCVCM = CInt -> CInt -> PC -> TCVCM --
137type TMCMCVCM = CInt -> CInt -> PD -> TCMCVCM --
138type TCMCMCVCM = CInt -> CInt -> PC -> TCMCVCM --
139type TCMCM = CInt -> CInt -> PC -> TCM --
140type TVCM = CInt -> PD -> TCM --
141type TCMVCM = CInt -> CInt -> PC -> TVCM --
142type TCMCMVCM = CInt -> CInt -> PC -> TCMVCM --
143type TCMCMCM = CInt -> CInt -> PC -> TCMCM --
144type TCV = CInt -> PC -> IO CInt --
145type TCVCV = CInt -> PC -> TCV --
146type TCVCVCV = CInt -> PC -> TCVCV --
147type TCVV = CInt -> PC -> TV --
148type TQV = CInt -> PQ -> IO CInt --
149type TQVQV = CInt -> PQ -> TQV --
150type TQVQVQV = CInt -> PQ -> TQVQV --
151type TQVF = CInt -> PQ -> TF --
152type TQM = CInt -> CInt -> PQ -> IO CInt --
153type TQMQM = CInt -> CInt -> PQ -> TQM --
154type TQMQMQM = CInt -> CInt -> PQ -> TQMQM --
155type TCMCV = CInt -> CInt -> PC -> TCV --
156type TVCV = CInt -> PD -> TCV --
157type TCVM = CInt -> PC -> TM --
158type TMCVM = CInt -> CInt -> PD -> TCVM --
159type TMMCVM = CInt -> CInt -> PD -> TMCVM --
160 107
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
index 5da94ca..ffc5c20 100644
--- a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
+++ b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
@@ -924,8 +924,8 @@ int nlfit(int method, int f(int, double*, int, double*),
924 924
925#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK } 925#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }
926 926
927int random_vector(int seed, int code, RVEC(r)) { 927int random_vector_GSL(int seed, int code, RVEC(r)) {
928 DEBUGMSG("random_vector") 928 DEBUGMSG("random_vector_GSL")
929 static gsl_rng * gen = NULL; 929 static gsl_rng * gen = NULL;
930 if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);} 930 if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
931 gsl_rng_set (gen, seed); 931 gsl_rng_set (gen, seed);
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
index 0a82d3f..fa125a0 100644
--- a/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
@@ -22,11 +22,11 @@ module Numeric.LinearAlgebra.Random (
22 22
23import Numeric.GSL.LinearAlgebra 23import Numeric.GSL.LinearAlgebra
24import Data.Packed.Numeric 24import Data.Packed.Numeric
25import Numeric.LinearAlgebra.Algorithms 25import Numeric.LinearAlgebra.Base(Seed,RandDist(..),cholSH)
26import System.Random(randomIO) 26import System.Random(randomIO)
27 27
28 28
29type Seed = Int 29
30 30
31-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate 31-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
32-- Gaussian distribution. 32-- Gaussian distribution.