----------------------------------------------------------------------------- -- | -- Module : Numeric.Random.Distribution.Poisson -- Copyright : (c) Matthew Donadio 2003 -- License : GPL -- -- Maintainer : m.p.donadio@ieee.org -- 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 testHead n lambda = mean \$ map fromIntegral \$ map (head . poisson lambda . randomRs (0,1) . mkStdGen) [1..n]