{-# 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 :: 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
uniformSample :: Seed
-> Int
-> [(Double,Double)]
-> Matrix Double
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
randm :: RandDist
-> Int
-> Int
-> 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 (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))
rand :: Int -> Int -> IO (Matrix Double)
rand :: Int -> Int -> IO (Matrix Double)
rand = RandDist -> Int -> Int -> IO (Matrix Double)
randm RandDist
Uniform
randn :: Int -> Int -> IO (Matrix Double)
randn :: Int -> Int -> IO (Matrix Double)
randn = RandDist -> Int -> Int -> IO (Matrix Double)
randm RandDist
Gaussian