{-# LANGUAGE BangPatterns #-}
module QuantLib.Stochastic.Random
( BoxMuller
, mkNormalGen
, NormalGenerator (..)
, InverseNormal
, mkInverseNormal
) where
import QuantLib.Math.InverseNormal
import QuantLib.Stochastic.PureMT
class RandomGenerator a where
create :: IO a
next :: a -> (Double, a)
splitWithSeed :: Integer -> a -> (a, a)
instance RandomGenerator PureMT where
create :: IO PureMT
create = IO PureMT
newPureMT
next :: PureMT -> (Double, PureMT)
next = PureMT -> (Double, PureMT)
randomDouble
splitWithSeed :: Integer -> PureMT -> (PureMT, PureMT)
splitWithSeed = Integer -> PureMT -> (PureMT, PureMT)
splitMTwithSeed
data BoxMuller a = BoxMuller {
BoxMuller a -> Bool
bmFirst :: Bool,
BoxMuller a -> Double
bmSecondValue :: Double,
BoxMuller a -> a
bmRng :: a
}
mkNormalGen :: RandomGenerator a => IO (BoxMuller a)
mkNormalGen :: IO (BoxMuller a)
mkNormalGen = do
a
rng <- IO a
forall a. RandomGenerator a => IO a
create
BoxMuller a -> IO (BoxMuller a)
forall (m :: * -> *) a. Monad m => a -> m a
return (BoxMuller a -> IO (BoxMuller a))
-> BoxMuller a -> IO (BoxMuller a)
forall a b. (a -> b) -> a -> b
$! a -> BoxMuller a
forall a. RandomGenerator a => a -> BoxMuller a
createNormalGen a
rng
createNormalGen :: RandomGenerator a => a -> BoxMuller a
createNormalGen :: a -> BoxMuller a
createNormalGen a
r = BoxMuller :: forall a. Bool -> Double -> a -> BoxMuller a
BoxMuller {
bmFirst :: Bool
bmFirst = Bool
True,
bmSecondValue :: Double
bmSecondValue = Double
0.0,
bmRng :: a
bmRng = a
r
}
class NormalGenerator a where
ngGetNext :: a -> (Double, a)
ngMkNew :: a -> IO a
ngSplit :: a -> (a, a)
ngSplit = Integer -> a -> (a, a)
forall a. NormalGenerator a => Integer -> a -> (a, a)
ngSplitWithSeed Integer
1
ngSplitWithSeed :: Integer -> a -> (a, a)
instance RandomGenerator a => NormalGenerator (BoxMuller a) where
ngMkNew :: BoxMuller a -> IO (BoxMuller a)
ngMkNew BoxMuller a
_ = IO (BoxMuller a)
forall a. RandomGenerator a => IO (BoxMuller a)
mkNormalGen
ngGetNext :: BoxMuller a -> (Double, BoxMuller a)
ngGetNext = BoxMuller a -> (Double, BoxMuller a)
forall a. RandomGenerator a => BoxMuller a -> (Double, BoxMuller a)
boxMullerGetNext
ngSplitWithSeed :: Integer -> BoxMuller a -> (BoxMuller a, BoxMuller a)
ngSplitWithSeed Integer
seed BoxMuller a
x = (BoxMuller a
x { bmRng :: a
bmRng = a
rng1 }, BoxMuller a
x { bmRng :: a
bmRng = a
rng2 })
where
(a
rng1, a
rng2) = Integer -> a -> (a, a)
forall a. RandomGenerator a => Integer -> a -> (a, a)
splitWithSeed Integer
seed (BoxMuller a -> a
forall a. BoxMuller a -> a
bmRng BoxMuller a
x)
boxMullerGetNext :: RandomGenerator a => BoxMuller a -> (Double, BoxMuller a)
boxMullerGetNext :: BoxMuller a -> (Double, BoxMuller a)
boxMullerGetNext (BoxMuller Bool
True Double
_ a
rng) = (Double
s1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ratio, BoxMuller :: forall a. Bool -> Double -> a -> BoxMuller a
BoxMuller {
bmFirst :: Bool
bmFirst = Bool
False,
bmSecondValue :: Double
bmSecondValue = Double
s2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
ratio,
bmRng :: a
bmRng = a
g2
})
where
(Double
x1, a
g1) = a -> (Double, a)
forall a. RandomGenerator a => a -> (Double, a)
next a
rng
(Double
x2, a
g2) = a -> (Double, a)
forall a. RandomGenerator a => a -> (Double, a)
next a
g1
(Double
r, Double
s1, Double
s2) = (Double, Double, Double)
getRs
ratio :: Double
ratio = Double -> Double
boxMullerRatio Double
r
getRs :: (Double, Double, Double)
getRs =
let
as1 :: Double
as1 = Double
2.0Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1.0
as2 :: Double
as2 = Double
2.0Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1.0
ar :: Double
ar = Double
s1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
s1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
s2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
s2
in
if Double
rDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>=Double
1.0 Bool -> Bool -> Bool
|| Double
rDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<=Double
0.0 then (Double, Double, Double)
getRs else (Double
ar, Double
as1, Double
as2)
boxMullerGetNext (BoxMuller Bool
False !Double
s !a
r) = (Double
s, Bool -> Double -> a -> BoxMuller a
forall a. Bool -> Double -> a -> BoxMuller a
BoxMuller Bool
True Double
s a
r)
{-# ANN boxMullerRatio "NoHerbie" #-}
boxMullerRatio :: Double -> Double
boxMullerRatio :: Double -> Double
boxMullerRatio Double
r = Double -> Double
forall a. Floating a => a -> a
sqrt (-Double
2.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
r Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
r)
newtype InverseNormal a = InverseNormal a
mkInverseNormal :: RandomGenerator a => IO (InverseNormal a)
mkInverseNormal :: IO (InverseNormal a)
mkInverseNormal = do
a
rng <- IO a
forall a. RandomGenerator a => IO a
create
InverseNormal a -> IO (InverseNormal a)
forall (m :: * -> *) a. Monad m => a -> m a
return (InverseNormal a -> IO (InverseNormal a))
-> InverseNormal a -> IO (InverseNormal a)
forall a b. (a -> b) -> a -> b
$! a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
rng
instance RandomGenerator a => NormalGenerator (InverseNormal a) where
ngMkNew :: InverseNormal a -> IO (InverseNormal a)
ngMkNew InverseNormal a
_ = IO (InverseNormal a)
forall a. RandomGenerator a => IO (InverseNormal a)
mkInverseNormal
ngGetNext :: InverseNormal a -> (Double, InverseNormal a)
ngGetNext (InverseNormal a
rng) = (Double -> Double
inverseNormal Double
x, a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
newRng)
where (Double
x, a
newRng) = a -> (Double, a)
forall a. RandomGenerator a => a -> (Double, a)
next a
rng
ngSplitWithSeed :: Integer -> InverseNormal a -> (InverseNormal a, InverseNormal a)
ngSplitWithSeed Integer
seed (InverseNormal a
x) = (a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
x1, a -> InverseNormal a
forall a. a -> InverseNormal a
InverseNormal a
x2)
where
(a
x1, a
x2) = Integer -> a -> (a, a)
forall a. RandomGenerator a => Integer -> a -> (a, a)
splitWithSeed Integer
seed a
x