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)
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_clt :: Int
-> (Double,Double)
-> [Double]
-> [Double]
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_bm :: (Double,Double)
-> [Double]
-> [Double]
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)
normal_ar :: (Double,Double)
-> [Double]
-> [Double]
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"
normal_r :: (Double,Double)
-> [Double]
-> [Double]
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
normal' [a]
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"normal_r: infinite list of random variables expected"