bellrand :: (RandomGen g, Floating a, Random a) => a -> a -> g -> (a, a, g)
-- Produces two randomly generated numbers from a normal distribution with
-- the given mean and standard deviation
bellrand mean stddev g =
let (u, g') = randomR (-1, 1) g
(v, g'') = randomR (-1, 1) g'
r = u * u + v * v
f = sqrt (negate 2 * log r / r)
in if 0 < r && r < 1 then (mean + stddev * u * f, mean + stddev * v * f, g'')
else bellrand mean stddev g''
bellrands :: (RandomGen g, Floating a, Random a) => a -> a -> g -> [a]
-- Like `bellrand`, but generates an infinite list
bellrands mean stddev g = concatMap (\(a, b, _) -> [a, b]) $ tail
$ iterate (\(_, _, h) -> bellrand mean stddev h) (0, 0, g)