```-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Random.Distribution.Poisson
--
-- Stability   :  experimental
-- Portability :  portable
--
-- UNTESTED
--
-- Module for transforming a list of uniform random variables
-- into a list of Poisson random variables.
--
-- Reference: Ross
--     Donald E. Knuth (1969). Seminumerical Algorithms, The Art of Computer Programming, Volume 2
--
----------------------------------------------------------------------------

module Numeric.Random.Distribution.Poisson (poisson, test, testHead) where

import Numeric.Statistics.Moment (mean)

import Data.List (mapAccumL)
import System.Random (randomRs, mkStdGen)

-- * Functions

{- |
Generates a list of poisson random variables from a list of uniforms.
-}

poisson :: Double    -- ^ lambda - expectation value, should be non-negative.
-> [Double]  -- ^ uniformly distributed values from the interval [0,1]
-> [Int]     -- ^ Poisson distributed outputs

poisson lambda (u:us) =
let e = exp (-lambda)
{- 'group' cannot replace segmentAfter here,
because it merges adjacent False values. -}
in  map (length . tail) . segmentAfter not . snd \$
mapAccumL
(\p ui ->
let b = p >= e
in  (if b then p*ui else ui, b))
u us
poisson _ [] =
error "poisson: list of uniformly distributed values must not be empty"

{- |
Split after every element that satisfies the predicate.

A candidate for a Utility module.
-}
segmentAfter :: (a -> Bool) -> [a] -> [[a]]
segmentAfter p =
foldr (\ x ~yt@(y:ys) -> if p x then [x]:yt else (x:y):ys) [[]]

{-
The expectation value,
and thus the mean of a sequence of Poisson distributed values,
approximates lambda.
-}

test :: Int -> Double -> Double
test n lambda =
mean \$ map fromIntegral \$
take n \$ poisson lambda \$
randomRs (0,1) \$ mkStdGen 1

{-
Only test the leading number of several Poisson lists.
-}
testHead :: Int -> Double -> Double