```-- Algorithms taken from Dodge and Jerse's Computer Music: Synthesis,
-- Composition, and Performance, Chapter 11.

module System.Random.Distributions (
-- * Random Distributions
linear, exponential, bilExp, gaussian, cauchy, poisson, frequency

-- * Utility Functions
, rands

) where

import System.Random

{- | Given a random number generator, generates a linearly distributed
random variable between 0 and 1.  Returns the random value together
with a new random number generator.  The probability density function
is given by

> f(x) = 2(1-x)  0 <= x <= 1
>      = 0       otherwise

-}
linear :: (RandomGen g, Floating a, Random a, Ord a) => g -> (a,g)
linear g0 =
let (r1, g1) = randomR (0, 1) g0
(r2, g2) = randomR (0, 1) g1
in (min r1 r2, g2)

{- | Takes a random number generator and produces another one
that avoids generating the given number.
-}
avoid :: (Random a, Eq a, RandomGen g) => a -> (g -> (a,g)) -> g -> (a,g)
avoid x f g = if r == x then avoid x f g' else (r,g')
where (r,g') = f g

{- | Generates an exponentially distributed random variable given a
generating a small number.  The mean of the distribution is
1/lambda.  The range of the generated number is [0,inf] although
the chance of getting a very large number is very small.

The probability density function is given by

> f(x) = lambda e^(-lambda * x)
-}
exponential :: (RandomGen g, Floating a, Random a, Eq a) =>
a  -- ^ horizontal spread of the function.
-> g  -- ^ a random number generator.
-> (a,g)
exponential lambda g0 = (-log r1 / lambda, g1)
where (r1, g1) = avoid 0 random g0

{- | Generates a random number with a bilateral exponential distribution.
Similar to exponential, but the mean of the distribution is 0 and
50% of the results fall between (-1/lambda, 1/lambda).

-}
bilExp :: (Floating a, Ord a, Random a, RandomGen g) =>
a    -- ^ horizontal spread of the function.
-> g    -- ^ a random number generator.
-> (a,g)
bilExp lambda g0 =
let (r', g1) = avoid 0 random g0
r = 2 * r'
u = if r > 1 then 2 - r else r
in (signum (1 - r) * log u / lambda, g1)

{- | Generates a random number with a Gaussian distribution.
-}
gaussian :: (Floating a, Random a, RandomGen g) =>
a     -- ^ standard deviation.
-> a     -- ^ mean.
-> g     -- ^ a random number generator.
-> (a,g)
gaussian stddev center g0 =
let n = 12
s = sum \$ take n \$ randoms g0
in (stddev * (s - fromIntegral n / 2) + center, fst (split g0))

{- | Generates a Cauchy-distributed random variable.
The distribution is symmetric with a mean of 0.

-}
cauchy :: (Floating a, Random a, RandomGen g, Eq a) =>
a   -- ^ alpha (density).
-> g   -- ^ a random number generator.
-> (a,g)
cauchy density g0 = (density * tan (u * pi), g1)
where (u, g1) = avoid 0.5 random g0

{- | Generates a Poisson-distributed random variable.
The given parameter lambda is the mean of the distribution.
If lambda is an integer, the probability that the result j=lambda-1
will be as great as that of j=lambda.  The Poisson distribution
is discrete. The returned value will be a non-negative
integer.

-}
poisson :: (Num t, Ord a, Floating a, RandomGen g, Random a) =>
a -> g -> (t, g)
poisson lambda g0 = (k 0 us, g1)
where v   = exp (-lambda)
us  = scanl1 (*) (randoms g0)
g1  = fst (split g0)
k n (u:us)
| u >= v     = k (n+1) us
| otherwise  = n
k _ [] = error "System.Random.Distributions.poisson: randoms did not return an infinite list"

{- | Given a list of weight-value pairs, generates a value randomly picked
from the list, weighting the probability of choosing each value by the
weight given.

-}
frequency :: (Floating w, Ord w, Random w, RandomGen g)
=> [(w, a)] -> g -> (a,g)
frequency xs g0 = (pick r xs, g1)
where (r, g1) = randomR (0, tot) g0
tot = sum (map fst xs)
pick n ((w,a):xs)
| n <= w    = a
| otherwise = pick (n-w) xs
pick _ [] = error "System.Random.Distributions.frequency: The impossible happened"

{- | Given a function generating a random number variable and a random
number generator, produces an infinite list of random values
generated from the given function.

-}
rands :: (RandomGen g, Random a) =>
(g -> (a,g)) -> g -> [a]
rands f g = x : rands f g' where (x,g') = f g

```