{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ExtendedDefaultRules #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}

-- | simulation to support testing of Mealy's using mwc-probability
module Data.Mealy.Simulate
  ( rvs,
    rvsp,
    create,
  )
where

import Control.Monad.Primitive (PrimState)
import NumHask.Prelude hiding (fold)
import System.Random.MWC.Probability

-- $setup
-- >>> :set -XFlexibleContexts
-- >>> import Data.Mealy
-- >>> gen <- create

-- | rvs creates a list of standard normal random variates.
--
-- >>> import Data.Mealy
-- >>> import Data.Mealy.Simulate
-- >>> gen <- create
-- >>> rvs gen 3
-- [1.8005943761746166e-2,0.36444481359059255,-1.2939898115295387]
--
-- >>> rs <- rvs gen 10000
-- >>> fold (ma 1) rs
-- 1.29805301109162e-2
--
-- >>> fold (std 1) rs
-- 1.0126527176272948
rvs :: Gen (PrimState IO) -> Int -> IO [Double]
rvs :: Gen (PrimState IO) -> Int -> IO [Double]
rvs Gen (PrimState IO)
gen Int
n = forall (m :: * -> *) a.
PrimMonad m =>
Int -> Prob m a -> Gen (PrimState m) -> m [a]
samples Int
n forall (m :: * -> *). PrimMonad m => Prob m Double
standardNormal Gen (PrimState IO)
gen

-- | rvsPair generates a list of correlated random variate tuples
--
-- >>> rvsp gen 3 0.8
-- [(1.8005943761746166e-2,7.074509906249835e-2),(0.36444481359059255,-0.7073208451897444),(-1.2939898115295387,-0.643930709405127)]
--
-- >>> rsp <- rvsp gen 10000 0.8
-- >>> fold (corr (ma 1) (std 1)) rsp
-- 0.8050112742986588
rvsp :: Gen (PrimState IO) -> Int -> Double -> IO [(Double, Double)]
rvsp :: Gen (PrimState IO) -> Int -> Double -> IO [(Double, Double)]
rvsp Gen (PrimState IO)
gen Int
n Double
c = do
  [Double]
s0 <- Gen (PrimState IO) -> Int -> IO [Double]
rvs Gen (PrimState IO)
gen Int
n
  [Double]
s1 <- Gen (PrimState IO) -> Int -> IO [Double]
rvs Gen (PrimState IO)
gen Int
n
  let s1' :: [Double]
s1' = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
y -> Double
c forall a. Multiplicative a => a -> a -> a
* Double
x forall a. Additive a => a -> a -> a
+ forall a. ExpField a => a -> a
sqrt (Double
1 forall a. Subtractive a => a -> a -> a
- Double
c forall a. Multiplicative a => a -> a -> a
* Double
c) forall a. Multiplicative a => a -> a -> a
* Double
y) [Double]
s0 [Double]
s1
  forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ forall a b. [a] -> [b] -> [(a, b)]
zip [Double]
s0 [Double]
s1'