diff --git a/README.markdown b/README.markdown
--- a/README.markdown
+++ b/README.markdown
@@ -14,21 +14,6 @@
 libraries and applications that use this library using a high level of
 optimisation.
 
-Suggested GHC options:
-
-    -O -funbox-strict-fields
-
-To illustrate, here are the times (in seconds) to generate and sum 250
-million random Word32 values, on a laptop with a 2.4GHz Core2 Duo
-P8600 processor, running Fedora 11 and GHC 6.10.3:
-
-    no flags   200+
-    -O           1.249
-    -O -fvia-C   0.991
-
-As the numbers above suggest, compiling without optimisation will
-yield unacceptable performance.
-
 
 # Get involved!
 
diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs
--- a/Statistics/Distribution.hs
+++ b/Statistics/Distribution.hs
@@ -24,6 +24,7 @@
       -- ** Random number generation
     , ContGen(..)
     , DiscreteGen(..)
+    , genContinous
       -- * Helper functions
     , findRoot
     , sumProbabilities
@@ -123,7 +124,13 @@
 class (DiscreteDistr d, ContGen d) => DiscreteGen d where
   genDiscreteVar :: PrimMonad m => d -> Gen (PrimState m) -> m Int
 
-
+-- | Generate variates from continous distribution using inverse
+--   transform rule.
+genContinous :: (ContDistr d, PrimMonad m) => d -> Gen (PrimState m) -> m Double
+genContinous d gen = do
+  x <- uniform gen
+  return $! quantile d x
+{-# INLINE genContinous #-}
 
 data P = P {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
diff --git a/Statistics/Distribution/Beta.hs b/Statistics/Distribution/Beta.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Distribution/Beta.hs
@@ -0,0 +1,98 @@
+{-# LANGUAGE DeriveDataTypeable #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Statistics.Distribution.Beta
+-- Copyright   :  (C) 2012 Edward Kmett,
+-- License     :  BSD-style (see the file LICENSE)
+--
+-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
+-- Stability   :  provisional
+-- Portability :  DeriveDataTypeable
+--
+----------------------------------------------------------------------------
+module Statistics.Distribution.Beta
+  ( BetaDistribution
+    -- * Constructor
+  , betaDistr
+  , improperBetaDistr
+    -- * Accessors
+  , bdAlpha
+  , bdBeta
+  ) where
+
+import Numeric.SpecFunctions           (incompleteBeta, invIncompleteBeta, logBeta)
+import Numeric.MathFunctions.Constants (m_NaN)
+import qualified Statistics.Distribution as D
+import Data.Typeable
+
+-- | The beta distribution
+data BetaDistribution = BD
+ { bdAlpha :: {-# UNPACK #-} !Double
+   -- ^ Alpha shape parameter
+ , bdBeta  :: {-# UNPACK #-} !Double
+   -- ^ Beta shape parameter
+ } deriving (Eq,Read,Show,Typeable)
+
+-- | Create beta distribution. Both shape parameters must be positive.
+betaDistr :: Double             -- ^ Shape parameter alpha
+          -> Double             -- ^ Shape parameter beta
+          -> BetaDistribution
+betaDistr a b
+  | a > 0 && b > 0 = improperBetaDistr a b
+  | otherwise      =
+      error $  "Statistics.Distribution.Beta.betaDistr: "
+            ++ "shape parameters must be positive. Got a = "
+            ++ show a
+            ++ " b = "
+            ++ show b
+{-# INLINE betaDistr #-}
+
+-- | Create beta distribution. This construtor doesn't check parameters.
+improperBetaDistr :: Double             -- ^ Shape parameter alpha
+                  -> Double             -- ^ Shape parameter beta
+                  -> BetaDistribution
+improperBetaDistr = BD
+{-# INLINE improperBetaDistr #-}
+
+instance D.Distribution BetaDistribution where
+  cumulative (BD a b) x
+    | x <= 0    = 0
+    | x >= 1    = 1
+    | otherwise = incompleteBeta a b x
+  {-# INLINE cumulative #-}
+
+instance D.Mean BetaDistribution where
+  mean (BD a b) = a / (a + b)
+  {-# INLINE mean #-}
+
+instance D.MaybeMean BetaDistribution where
+  maybeMean = Just . D.mean
+  {-# INLINE maybeMean #-}
+
+instance D.Variance BetaDistribution where
+  variance (BD a b) = a*b / (apb*apb*(apb+1))
+    where apb = a + b
+  {-# INLINE variance #-}
+
+instance D.MaybeVariance BetaDistribution where
+  maybeVariance = Just . D.variance
+  {-# INLINE maybeVariance #-}
+
+instance D.ContDistr BetaDistribution where
+  density (BD a b) x
+   | a <= 0 || b <= 0 = m_NaN
+   | x <= 0 = 0
+   | x >= 1 = 0
+   | otherwise = exp $ (a-1)*log x + (b-1)*log (1-x) - logBeta a b
+  {-# INLINE density #-}
+
+  quantile (BD a b) p
+    | p == 0         = 0
+    | p == 1         = 1
+    | p > 0 && p < 1 = invIncompleteBeta a b p
+    | otherwise      =
+        error $ "Statistics.Distribution.Gamma.quantile: p must be in [0,1] range. Got: "++show p
+  {-# INLINE quantile #-}
+
+instance D.ContGen BetaDistribution where
+  genContVar = D.genContinous
diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs
--- a/Statistics/Distribution/Binomial.hs
+++ b/Statistics/Distribution/Binomial.hs
@@ -85,13 +85,13 @@
 {-# INLINE variance #-}
 
 -- | Construct binomial distribution. Number of trials must be
---   positive and probability must be in [0,1] range
+--   non-negative and probability must be in [0,1] range
 binomial :: Int                 -- ^ Number of trials.
          -> Double              -- ^ Probability.
          -> BinomialDistribution
 binomial n p 
-  | n <= 0         = 
-    error $ msg ++ "number of trials must be positive. Got " ++ show n
+  | n < 0          =
+    error $ msg ++ "number of trials must be non-negative. Got " ++ show n
   | p < 0 || p > 1 = 
     error $ msg++"probability must be in [0,1] range. Got " ++ show p
   | otherwise      = BD n p
diff --git a/Statistics/Distribution/CauchyLorentz.hs b/Statistics/Distribution/CauchyLorentz.hs
--- a/Statistics/Distribution/CauchyLorentz.hs
+++ b/Statistics/Distribution/CauchyLorentz.hs
@@ -63,3 +63,6 @@
     | p == 1         =  1 / 0
     | otherwise      =
       error $ "Statistics.Distribution.CauchyLorentz..quantile: p must be in [0,1] range. Got: "++show p
+
+instance D.ContGen CauchyDistribution where
+  genContVar = D.genContinous
diff --git a/Statistics/Distribution/ChiSquared.hs b/Statistics/Distribution/ChiSquared.hs
--- a/Statistics/Distribution/ChiSquared.hs
+++ b/Statistics/Distribution/ChiSquared.hs
@@ -21,7 +21,8 @@
 import Data.Typeable         (Typeable)
 import Numeric.SpecFunctions (incompleteGamma,invIncompleteGamma,logGamma)
 
-import qualified Statistics.Distribution as D
+import qualified Statistics.Distribution         as D
+import qualified System.Random.MWC.Distributions as MWC
 
 
 -- | Chi-squared distribution
@@ -63,6 +64,10 @@
 instance D.MaybeVariance ChiSquared where
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
+
+instance D.ContGen ChiSquared where
+    genContVar (ChiSquared n) = MWC.chiSquare n
+
 
 cumulative :: ChiSquared -> Double -> Double
 cumulative chi x
diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs
--- a/Statistics/Distribution/Exponential.hs
+++ b/Statistics/Distribution/Exponential.hs
@@ -24,8 +24,9 @@
     ) where
 
 import Data.Typeable (Typeable)
-import qualified Statistics.Distribution as D
-import qualified Statistics.Sample as S
+import qualified Statistics.Distribution         as D
+import qualified Statistics.Sample               as S
+import qualified System.Random.MWC.Distributions as MWC
 import Statistics.Types (Sample)
 
 newtype ExponentialDistribution = ED {
@@ -54,6 +55,9 @@
 instance D.MaybeVariance ExponentialDistribution where
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
+
+instance D.ContGen ExponentialDistribution where
+  genContVar = MWC.exponential . edLambda
 
 cumulative :: ExponentialDistribution -> Double -> Double
 cumulative (ED l) x | x <= 0    = 0
diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs
--- a/Statistics/Distribution/FDistribution.hs
+++ b/Statistics/Distribution/FDistribution.hs
@@ -22,7 +22,7 @@
 
 
 
--- | Student-T distribution
+-- | F distribution
 data FDistribution = F { fDistributionNDF1 :: {-# UNPACK #-} !Double
                        , fDistributionNDF2 :: {-# UNPACK #-} !Double
                        , _pdfFactor        :: {-# UNPACK #-} !Double
@@ -74,6 +74,9 @@
   maybeStdDev (F n m _) 
     | m > 4     = Just $ 2 * sqr m * (m + n - 2) / (n * sqr (m - 2) * (m - 4))
     | otherwise = Nothing
+
+instance D.ContGen FDistribution where
+  genContVar = D.genContinous
 
 sqr :: Double -> Double
 sqr x = x * x
diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs
--- a/Statistics/Distribution/Gamma.hs
+++ b/Statistics/Distribution/Gamma.hs
@@ -19,6 +19,7 @@
       GammaDistribution
     -- * Constructors
     , gammaDistr
+    , improperGammaDistr
     -- * Accessors
     , gdShape
     , gdScale
@@ -27,8 +28,9 @@
 import Data.Typeable (Typeable)
 import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN)
 import Numeric.SpecFunctions           (incompleteGamma, invIncompleteGamma)
-import Statistics.Distribution.Poisson.Internal as Poisson
-import qualified Statistics.Distribution        as D
+import Statistics.Distribution.Poisson.Internal  as Poisson
+import qualified Statistics.Distribution         as D
+import qualified System.Random.MWC.Distributions as MWC
 
 -- | The gamma distribution.
 data GammaDistribution = GD {
@@ -44,10 +46,18 @@
 gammaDistr k theta
   | k     <= 0 = error $ msg ++ "shape must be positive. Got " ++ show k
   | theta <= 0 = error $ msg ++ "scale must be positive. Got " ++ show theta
-  | otherwise  = GD k theta
+  | otherwise  = improperGammaDistr k theta
     where msg = "Statistics.Distribution.Gamma.gammaDistr: "
 {-# INLINE gammaDistr #-}
 
+-- | Create gamma distribution. This constructor do not check whether
+--   parameters are valid
+improperGammaDistr :: Double            -- ^ Shape parameter. /k/
+                   -> Double            -- ^ Scale parameter, &#977;.
+                   -> GammaDistribution
+improperGammaDistr = GD
+{-# INLINE improperGammaDistr #-}
+
 instance D.Distribution GammaDistribution where
     cumulative = cumulative
 
@@ -70,6 +80,8 @@
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
 
+instance D.ContGen GammaDistribution where
+    genContVar (GD a l) = MWC.gamma a l
 
 
 density :: GammaDistribution -> Double -> Double
diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs
--- a/Statistics/Distribution/Geometric.hs
+++ b/Statistics/Distribution/Geometric.hs
@@ -59,10 +59,9 @@
 geometric :: Double                -- ^ Success rate
           -> GeometricDistribution
 geometric x
-  | x < 0 || x > 1 = 
+  | x >= 0 && x <= 1 = GD x
+  | otherwise        =
     error $ "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x
-  | otherwise      =
-    GD x
 {-# INLINE geometric #-}
 
 probability :: GeometricDistribution -> Int -> Double
diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs
--- a/Statistics/Distribution/Normal.hs
+++ b/Statistics/Distribution/Normal.hs
@@ -57,8 +57,8 @@
     stdDev = stdDev
 
 instance D.ContGen NormalDistribution where
-    genContVar d gen = do x <- MWC.standard gen
-                          return $! stdDev d * (x - mean d)
+    genContVar d = MWC.normal (mean d) (stdDev d)
+    {-# INLINE genContVar #-}
 
 -- | Standard normal distribution with mean equal to 0 and variance equal to 1
 standard :: NormalDistribution
diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs
--- a/Statistics/Distribution/Poisson.hs
+++ b/Statistics/Distribution/Poisson.hs
@@ -63,9 +63,9 @@
 -- | Create Poisson distribution.
 poisson :: Double -> PoissonDistribution
 poisson l
-  | l <= 0    = error $ "Statistics.Distribution.Poisson.poisson:\
-                        \ lambda must be positive. Got " ++ show l
-  | otherwise = PD l
+  | l >=  0   = PD l
+  | otherwise = error $ "Statistics.Distribution.Poisson.poisson:\
+                        \ lambda must be non-negative. Got " ++ show l
 {-# INLINE poisson #-}
 
 -- $references
diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs
--- a/Statistics/Distribution/StudentT.hs
+++ b/Statistics/Distribution/StudentT.hs
@@ -69,3 +69,6 @@
 instance D.MaybeVariance StudentT where
   maybeStdDev (StudentT ndf) | ndf > 2   = Just $ ndf / (ndf - 2)
                              | otherwise = Nothing
+
+instance D.ContGen StudentT where
+  genContVar = D.genContinous
diff --git a/Statistics/Distribution/Uniform.hs b/Statistics/Distribution/Uniform.hs
--- a/Statistics/Distribution/Uniform.hs
+++ b/Statistics/Distribution/Uniform.hs
@@ -9,19 +9,26 @@
 -- Portability : portable
 --
 -- Variate distributed uniformly in the interval.
-module Statistics.Distribution.Uniform (
-    UniformDistribution
-  , uniformDistr
-  ) where
+module Statistics.Distribution.Uniform
+    (
+      UniformDistribution
+    -- * Constructors
+    , uniformDistr
+    -- ** Accessors
+    , uniformA
+    , uniformB
+    ) where
 
 import Data.Typeable (Typeable)
 import qualified Statistics.Distribution as D
 import qualified System.Random.MWC       as MWC
 
 
--- | Uniform distribution
-data UniformDistribution = UniformDistribution {-# UNPACK #-} !Double {-# UNPACK #-} !Double
-                           deriving (Eq,Show,Read,Typeable)
+-- | Uniform distribution from A to B
+data UniformDistribution = UniformDistribution {
+      uniformA :: {-# UNPACK #-} !Double -- ^ Low boundary of distribution
+    , uniformB :: {-# UNPACK #-} !Double -- ^ Upper boundary of distribution
+    } deriving (Eq, Read, Show, Typeable)
 
 -- | Create uniform distribution.
 uniformDistr :: Double -> Double -> UniformDistribution
diff --git a/Statistics/Sample/Histogram.hs b/Statistics/Sample/Histogram.hs
--- a/Statistics/Sample/Histogram.hs
+++ b/Statistics/Sample/Histogram.hs
@@ -19,6 +19,7 @@
     , range
     ) where
 
+import Numeric.MathFunctions.Constants (m_epsilon)
 import Statistics.Function (minMax)
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Generic.Mutable as GM
@@ -72,7 +73,7 @@
          GM.write bins b . (+1) =<< GM.read bins b
          go (i+1)
        len = G.length xs
-       d = (hi - lo) / fromIntegral numBins
+       d = ((hi - lo) * (1 + realToFrac m_epsilon)) / fromIntegral numBins
 {-# INLINE histogram_ #-}
 
 -- | /O(n)/ Compute decent defaults for the lower and upper bounds of
diff --git a/Statistics/Sample/KernelDensity.hs b/Statistics/Sample/KernelDensity.hs
--- a/Statistics/Sample/KernelDensity.hs
+++ b/Statistics/Sample/KernelDensity.hs
@@ -52,7 +52,8 @@
 kde n0 xs = kde_ n0 (lo - range / 10) (hi + range / 10) xs
   where
     (lo,hi) = minMax xs
-    range = hi - lo
+    range   | U.length xs <= 1 = 1       -- Unreasonable guess
+            | otherwise        = hi - lo
 
 -- | Gaussian kernel density estimator for one-dimensional data, using
 -- the method of Botev et al.
@@ -72,19 +73,20 @@
      -- ^ Upper bound (@max@) of the mesh range.
      -> U.Vector Double -> (U.Vector Double, U.Vector Double)
 kde_ n0 min max xs
+  | U.null xs = error "Statistics.KernelDensity.kde: empty sample"
   | n0 < 1    = error "Statistics.KernelDensity.kde: invalid number of points"
   | otherwise = (mesh, density)
   where
     mesh = G.generate ni $ \z -> min + (d * fromIntegral z)
         where d = r / (n-1)
-    density = G.map (/r) . idct $ G.zipWith f a (G.enumFromTo 0 (n-1))
+    density = G.map (/(2 * r)) . idct $ G.zipWith f a (G.enumFromTo 0 (n-1))
       where f b z = b * exp (sqr z * sqr pi * t_star * (-0.5))
-    !n = fromIntegral ni
+    !n  = fromIntegral ni
     !ni = nextHighestPowerOfTwo n0
-    !r = max - min
-    a = dct . G.map (/ G.sum h) $ h
-      where h = G.map (/ len) $ histogram_ ni min max xs
-    !len = fromIntegral (G.length xs)
+    !r  = max - min
+    a   = dct . G.map (/ G.sum h) $ h
+        where h = G.map (/ len) $ histogram_ ni min max xs
+    !len    = fromIntegral (G.length xs)
     !t_star = fromRoot (0.28 * len ** (-0.4)) . ridders 1e-14 (0,0.1) $ \x ->
               x - (len * (2 * sqrt pi) * go 6 (f 7 x)) ** (-0.4)
       where
@@ -94,10 +96,9 @@
                 iv = G.map sqr $ G.enumFromTo 1 (n-1)
         go s !h | s == 1    = h
                 | otherwise = go (s-1) (f s time)
-          where time = (2 * const * k0 / len / h) ** (2 / (3 + 2 * s))
+          where time  = (2 * const * k0 / len / h) ** (2 / (3 + 2 * s))
                 const = (1 + 0.5 ** (s+0.5)) / 3
-                k0 = U.product (G.enumFromThenTo 1 3 (2*s-1)) / m_sqrt_2_pi
-    _bandwidth = sqrt t_star * r
+                k0    = U.product (G.enumFromThenTo 1 3 (2*s-1)) / m_sqrt_2_pi
     sqr x = x * x
 
 -- $references
diff --git a/Statistics/Test/WilcoxonT.hs b/Statistics/Test/WilcoxonT.hs
--- a/Statistics/Test/WilcoxonT.hs
+++ b/Statistics/Test/WilcoxonT.hs
@@ -74,7 +74,7 @@
 -- all the coefficients by r down the list.
 --
 -- This list will be processed lazily from the head.
-coefficients :: Int -> [Int]
+coefficients :: Int -> [Integer]
 coefficients 1 = [1, 1] -- 1 + x
 coefficients r = let coeffs = coefficients (r-1)
                      (firstR, rest) = splitAt r coeffs
@@ -86,7 +86,10 @@
 
 -- This list will be processed lazily from the head.
 summedCoefficients :: Int -> [Double]
-summedCoefficients = map fromIntegral . scanl1 (+) . coefficients
+summedCoefficients n
+  | n < 1     = error "Statistics.Test.WilcoxonT.summedCoefficients: nonpositive sample size"
+  | n > 1023  = error "Statistics.Test.WilcoxonT.summedCoefficients: sample is too large (see bug #18)"
+  | otherwise = map fromIntegral $ scanl1 (+) $ coefficients n
 
 -- | Tests whether a given result from a Wilcoxon signed-rank matched-pairs test
 -- is significant at the given level.
diff --git a/Statistics/Transform.hs b/Statistics/Transform.hs
--- a/Statistics/Transform.hs
+++ b/Statistics/Transform.hs
@@ -43,37 +43,50 @@
 
 -- | Discrete cosine transform (DCT-II).
 dct :: U.Vector Double -> U.Vector Double
-dct = dct_ . G.map (:+0)
+dct = dctWorker . G.map (:+0)
 
--- | Discrete cosine transform, with complex coefficients (DCT-II).
+-- | Discrete cosine transform (DCT-II). Only real part of vector is
+--   transformed, imaginary part is ignored.
 dct_ :: U.Vector CD -> U.Vector Double
-dct_ xs = G.map realPart $ G.zipWith (*) weights (fft interleaved)
+dct_ = dctWorker . G.map (\(i :+ _) -> i :+ 0)
+
+dctWorker :: U.Vector CD -> U.Vector Double
+dctWorker xs
+  = G.map realPart $ G.zipWith (*) weights (fft interleaved)
   where
     interleaved = G.backpermute xs $ G.enumFromThenTo 0 2 (len-2) G.++
                                      G.enumFromThenTo (len-1) (len-3) 1
-    weights = G.cons 1 . G.generate (len-1) $ \x ->
+    weights = G.cons 2 . G.generate (len-1) $ \x ->
               2 * exp ((0:+(-1))*fi (x+1)*pi/(2*n))
       where n = fi len
     len = G.length xs
 
+
+
 -- | Inverse discrete cosine transform (DCT-III). It's inverse of
 -- 'dct' only up to scale parameter:
 --
 -- > (idct . dct) x = (* lenngth x)
 idct :: U.Vector Double -> U.Vector Double
-idct = idct_ . G.map (:+0)
+idct = idctWorker . G.map (:+0)
 
--- | Inverse discrete cosine transform, with complex coefficients
--- (DCT-III).
+-- | Inverse discrete cosine transform (DCT-III). Only real part of vector is
+--   transformed, imaginary part is ignored.
 idct_ :: U.Vector CD -> U.Vector Double
-idct_ xs = G.generate len interleave
+idct_ = idctWorker . G.map (\(i :+ _) -> i :+ 0)
+
+idctWorker :: U.Vector CD -> U.Vector Double
+idctWorker xs = G.generate len interleave
   where
     interleave z | even z    = vals `G.unsafeIndex` halve z
                  | otherwise = vals `G.unsafeIndex` (len - halve z - 1)
     vals = G.map realPart . ifft $ G.zipWith (*) weights xs
-    weights = G.generate len $ \x -> n * exp ((0:+1)*fi x*pi/(2*n))
+    weights 
+      = G.cons n
+      $ G.generate (len - 1) $ \x -> 2 * n * exp ((0:+1) * fi (x+1) * pi/(2*n))
       where n = fi len
     len = G.length xs
+
 
 -- | Inverse fast Fourier transform.
 ifft :: U.Vector CD -> U.Vector CD
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.10.1.0
+version:        0.10.2.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -22,6 +22,26 @@
   * Common statistical tests for significant differences between
     samples.
   .
+  Changes in 0.10.2.0
+  .
+  * Bugs in DCT and IDCT are fixed.
+  .
+  * Accesors for uniform distribution are added.
+  .
+  * 'ContGen' instances for all continous distribtuions are added.
+  .
+  * Beta distribution is added.
+  .
+  * Constructor for improper gamma distribtuion is added.
+  .
+  * Binomial distribution allows zero trials.
+  .
+  * Poisson distribution now accept zero parameter.
+  .
+  * Integer overflow in caculation of Wilcoxon-T test is fixed.
+  .
+  * Bug in 'ContGen' instance for normal distribution is fixed.
+  .
   Changes in 0.10.1.0
   .
   * Kolmogorov-Smirnov nonparametric test added.
@@ -135,6 +155,7 @@
     Statistics.Autocorrelation
     Statistics.Constants
     Statistics.Distribution
+    Statistics.Distribution.Beta
     Statistics.Distribution.Binomial
     Statistics.Distribution.CauchyLorentz
     Statistics.Distribution.ChiSquared
@@ -199,10 +220,11 @@
   other-modules:
     Tests.Distribution
     Tests.Helpers
-    Tests.Math
-    Tests.Math.Tables
+    Tests.Function
     Tests.NonparametricTest
+    Tests.NonparametricTest.Table
     Tests.Transform
+    Tests.KDE
 
   ghc-options:
     -Wall -threaded -rtsopts
@@ -215,6 +237,7 @@
     test-framework,
     test-framework-quickcheck2,
     test-framework-hunit,
+    math-functions,
     statistics,
     primitive,
     vector,
diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/Distribution.hs
@@ -0,0 +1,284 @@
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# LANGUAGE ScopedTypeVariables  #-}
+-- Required for Param
+{-# LANGUAGE FlexibleInstances    #-}
+{-# LANGUAGE OverlappingInstances #-}
+{-# LANGUAGE ViewPatterns #-}
+module Tests.Distribution (
+    distributionTests
+  ) where
+
+import Control.Applicative
+import Control.Exception
+
+import Data.List     (find)
+import Data.Typeable (Typeable)
+
+import qualified Numeric.IEEE    as IEEE
+
+import Test.Framework                       (Test,testGroup)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.QuickCheck         as QC
+import Test.QuickCheck.Monadic as QC
+import Text.Printf
+
+import Statistics.Distribution
+import Statistics.Distribution.Beta
+import Statistics.Distribution.Binomial
+import Statistics.Distribution.ChiSquared
+import Statistics.Distribution.CauchyLorentz
+import Statistics.Distribution.Exponential
+import Statistics.Distribution.FDistribution
+import Statistics.Distribution.Gamma
+import Statistics.Distribution.Geometric
+import Statistics.Distribution.Hypergeometric
+import Statistics.Distribution.Normal
+import Statistics.Distribution.Poisson
+import Statistics.Distribution.StudentT
+import Statistics.Distribution.Uniform
+
+import Prelude hiding (catch)
+
+import Tests.Helpers
+
+
+-- | Tests for all distributions
+distributionTests :: Test
+distributionTests = testGroup "Tests for all distributions"
+  [ contDistrTests (T :: T BetaDistribution        )
+  , contDistrTests (T :: T CauchyDistribution      )
+  , contDistrTests (T :: T ChiSquared              )
+  , contDistrTests (T :: T ExponentialDistribution )
+  , contDistrTests (T :: T GammaDistribution       )
+  , contDistrTests (T :: T NormalDistribution      )
+  , contDistrTests (T :: T UniformDistribution     )
+  , contDistrTests (T :: T StudentT                )
+  , contDistrTests (T :: T FDistribution           )
+
+  , discreteDistrTests (T :: T BinomialDistribution       )
+  , discreteDistrTests (T :: T GeometricDistribution      )
+  , discreteDistrTests (T :: T HypergeometricDistribution )
+  , discreteDistrTests (T :: T PoissonDistribution        )
+
+  , unitTests
+  ]
+
+----------------------------------------------------------------
+-- Tests
+----------------------------------------------------------------
+
+-- Tests for continous distribution
+contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
+contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
+  cdfTests t ++
+  [ testProperty "PDF sanity"              $ pdfSanityCheck   t
+  , testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t
+  , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t
+  ]
+
+-- Tests for discrete distribution
+discreteDistrTests :: (Param d, DiscreteDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
+discreteDistrTests t = testGroup ("Tests for: " ++ typeName t) $
+  cdfTests t ++
+  [ testProperty "Prob. sanity"         $ probSanityCheck       t
+  , testProperty "CDF is sum of prob."  $ discreteCDFcorrect    t
+  ]
+
+-- Tests for distributions which have CDF
+cdfTests :: (Param d, Distribution d, QC.Arbitrary d, Show d) => T d -> [Test]
+cdfTests t =
+  [ testProperty "C.D.F. sanity"        $ cdfSanityCheck         t
+  , testProperty "CDF limit at +∞"      $ cdfLimitAtPosInfinity  t
+  , testProperty "CDF limit at -∞"      $ cdfLimitAtNegInfinity  t
+  , testProperty "CDF is nondecreasing" $ cdfIsNondecreasing     t
+  , testProperty "1-CDF is correct"     $ cdfComplementIsCorrect t
+  ]
+----------------------------------------------------------------
+
+-- CDF is in [0,1] range
+cdfSanityCheck :: (Distribution d) => T d -> d -> Double -> Bool
+cdfSanityCheck _ d x = c >= 0 && c <= 1 
+  where c = cumulative d x
+
+-- CDF never decreases
+cdfIsNondecreasing :: (Distribution d) => T d -> d -> Double -> Double -> Bool
+cdfIsNondecreasing _ d = monotonicallyIncreasesIEEE $ cumulative d
+
+-- CDF limit at +∞ is 1
+cdfLimitAtPosInfinity :: (Param d, Distribution d) => T d -> d -> Property
+cdfLimitAtPosInfinity _ d =
+  okForInfLimit d ==> printTestCase ("Last elements: " ++ show (drop 990 probs))
+                    $ Just 1.0 == (find (>=1) probs)
+  where
+    probs = take 1000 $ map (cumulative d) $ iterate (*1.4) 1
+
+-- CDF limit at -∞ is 0
+cdfLimitAtNegInfinity :: (Param d, Distribution d) => T d -> d -> Property
+cdfLimitAtNegInfinity _ d =
+  okForInfLimit d ==> printTestCase ("Last elements: " ++ show (drop 990 probs))
+                    $ case find (< IEEE.epsilon) probs of
+                        Nothing -> False
+                        Just p  -> p >= 0
+  where
+    probs = take 1000 $ map (cumulative d) $ iterate (*1.4) (-1)
+
+-- CDF's complement is implemented correctly
+cdfComplementIsCorrect :: (Distribution d) => T d -> d -> Double -> Bool
+cdfComplementIsCorrect _ d x = (eq 1e-14) 1 (cumulative d x + complCumulative d x)
+
+
+-- PDF is positive
+pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool
+pdfSanityCheck _ d x = p >= 0
+  where p = density d x
+
+-- Quantile is inverse of CDF
+quantileIsInvCDF :: (Param d, ContDistr d) => T d -> d -> Double -> Property
+quantileIsInvCDF _ d (snd . properFraction -> p) =
+  p > 0 && p < 1  ==> ( printTestCase (printf "Quantile     = %g" q )
+                      $ printTestCase (printf "Probability  = %g" p )
+                      $ printTestCase (printf "Probability' = %g" p')
+                      $ printTestCase (printf "Error        = %e" (abs $ p - p'))
+                      $ abs (p - p') < invQuantilePrec d
+                      )
+  where
+    q  = quantile   d p
+    p' = cumulative d q
+
+-- Test that quantile fails if p<0 or p>1
+quantileShouldFail :: (ContDistr d) => T d -> d -> Double -> Property
+quantileShouldFail _ d p =
+  p < 0 || p > 1 ==> QC.monadicIO $ do r <- QC.run $ catch
+                                              (do { return $! quantile d p; return False })
+                                              (\(e :: SomeException) -> return True)
+                                       QC.assert r
+
+
+-- Probability is in [0,1] range
+probSanityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Bool
+probSanityCheck _ d x = p >= 0 && p <= 1 
+  where p = probability d x
+
+-- Check that discrete CDF is correct
+discreteCDFcorrect :: (DiscreteDistr d) => T d -> d -> Int -> Int -> Property
+discreteCDFcorrect _ d a b
+  = printTestCase (printf "CDF = %g" p1)
+  $ printTestCase (printf "Sum = %g" p2)
+  $ printTestCase (printf "Δ   = %g" (abs (p1 - p2)))
+  $ abs (p1 - p2) < 3e-10
+  -- Avoid too large differeneces. Otherwise there is to much to sum
+  --
+  -- Absolute difference is used guard againist precision loss when
+  -- close values of CDF are subtracted
+  where
+    n  = min a b
+    m  = n + (abs (a - b) `mod` 100)
+    p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)
+    p2 = sum $ map (probability d) [n .. m]
+
+
+    
+----------------------------------------------------------------
+-- Arbitrary instances for ditributions
+----------------------------------------------------------------
+
+instance QC.Arbitrary BinomialDistribution where
+  arbitrary = binomial <$> QC.choose (1,100) <*> QC.choose (0,1)
+instance QC.Arbitrary ExponentialDistribution where
+  arbitrary = exponential <$> QC.choose (0,100)
+instance QC.Arbitrary GammaDistribution where
+  arbitrary = gammaDistr <$> QC.choose (0.1,10) <*> QC.choose (0.1,10)
+instance QC.Arbitrary BetaDistribution where
+  arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,10)
+instance QC.Arbitrary GeometricDistribution where
+  arbitrary = geometric <$> QC.choose (0,1)
+instance QC.Arbitrary HypergeometricDistribution where
+  arbitrary = do l <- QC.choose (1,20)
+                 m <- QC.choose (0,l)
+                 k <- QC.choose (1,l)
+                 return $ hypergeometric m l k
+instance QC.Arbitrary NormalDistribution where
+  arbitrary = normalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-3, 1e3)
+instance QC.Arbitrary PoissonDistribution where
+  arbitrary = poisson <$> QC.choose (0,1)
+instance QC.Arbitrary ChiSquared where
+  arbitrary = chiSquared <$> QC.choose (1,100)
+instance QC.Arbitrary UniformDistribution where
+  arbitrary = do a <- QC.arbitrary
+                 b <- QC.arbitrary `suchThat` (/= a)
+                 return $ uniformDistr a b
+instance QC.Arbitrary CauchyDistribution where
+  arbitrary = cauchyDistribution
+                <$> arbitrary
+                <*> ((abs <$> arbitrary) `suchThat` (> 0))
+instance QC.Arbitrary StudentT where
+  arbitrary = studentT <$> ((abs <$> arbitrary) `suchThat` (>0))
+instance QC.Arbitrary FDistribution where
+  arbitrary =  fDistribution 
+           <$> ((abs <$> arbitrary) `suchThat` (>0))
+           <*> ((abs <$> arbitrary) `suchThat` (>0))
+
+
+
+-- Parameters for distribution testing. Some distribution require
+-- relaxing parameters a bit
+class Param a where
+  -- Precision for quantileIsInvCDF
+  invQuantilePrec :: a -> Double
+  invQuantilePrec _ = 1e-14
+  -- Distribution is OK for testing limits
+  okForInfLimit :: a -> Bool
+  okForInfLimit _ = True
+
+
+instance Param a
+
+instance Param StudentT where
+  invQuantilePrec _ = 1e-13
+  okForInfLimit   d = studentTndf d > 0.75
+
+instance Param FDistribution where
+  invQuantilePrec _ = 1e-12
+
+
+
+----------------------------------------------------------------
+-- Unit tests
+----------------------------------------------------------------
+
+unitTests :: Test
+unitTests = testGroup "Unit tests"
+  [ testAssertion "density (gammaDistr 150 1/150) 1 == 4.883311" $
+      4.883311418525483 =~ (density (gammaDistr 150 (1/150)) 1)
+    -- Student-T
+  , testStudentPDF 0.3  1.34  0.0648215  -- PDF
+  , testStudentPDF 1    0.42  0.27058
+  , testStudentPDF 4.4  0.33  0.352994
+  , testStudentCDF 0.3  3.34  0.757146   -- CDF
+  , testStudentCDF 1    0.42  0.626569
+  , testStudentCDF 4.4  0.33  0.621739
+    -- F-distribution
+  , testFdistrPDF  1  3   3     (1/(6 * pi)) -- PDF
+  , testFdistrPDF  2  2   1.2   0.206612
+  , testFdistrPDF  10 12  8     0.000385613179281892790166
+  , testFdistrCDF  1  3   3     0.81830988618379067153 -- CDF
+  , testFdistrCDF  2  2   1.2   0.545455
+  , testFdistrCDF  10 12  8     0.99935509863451408041
+  ]
+  where
+    -- Student-T
+    testStudentPDF ndf x exact
+      = testAssertion (printf "density (studentT %f) %f ≈ %f" ndf x exact)
+      $ eq 1e-5  exact  (density (studentT ndf) x)
+    testStudentCDF ndf x exact
+      = testAssertion (printf "cumulative (studentT %f) %f ≈ %f" ndf x exact)
+      $ eq 1e-5  exact  (cumulative (studentT ndf) x)
+    -- F-distribution
+    testFdistrPDF n m x exact
+      = testAssertion (printf "density (fDistribution %i %i) %f ≈ %f [got %f]" n m x exact d)
+      $ eq 1e-5  exact d
+      where d = density (fDistribution n m) x
+    testFdistrCDF n m x exact
+      = testAssertion (printf "cumulative (fDistribution %i %i) %f ≈ %f [got %f]" n m x exact d)
+      $ eq 1e-5  exact d
+      where d = cumulative (fDistribution n m) x
diff --git a/tests/Tests/Function.hs b/tests/Tests/Function.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/Function.hs
@@ -0,0 +1,26 @@
+module Tests.Function ( tests ) where
+
+import qualified Data.Vector.Unboxed as U
+import           Data.Vector.Unboxed   ((!))
+
+import Test.QuickCheck
+import Test.Framework
+import Test.Framework.Providers.QuickCheck2
+
+import Statistics.Function
+
+
+
+tests :: Test
+tests = testGroup "S.Function"
+  [ testProperty "Sort is sort" p_sort
+  ]
+
+
+p_sort :: [Double] -> Property
+p_sort xs =
+  not (null xs) ==> U.all (uncurry (<=)) (U.zip v $ U.tail v)
+    where
+      v = sort $ U.fromList xs
+
+      
diff --git a/tests/Tests/Helpers.hs b/tests/Tests/Helpers.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/Helpers.hs
@@ -0,0 +1,100 @@
+-- | Helpers for testing
+module Tests.Helpers (
+    -- * helpers
+    T(..)
+  , typeName
+  , eq
+  , eqC
+  , (=~)
+    -- * Generic QC tests
+  , monotonicallyIncreases
+  , monotonicallyIncreasesIEEE
+    -- * HUnit helpers
+  , testAssertion
+  , testEquality
+  ) where
+
+import Data.Complex
+import Data.Typeable
+
+import qualified Numeric.IEEE    as IEEE
+
+import qualified Test.HUnit      as HU
+import Test.Framework
+import Test.Framework.Providers.HUnit
+
+import Numeric.MathFunctions.Constants
+
+
+
+----------------------------------------------------------------
+-- Helpers
+----------------------------------------------------------------
+
+-- | Phantom typed value used to select right instance in QC tests
+data T a = T
+
+-- | String representation of type name
+typeName :: Typeable a => T a -> String
+typeName = show . typeOf . typeParam
+  where
+    typeParam :: T a -> a
+    typeParam _ = undefined
+
+-- | Approximate equality for 'Double'. Doesn't work well for numbers
+--   which are almost zero.
+eq :: Double                    -- ^ Relative error
+   -> Double -> Double -> Bool
+eq eps a b 
+  | a == 0 && b == 0 = True
+  | otherwise        = abs (a - b) <= eps * max (abs a) (abs b)
+
+-- | Approximate equality for 'Complex Double'
+eqC :: Double                   -- ^ Relative error
+    -> Complex Double
+    -> Complex Double
+    -> Bool
+eqC eps a@(ar :+ ai) b@(br :+ bi)
+  | a == 0 && b == 0 = True
+  | otherwise        = abs (ar - br) <= eps * d
+                    && abs (ai - bi) <= eps * d
+  where
+    d = max (realPart $ abs a) (realPart $ abs b)
+
+
+-- | Approximately equal up to 1 ulp
+(=~) :: Double -> Double -> Bool
+(=~) = eq m_epsilon
+
+
+----------------------------------------------------------------
+-- Generic QC
+----------------------------------------------------------------
+
+-- Check that function is nondecreasing
+monotonicallyIncreases :: (Ord a, Ord b) => (a -> b) -> a -> a -> Bool
+monotonicallyIncreases f x1 x2 = f (min x1 x2) <= f (max x1 x2)
+
+-- Check that function is nondecreasing taking rounding errors into
+-- account.
+--
+-- In fact funstion is allowed to decrease less than one ulp in order
+-- to guard againist problems with excess precision. On x86 FPU works
+-- with 80-bit numbers but doubles are 64-bit so rounding happens
+-- whenever values are moved from registers to memory
+monotonicallyIncreasesIEEE :: (Ord a, IEEE.IEEE b)  => (a -> b) -> a -> a -> Bool
+monotonicallyIncreasesIEEE f x1 x2 =
+  y1 <= y2 || (y1 - y2) < y2 * IEEE.epsilon
+  where
+    y1 = f (min x1 x2)
+    y2 = f (max x1 x2)
+
+----------------------------------------------------------------
+-- HUnit helpers
+----------------------------------------------------------------
+
+testAssertion :: String -> Bool -> Test
+testAssertion str cont = testCase str $ HU.assertBool str cont
+
+testEquality :: (Show a, Eq a) => String -> a -> a -> Test
+testEquality msg a b = testCase msg $ HU.assertEqual msg a b
diff --git a/tests/Tests/KDE.hs b/tests/Tests/KDE.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/KDE.hs
@@ -0,0 +1,45 @@
+-- | Tests for Kernel density estimates.
+module Tests.KDE ( 
+  tests 
+  )where
+
+import Data.Vector.Unboxed ((!))
+import qualified Data.Vector.Unboxed as U
+
+import Test.Framework                       (Test, testGroup)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.QuickCheck                      
+import Text.Printf
+
+import Statistics.Sample.KernelDensity
+
+
+
+tests :: Test
+tests = testGroup "KDE"
+  [ testProperty "integral(PDF) == 1" t_densityIsPDF
+  ]
+
+t_densityIsPDF :: [Double] -> Property
+t_densityIsPDF vec 
+  = not (null vec) ==> test
+  where
+    (xs,ys)  = kde 4096 (U.fromList vec)
+    step     = (xs ! 1) - (xs ! 0)
+    integral = integratePDF step ys
+    --
+    test = printTestCase (printf "Integral %f" integral)
+         $ abs (1 - integral) <= 1e-3
+
+          
+
+integratePDF :: Double -> U.Vector Double -> Double
+integratePDF step vec 
+  = step * U.sum (U.zipWith (*) vec weights)
+  where
+    n       = U.length vec
+    weights = U.generate n go
+      where
+        go i | i == 0    = 0.5
+             | i == n-1  = 0.5
+             | otherwise = 1
diff --git a/tests/Tests/NonparametricTest.hs b/tests/Tests/NonparametricTest.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/NonparametricTest.hs
@@ -0,0 +1,233 @@
+-- Tests for Statistics.Test.NonParametric
+module Tests.NonparametricTest (
+  nonparametricTests
+  ) where
+
+
+import qualified Data.Vector.Unboxed as U
+import Test.HUnit                     (assertEqual)
+import Test.Framework
+import Test.Framework.Providers.HUnit
+
+import Statistics.Test.MannWhitneyU
+import Statistics.Test.WilcoxonT
+
+import Tests.Helpers
+import Tests.NonparametricTest.Table
+
+import Statistics.Test.KolmogorovSmirnov
+import Statistics.Distribution.Normal    (standard)
+
+
+
+nonparametricTests :: Test
+nonparametricTests = testGroup "Nonparametric tests"
+                   $ concat [ mannWhitneyTests
+                            , wilcoxonSumTests
+                            , wilcoxonPairTests
+                            , kolmogorovSmirnovDTest
+                            ]
+
+----------------------------------------------------------------
+
+mannWhitneyTests :: [Test]
+mannWhitneyTests = zipWith test [(0::Int)..] testData ++
+  [ testEquality "Mann-Whitney U Critical Values, m=1"
+      (replicate (20*3) Nothing)
+      [mannWhitneyUCriticalValue (1,x) p | x <- [1..20], p <- [0.005,0.01,0.025]]
+  , testEquality "Mann-Whitney U Critical Values, m=2, p=0.025"
+      (replicate 7 Nothing ++ map Just [0,0,0,0,1,1,1,1,1,2,2,2,2])
+      [mannWhitneyUCriticalValue (2,x) 0.025 | x <- [1..20]]
+  , testEquality "Mann-Whitney U Critical Values, m=6, p=0.05"
+      (replicate 1 Nothing ++ map Just [0, 2,3,5,7,8,10,12,14,16,17,19,21,23,25,26,28,30,32])
+      [mannWhitneyUCriticalValue (6,x) 0.05 | x <- [1..20]]
+  , testEquality "Mann-Whitney U Critical Values, m=20, p=0.025"
+      (replicate 1 Nothing ++ map Just [2,8,14,20,27,34,41,48,55,62,69,76,83,90,98,105,112,119,127])
+      [mannWhitneyUCriticalValue (20,x) 0.025 | x <- [1..20]]
+  ]
+  where
+    test n (a, b, c, d)
+      = testCase "Mann-Whitney" $ do
+          assertEqual ("Mann-Whitney U "     ++ show n) c us
+          assertEqual ("Mann-Whitney U Sig " ++ show n) d ss
+      where
+        us = mannWhitneyU (U.fromList a) (U.fromList b)
+        ss = mannWhitneyUSignificant TwoTailed (length a, length b) 0.05 us
+    -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
+    testData :: [([Double], [Double], (Double, Double), Maybe TestResult)]
+    testData = [ ( [3,4,2,6,2,5]
+                 , [9,7,5,10,6,8]
+                 , (2, 34)
+                 , Just Significant
+                 )
+               , ( [540,480,600,590,605]
+                 , [760,890,1105,595,940]
+                 , (2, 23)
+                 , Just Significant
+                 )
+               , ( [19,22,16,29,24]
+                 , [20,11,17,12]
+                 , (17, 3)
+                 , Just NotSignificant
+                 )
+               , ( [126,148,85,61, 179,93, 45,189,85,93]
+                 , [194,128,69,135,171,149,89,248,79,137]
+                 , (35,65)
+                 , Just NotSignificant
+                 )
+               , ( [1..30]
+                 , [1..30]
+                 , (450,450)
+                 , Just NotSignificant
+                 )
+               , ( [1 .. 30]
+                 , [11.5 .. 40 ]
+                 , (190.0,710.0)
+                 , Just Significant
+                 )
+               ]
+
+wilcoxonSumTests :: [Test]
+wilcoxonSumTests = zipWith test [(0::Int)..] testData
+  where
+    test n (a, b, c) = testCase "Wilcoxon Sum"
+                     $ assertEqual ("Wilcoxon Sum " ++ show n) c (wilcoxonRankSums (U.fromList a) (U.fromList b))
+    -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
+    testData :: [([Double], [Double], (Double, Double))]
+    testData = [ ( [8.50,9.48,8.65,8.16,8.83,7.76,8.63]
+                 , [8.27,8.20,8.25,8.14,9.00,8.10,7.20,8.32,7.70]
+                 , (75, 61)
+                 )
+               , ( [0.45,0.50,0.61,0.63,0.75,0.85,0.93]
+                 , [0.44,0.45,0.52,0.53,0.56,0.58,0.58,0.65,0.79]
+                 , (71.5, 64.5)
+                 )
+               ]
+
+wilcoxonPairTests :: [Test]
+wilcoxonPairTests = zipWith test [(0::Int)..] testData ++
+  -- Taken from the Mitic paper:
+  [ testAssertion "Sig 16, 35" (to4dp 0.0467 $ wilcoxonMatchedPairSignificance 16 35)
+  , testAssertion "Sig 16, 36" (to4dp 0.0523 $ wilcoxonMatchedPairSignificance 16 36)
+  , testEquality   "Wilcoxon critical values, p=0.05"
+      (replicate 4 Nothing ++ map Just [0,2,3,5,8,10,13,17,21,25,30,35,41,47,53,60,67,75,83,91,100,110,119])
+      [wilcoxonMatchedPairCriticalValue x 0.05 | x <- [1..27]]
+  , testEquality "Wilcoxon critical values, p=0.025"
+      (replicate 5 Nothing ++ map Just [0,2,3,5,8,10,13,17,21,25,29,34,40,46,52,58,65,73,81,89,98,107])
+      [wilcoxonMatchedPairCriticalValue x 0.025 | x <- [1..27]]
+  , testEquality "Wilcoxon critical values, p=0.01"
+      (replicate 6 Nothing ++ map Just [0,1,3,5,7,9,12,15,19,23,27,32,37,43,49,55,62,69,76,84,92])
+      [wilcoxonMatchedPairCriticalValue x 0.01 | x <- [1..27]]
+  , testEquality "Wilcoxon critical values, p=0.005"
+      (replicate 7 Nothing ++ map Just [0,1,3,5,7,9,12,15,19,23,27,32,37,42,48,54,61,68,75,83])
+      [wilcoxonMatchedPairCriticalValue x 0.005 | x <- [1..27]]
+  ]
+  where
+    test n (a, b, c) = testEquality ("Wilcoxon Paired " ++ show n) c res
+      where res = (wilcoxonMatchedPairSignedRank (U.fromList a) (U.fromList b))
+
+    -- List of (Sample A, Sample B, (Positive Rank, Negative Rank))
+    testData :: [([Double], [Double], (Double, Double))]
+    testData = [ ([1..10], [1..10], (0, 0     ))
+               , ([1..5],  [6..10], (0, 5*(-3)))
+               -- Worked example from the Internet:
+               , ( [125,115,130,140,140,115,140,125,140,135]
+                 , [110,122,125,120,140,124,123,137,135,145]
+                 , ( sum $ filter (> 0) [7,-3,1.5,9,0,-4,8,-6,1.5,-5]
+                   , sum $ filter (< 0) [7,-3,1.5,9,0,-4,8,-6,1.5,-5]
+                   )
+                 )
+               -- Worked examples from books/papers:
+               , ( [2.4,1.9,2.3,1.9,2.4,2.5]
+                 , [2.0,2.1,2.0,2.0,1.8,2.0]
+                 , (18, -3)
+                 )
+               , ( [130,170,125,170,130,130,145,160]
+                 , [120,163,120,135,143,136,144,120]
+                 , (27, -9)
+                 )
+               , ( [540,580,600,680,430,740,600,690,605,520]
+                 , [760,710,1105,880,500,990,1050,640,595,520]
+                 , (3, -42)
+                 )
+               ]
+    to4dp tgt x = x >= tgt - 0.00005 && x < tgt + 0.00005
+
+
+
+----------------------------------------------------------------
+-- K-S test
+----------------------------------------------------------------
+
+
+kolmogorovSmirnovDTest :: [Test]
+kolmogorovSmirnovDTest =
+  [ testAssertion "K-S D statistics" $
+    and [ eq 1e-6 (kolmogorovSmirnovD standard (toU sample)) reference
+        | (reference,sample) <- tableKSD
+        ]
+  , testAssertion "K-S 2-sample statistics" $
+    and [ eq 1e-6 (kolmogorovSmirnov2D (toU xs) (toU ys)) reference
+        | (reference,xs,ys) <- tableKS2D
+        ]
+  , testAssertion "K-S probability" $
+    and [ eq 1e-14 (kolmogorovSmirnovProbability n d) p
+        | (d,n,p) <- testData
+        ]
+  ]
+  where
+    toU = U.fromList
+    -- Test data for the calculation of cumulative probability 
+    -- P(D[n] < d).
+    -- 
+    -- Test data is:
+    --    (D[n], n, p)
+    -- Table is generated using sample program from paper
+    testData :: [(Double,Int,Double)]
+    testData = 
+      [ (0.09           ,    3, 0                   )
+      , (0.2            ,    3, 0.00177777777777778 )
+      , (0.301          ,    3, 0.116357025777778   )
+      , (0.392          ,    3, 0.383127210666667   )
+      , (0.5003         ,    3, 0.667366306558667   )
+      , (0.604          ,    3, 0.861569877333333   )
+      , (0.699          ,    3, 0.945458198         )
+      , (0.802          ,    3, 0.984475216         )
+      , (0.9            ,    3, 0.998               )
+      , (0.09           ,    5, 0                   )
+      , (0.2            ,    5, 0.0384              )
+      , (0.301          ,    5, 0.33993786080016    )
+      , (0.392          ,    5, 0.66931908083712    )
+      , (0.5003         ,    5, 0.888397260183794   )
+      , (0.604          ,    5, 0.971609957879808   )
+      , (0.699          ,    5, 0.994331075994008   )
+      , (0.802          ,    5, 0.999391366368064   )
+      , (0.9            ,    5, 0.99998             )
+      , (0.09           ,    8, 3.37615237575e-06   )
+      , (0.2            ,    8, 0.151622071801758   )
+      , (0.301          ,    8, 0.613891042670582   )
+      , (0.392          ,    8, 0.871491561427005   )
+      , (0.5003         ,    8, 0.977534089199071   )
+      , (0.604          ,    8, 0.997473116268255   )
+      , (0.699          ,    8, 0.999806082005123   )
+      , (0.802          ,    8, 0.999995133786947   )
+      , (0.9            ,    8, 0.99999998          )
+      , (0.09           ,   10, 3.89639433093119e-05)
+      , (0.2            ,   10, 0.25128096          )
+      , (0.301          ,   10, 0.732913126355935   )
+      , (0.392          ,   10, 0.932185254518767   )
+      , (0.5003         ,   10, 0.992276179340446   )
+      , (0.604          ,   10, 0.999495533516769   )
+      , (0.699          ,   10, 0.999979691783985   )
+      , (0.802          ,   10, 0.999999801409237   )
+      , (0.09           ,   20, 0.00794502217168886 )
+      , (0.2            ,   20, 0.647279826376584   )
+      , (0.301          ,   20, 0.958017466965765   )
+      , (0.392          ,   20, 0.997206424843499   )
+      , (0.5003         ,   20, 0.999962641414228   )
+      , (0.09           ,   30, 0.0498147538075168  )
+      , (0.2            ,   30, 0.842030838984526   )
+      , (0.301          ,   30, 0.993403560017612   )
+      , (0.392          ,   30, 0.99988478803318    )
+      , (0.09           ,  100, 0.629367974413669   )
+      ]
diff --git a/tests/Tests/NonparametricTest/Table.hs b/tests/Tests/NonparametricTest/Table.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/NonparametricTest/Table.hs
@@ -0,0 +1,36 @@
+module Tests.NonparametricTest.Table where
+
+-- Table for Kolmogorov-Smirnov statistics for standard normal
+-- distribution. Generated using R.
+--
+-- First element of tuple is D second is sample for which it was
+-- calculated. 
+tableKSD :: [(Double,[Double])]
+tableKSD = 
+  [ (0.2012078,[1.360645,-0.3151904,-1.245443,0.1741977,-0.1421206,-1.798246,1.171594,-1.335844,-5.050093e-2,1.030063,-1.849005,0.6491455,-0.7028004])
+  , (0.2569956,[0.3884734,-1.227821,-0.4166262,0.429118,-0.9280124,0.8025867,-0.6703089,-0.2124872,0.1224496,0.1087734,-4.285284e-2,-1.039936,-0.7071956])
+  , (0.1960356,[-1.814745,-0.6327167,0.7082493,0.6264716,1.02061,-0.4094635,0.821026,-0.4255047,-0.4820728,-0.2239833,0.648517,1.114283,0.3610216])
+  , (0.2095746,[0.187011,0.1805498,0.4448389,0.6065506,0.2308673,0.5292549,-1.489902,-1.455191,0.5449396,-0.1436403,-0.7977073,-0.2693545,0.8260888,-1.474473,-2.158696e-2,-0.1455387])
+  , (0.1922603,[0.5772317,-1.255561,1.605823,0.4923361,0.2470848,1.176101,-0.3767689,-0.6896885,0.4509345,-0.5048447,0.9436534,1.025599,0.2998393,-3.415219e-2,1.264315,-1.44433,-1.646449e-2])
+  , (0.2173401,[1.812807,-0.8687497,-0.5710508,1.003647,1.142621,0.6546577,-6.083323e-3,1.628574e-2,1.067499,-1.953143,-0.6060077,1.90859,-0.7480553,0.6715162,-0.928759,1.862,1.604621,-0.2171044,-0.1835918])
+  , (0.2510541,[-0.4769572,1.062319,0.9952284,1.198086,1.015589,-0.4154523,-0.6711762,1.202902,0.2217098,5.381759e-2,0.6679715,0.2551287,-0.1371492])
+  , (0.1996022,[1.158607,-0.7354863,1.526559,-0.7855418,-2.82999,-0.6045106,-0.1830228,0.3306812,-0.819657,-1.223715,0.2536423,-0.4155781,1.447042])
+  , (0.2284761,[1.239965,0.8187093,0.5199788,1.172072,0.748259,1.869376e-2,0.1625921,-1.712065,0.7043582,-1.702702,-0.4792806,-0.1023351,0.1187189])
+  , (0.2337866,[0.9417261,-0.1024297,-0.7354359,1.099991,0.801984,-0.3745397,-1.749564,1.795771,1.099963,-0.605557,-2.035897,1.893603,-0.3468928,-0.2593938,2.100988,0.9665698,0.8757091,0.7696328,0.8730729,-0.3990352,2.04361,-0.4617864,-0.155021,2.15774,0.2687795,-0.9853512,-0.3264898,1.260026,4.267695,-0.5571145,0.6307067,-0.1691405,-1.730686])
+  , (0.3389167,[2.025542,-1.542641,-1.090238,3.99027,9.949129e-2,-0.8974433,-2.508418,6.390346,-2.675515,1.154459,1.688072,2.220727,-0.4743102])
+  , (0.4920231,[0.5192906,-3.260813,-1.245185,1.693084,3.561318,4.058924,2.27063,0.9446943,4.794159,-3.423733,0.8240817,0.644059,0.900175,1.932513,1.024586,2.82823,2.072192,-0.353231,-0.4319673,1.505952,1.0199,4.555054,2.364929,5.531467,3.279415,3.19821,2.726925,1.680027,-0.9041334,-0.8246765,-1.343979,8.454955,1.354581])
+  , (0.6727408,[-6.705672,-3.193988,-4.612611,-3.207994,-5.070172,-6.141169,-0.397149,-4.093359,-1.204801,-3.986585,-2.724662,0.9868107,-6.295266,-5.95839,-6.35114,-1.679555,-2.635889,-4.050329,1.557428,-2.548465,-0.9073924,-1.502018,-4.535688,-4.158818,-8.833434,-5.944697,-1.569672,-4.70399,-7.832059,-4.093708,-8.393417,-2.085432,-7.06495,-0.4230419,-3.046822,-3.23895,-0.9265873,-9.227822,3.293713,-5.593577,-5.942398,-4.358421,2.660044,-4.301572,-1.258879,0.1499903,3.572833,-3.19844,0.8652432,-0.3025793,-1.576673,-7.666265,-6.751463,-1.398944,-2.690656,-1.429654,7.508364e-2,0.7998344,-3.562074,-1.021431,1.342968,2.110244,-7.561497,-2.372083,-3.649193,-5.7723,-1.068083,0.7537809,-4.569546,-1.198005,-5.638384,-1.227226,-1.195852,-1.118175,-9.130527,0.9675821,-2.497391,0.5988562,-1.965783,-4.25741,-6.547006,-1.459294,-2.380556,-3.977307,-7.809006,-4.276819,-4.028746,-9.055546,-3.599239,-1.470512,-8.253329,-1.351687,-4.269324,-6.140353,-6.30808,-1.834091,-3.135146,-9.391791,3.117815,-5.554733,-2.556769,-3.287376,-2.064013,-5.741995,-5.047918,-4.808841,-1.488526,-0.2351115,-5.760833,-2.722929,-7.012353,2.281171,-3.890514,-1.516824,-1.41011,-1.828457,-5.561244,-3.472142,-10.16919,-0.4369042,-5.698953,-4.587462,-4.897086])
+  ]
+
+-- Table for 2-sample Kolmogorov-Smirnov statistics. Generated using R
+--
+-- First element is D, second and third are samples
+tableKS2D :: [(Double,[Double],[Double])]
+tableKS2D =
+  [ (0.2820513,[-0.4212928,2.146532,0.7585263,-0.5086105,-0.7725486,6.235548e-2,-0.1849861,0.861972,-0.1958534,-3.379697e-2,-1.316854,0.6701269],[0.4957582,0.4241167,0.9822869,0.4504248,-0.1749617,1.178098,-1.117222,-0.859273,0.3073736,0.4344583,-0.4761338,-1.332374,1.487291])
+  , (0.2820513,[-0.712252,0.7990333,-0.7968473,1.443609,1.163096,-1.349071,-0.1553941,-2.003104,-0.3400618,-0.7019282,0.183293,-0.2352167],[-0.4622455,-0.8132221,0.1161614,-1.472115e-2,1.001454,-6.557789e-2,-0.2531216,-1.032432,0.4105478,1.749614,0.9722899,5.850942e-2,-0.3352746])
+  , (0.2564103,[0.3509882,-0.2982833,1.314731,1.264223,-0.8156374,0.3734029,-3.288915e-2,0.6766016,0.9786335,0.1079949,-0.4211722,1.58656],[0.8024675,7.464538e-2,0.2739861,-2.334255e-2,0.5611802,0.6683374,0.4358206,0.349843,1.207834,1.402578,-0.4049183,0.4286042,1.665129])
+  , (0.1833333,[1.376196,9.926384e-2,2.199292,-2.04993,0.5585353,-0.4812132,0.1041527,2.084774,0.71194,-1.398245,-4.458574e-2,1.484945,-1.473182,1.020076,-0.7019646,0.2182066,-1.702963,-0.3522622,-0.8129267,-0.6338972],[-1.020371,0.3323861,1.513288,0.1958708,-1.0723,5.323446e-2,-0.9993713,-0.7046356,-0.6781067,-0.4471603,1.512042,-0.2650665,-4.765228e-2,-1.501205,1.228664,0.5332935,-0.2960315,-0.1509683])
+  , (0.5666667,[0.7145305,0.1255674,2.001531,0.1419216],[2.113474,-0.3352839,-0.4962429,-1.386079,0.6404667,-0.7145304,0.1084008,-0.9821421,-2.270472,-1.003846,-0.5644588,2.699695,-1.296494,-0.1538839,1.319094,-1.127544,0.3568889,0.2004726,-1.313291,0.3581084,0.3313498,0.9336278,0.9850203,-1.309506,1.170459,-0.7517466,-1.771269,0.7156381,-1.129691,0.877729])
+  , (0.5,[0.6950626,0.1643805,-0.3102472,0.4810762,0.1844602,1.338836,-0.8083386,-0.5482141,0.9532421,-0.2644837],[7.527945,-1.95654,1.513725,-1.318431,2.453895,0.2078194,0.7371092,2.834245,-2.134794,3.938259])
+  ] 
diff --git a/tests/Tests/Transform.hs b/tests/Tests/Transform.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/Transform.hs
@@ -0,0 +1,145 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts  #-}
+{-# LANGUAGE ViewPatterns      #-}
+module Tests.Transform
+    (
+      tests
+    ) where
+
+import Data.Bits             ((.&.), shiftL)
+import Data.Complex          (Complex((:+)))
+import Data.Functor          ((<$>))
+import Statistics.Function   (within)
+import Statistics.Transform
+
+import Test.Framework                       (Test, testGroup)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.QuickCheck                      (Positive(..),Property,Arbitrary(..),Gen,choose,vectorOf,
+                                             printTestCase, quickCheck)
+
+import Text.Printf
+
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
+
+import Tests.Helpers
+
+
+
+tests :: Test
+tests = testGroup "fft" [
+          testProperty "t_impulse"        t_impulse
+        , testProperty "t_impulse_offset" t_impulse_offset
+        , testProperty "ifft . fft = id"  (t_fftInverse $ ifft . fft)
+        , testProperty "fft . ifft = id"  (t_fftInverse $ fft . ifft)
+        , testProperty "idct . dct = id [up to scale]"
+            (t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))
+        , testProperty "dct . idct = id [up to scale]"
+            (t_fftInverse (\v -> U.map (/ (2 * fromIntegral (U.length v))) $ idct $ dct v))
+          -- Exact small size DCT
+          -- 2
+        , testDCT [1,0] $ map (*2) [1, cos (pi/4)   ]
+        , testDCT [0,1] $ map (*2) [1, cos (3*pi/4) ]
+          -- 4
+        , testDCT [1,0,0,0] $ map (*2) [1, cos(  pi/8), cos( 2*pi/8), cos( 3*pi/8)]
+        , testDCT [0,1,0,0] $ map (*2) [1, cos(3*pi/8), cos( 6*pi/8), cos( 9*pi/8)]
+        , testDCT [0,0,1,0] $ map (*2) [1, cos(5*pi/8), cos(10*pi/8), cos(15*pi/8)]
+        , testDCT [0,0,0,1] $ map (*2) [1, cos(7*pi/8), cos(14*pi/8), cos(21*pi/8)]
+          -- Exact small size IDCT
+          -- 2
+        , testIDCT [1,0]            [1,         1          ]
+        , testIDCT [0,1] $ map (*2) [cos(pi/4), cos(3*pi/4)]
+          -- 4
+        , testIDCT [1,0,0,0]            [1,            1,            1,            1            ]
+        , testIDCT [0,1,0,0] $ map (*2) [cos(   pi/8), cos( 3*pi/8), cos( 5*pi/8), cos( 7*pi/8) ]
+        , testIDCT [0,0,1,0] $ map (*2) [cos( 2*pi/8), cos( 6*pi/8), cos(10*pi/8), cos(14*pi/8) ]
+        , testIDCT [0,0,0,1] $ map (*2) [cos( 3*pi/8), cos( 9*pi/8), cos(15*pi/8), cos(21*pi/8) ]
+        ]
+
+-- A single real-valued impulse at the beginning of an otherwise zero
+-- vector should be replicated in every real component of the result,
+-- and all the imaginary components should be zero.
+t_impulse :: Double -> Positive Int -> Bool
+t_impulse k (Positive m) = G.all (c_near i) (fft v)
+  where v = i `G.cons` G.replicate (n-1) 0
+        i = k :+ 0
+        n = 1 `shiftL` (m .&. 6)
+
+-- If a real-valued impulse is offset from the beginning of an
+-- otherwise zero vector, the sum-of-squares of each component of the
+-- result should equal the square of the impulse.
+t_impulse_offset :: Double -> Positive Int -> Positive Int -> Bool
+t_impulse_offset k (Positive x) (Positive m) = G.all ok (fft v)
+  where v = G.concat [G.replicate xn 0, G.singleton i, G.replicate (n-xn-1) 0]
+        ok (re :+ im) = within ulps (re*re + im*im) (k*k)
+        i  = k :+ 0
+        xn = x `rem` n
+        n  = 1 `shiftL` (m .&. 6)
+
+-- Test that (ifft . fft ≈ id)
+--
+-- Approximate equality here is tricky. Smaller values of vector tend
+-- to have large relative error. Thus we should test that vectors as
+-- whole are approximate equal.
+t_fftInverse :: (HasNorm (U.Vector a), U.Unbox a, Num a, Show a, Arbitrary a)
+             => (U.Vector a -> U.Vector a) -> Property
+t_fftInverse roundtrip = do
+  x <- genFftVector
+  let n  = G.length x
+      x' = roundtrip x
+      d  = G.zipWith (-) x x'
+      nd = vectorNorm d
+      nx = vectorNorm x
+  id $ printTestCase "Original vector"
+     $ printTestCase (show x )
+     $ printTestCase "Transformed one"
+     $ printTestCase (show x')
+     $ printTestCase (printf "Length = %i" n)
+     $ printTestCase (printf "|x - x'| / |x| = %.6g" (nd / nx))
+     $ nd <= 3e-14 * nx
+
+-- Test discrete cosine transform
+testDCT :: [Double] -> [Double] -> Test
+testDCT (U.fromList -> vec) (U.fromList -> res)
+  = testAssertion ("DCT test for " ++ show vec)
+  $ vecEqual 3e-14 (dct vec) res
+
+-- Test inverse discrete cosine transform
+testIDCT :: [Double] -> [Double] -> Test
+testIDCT (U.fromList -> vec) (U.fromList -> res)
+  = testAssertion ("IDCT test for " ++ show vec)
+  $ vecEqual 3e-14 (idct vec) res
+
+
+
+----------------------------------------------------------------
+
+-- With an error tolerance of 8 ULPs, a million QuickCheck tests are
+-- likely to all succeed. With a tolerance of 7, we fail around the
+-- half million mark.
+ulps :: Int
+ulps = 8
+
+c_near :: CD -> CD -> Bool
+c_near (a :+ b) (c :+ d) = within ulps a c && within ulps b d
+
+-- Arbitrary vector for FFT od DCT
+genFftVector :: (U.Unbox a, Arbitrary a) => Gen (U.Vector a)
+genFftVector = do
+  n <- (2^)  <$> choose (1,9::Int)    -- Size of vector
+  G.fromList <$> vectorOf n arbitrary -- Vector to transform
+
+-- Ad-hoc type class for calculation of vector norm
+class HasNorm a where
+  vectorNorm :: a -> Double
+
+instance HasNorm (U.Vector Double) where
+  vectorNorm = sqrt . U.sum . U.map (\x -> x*x)
+
+instance HasNorm (U.Vector CD) where
+  vectorNorm = sqrt . U.sum . U.map (\(x :+ y) -> x*x + y*y)
+
+-- Approximate equality for vectors
+vecEqual :: Double -> U.Vector Double -> U.Vector Double -> Bool
+vecEqual ε v u
+  = vectorNorm (U.zipWith (-) v u) < ε * vectorNorm v
diff --git a/tests/tests.hs b/tests/tests.hs
new file mode 100644
--- /dev/null
+++ b/tests/tests.hs
@@ -0,0 +1,15 @@
+import Test.Framework       (defaultMain)
+
+import Tests.Distribution
+import Tests.NonparametricTest
+import qualified Tests.Transform
+import qualified Tests.Function
+import qualified Tests.KDE
+
+main :: IO ()
+main = defaultMain [ distributionTests 
+                   , nonparametricTests
+                   , Tests.Transform.tests
+                   , Tests.Function.tests
+                   , Tests.KDE.tests
+                   ]