-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Random.Distribution.Normal
-- Copyright   :  (c) Matthew Donadio 2003
-- License     :  GPL
--
-- Maintainer  :  m.p.donadio@ieee.org
-- Stability   :  experimental
-- Portability :  portable
--
-- Module for transforming a list of uniform random variables into a
-- list of normal random variables.
--
-----------------------------------------------------------------------------

-- TODO: The speedup from Ross for the A-R method

-- TODO: Marsaglia's ziggurat method

-- TODO: Leva' method

-- TODO: Ahrens-Dieter method

module Numeric.Random.Distribution.Normal (normal_clt, normal_bm,
                                           normal_ar, normal_r) where

import DSP.Basic (interleave, uninterleave, norm2sqr, toMaybe)
import Data.Maybe (mapMaybe)

-- * Functions

-- adjust takes a unit normal random variable and sets the mean and
-- variance to whatever is needed.

adjust :: (Double,Double) -> Double -> Double
adjust :: (Double, Double) -> Double -> Double
adjust (Double
mu,Double
sigma) Double
x = Double
mu forall a. Num a => a -> a -> a
+ Double
sigma forall a. Num a => a -> a -> a
* Double
x

-- | Normal random variables via the Central Limit Theorm (not explicity
-- given, but see Ross)
--
-- If mu=0 and sigma=1, then this will generate numbers in the range
-- [-n/2,n/2]

normal_clt :: Int             -- ^ Number of uniforms to sum
           -> (Double,Double) -- ^ (mu,sigma)
           -> [Double]        -- ^ U
           -> [Double]        -- ^ X

normal_clt :: Int -> (Double, Double) -> [Double] -> [Double]
normal_clt Int
n (Double, Double)
muSigma [Double]
u = forall a b. (a -> b) -> [a] -> [b]
map ((Double, Double) -> Double -> Double
adjust (Double, Double)
muSigma) forall a b. (a -> b) -> a -> b
$ [Double] -> [Double]
normal' [Double]
u
    where normal' :: [Double] -> [Double]
normal' [Double]
us = Double
var_adj forall a. Num a => a -> a -> a
* ((forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [a]
take Int
n [Double]
us) forall a. Num a => a -> a -> a
- Double
mean_adj) forall a. a -> [a] -> [a]
: ([Double] -> [Double]
normal' forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [a]
drop Int
n [Double]
us)
          var_adj :: Double
var_adj  = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ Double
12 forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
          mean_adj :: Double
mean_adj = forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n forall a. Fractional a => a -> a -> a
/ Double
2

-- | Normal random variables via the Box-Mueller Polar Method (Ross, pp
-- 450--452)
--
-- If mu=0 and sigma=1, then this will generate numbers in the range
-- [-8.57,8.57] assuing that the uniform RNG is really giving full
-- precision for doubles.

normal_bm :: (Double,Double) -- ^ (mu,sigma)
          -> [Double]        -- ^ U
          -> [Double]        -- ^ X

normal_bm :: (Double, Double) -> [Double] -> [Double]
normal_bm (Double, Double)
muSigma =
   forall a b. (a -> b) -> [a] -> [b]
map ((Double, Double) -> Double -> Double
adjust (Double, Double)
muSigma) forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall a. [a] -> [a] -> [a]
interleave forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. [(a, b)] -> ([a], [b])
unzip forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe forall a. (Floating a, Ord a) => (a, a) -> Maybe (a, a)
normalDist forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry forall a b. [a] -> [b] -> [(a, b)]
zip forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. [a] -> ([a], [a])
uninterleave forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a -> b) -> [a] -> [b]
map (\Double
u -> Double
2forall a. Num a => a -> a -> a
*Double
uforall a. Num a => a -> a -> a
-Double
1)

normalDist :: (Floating a, Ord a) => (a,a) -> Maybe (a,a)
normalDist :: forall a. (Floating a, Ord a) => (a, a) -> Maybe (a, a)
normalDist z :: (a, a)
z@(a
x,a
y) =
   let norm2 :: a
norm2 = forall a. Num a => (a, a) -> a
norm2sqr (a, a)
z
       p :: a
p = forall a. Floating a => a -> a
sqrt (-a
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log a
norm2) forall a. Fractional a => a -> a -> a
/ a
norm2
   in  forall a. Bool -> a -> Maybe a
toMaybe (a
norm2forall a. Ord a => a -> a -> Bool
<=a
1) (a
pforall a. Num a => a -> a -> a
*a
x, a
pforall a. Num a => a -> a -> a
*a
y)


-- | Acceptance-Rejection Method (Ross, pp 448--450)
--
-- If mu=0 and sigma=1, then this will generate numbers in the range
-- [-36.74,36.74] assuming that the uniform RNG is really giving full
-- precision for doubles.

normal_ar :: (Double,Double) -- ^ (mu,sigma)
          -> [Double]        -- ^ U
          -> [Double]        -- ^ X

normal_ar :: (Double, Double) -> [Double] -> [Double]
normal_ar (Double, Double)
muSigma [Double]
u = forall a b. (a -> b) -> [a] -> [b]
map ((Double, Double) -> Double -> Double
adjust (Double, Double)
muSigma) forall a b. (a -> b) -> a -> b
$ forall {a}. (Floating a, Ord a) => [a] -> [a]
normal' [Double]
u
    where normal' :: [a] -> [a]
normal' (a
u1:a
u2:a
u3:[a]
us) | a
y forall a. Ord a => a -> a -> Bool
> a
0     = a
z forall a. a -> [a] -> [a]
: [a] -> [a]
normal' [a]
us
                                | Bool
otherwise = [a] -> [a]
normal' (a
u3forall a. a -> [a] -> [a]
:[a]
us)
              where y1 :: a
y1 = -forall a. Floating a => a -> a
log a
u1
                    y2 :: a
y2 = -forall a. Floating a => a -> a
log a
u2
                    y :: a
y  = a
y2 forall a. Num a => a -> a -> a
- (a
y1 forall a. Num a => a -> a -> a
- a
1)forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int) forall a. Fractional a => a -> a -> a
/ a
2
                    z :: a
z = if a
u3 forall a. Ord a => a -> a -> Bool
<= a
0.5 then a
y1 else -a
y1
          normal' [a]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"normal_ar: infinite list of random variables expected"

-- | Ratio Method (Kinderman-Monahan) (Knuth, v2, 2ed, pp 125--127)
--
-- If mu=0 and sigma=1, then this will generate numbers in the range
-- [-1e15,1e15] (?) assuming that the uniform RNG is really giving full
-- precision for doubles.

normal_r :: (Double,Double) -- ^ (mu,sigma)
         -> [Double]        -- ^ U
         -> [Double]        -- ^ X

normal_r :: (Double, Double) -> [Double] -> [Double]
normal_r (Double, Double)
muSigma = forall a b. (a -> b) -> [a] -> [b]
map ((Double, Double) -> Double -> Double
adjust (Double, Double)
muSigma) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall {a}. (Ord a, Floating a) => [a] -> [a]
normal'
    where normal' :: [a] -> [a]
normal' (a
u:a
v:[a]
us) | a
xforall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int) forall a. Ord a => a -> a -> Bool
<= -a
4 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log a
u = a
x forall a. a -> [a] -> [a]
: [a] -> [a]
normal' [a]
us
                           | Bool
otherwise                = [a] -> [a]
normal' [a]
us
              where x :: a
x = a
a forall a. Num a => a -> a -> a
* (a
v forall a. Num a => a -> a -> a
- a
0.5) forall a. Fractional a => a -> a -> a
/ a
u
                    a :: a
a = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ a
8 forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp a
1 -- 1.71552776992141359295
          normal' [a]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"normal_r: infinite list of random variables expected"