-- | Yet another simple module for implementing statistics stuff.

module Statistics (normal, normals, stdnormal, stdnormals, module System.Random) where
import System.Random

stdnormal :: StdGen -> (Double,StdGen)
stdnormal = normal 0 1

stdnormals :: StdGen -> [Double]
stdnormals = normals 0 1

normal :: Double -> Double -> StdGen -> (Double,StdGen)
normal mu sigma = \g -> let (x,g') = randomR (0,1) g in (invcumnorm mu sigma x,g') 

normals :: Double -> Double -> StdGen -> [Double]
normals mu sigma = \g -> let (x,g') = normal mu sigma g in x : normals mu sigma g'

lognormal :: Double -> Double -> StdGen -> Double
lognormal mu sigma = undefined

-- support

invcumnorm mu sigma z = mu + search (-limit*sigma) (limit*sigma)
    where search a b = let c = (a+b)/2
                           cn = cumnorm 0 sigma c
                       in if abs (z - cn) < 10*epsilon || abs (a-b) < epsilon then c
                            else if cn > z then search a c
                                 else search c b

cumstdnorm :: Double -> Double
cumstdnorm x = 0.5*(1+erf (x/sqrt 2))

cumnorm :: Double -> Double -> Double -> Double
cumnorm mu sigma x = 0.5*(1+erf((x-mu)/(sigma*sqrt 2)))

-- taylor expansion, see wikipedia "error function".  Tested within the range (-limit..limit)
erf :: Double -> Double
erf x | x>limit         = 1
      | x< negate limit = 0
      | otherwise   = (2/sqrt pi)*sum (reverse $ takeWhile ((>=epsilon).abs) [ ((-1)**n*x**(2*n+1)) / (fac n*(2*n+1)) | n <- [0..]])

epsilon = 0.0000000001
limit = 4.4 :: Double
fac x = product [2..x]