random-variates 0.1.0.0 → 0.1.1.0

raw patch · 21 files changed

+848/−367 lines, 21 filesdep +HUnitdep +directorydep +erfdep ~basedep ~containersPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: HUnit, directory, erf, mtl, random-variates

Dependency ranges changed: base, containers

API changes (from Hackage documentation)

Files

• random-variates.cabal +36/−16
• src/Helpers.hs +0/−56
• src/Stochastic/Analysis.hs +58/−0
• src/Stochastic/Bernoulli.hs +0/−20
• src/Stochastic/Binomial.hs +0/−20
• src/Stochastic/Distribution.hs +0/−26
• src/Stochastic/Distribution/Continuous.hs +18/−0
• src/Stochastic/Distribution/Discrete.hs +16/−0
• src/Stochastic/Distributions.hs +54/−71
• src/Stochastic/Distributions/Continuous.hs +90/−0
• src/Stochastic/Distributions/Discrete.hs +171/−0
• src/Stochastic/Exponential.hs +0/−19
• src/Stochastic/Generator.hs +45/−0
• src/Stochastic/Geometric.hs +0/−23
• src/Stochastic/Normal.hs +0/−25
• src/Stochastic/Poisson.hs +0/−24
• src/Stochastic/Tools.hs +182/−0
• src/Stochastic/Uniform.hs +0/−22
• src/Stochastic/ZipF.hs +0/−45
• tests/unit.hs +149/−0
• tests/vis.hs +29/−0

@@ -2,9 +2,9 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                random-variates-version:             0.1.0.0+version:             0.1.1.0 synopsis:            "Uniform RNG => Non-Uniform RNGs"-description:         "Collection of transforms uniform random number generators (RNGs) into any of a dozen common RNGs. Each presenting several common interfaces"   +description:         "Collection of transforms uniform random number generators (RNGs) into any of a dozen common RNGs. Each presenting several common interfaces. Additionally Empirical distributions can be sampled from and tested (chi-squared) against theoretical distributions."    license:             MIT license-file:        LICENSE author:              Keynan James Pratt <keynan.pratt@gmail.com>@@ -21,23 +21,43 @@   location:            git@bitbucket.org:kpratt/random-variate.git      library-  exposed-modules:     Stochastic.Distribution-                       Stochastic.Distributions-                       Stochastic.Uniform+  exposed-modules:+                  Stochastic.Generator+                  Stochastic.Distribution.Continuous+                  Stochastic.Distribution.Discrete+                  Stochastic.Distributions+                  Stochastic.Distributions.Continuous+                  Stochastic.Distributions.Discrete+                  Stochastic.Analysis+                  Stochastic.Tools -  other-modules:       Stochastic.Bernoulli-                       Stochastic.Binomial-                       Stochastic.Exponential-                       Stochastic.Geometric-                       Stochastic.Normal-                       Stochastic.Poisson-                       Stochastic.ZipF-                       Helpers+  other-modules:    build-depends:       -                       base >=4.6 && <5.0-                       , reinterpret-cast >= 0.1.0-                       , containers >= 0.5.0.0+                         base >=4.6 && <5.0+                       , containers >= 0.5.7.0                        , lens >=4.13                        , random >=1.1+                       , reinterpret-cast >= 0.1.0+                       , mtl >= 2.2+                       , erf >= 2.0   hs-source-dirs:      src   default-language:    Haskell2010++Test-Suite vis+  type:                exitcode-stdio-1.0+  main-is:             tests/vis.hs+  build-depends:       base+                     , directory       >= 1.2.0.1+                     , random-variates >=0.1+  default-language:    Haskell2010+++Test-Suite units+  type:                exitcode-stdio-1.0+  main-is:             tests/unit.hs+  build-depends:       HUnit >= 1.3+                     , base+                     , random-variates >=0.1+  default-language:    Haskell2010++  

@@ -1,56 +0,0 @@-module Helpers(maybeHead,-               headOrElse,-               statefully,-               mapTuple,-               statefullyTakeWhile,-               histogram) where--import Data.Maybe-import qualified Data.Map as Map--maybeHead :: [a] -> Maybe a-maybeHead []     = Nothing-maybeHead (x:xs) = Just x--headOrElse :: a -> [a] -> a-headOrElse d ls = fromMaybe d (maybeHead ls)--statefully :: (g -> (a, g)) -> Int -> g -> ([a], g)-statefully f n g0 = case n of-    0 -> ([], g0)-    x -> (r:rest, g2)-      where-        (r, g1) = f g0-        (rest, g2) = statefully f (n-1) g1--statefullyTakeWhile :: (g -> (a, g)) ->-                       ([a] -> Bool) ->-                       g ->-                       ([a], g)-statefullyTakeWhile f p g0 = r ([], g0)-  where-    r (list, g1)-      | p list    = (list, g1)-      | otherwise = mapTuple (\l -> l:list) (id) (f g1)---mapTuple :: (a -> b) -> (c -> d) -> (a, c) -> (b, d)    -mapTuple f g (x, y) = (f x, g y)--histogram :: Double -> [Double] -> [Int]-histogram precision ls = map lookup0 [0..limit]-  where-    limit = truncate $ (1.0 / precision)-    divs = map (\x -> x / precision) ls-    ints = map (truncate) divs-    m   = foldl-          (\b a ->-            Map.insert a ((fromMaybe 0 (Map.lookup a b))+1) b)-          Map.empty-          ints-    lookup0 = (\x -> fromMaybe 0 $ Map.lookup x m)-    -toDbl = fromInteger . toInteger---

@@ -0,0 +1,58 @@+module Stochastic.Analysis(chiSquaredTest, discreteChiSquaredTest) where++import qualified Stochastic.Distributions.Continuous as C+import qualified Stochastic.Distributions.Discrete as D++import Data.List(sort)+import Stochastic.Distributions+import Stochastic.Tools++chiCritD :: (C.ContinuousDistribution g)+            => g+            -> Empirical+            -> Double+            -> Double+chiCritD theory emp =+  chiCrit (degreesOfFreedom emp) (C.degreesOfFreedom theory)++-- k = number of intervals in the empirical histogram+-- s = number of parameters to the theoretical distribution+chiCrit :: Int -> Int -> Double -> Double+chiCrit k s alpha = C.cdf' chi (1-alpha)+  where+    df = k - s - 1+    chi = C.ChiSquared df (stdBase 42)++-- empirical distributions should be+-- given as the second argument+chiSquaredTest :: (C.ContinuousDistribution g)+                  => g+                  -> Empirical+                  -> [Double]+                  -> Double+chiSquaredTest c d sampleAt = if (isNaN final) then (error $ "frog") else final+  where+    final = sum $ fmap f sampleAt+    f x = let o =  cdf' d $ C.cdf c x in+           let e = x in+           let ret = ((e-o)**2) / e in+           if (isNaN ret)+           then error $ (show e) ++ " " ++ (show o) ++ " " ++ (show ret) ++ " " ++ (show (0/0))+           else ret+++discreteChiSquaredTest :: (D.DiscreteDistribution g)+                          => g+                          -> Empirical+                          -> [Int]+                          -> Double+discreteChiSquaredTest c d sampleAt = sum $ fmap f sampleAt+  where+    f :: Int -> Double+    f x = let e = D.cdf' c $ cdf d $ toDbl x in+           let o = x in+           toDbl ((e-o)^2) / toDbl e+    toDbl = fromInteger . toInteger+++

@@ -1,20 +0,0 @@-module Stochastic.Bernoulli(Bernoulli, mkBernoulli) where--import Stochastic.Distribution-import Stochastic.Uniform-import Helpers--data Bernoulli = Bernoulli Double Uniform--mkBernoulli :: Uniform -> Double -> Bernoulli-mkBernoulli base p = Bernoulli p base--instance DiscreteDistribution Bernoulli where-  randIntIn (a, b) (Bernoulli p g0) =-    mapTuple-    (\x -> a + ((floor $ x + (1-p)) * (b-a)) )-    (Bernoulli p)-    (randDouble g0)-  randInt g0 = randIntIn (0, 1) g0--toDbl = fromInteger . toInteger

@@ -1,20 +0,0 @@-module Stochastic.Binomial(Binomial, mkBinomial) where--import Stochastic.Distribution-import Stochastic.Bernoulli-import Helpers--data Binomial = Binomial Int Bernoulli--mkBinomial :: Bernoulli -> Int -> Binomial-mkBinomial base n = Binomial n base--instance DiscreteDistribution Binomial where-  randIntIn (a, b) (Binomial n g0) =-    mapTuple-    (\xs -> (a - 1) + (sum xs))-    (Binomial n)-    (randInts (b - a + 1) g0)-  randInt (Binomial n g0) = randIntIn (0, n) (Binomial n g0)--toDbl = fromInteger . toInteger

@@ -1,26 +0,0 @@-{-# LANGUAGE DefaultSignatures #-}--module Stochastic.Distribution(ContinuousDistribution(..), DiscreteDistribution(..)) where--import Helpers(histogram, statefully)--class DiscreteDistribution g where-  randInt :: g -> (Int, g)-  randInt g = randIntIn (0, maxBound::Int) g-  randInts :: Int -> g -> ([Int], g)-  randInts n g0 = statefully (randInt) n g0-  randIntIn :: (Int, Int) -> g -> (Int, g)-  default randIntIn :: (ContinuousDistribution g) =>-                       (Int, Int) -> g -> (Int, g)-  randIntIn (a, b) g0 = ((ceiling (toDbl (b - a + 1) * d)) + (a-1), g1)-    where-      (d, g1) = randDouble g0-  --class ContinuousDistribution g where-  randDouble :: g -> (Double, g)-  randDoubles :: Int -> g -> ([Double], g)-  randDoubles n g0 = statefully (randDouble) n g0---toDbl = fromInteger . toInteger

@@ -0,0 +1,18 @@+module Stochastic.Distribution.Continuous where++import Stochastic.Tools++class ContinuousDistribution g where+  rand  :: g -> (Double, g)+  +  rands :: Int -> g -> ([Double], g)+  rands n g0 = statefully (rand) n g0++  cdf  :: g -> Double -> Double+  cdf' :: g -> Double -> Double+  +  pdf  :: g -> Double -> Double -> Double+  pdf g a b = (cdf g b) - (cdf g a)++  degreesOfFreedom :: g -> Int+  

@@ -0,0 +1,16 @@+module Stochastic.Distribution.Discrete where++import Stochastic.Tools++class DiscreteDistribution g where+  rand :: g -> (Int, g)++  rands :: Int -> g -> ([Int], g)+  rands n g0 = statefully (rand) n g0++  cdf  :: g -> Int -> Double+  +  cdf' :: g -> Double -> Int+  +  pmf  :: g -> Int -> Double+

@@ -1,89 +1,72 @@+{-# LANGUAGE TypeFamilies #-} module Stochastic.Distributions(-  ContinuousDistribution(..)-  ,DiscreteDistribution(..)-  ,Distributions(..)-  ,mkDistributions-  ,stdDistributions-  ,liftD-  ,liftC-  ,liftDN-  ,liftCN---  ,plot+  UniformBase(rDouble)+  ,stdBase+  ,Empirical(..)+  ,mkEmpirical   ) where -import Helpers(histogram, statefully, mapTuple)-import Stochastic.Distribution+import System.Random+import Control.Monad.State.Lazy+import Stochastic.Tools -import qualified Stochastic.Uniform     as Uni-import qualified Stochastic.ZipF        as Zipf-import qualified Stochastic.Geometric   as Geo-import qualified Stochastic.Exponential as Exp-import qualified Stochastic.Poisson     as Poi-import qualified Stochastic.Normal      as Nor-import qualified Stochastic.Bernoulli   as Ber-import qualified Stochastic.Binomial    as Bin+data UniformBase = UniformBase {+  rDouble :: (Double, UniformBase)+} -data Distributions = Distributions {-  mkUniform    :: Int                     -> Uni.Uniform-  ,mkExp       :: Int -> Double           -> Exp.Exponential-  ,mkNormal    :: Int -> Double -> Double -> Nor.Normal -  ,mkZipF      :: Int -> Int    -> Double -> Zipf.ZipF-  ,mkGeometric :: Int -> Double           -> Geo.Geometric-  ,mkPoisson   :: Int -> Double           -> Poi.Poisson-  ,mkBernoulli :: Int -> Double           -> Ber.Bernoulli-  ,mkBinomial  :: Int -> Double -> Int    -> Bin.Binomial-}+stdGen2Uni gen = UniformBase {+  rDouble = mapTuple+            (id)+            (stdGen2Uni)+            (randomR (0,1) gen)+  } -stdDistributions = mkDistributions Uni.stdUniform+stdBase s = stdGen2Uni (mkStdGen s) -mkDistributions uniform = Distributions {-  mkUniform     = uniform-  , mkZipF      = apiMkZipF      (uniform)-  , mkGeometric = apiMkGeometric (uniform)-  , mkExp       = apiMkExp       (uniform)-  , mkPoisson   = apiMkPoisson   (uniform)-  , mkNormal    = apiMkNormal    (uniform)-  , mkBernoulli = apiMkBernoulli (uniform)-  , mkBinomial  = apiMkBinomial  (uniform)-} -apiMkZipF      mkUni seed n slope = Zipf.mkZipF      (mkUni seed) n slope-apiMkGeometric mkUni seed p       = Geo.mkGeometric  (mkUni seed) p-apiMkExp       mkUni seed y       = Exp.mkExp        (mkUni seed) y-apiMkPoisson   mkUni seed y       = Poi.mkPoisson $  apiMkExp (mkUni) seed y-apiMkNormal    mkUni seed m d     = Nor.mkNormal     (mkUni seed) m d-apiMkBernoulli mkUni seed d       = Ber.mkBernoulli  (mkUni seed) d-apiMkBinomial  mkUni seed d n     = Bin.mkBinomial   (apiMkBernoulli (mkUni) seed d) n+data Empirical = Empirical {+  degreesOfFreedom :: Int,+  cdf  :: Double -> Double,+  cdf' :: Double -> Double+  } -liftD :: (DiscreteDistribution g) => (Int, Int) -> (Int -> a) -> (g -> (a, g))-liftD range f g0 = ((f r), g1)+empiricalCDF' hist x = u - ((c - x)/s)   where-    (r, g1) = (randIntIn range g0)+    interval = head $ filter (\y -> cum_frequence y > x) $ hist+    u = upper_bound interval+    s = slope interval+    c = cum_frequence interval -liftDN :: (DiscreteDistribution g) => [(Int, Int)] -> ([Int] -> a) -> (g -> (a, g))-liftDN ranges f g0 = mapTuple (f) (id) (h ranges g0)+  +                 +empiricalCDF hist x+  | x <  (lower_bound $ head hist) = 0+  | x >= (upper_bound $ head $ reverse hist) = 1+  | otherwise =+    f $ filter (\y -> lower_bound y <= x  && upper_bound y >= x) hist   where-    h [] g1 = ([], g1)-    h (r:rs) g1 = let (s, g2) = randIntIn r g1 in mapTuple (\x -> s:x) (id) (h rs g2)+    f [] =+      error $ "x "++(show x)++" not in histogram range\n" +++      (foldr (\h str -> (show (lower_bound h, upper_bound h, frequency h)) ++ "\n" ++ str) "" hist)+    f (y:ys) =+      let part = (x - (lower_bound y)) in+      let step = part * slope y        in+      (cum_frequence y) + (step * (rel_frequence y))+     -liftC :: (ContinuousDistribution g) => (Double -> a) -> (g -> (a, g))-liftC f g0 = ((f r), g1)++mkEmpirical :: [Double] -> Empirical+mkEmpirical samples = Empirical {+  degreesOfFreedom = length h,+  cdf  = empiricalCDF h,+  cdf' = empiricalCDF' h+  }   where-    (r, g1) = (randDouble g0)+    h = fIHistogram samples+    count :: Double+    count = fromInteger . toInteger $ length samples -liftCN :: (ContinuousDistribution g) => Int -> ([Double] -> a) -> (g -> (a, g))-liftCN n f g0 = mapTuple (f) (id) (randDoubles n g0) -plot :: DiscreteDistribution g => g -> Int -> Int -> [Int]-plot g0 interval samples = []-{--plot g n samples = map (truncate . (+0.5) . (*100))-  $ map-  (\x ->-    (toDbl x)/(toDbl samples))-  $ histogram (1.0 / (toDbl n))-  (fst (randInts g samples))--}-toDbl = fromInteger . toInteger+ 

@@ -0,0 +1,90 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeFamilies #-}+module Stochastic.Distributions.Continuous(+  mkUniform+  ,mkExp+  ,mkNormal+  ,mkEmpirical+  ,Dist(..)+  ,ContinuousDistribution(..)+  ) where++import Data.Maybe+import Control.Monad.State.Lazy+--import Stochastic.Analysis+import Stochastic.Generator++import Stochastic.Distributions(UniformBase, rDouble)+import qualified Stochastic.Distributions as B(cdf, mkEmpirical, Empirical)+import Stochastic.Distribution.Continuous+import Stochastic.Tools+import Data.Number.Erf++instance ContinuousDistribution UniformBase where+  rand uni = rDouble uni+  cdf  _ x = x+  cdf' _ p = p+  degreesOfFreedom _ = 0+++instance Generator Dist where+  type (From Dist) = Double+  nextG = state $ \ g0 -> rand g0++instance Generator UniformBase where+  type (From UniformBase) = Double+  nextG = state $ \ g0 -> rDouble g0+++data Dist =+  Uniform UniformBase+  | Exponential Double UniformBase+  | Normal Double Double (Maybe Double) UniformBase+  | ChiSquared Int UniformBase+  | Empirical B.Empirical UniformBase+    -- empirical points, lo, [(point, mass)]++mkEmpirical :: UniformBase -> [Double] -> Dist+mkEmpirical base samples = Empirical (B.mkEmpirical samples) base++mkExp :: UniformBase -> Double -> Dist+mkExp base y = Exponential y base+mkNormal :: UniformBase -> Double -> Double -> Dist+mkNormal uni mean dev = Normal mean dev Nothing uni+mkUniform :: UniformBase -> Dist+mkUniform uni = Uniform uni++instance ContinuousDistribution Dist where+  rand (Uniform uni) = mapTuple (id) (Uniform) (rand uni)+  rand (Exponential y u) =+    mapTuple (\x -> (-1.0/y) * (log $ x)) (Exponential y) (rand u)+  rand (Normal mean dev m uni) = f m+    where+      f (Just x) = (x, (Normal mean dev Nothing uni'))+      f Nothing  = (y, (Normal mean dev (Just z) uni'))+      ([u1, u2], uni') = rands 2 uni+      from_u g = mean + dev * (sqrt (-2 * (log u1))) * ( g (2 * pi * u2) )+      y = from_u (sin)+      z = from_u (cos)++  cdf  (Uniform _) x = x+  cdf  (Exponential y _) x = 1 - (1 / (exp (y*x)))+  cdf  (Normal u s _ _) x =+    0.5 * (1 + (erf ((x-u)/(s * (sqrt 2))) ))+  cdf (ChiSquared k _) x = (1/(gamma (kd/2))) * lig+    where+      kd = fromInteger $ toInteger k+      lig = lower_incomplete_gamma (kd /2) (x/2)+  cdf (Empirical b _) x = B.cdf b x+         ++  cdf' (Uniform _) p = p+  cdf' (Exponential y _) p = -(log (1-p)) / y+  cdf' (Normal u s _ _) p =+    u + (s * (sqrt 2) * (inverf(2*p-1)))++  degreesOfFreedom (Uniform _) = 0+  degreesOfFreedom (Exponential _ _) = 1+  degreesOfFreedom (Normal _ _ _ _) = 2++

@@ -0,0 +1,171 @@+{-# LANGUAGE TypeFamilies #-}+module Stochastic.Distributions.Discrete(+  mkBinomial+  ,mkBernoulli+  ,mkPoisson+  ,mkZipF+  ,mkGeometric+  ,Stochastic.Distributions.Discrete.Dist(..)+  ,DiscreteDistribution(..)+  ) where++import Data.Maybe+import Control.Monad.State.Lazy+import Stochastic.Tools+import Stochastic.Generator(Generator(..), foldWhile)+import Stochastic.Distributions+import Stochastic.Distribution.Discrete+import qualified Stochastic.Distributions.Continuous as C++instance Generator Dist where+  type (From Dist) = Int+  nextG = state $ \ g0 -> rand g0++data Dist =+  Uniform Int Int UniformBase+  | Poisson C.Dist+  | Geometric Double UniformBase+  | Bernoulli Double UniformBase+  | Binomial Int Double DiscreteCache UniformBase+  | ZipF Int Double DiscreteCache UniformBase++mkBinomial :: UniformBase -> Double -> Int -> Dist+mkBinomial base p n = Binomial n p cache base+  where+    nd = toDbl n+    cache :: DiscreteCache+    cache = foldr (\(w,x,y,z) l -> (w,y,z):l) [] (create n)+    create :: Int -> [(Int, Double, Double, Double)]+    create 0 = let pmfk = (nd * (log (1 - p)) )+               in [(0, pmfk, exp pmfk, exp pmfk)]+    create k = (k, lpmfk, pmfk, cdfk) : sub+      where+        sub = create (k-1)+        (j, lpmfj, pmfj, cdfj) = head sub+        pmfk = exp lpmfk+        kd = toDbl k+        lpmfk = lpmfj ++                (log p) - (log (1- p)) -+                (log (kd)) + (log (nd-kd+1))+        cdfk = pmfk + cdfj+        +-- k, pmf k, cdf k+-- do in log space+type DiscreteCache = [(Int, Double, Double)]++mkUniform :: UniformBase -> Int -> Int -> Dist+mkUniform base a b = Uniform a b base++mkBernoulli :: UniformBase -> Double -> Dist+mkBernoulli base p = Bernoulli p base++mkPoisson :: UniformBase -> Double -> Dist+mkPoisson base y = Poisson (C.mkExp base y)++mkGeometric :: UniformBase -> Double -> Dist+mkGeometric base p = Geometric p base++mkZipF :: UniformBase -> Int -> Double -> Dist+mkZipF base n slope = ZipF n slope cache base+  where+    hns = sum $ take n $ harmonics slope+    cache :: DiscreteCache+    cache = create n+    create :: Int -> DiscreteCache+    create 1 = let pmfk = 1 / hns in [(1,pmfk,pmfk)]+    create k =+      (k, (1/(((toDbl k)**slope) * hns)), cdfj + pmfj) : sub+      where+        sub             = create (k-1)+        (j, pmfj, cdfj) = head sub++instance DiscreteDistribution Dist where+  rand (Binomial n p cache g0) =+    mapTuple+    (\u -> +      length $ filter (pred u) cache)+    (Binomial n p cache)+    (C.rand g0)+    where pred u (k, pmf, cdf) = cdf < u+  rand (Geometric p g0) =+    mapTuple+    (\u -> ceiling $ (log u) / (log (1-p)))+    (Geometric p)+    (C.rand g0)+  rand (Poisson g0@(C.Exponential y u)) =+    mapTuple+    (\x -> length x)+    (Poisson)+    (runState (foldWhile (+) 0.0 (<1.0)) g0)+  rand (Bernoulli p g0) =+    mapTuple+    (\x -> if (x >= p) then 1 else 0)+    (Bernoulli p)+    (C.rand g0)+  rand (ZipF n slope cache u0) = +    mapTuple+    (\u ->+      length $ filter (pred u) cache)+    (ZipF n slope cache)+    (C.rand u0)+    where pred u (k, pmf, cdf) = cdf < u+  rand (Uniform a b g0) =+    mapTuple+    (\x -> truncate (toDbl (b - a) * x + toDbl a))+    (Uniform a b)+    (C.rand g0)++  cdf (Poisson (C.Exponential y _)) x =+    (1/(exp y)) * (sum [ (y ** (toDbl i)) / (fromInteger $ fac i) | i <- [0..x]])+  cdf (Geometric p _) x = 1 - (1-p)^x+  cdf (Bernoulli p _) x+    | x < 0     = 0+    | x >= 1    = 1+    | otherwise = p+  cdf (Binomial n p cache _) x = r+    where+      (_, _, r) = head $ filter (\(w,_,_) -> (w==x)) cache+  cdf (ZipF n s cache _) x = r+    where+      (_, _, r) =+        head $ filter (\(k, _, _) -> x == k) cache +  cdf (Uniform a b _) x = toDbl (x-a) / toDbl (b-a)++  cdf' g@(Poisson (C.Exponential y _)) x =+    sum . fst $ runState (fold) [1..]+    where+      reduce :: Int -> Double -> Double+      reduce = (\y p -> (pmf g y) + p)+      fold :: State [Int] [From [Int]]+      fold = foldWhile (reduce) 0 (<x)+  cdf' (Geometric p _) x =+    ceiling $ (log (1-x)) / (log (1-p))+  cdf' (Bernoulli p _) x+    | x > p     = 1+    | otherwise = 0+  cdf' (Binomial n p cache _) x = r+    where+      (r, _, _) =+        head $ reverse $ filter (\(_, _, z) -> z > x) cache+  cdf' (ZipF n s cache _) x = r+    where+      (r, _, _) =+        head $ reverse $ filter (\(_, _, z) -> z > x) cache+  cdf' (Uniform a b _) x = truncate $ toDbl (b-a) * x + toDbl a+  +  pmf (Poisson (C.Exponential y _)) x =+    (y^x) / ( exp y * (fromInteger $ fac x) )+  pmf (Geometric p _) x = 1 - (1-p)^(x-1)+  pmf (Bernoulli p _) 0 = 1-p+  pmf (Bernoulli p _) 1 = p+  pmf (Binomial n p cache _) x = r+    where+      (_, r, _) = head $ filter (\(w,_,_) -> (w==x)) cache+  pmf (ZipF n s cache _) k = r+    where+      (_, r, _) = head $ filter (\(w,_,_) -> (w==k)) cache+  pmf (Uniform a b _) x = toDbl x/toDbl (b-a)+++toDbl = fromInteger . toInteger+

@@ -1,19 +0,0 @@-module Stochastic.Exponential(Exponential, mkExp) where--import Stochastic.Distribution-import Stochastic.Uniform-import Data.Word-import Data.ReinterpretCast-import Helpers--data Exponential = Exponential Double Uniform--viaWord :: Double -> Int-viaWord w =  fromInteger $ toInteger $ doubleToWord w--mkExp :: Uniform -> Double -> Exponential-mkExp base y = Exponential y base--instance ContinuousDistribution Exponential where-  randDouble (Exponential y u) =-    mapTuple (\x -> (-1.0/y) * (log $ x)) (Exponential y) (randDouble u)

@@ -0,0 +1,45 @@+{-# LANGUAGE TypeFamilies #-}+module Stochastic.Generator where++import Control.Monad.State.Lazy++class Generator g where+  type From g+  nextG   :: State g (From g)+  nextN  :: Int -> State g [(From g)]+  nextN 0 = state $ \g0 -> ([], g0)+  nextN n =+    do+      x  <- nextG+      xs <- nextN (n-1)+      return (x:xs)++instance Generator [a] where+  type From [a] = a+  nextG = state $ \ g0 -> (head g0, tail g0)++foldWhile :: Generator g+         => (From g -> a -> a)+         -> a+         -> (a -> Bool)+         -> State g [From g]+foldWhile f z p =+  do+   x  <- nextG+   let y = f x z+   xs <- if (p y)+         then foldWhile f y p +         else return []+   return (x:xs)++while :: Generator g => ((From g) -> Bool) -> State g [From g]+while p =+  do+   x  <- nextG+   xs <- if (p x)+         then while p+         else return []+   return (x:xs)+++

@@ -1,23 +0,0 @@-module Stochastic.Geometric (mkGeometric, Geometric) where--import Stochastic.Distribution-import Stochastic.Uniform-import Helpers--data Geometric = Geometric Double Uniform--mkGeometric :: Uniform -> Double -> Geometric-mkGeometric base p = Geometric p base--instance DiscreteDistribution Geometric where-  randIntIn (a, b) (Geometric p g0) = mapTuple-               (\u -> trunc $ invert u)-               (Geometric p)-               (randDouble g0)-    where-      invert u = ceiling $ (log u) / (log (1-p))-      trunc x-        | x > (b-a) = b-        | otherwise = a+x--toDbl = fromInteger . toInteger

@@ -1,25 +0,0 @@-module Stochastic.Normal(mkNormal, Normal) where--import Data.Maybe-import Stochastic.Distribution-import Stochastic.Uniform-import Helpers--data Normal = Normal Double Double (Maybe Double) Uniform--mkNormal :: Uniform -> Double -> Double -> Normal-mkNormal uni mean dev = Normal mean dev Nothing uni--toDbl :: Int -> Double-toDbl = fromInteger . toInteger--instance ContinuousDistribution Normal where-  randDouble (Normal mean dev m uni) = f m-    where-      f (Just x) = (x, (Normal mean dev Nothing uni))-      f Nothing  = (y, (Normal mean dev (Just z) uni))-      ([u1, u2], uni') = randDoubles 2 uni-      from_u g = mean + dev * (sqrt (-2 * (log u1))) * ( g (2 * pi * u2) )-      y = from_u (sin)-      z = from_u (cos)-    

@@ -1,24 +0,0 @@-module Stochastic.Poisson(mkPoisson, Poisson) where--import Stochastic.Distribution-import Stochastic.Exponential-import Helpers--data Poisson = Poisson Exponential--mkPoisson :: Exponential -> Poisson-mkPoisson base = Poisson base--toDbl :: Int -> Double-toDbl = fromInteger . toInteger--instance DiscreteDistribution Poisson where-  randIntIn (a, b) (Poisson g0) = mapTuple-                                  (\x -> min (x+a-1) b)-                                  (Poisson)-                                  (f 0 0 g0)-    where-      f x s g1-        | s > 1     = (x-1, g1)-        | otherwise = f (x+1) (s+y) g2-          where (y, g2) = (randDouble g1)

@@ -0,0 +1,182 @@+{-# LANGUAGE TypeFamilies #-}+module Stochastic.Tools where++import Control.Monad.State.Lazy+import Data.Maybe+import qualified Data.Map as Map+import Data.List(sort)++maybeHead :: [a] -> Maybe a+maybeHead []     = Nothing+maybeHead (x:xs) = Just x++headOrElse :: a -> [a] -> a+headOrElse d ls = fromMaybe d (maybeHead ls)++mapTuple :: (a -> b) -> (c -> d) -> (a, c) -> (b, d)    +mapTuple f g (x, y) = (f x, g y)++histogram :: Double -> [Double] -> [Int]+histogram precision ls = map lookup0 [0..limit]+  where+    limit = truncate $ (1.0 / precision)+    divs = map (\x -> x / precision) ls+    ints = map (truncate) divs+    m   = foldl+          (\b a ->+            Map.insert a ((fromMaybe 0 (Map.lookup a b))+1) b)+          Map.empty+          ints+    lookup0 = (\x -> fromMaybe 0 $ Map.lookup x m)++comb :: Int -> Int -> Double+comb n k =+  (fromInteger $ foldr (*) 1 ls) / (fromInteger $ fac (n-k))+  where+    ls = [(toInteger k)..(toInteger n)]+++gamma :: Double -> Double+gamma = stirlingsApprox++lower_incomplete_gamma :: Double -> Double -> Double+lower_incomplete_gamma = estimate_lower_gamma 10++estimate_lower_gamma n s x = (x**s)*(gamma s) * (exp (-x)) * approx+  where+    approx = sum $ fmap term [0..n]+    term k = (x**k) / (gamma (s + k + 1))++stirlingsApprox n +  | n < 0.5   = stirlingsApprox (1-n)+  | otherwise = (sqrt (2*pi*n)) * (exp $ n * ((log n) - 1))++fac :: Int -> Integer+fac n = head $ drop (n-1) factorial++factorial :: [Integer]+factorial = 1 : 1 : (tail (zipWith (*) (factorial) [1..]))++fib :: Int -> Integer+fib n = head $ drop (n-1) fibinacci++fibinacci :: [Integer]+fibinacci = 1 : 1 : (zipWith (+) fibinacci (tail fibinacci))++harmonics :: Double -> [Double]+harmonics s = (h 1)+  where+    h n = (1.0/(n**s)) : h (n+1)++statefully :: (g -> (a, g)) -> Int -> g -> ([a], g)+statefully f n g0 = case n of+    0 -> ([], g0)+    x -> (r:rest, g2)+      where+        (r, g1) = f g0+        (rest, g2) = statefully f (n-1) g1++statefullyTakeWhile :: (g -> (a, g)) ->+                       ([a] -> Bool) ->+                       g ->+                       ([a], g)+statefullyTakeWhile f p g0 = r ([], g0)+  where+    r (list, g1)+      | p list    = (list, g1)+      | otherwise = mapTuple (\l -> l:list) (id) (f g1)++type Histogram = [Datagram]+data Datagram = Datagram {+  lower_bound :: Double,+  upper_bound :: Double,+  frequency :: Int,+  rel_frequence :: Double,+  cum_frequence :: Double,+  slope :: Double+  } deriving (Show, Eq)++++--peice wise linear +fIHistogram :: [Double] -> Histogram+fIHistogram list =+  datagramFromRaw len iSize $ bin (lowerOf iSize) list+  where+    toDbl = fromInteger . toInteger+    len = length list+    iSize :: Double+    iSize = (maxf list -+             minf list) /+            iCount+    iCount = sqrt (toDbl $ len)++maxf (x:xs) = foldl (max) x xs+minf (x:xs) = foldl (min) x xs++datagramFromRaw :: Int -> Double -> [(Double, Int)] -> [Datagram]+datagramFromRaw count iSize = accMap f z+  where+    f acc (x, c)  = Datagram {+      lower_bound = x,+      upper_bound = x + iSize,+      frequency = c,+      rel_frequence = rf,+      cum_frequence = rf + (cum_frequence acc),+      slope = rf / iSize+      }+      where rf = (toDbl c / toDbl count)+    z (x, c) = Datagram {+      lower_bound = x,+      upper_bound = x + iSize,+      frequency = c,+      rel_frequence = rf,+      cum_frequence = rf,+      slope = rf / iSize+      }+      where rf = (toDbl c / toDbl count)+    toDbl = fromInteger . toInteger++accMap :: (b -> a -> b) -> (a -> b) -> [a] -> [b]+accMap f z ls = g (z $ head ls) (tail ls)+  where+    g prev []     = []+    g prev (x:xs) = new : (g new xs)+      where new = (f prev x)+{-+Takes the set +w : x : y : z : []+and produces the set +(f w x) : (f x y) : (f y z) : []+-}    +biFold :: (a -> a -> b) -> [a] -> [b]+biFold f = g+  where+    g []  = []+    g [x] = []+    g (x:y:ys) = (f x y) : (g (y:ys))++dependentMap :: (a -> a -> b) -> (a -> b) -> [a] -> [b]+dependentMap f z xs = (z $ head xs) : (biFold f xs)+    +    +bin :: (Double -> Double) -> [Double] -> [(Double, Int)]+bin lowerOf ls = group $ map (lowerOf) ls++lowerOf :: Double -> Double -> Double+lowerOf iSize v = iSize * (fromInteger $ truncate (v / iSize))+++    +group :: (Ord a) => [a] -> [(a, Int)]+group = groupSeq . sort+    +groupSeq :: (Eq a) => [a] -> [(a, Int)]+groupSeq list = foldr g [] list+  where+    g x [] = [(x, 1)]+    g x ((y,c):ys)+      | x == y    = (y, c+1):ys+      | otherwise = (x, 1):(y,c):ys++

@@ -1,22 +0,0 @@-module Stochastic.Uniform (Uniform(..)-                          , stdUniform) where--import Helpers-import System.Random-import Stochastic.Distribution---instance ContinuousDistribution Uniform where-  randDouble uni = rDouble uni --data Uniform = Uniform {-  rDouble :: (Double, Uniform)-}--stdGen2Uni gen = Uniform {-  rDouble = mapTuple (id) (stdGen2Uni) (randomR (0,1) gen)-  }--stdUniform s = stdGen2Uni (mkStdGen s)--

@@ -1,45 +0,0 @@-module Stochastic.ZipF (mkZipF, ZipF) where--import Helpers-import Stochastic.Uniform-import Data.Maybe-import Stochastic.Distribution--data ZipF = ZipF Int Double Uniform--mkZipF :: Uniform -> Int -> Double -> ZipF-mkZipF base n slope = ZipF n slope base---- index, number, cumlative-harmonics :: Double -> [(Double, Double, Double)]-harmonics s = (h 1 0)-  where-    h n acc = (n, v, v+acc) : h (n+1) (v+acc)-      where-        v = (1.0/(n**s))--toDbl = fromInteger . toInteger--f n s d = h2 $ h1 hs-  where-    hs      = harmonics s-    mx      = _3 $ head $ drop (n-1) $ take n $ hs-    h1 xs    = headOrElse (toDbl n, 0, 0) $-              (dropWhile (\(i, v, c) -> c < (d * mx))) (take (n) xs)-    h2 (x, _, _) = truncate x-  -g = f 10 1--instance DiscreteDistribution ZipF where-  randIntIn (a, b) (ZipF n slope u0) = (f n slope d, ZipF n slope u1)-    where-      (d, u1) = randDouble u0--_1 :: (a, b, c) -> a-_1 (x, y, z) = x--_3 :: (a, b, c) -> c-_3 (x, y, z) = z---

@@ -0,0 +1,149 @@+import Stochastic.Distributions+import Stochastic.Analysis+import Stochastic.Tools+import qualified Stochastic.Distributions.Continuous as C+import qualified Stochastic.Distributions.Discrete as D+import Test.HUnit+import Data.Monoid+import Control.Monad+--import Utils+import Control.Monad (unless)++tests = TestList [+  TestLabel "Binomial pmf" bpmf1+  ,TestLabel "ZipF pmf" zipFpmf1+  ,TestLabel "Normal cdf 0" normalCDF0+  ,TestLabel "Normal cdf 1" normalCDF1+  ,TestLabel "Normal cdf 2" normalCDF2+  ,TestLabel "Normal cdf 3" normalCDF3+  ,TestLabel "Normal ChiSquared" chinorm+  ,TestLabel "Exponential ChiSquared" chiexp+  ,TestLabel "Uniform ChiSquared" chiuni+  ,TestLabel "Poisson ChiSquared" chiPoisson+  ,TestLabel "Test of grouping" groupTest+--  ,TestLabel "Test histogram" binTest+  ]+++--+-- Discrete Tests+--+bpmf1 = TestCase (assertWithinDelta "binomial's pmf must sum to 1" (1e-15) e a)+  where+    b = D.mkBinomial (stdBase 42) 0.4 10+    e = 1+    a = sum [ D.pmf b k | k <- [0..10] ]++zipFpmf1 = TestCase (assertWithinDelta "ZipF's pmf must sum to 1" (1e-15) e a)+  where+    b = D.mkZipF (stdBase 42) 10 1+    e = 1+    a = sum [ D.pmf b k | k <- [1..10] ]++--+-- Continuous Tests+--++stdNormal = C.mkNormal (stdBase 42) 0 1++normalCDF0 = TestCase (assertWithinDelta+                       "Normal CDF at 0 standard deviations"+                       (1e-15) 0.5    $ C.cdf stdNormal 0)+normalCDF1 = TestCase (assertWithinDelta+                       "Normal CDF at 1 standard deviations"+                       (1e-4)  0.8413 $ C.cdf stdNormal 1)+normalCDF2 = TestCase (assertWithinDelta+                       "Normal CDF at 2 standard deviations"+                       (1e-4)  0.9772 $ C.cdf stdNormal 2)+normalCDF3 = TestCase (assertWithinDelta+                       "Normal CDF at 0 standard deviations"+                       (1e-4)  0.9987 $ C.cdf stdNormal 3)++--+-- Stochastic Tests+--++chiTestRandom :: C.Dist+                 -> Double+chiTestRandom g = chiSquaredTest g (mkEmpirical samples) bounds+  where+    (samples, _) = C.rands 1000 g+    hist = fIHistogram samples+    bounds = fmap (lower_bound) hist++chinorm = TestCase (+  assertWithinDelta+  "Sample of the normal distribution passes the ChiSquaredTest with HIGH confidence"+  (2e-1) 0 (chiTestRandom stdNormal))++chiuni = TestCase (+  assertWithinDelta+  "Sample of the uniform distribution passes the ChiSquaredTest with HIGH confidence"+  (1e-1) 0 (chiTestRandom $ C.mkUniform $ stdBase 42))++chiexp = TestCase (+  assertWithinDelta+  "Sample of the exponential distribution passes the ChiSquaredTest with HIGH confidence"+  (2e-1) 0 (chiTestRandom $ C.mkExp (stdBase 42) 1))++groupTest = TestCase+            (assertEqual+             "groups where non optimal"+             ([(1,2), (2,3), (3,1)])+             (group [2,1,1,2,2,3])+            )++chiPoisson = TestCase (+  assertWithinDelta+  "Sample of the poisson distribution passes the ChiSquaredTest with HIGH confidence"+  (2e-1) 0 (chiTestDiscreteRandom $ D.mkPoisson (stdBase 42) 1))++chiTestDiscreteRandom :: D.Dist+                      -> Double+chiTestDiscreteRandom g = discreteChiSquaredTest g (mkEmpirical samples) bounds+  where+    samples = fmap (toDbl) $ fst $ D.rands 1000 g+    hist = fIHistogram samples+    bounds = fmap (fromDbl . lower_bound) hist+    toDbl = fromInteger . toInteger+    fromDbl = fromInteger . truncate++            +{-+binTest = TestCase+          (assertEqual+           "bins had holes"+           ([])+           (fmap (\x -> (lower_bound x, upper_bound x, frequency x))+            $ fIHistogram samples)+          )+  where+    (samples, _) = C.rands 1000 stdNormal+-}++{-+ let norm = mkNormal (stdBase 42) 0 1 in                                let (samples,_) = rands 1000 norm in                                   let emp = mkEmpirical samples in                                       let hist = fIHistogram samples in                                      let bounds = fmap (lower_bound) hist in                                let chi = chiSquaredTest norm emp bounds in                            let imp = fmap (\x -> (lower_bound x, frequency x, C.cdf norm (lower_bound x), C.cdf emp (lower_bound x)) ) hist in mapM (putStrLn.show) imp-}++--+-- Custom Assertions+--++assertWithinDelta ::   String  -- ^ The message prefix+               -> Double  -- ^ The maximum difference between expected and actual+               -> Double  -- ^ The expected value+               -> Double  -- ^ The actual value+               -> Assertion+assertWithinDelta preface delta expected actual = +  unless passed (assertFailure msg)+  where+    passed = (abs (expected - actual) < delta)+    msg = (if null preface then "" else preface ++ "\n") +++          "expected: " ++ show expected ++ "\n but got: " ++ show actual+  +++main :: IO ()+main =+  do+    counts <- runTestTT tests+    return ()

@@ -0,0 +1,29 @@+import Stochastic.Generator+import Stochastic.Distribution.Continuous+import Stochastic.Distributions+import Stochastic.Distributions.Continuous+import System.Directory+import System.IO++dumpSamples :: (String, [Double]) -> IO ()+dumpSamples  (gen_name, samples) = do+  _ <- createDirectoryIfMissing True "dist/samples"+  _ <- withFile+       ("dist/samples/" ++ gen_name ++ ".dat")+       WriteMode+       (\handle -> mapM_+                   (hPrint handle)+                   samples+       )+  return ()++cds seed sample_size = [+  ("uniform",          fst $ rands sample_size $ mkUniform (stdBase seed)    )+  ,("exponential",     fst $ rands sample_size $ mkExp     (stdBase seed)   1)+  ,("normal",          fst $ rands sample_size $ mkNormal  (stdBase seed) 0 1)+  ]++main :: IO ()+main =+  let gens = cds 42 100000 in+  mapM_ (dumpSamples) gens