summaryrefslogtreecommitdiff
path: root/packages/gsl/src/Numeric/GSL/Random.hs
blob: f1f49e5bc41293c8065e372de789372c1310e663 (plain)
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
-----------------------------------------------------------------------------
-- |
-- 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(cholSH)
import Numeric.Container hiding (
    randomVector,
    gaussianSample,
    uniformSample
    )
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 = dim med
    meds = konst 1 n `outer` med
    rs = reshape c $ randomVector seed Gaussian (c * n)
    m = rs `mXm` cholSH cov `add` 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 = dim a
    dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
    am = konst 1 n `outer` a
    m = fromColumns (zipWith scale cs dat) `add` 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