{-# LANGUAGE CPP #-}
-----------------------------------------------------------------------------
-- |
-- 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)
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

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
               -> Herm   Double -- ^ covariance matrix
               -> Matrix Double -- ^ result
gaussianSample :: Int -> Int -> Vector Double -> Herm Double -> Matrix Double
gaussianSample Int
seed Int
n Vector Double
med Herm Double
cov = Matrix Double
m where
    c :: IndexOf Vector
c = Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Vector Double
med
    meds :: Matrix Double
meds = Double -> Int -> Vector Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst Double
1 Int
n Vector Double -> Vector Double -> Matrix Double
forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
med
    rs :: Matrix Double
rs = Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
reshape Int
c (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Int -> RandDist -> Int -> Vector Double
randomVector Int
seed RandDist
Gaussian (Int
c Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n)
    m :: Matrix Double
m = Matrix Double
rs Matrix Double -> Matrix Double -> Matrix Double
forall t. Numeric t => Matrix t -> Matrix t -> Matrix t
<> Herm Double -> Matrix Double
forall t. Field t => Herm t -> Matrix t
chol Herm Double
cov Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ Matrix Double
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 :: Int -> Int -> [(Double, Double)] -> Matrix Double
uniformSample Int
seed Int
n [(Double, Double)]
rgs = Matrix Double
m where
    ([Double]
as,[Double]
bs) = [(Double, Double)] -> ([Double], [Double])
forall a b. [(a, b)] -> ([a], [b])
unzip [(Double, Double)]
rgs
    a :: Vector Double
a = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList [Double]
as
    cs :: [Double]
cs = (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
subtract [Double]
as [Double]
bs
    d :: IndexOf Vector
d = Vector Double -> IndexOf Vector
forall (c :: * -> *) t. Container c t => c t -> IndexOf c
size Vector Double
a
    dat :: [Vector Double]
dat = Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix Double -> [Vector Double])
-> Matrix Double -> [Vector Double]
forall a b. (a -> b) -> a -> b
$ Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
reshape Int
n (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Int -> RandDist -> Int -> Vector Double
randomVector Int
seed RandDist
Uniform (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d)
    am :: Matrix Double
am = Double -> Int -> Vector Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst Double
1 Int
n Vector Double -> Vector Double -> Matrix Double
forall t. Product t => Vector t -> Vector t -> Matrix t
`outer` Vector Double
a
    m :: Matrix Double
m = [Vector Double] -> Matrix Double
forall t. Element t => [Vector t] -> Matrix t
fromColumns ((Double -> Vector Double -> Vector Double)
-> [Double] -> [Vector Double] -> [Vector Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Vector Double -> Vector Double
forall t (c :: * -> *). Linear t c => t -> c t -> c t
scale [Double]
cs [Vector Double]
dat) Matrix Double -> Matrix Double -> Matrix Double
forall a. Num a => a -> a -> a
+ Matrix Double
am

-- | pseudorandom matrix with uniform elements between 0 and 1
randm :: RandDist
     -> Int -- ^ rows
     -> Int -- ^ columns
     -> IO (Matrix Double)
randm :: RandDist -> Int -> Int -> IO (Matrix Double)
randm RandDist
d Int
r Int
c = do
    Int
seed <- IO Int
forall a (m :: * -> *). (Random a, MonadIO m) => m a
randomIO
    Matrix Double -> IO (Matrix Double)
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
reshape Int
c (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Int -> RandDist -> Int -> Vector Double
randomVector Int
seed RandDist
d (Int
rInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
c))

-- | pseudorandom matrix with uniform elements between 0 and 1
rand :: Int -> Int -> IO (Matrix Double)
rand :: Int -> Int -> IO (Matrix Double)
rand = RandDist -> Int -> Int -> IO (Matrix Double)
randm RandDist
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 :: Int -> Int -> IO (Matrix Double)
randn = RandDist -> Int -> Int -> IO (Matrix Double)
randm RandDist
Gaussian