{-# LANGUAGE CPP #-}
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
gaussianSample :: Seed
-> Int
-> Vector Double
-> Herm Double
-> Matrix Double
gaussianSample :: Seed -> Seed -> Vector Double -> Herm Double -> Matrix Double
gaussianSample seed :: Seed
seed n :: Seed
n med :: Vector Double
med cov :: 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 -> Seed -> Vector Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst 1 Seed
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 = Seed -> Vector Double -> Matrix Double
forall t. Storable t => Seed -> Vector t -> Matrix t
reshape Seed
c (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Seed -> RandDist -> Seed -> Vector Double
randomVector Seed
seed RandDist
Gaussian (Seed
c Seed -> Seed -> Seed
forall a. Num a => a -> a -> a
* Seed
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
uniformSample :: Seed
-> Int
-> [(Double,Double)]
-> Matrix Double
uniformSample :: Seed -> Seed -> [(Double, Double)] -> Matrix Double
uniformSample seed :: Seed
seed n :: Seed
n rgs :: [(Double, Double)]
rgs = Matrix Double
m where
(as :: [Double]
as,bs :: [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
$ Seed -> Vector Double -> Matrix Double
forall t. Storable t => Seed -> Vector t -> Matrix t
reshape Seed
n (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Seed -> RandDist -> Seed -> Vector Double
randomVector Seed
seed RandDist
Uniform (Seed
nSeed -> Seed -> Seed
forall a. Num a => a -> a -> a
*Seed
d)
am :: Matrix Double
am = Double -> Seed -> Vector Double
forall e d (c :: * -> *). Konst e d c => e -> d -> c e
konst 1 Seed
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
randm :: RandDist
-> Int
-> Int
-> IO (Matrix Double)
randm :: RandDist -> Seed -> Seed -> IO (Matrix Double)
randm d :: RandDist
d r :: Seed
r c :: Seed
c = do
Seed
seed <- IO Seed
forall a. Random a => IO a
randomIO
Matrix Double -> IO (Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Seed -> Vector Double -> Matrix Double
forall t. Storable t => Seed -> Vector t -> Matrix t
reshape Seed
c (Vector Double -> Matrix Double) -> Vector Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$ Seed -> RandDist -> Seed -> Vector Double
randomVector Seed
seed RandDist
d (Seed
rSeed -> Seed -> Seed
forall a. Num a => a -> a -> a
*Seed
c))
rand :: Int -> Int -> IO (Matrix Double)
rand :: Seed -> Seed -> IO (Matrix Double)
rand = RandDist -> Seed -> Seed -> IO (Matrix Double)
randm RandDist
Uniform
randn :: Int -> Int -> IO (Matrix Double)
randn :: Seed -> Seed -> IO (Matrix Double)
randn = RandDist -> Seed -> Seed -> IO (Matrix Double)
randm RandDist
Gaussian