random-variates 0.1.3.0 → 0.1.4.0

raw patch · 8 files changed

+132/−77 lines, 8 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• random-variates.cabal +2/−1
• src/Stochastic/Analysis.hs +12/−9
• src/Stochastic/Distributions.hs +3/−11
• src/Stochastic/Distributions/Continuous.hs +44/−20
• src/Stochastic/Distributions/Discrete.hs +34/−25
• src/Stochastic/Generator.hs +2/−2
• src/Stochastic/Tools.hs +14/−6
• tests/unit.hs +21/−3

@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                random-variates-version:             0.1.3.0+version:             0.1.4.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. Additionally Empirical distributions can be sampled from and tested (chi-squared) against theoretical distributions."    license:             MIT@@ -68,6 +68,7 @@   build-depends:       HUnit >= 1.3                      , base                      , random-variates >=0.1+                     , random >=1.1   default-language:    Haskell2010     

@@ -33,12 +33,11 @@ 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+    f x =  let o = cdf d x in+           let e = C.cdf c x in+           if (e == 0)+           then ((e-o)**2)+           else ((e-o)**2) / e   discreteChiSquaredTest :: (D.DiscreteDistribution g)@@ -49,9 +48,13 @@ 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+    f 0 = 0+    f x =  let e = D.cdf c x in+           let o = cdf d $ toDbl x in+           if (e == 0)+           then ((e-o)**2)+           else ((e-o)**2) / e+    toDbl :: Int -> Double      toDbl = fromInteger . toInteger  

@@ -36,14 +36,10 @@   let words = [w1,w2,w3,w4,w5,w6,w7,w8] :: String   return $ runGet getWord64host $ LBS.fromStrict $ B.pack words -mk g = UniformBase {-  rDouble = mapTuple (id) (mk) (random g)-  }--stdBase :: Integer -> UniformBase-stdBase s = mk $ xorshift128plus s+stdBase :: Integer -> UniformRandom+stdBase s = xorshift128plus s -seededBase :: IO UniformBase+seededBase :: IO UniformRandom seededBase = do   word <- withBinaryFile "/dev/random/" ReadMode (readWord64)   let seed = toInteger word@@ -61,8 +57,6 @@     u = upper_bound interval     s = slope interval     c = cum_frequence interval--                     empiricalCDF hist x   | x <  (lower_bound $ head hist) = 0@@ -78,8 +72,6 @@       let step = part * slope y        in       (cum_frequence y) + (step * (rel_frequence y))     -- mkEmpirical :: [Double] -> Empirical mkEmpirical samples = Empirical {   degreesOfFreedom = length h,

@@ -1,3 +1,5 @@+{-# LANGUAGE ExistentialQuantification #-}+ module Stochastic.Distributions.Continuous(   mkUniform   ,mkExp@@ -12,46 +14,68 @@ --import Stochastic.Analysis import Stochastic.Generator -import Stochastic.Distributions(UniformBase, rDouble)+import Stochastic.Distributions(stdBase) import qualified Stochastic.Distributions as B(cdf, mkEmpirical, Empirical) import Stochastic.Distribution.Continuous import Stochastic.Tools import Data.Number.Erf+import System.Random -instance ContinuousDistribution UniformBase where-  rand uni = rDouble uni-  cdf  _ x = x-  cdf' _ p = p-  degreesOfFreedom _ = 0+data Dist = forall a . RandomGen a => Uniform a+          | forall a . RandomGen a => Exponential Double a+          | forall a . RandomGen a => Normal Double Double (Maybe Double) a+          | forall a . RandomGen a => ChiSquared Int a+          | forall a . RandomGen a => Empirical B.Empirical a+    -- empirical points, lo, [(point, mass)]+instance RandomGen Dist where+  next g@(Uniform uni) = mapTuple id Uniform $ next g+  next g@(Exponential y _) =+    let (x, g') = rand g in+    let (_, scale) = genRange g in+    let x' = truncate ((fromIntegral scale) * x) in+    if (x' < scale)+    then (x', g')+    else next g'+  next g@(Normal mean dev _ _) =+    let (x, g') = rand g in+    (truncate x, g') +  genRange g@(Uniform _) = (minBound :: Int, maxBound :: Int)+  genRange g@(Exponential _ _) = (0, maxBound `div` 4096 :: Int)+  genRange g@(Normal mean dev _ _) =+    let pm = (dev * 6.66) in+    (ceiling $ mean - pm, ceiling $ mean + pm) -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 :: forall a . RandomGen a => a -> [Double] -> Dist mkEmpirical base samples = Empirical (B.mkEmpirical samples) base -mkExp :: UniformBase -> Double -> Dist+mkExp :: forall a . RandomGen a => a ->  Double -> Dist mkExp base y = Exponential y base-mkNormal :: UniformBase -> Double -> Double -> Dist++mkNormal :: forall a . RandomGen a => a -> Double -> Double -> Dist mkNormal uni mean dev = Normal mean dev Nothing uni-mkUniform :: UniformBase -> Dist++mkUniform :: forall a . RandomGen a => a -> Dist mkUniform uni = Uniform uni +intWordDbl :: Int -> Double+intWordDbl x = fromRational $ toRational ((fromInteger $ toInteger x) :: Word)++randomN :: forall a . forall b . (RandomGen a, Random b) => Int -> a -> ([b], a)+randomN n = genTake (random) n+  + instance ContinuousDistribution Dist where-  rand (Uniform uni) = mapTuple (id) (Uniform) (rand uni)+  rand (Uniform uni) = mapTuple (id) (Uniform) (random uni)   rand (Exponential y u) =-    mapTuple (\x -> (-1.0/y) * (log $ x)) (Exponential y) (rand u)+    mapTuple ((\x -> -(log $ x) / y)) (Exponential y) (random 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+      (vs, uni') = randomN 2 uni+      [u1, u2] = map (id) vs       from_u g = mean + dev * (sqrt (-2 * (log u1))) * ( g (2 * pi * u2) )       y = from_u (sin)       z = from_u (cos)

@@ -1,3 +1,5 @@+{-# LANGUAGE ExistentialQuantification #-}+ module Stochastic.Distributions.Discrete(   mkBinomial   ,mkBernoulli@@ -14,17 +16,21 @@ import Stochastic.Generator(foldGenWhile) import Stochastic.Distributions import Stochastic.Distribution.Discrete+import Stochastic.Generator+import System.Random import qualified Stochastic.Distributions.Continuous as C -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+data Dist = forall a . RandomGen a => Uniform Int Int a+          | Poisson C.Dist+          | forall a . RandomGen a => Geometric Double a+          | forall a . RandomGen a => Bernoulli Double a+          | forall a . RandomGen a => Binomial Int Double DiscreteCache a+          | forall a . RandomGen a => ZipF Int Double DiscreteCache a -mkBinomial :: UniformBase -> Double -> Int -> Dist+instance RandomGen Dist where+  next g = rand g++mkBinomial :: forall a . RandomGen a => a -> Double -> Int -> Dist mkBinomial base p n = Binomial n p cache base   where     nd = toDbl n@@ -48,19 +54,19 @@ -- do in log space type DiscreteCache = [(Int, Double, Double)] -mkUniform :: UniformBase -> Int -> Int -> Dist+mkUniform :: forall a . RandomGen a => a -> Int -> Int -> Dist mkUniform base a b = Uniform a b base -mkBernoulli :: UniformBase -> Double -> Dist+mkBernoulli :: forall a . RandomGen a => a -> Double -> Dist mkBernoulli base p = Bernoulli p base -mkPoisson :: UniformBase -> Double -> Dist+mkPoisson :: forall a . RandomGen a => a -> Double -> Dist mkPoisson base y = Poisson (C.mkExp base y) -mkGeometric :: UniformBase -> Double -> Dist+mkGeometric :: forall a . RandomGen a => a -> Double -> Dist mkGeometric base p = Geometric p base -mkZipF :: UniformBase -> Int -> Double -> Dist+mkZipF :: forall a . RandomGen a => a -> Int -> Double -> Dist mkZipF base n slope = ZipF n slope cache base   where     hns = sum $ take n $ harmonics slope@@ -74,41 +80,44 @@         sub             = create (k-1)         (j, pmfj, cdfj) = head sub +nextN :: forall a . RandomGen a => Int -> a -> ([Int], a)+nextN n = genTake (next) n+ instance DiscreteDistribution Dist where   rand (Binomial n p cache g0) =     mapTuple     (\u -> -      length $ filter (pred u) cache)+      length $ filter (pred u) cache  )     (Binomial n p cache)-    (C.rand g0)+    (random g0)     where pred u (k, pmf, cdf) = cdf < u   rand (Geometric p g0) =     mapTuple-    (\u -> ceiling $ (log u) / (log (1-p)))+    ((\u -> ceiling $ (log u) / (log (1-p))))     (Geometric p)-    (C.rand g0)+    (random g0)   rand (Poisson g0@(C.Exponential y u)) =     mapTuple     (\x -> length x)     (Poisson)-    ((foldGenWhile (C.rand) (+) 0.0 (<1.0)) g0)+    ((foldGenWhile (C.rand) (+) (0.0) (<1.0)) g0)   rand (Bernoulli p g0) =     mapTuple-    (\x -> if (x >= p) then 1 else 0)+    (\x -> if (x >= p) then 1 else 0)      (Bernoulli p)-    (C.rand g0)+    (random g0)   rand (ZipF n slope cache u0) =      mapTuple     (\u ->-      length $ filter (pred u) cache)+      length $ filter (pred u) cache)      (ZipF n slope cache)-    (C.rand u0)+    (random 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)+    (random g0)    cdf (Poisson (C.Exponential y _)) x =     (1/(exp y)) * (sum [ (y ** (toDbl i)) / (fromInteger $ fac i) | i <- [0..x]])@@ -119,11 +128,11 @@     | otherwise = p   cdf (Binomial n p cache _) x = r     where-      (_, _, r) = head $ filter (\(w,_,_) -> (w==x)) cache+      (_, _, r) = fromMaybe (0, 0, 1) $ maybeHead $ filter (\(w,_,_) -> (w==x)) cache   cdf (ZipF n s cache _) x = r     where       (_, _, r) =-        head $ filter (\(k, _, _) -> x == k) cache +        fromMaybe (0, 0, 1) $ maybeHead $ filter (\(k, _, _) -> x == k) cache    cdf (Uniform a b _) x = toDbl (x-a) / toDbl (b-a)    cdf' g@(Poisson (C.Exponential y _)) x =

@@ -38,14 +38,14 @@             let (xs, g2) = h g1 in             if (p x) then (x:xs, g2) else ([], g1) -genTake :: (g -> (a,g)) -> Integer -> (g -> ([a], g))+genTake :: (Eq b, Num b) => (g -> (a,g)) -> b -> (g -> ([a], g)) genTake f 0 g0 = ([], g0) genTake f n g0 = ((x:xs), g2)   where     (x, g1) = f g0     (xs, g2) = genTake f (n-1) g1 -dropGen :: (g -> (a,g)) -> Integer -> g -> g+dropGen :: (Eq b, Num b) => (g -> (a,g)) -> b -> g -> g dropGen f = d   where     d 0 g0 = g0

@@ -52,11 +52,17 @@   | otherwise = (sqrt (2*pi*n)) * (exp $ n * ((log n) - 1))  fac :: Int -> Integer-fac n = head $ drop (n-1) factorial+fac n = head $ drop (n) factorial  factorial :: [Integer]-factorial = 1 : 1 : (tail (zipWith (*) (factorial) [1..]))+factorial = 1 : (acc (*) 1 [1..]) +acc :: (a -> a -> a) -> a -> [a] -> [a]+acc c z ls = f z ls+  where+    f y [] = []+    f y (x:xs) = (x `c` y) : (f (x `c` y) xs)+ fib :: Int -> Integer fib n = head $ drop (n-1) fibinacci @@ -94,9 +100,10 @@   rel_frequence :: Double,   cum_frequence :: Double,   slope :: Double-  } deriving (Show, Eq)-+  } deriving (Eq) +instance Show Datagram where+  show d = show (lower_bound d, upper_bound d, rel_frequence d)  --peice wise linear  fIHistogram :: [Double] -> Histogram@@ -114,6 +121,7 @@ maxf (x:xs) = foldl (max) x xs minf (x:xs) = foldl (min) x xs +-- TODO fix 'holes' datagramFromRaw :: Int -> Double -> [(Double, Int)] -> [Datagram] datagramFromRaw count iSize = accMap f z   where@@ -140,8 +148,8 @@ 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)+    g prev []     = [prev]+    g prev (x:xs) = prev : (g new xs)       where new = (f prev x) {- Takes the set 

@@ -1,5 +1,9 @@+module Main where+ import Stochastic.Distributions import Stochastic.Analysis+import Stochastic.Uniform+import System.Random import Stochastic.Tools import qualified Stochastic.Distributions.Continuous as C import qualified Stochastic.Distributions.Discrete as D@@ -20,8 +24,9 @@   ,TestLabel "Exponential ChiSquared" chiexp   ,TestLabel "Uniform ChiSquared" chiuni   ,TestLabel "Poisson ChiSquared" chiPoisson-  ,TestLabel "Test of grouping" groupTest+{-  ,TestLabel "Test of grouping" groupTest --  ,TestLabel "Test histogram" binTest+-}   ]  @@ -84,7 +89,7 @@ chiexp = TestCase (   assertWithinDelta   "Sample of the exponential distribution passes the ChiSquaredTest with HIGH confidence"-  (2e-1) 0 (chiTestRandom $ C.mkExp (stdBase 42) 1))+  (0.25) 0 (chiTestRandom $ C.mkExp (stdBase 42) 1))  groupTest = TestCase             (assertEqual@@ -96,7 +101,7 @@ chiPoisson = TestCase (   assertWithinDelta   "Sample of the poisson distribution passes the ChiSquaredTest with HIGH confidence"-  (2e-1) 0 (chiTestDiscreteRandom $ D.mkPoisson (stdBase 42) 1))+  (1.2) 0 (chiTestDiscreteRandom $ D.mkPoisson (stdBase 42) 1))  chiTestDiscreteRandom :: D.Dist                       -> Double@@ -123,6 +128,19 @@  {-  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-}++pleasantAPIExample :: [(Double, Int, Double)]+pleasantAPIExample =+  let g = stdBase 42 in+  let ([x, y, z], g') = nWayAllocate 20 3 g in+  let g1 = C.mkExp x 1 in+  let g2 = D.mkPoisson y 1 in+  let g3 = C.mkNormal z 0 1 in+  zip3 (randoms g1 :: [Double]) (randoms g2 :: [Int])  (randoms g3 :: [Double]) +      +      +      +  -- -- Custom Assertions