diff --git a/random-variates.cabal b/random-variates.cabal
--- a/random-variates.cabal
+++ b/random-variates.cabal
@@ -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
+
+  
diff --git a/src/Helpers.hs b/src/Helpers.hs
deleted file mode 100644
--- a/src/Helpers.hs
+++ /dev/null
@@ -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
-
-
-
diff --git a/src/Stochastic/Analysis.hs b/src/Stochastic/Analysis.hs
new file mode 100644
--- /dev/null
+++ b/src/Stochastic/Analysis.hs
@@ -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
+
+
+
diff --git a/src/Stochastic/Bernoulli.hs b/src/Stochastic/Bernoulli.hs
deleted file mode 100644
--- a/src/Stochastic/Bernoulli.hs
+++ /dev/null
@@ -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
diff --git a/src/Stochastic/Binomial.hs b/src/Stochastic/Binomial.hs
deleted file mode 100644
--- a/src/Stochastic/Binomial.hs
+++ /dev/null
@@ -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
diff --git a/src/Stochastic/Distribution.hs b/src/Stochastic/Distribution.hs
deleted file mode 100644
--- a/src/Stochastic/Distribution.hs
+++ /dev/null
@@ -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
diff --git a/src/Stochastic/Distribution/Continuous.hs b/src/Stochastic/Distribution/Continuous.hs
new file mode 100644
--- /dev/null
+++ b/src/Stochastic/Distribution/Continuous.hs
@@ -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
+  
diff --git a/src/Stochastic/Distribution/Discrete.hs b/src/Stochastic/Distribution/Discrete.hs
new file mode 100644
--- /dev/null
+++ b/src/Stochastic/Distribution/Discrete.hs
@@ -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
+
diff --git a/src/Stochastic/Distributions.hs b/src/Stochastic/Distributions.hs
--- a/src/Stochastic/Distributions.hs
+++ b/src/Stochastic/Distributions.hs
@@ -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
+
 
diff --git a/src/Stochastic/Distributions/Continuous.hs b/src/Stochastic/Distributions/Continuous.hs
new file mode 100644
--- /dev/null
+++ b/src/Stochastic/Distributions/Continuous.hs
@@ -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
+
+
diff --git a/src/Stochastic/Distributions/Discrete.hs b/src/Stochastic/Distributions/Discrete.hs
new file mode 100644
--- /dev/null
+++ b/src/Stochastic/Distributions/Discrete.hs
@@ -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
+
diff --git a/src/Stochastic/Exponential.hs b/src/Stochastic/Exponential.hs
deleted file mode 100644
--- a/src/Stochastic/Exponential.hs
+++ /dev/null
@@ -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)
diff --git a/src/Stochastic/Generator.hs b/src/Stochastic/Generator.hs
new file mode 100644
--- /dev/null
+++ b/src/Stochastic/Generator.hs
@@ -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)
+
+
+
diff --git a/src/Stochastic/Geometric.hs b/src/Stochastic/Geometric.hs
deleted file mode 100644
--- a/src/Stochastic/Geometric.hs
+++ /dev/null
@@ -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
diff --git a/src/Stochastic/Normal.hs b/src/Stochastic/Normal.hs
deleted file mode 100644
--- a/src/Stochastic/Normal.hs
+++ /dev/null
@@ -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)
-    
diff --git a/src/Stochastic/Poisson.hs b/src/Stochastic/Poisson.hs
deleted file mode 100644
--- a/src/Stochastic/Poisson.hs
+++ /dev/null
@@ -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)
diff --git a/src/Stochastic/Tools.hs b/src/Stochastic/Tools.hs
new file mode 100644
--- /dev/null
+++ b/src/Stochastic/Tools.hs
@@ -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
+
+
diff --git a/src/Stochastic/Uniform.hs b/src/Stochastic/Uniform.hs
deleted file mode 100644
--- a/src/Stochastic/Uniform.hs
+++ /dev/null
@@ -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)
-
-
diff --git a/src/Stochastic/ZipF.hs b/src/Stochastic/ZipF.hs
deleted file mode 100644
--- a/src/Stochastic/ZipF.hs
+++ /dev/null
@@ -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
-
-
-
diff --git a/tests/unit.hs b/tests/unit.hs
new file mode 100644
--- /dev/null
+++ b/tests/unit.hs
@@ -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 ()
diff --git a/tests/vis.hs b/tests/vis.hs
new file mode 100644
--- /dev/null
+++ b/tests/vis.hs
@@ -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