1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
-----------------------------------------------------------------------------
-- |
-- Module : Numeric.GSL.Random
-- Copyright : (c) Alberto Ruiz 2009-14
-- License : GPL
--
-- Maintainer : Alberto Ruiz
-- Stability : provisional
--
-- Random vectors and matrices.
--
-----------------------------------------------------------------------------
module Numeric.GSL.Random (
Seed,
RandDist(..),
randomVector,
gaussianSample,
uniformSample,
rand, randn
) where
import Numeric.GSL.Vector
import Numeric.LinearAlgebra.HMatrix hiding (
randomVector,
gaussianSample,
uniformSample,
Seed,
rand,
randn
)
import System.Random(randomIO)
type Seed = Int
-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
-- Gaussian distribution.
gaussianSample :: Seed
-> Int -- ^ number of rows
-> Vector Double -- ^ mean vector
-> Matrix Double -- ^ covariance matrix
-> Matrix Double -- ^ result
gaussianSample seed n med cov = m where
c = size med
meds = konst 1 n `outer` med
rs = reshape c $ randomVector seed Gaussian (c * n)
m = rs <> cholSH cov + meds
-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
-- uniform distribution.
uniformSample :: Seed
-> Int -- ^ number of rows
-> [(Double,Double)] -- ^ ranges for each column
-> Matrix Double -- ^ result
uniformSample seed n rgs = m where
(as,bs) = unzip rgs
a = fromList as
cs = zipWith subtract as bs
d = size a
dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
am = konst 1 n `outer` a
m = fromColumns (zipWith scale cs dat) + am
-- | pseudorandom matrix with uniform elements between 0 and 1
randm :: RandDist
-> Int -- ^ rows
-> Int -- ^ columns
-> IO (Matrix Double)
randm d r c = do
seed <- randomIO
return (reshape c $ randomVector seed d (r*c))
-- | pseudorandom matrix with uniform elements between 0 and 1
rand :: Int -> Int -> IO (Matrix Double)
rand = randm Uniform
{- | pseudorandom matrix with normal elements
>>> x <- randn 3 5
>>> disp 3 x
3x5
0.386 -1.141 0.491 -0.510 1.512
0.069 -0.919 1.022 -0.181 0.745
0.313 -0.670 -0.097 -1.575 -0.583
-}
randn :: Int -> Int -> IO (Matrix Double)
randn = randm Gaussian
|
