diff --git a/ChangeLog b/ChangeLog
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,18 @@
+-*- text -*-
+
+Changes in 0.10.6.0
+
+  * The Estimator type has become an algebraic data type.  This allows
+    the jackknife function to potentially use more efficient jackknife
+    implementations.
+
+  * jackknifeMean, jackknifeStdDev, jackknifeVariance,
+    jackknifeVarianceUnb: new functions.  These have O(n) cost instead
+    of the O(n^2) cost of the standard jackknife.
+
+  * The mean function has been renamed to welfordMean; a new
+    implementation of mean has better numerical accuracy in almost all
+    cases.
 
 Changes in 0.10.5.2
 
diff --git a/Statistics/Autocorrelation.hs b/Statistics/Autocorrelation.hs
--- a/Statistics/Autocorrelation.hs
+++ b/Statistics/Autocorrelation.hs
@@ -17,7 +17,9 @@
     , autocorrelation
     ) where
 
+import Prelude hiding (sum)
 import Statistics.Sample (mean)
+import Statistics.Sample.Internal (sum)
 import qualified Data.Vector.Generic as G
 
 -- | Compute the autocovariance of a sample, i.e. the covariance of
@@ -25,7 +27,7 @@
 autocovariance :: (G.Vector v Double, G.Vector v Int) => v Double -> v Double
 autocovariance a = G.map f . G.enumFromTo 0 $ l-2
   where
-    f k = G.sum (G.zipWith (*) (G.take (l-k) c) (G.slice k (l-k) c))
+    f k = sum (G.zipWith (*) (G.take (l-k) c) (G.slice k (l-k) c))
           / fromIntegral l
     c   = G.map (subtract (mean a)) a
     l   = G.length a
diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs
--- a/Statistics/Distribution.hs
+++ b/Statistics/Distribution.hs
@@ -32,14 +32,14 @@
     , sumProbabilities
     ) where
 
-import Control.Applicative     ((<$>), Applicative(..))
+import Control.Applicative ((<$>), Applicative(..))
 import Control.Monad.Primitive (PrimMonad,PrimState)
-
+import Prelude hiding (sum)
+import Statistics.Sample.Internal (sum)
+import System.Random.MWC (Gen, uniform)
 import qualified Data.Vector.Unboxed as U
-import System.Random.MWC
 
 
-
 -- | Type class common to all distributions. Only c.d.f. could be
 -- defined for both discrete and continous distributions.
 class Distribution d where
@@ -141,7 +141,7 @@
 class (Distribution d) => MaybeEntropy d where
   -- | Returns the entropy of a distribution, in nats, if such is defined.
   maybeEntropy :: d -> Maybe Double
-  
+
 -- | Type class for distributions with entropy, meaning Shannon
 --   entropy in the case of a discrete distribution, or differential
 --   entropy in the case of a continuous one.  If the distribution has
@@ -178,7 +178,7 @@
 -- bisection with the given guess as a starting point.  The upper and
 -- lower bounds specify the interval in which the probability
 -- distribution reaches the value /p/.
-findRoot :: ContDistr d => 
+findRoot :: ContDistr d =>
             d                   -- ^ Distribution
          -> Double              -- ^ Probability /p/
          -> Double              -- ^ Initial guess
@@ -207,7 +207,7 @@
 -- | Sum probabilities in inclusive interval.
 sumProbabilities :: DiscreteDistr d => d -> Int -> Int -> Double
 sumProbabilities d low hi =
-  -- Return value is forced to be less than 1 to guard againist roundoff errors. 
+  -- Return value is forced to be less than 1 to guard againist roundoff errors.
   -- ATTENTION! this check should be removed for testing or it could mask bugs.
-  min 1 . U.sum . U.map (probability d) $ U.enumFromTo low hi
+  min 1 . sum . U.map (probability d) $ U.enumFromTo low hi
 {-# INLINE sumProbabilities #-}
diff --git a/Statistics/Distribution/Beta.hs b/Statistics/Distribution/Beta.hs
--- a/Statistics/Distribution/Beta.hs
+++ b/Statistics/Distribution/Beta.hs
@@ -27,6 +27,8 @@
   incompleteBeta, invIncompleteBeta, logBeta, digamma)
 import Numeric.MathFunctions.Constants (m_NaN)
 import qualified Statistics.Distribution as D
+import Data.Binary (put, get)
+import Control.Applicative ((<$>), (<*>))
 
 -- | The beta distribution
 data BetaDistribution = BD
@@ -36,7 +38,9 @@
    -- ^ Beta shape parameter
  } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary BetaDistribution
+instance Binary BetaDistribution where
+    put (BD x y) = put x >> put y
+    get = BD <$> get <*> get
 
 -- | Create beta distribution. Both shape parameters must be positive.
 betaDistr :: Double             -- ^ Shape parameter alpha
@@ -85,12 +89,12 @@
 
 instance D.Entropy BetaDistribution where
   entropy (BD a b) =
-    logBeta a b 
+    logBeta a b
     - (a-1) * digamma a
     - (b-1) * digamma b
     + (a+b-2) * digamma (a+b)
   {-# INLINE entropy #-}
-    
+
 instance D.MaybeEntropy BetaDistribution where
   maybeEntropy = Just . D.entropy
   {-# INLINE maybeEntropy #-}
diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs
--- a/Statistics/Distribution/Binomial.hs
+++ b/Statistics/Distribution/Binomial.hs
@@ -30,6 +30,8 @@
 import qualified Statistics.Distribution.Poisson.Internal as I
 import Numeric.SpecFunctions (choose,incompleteBeta)
 import Numeric.MathFunctions.Constants (m_epsilon)
+import Data.Binary (put, get)
+import Control.Applicative ((<$>), (<*>))
 
 
 -- | The binomial distribution.
@@ -40,7 +42,9 @@
     -- ^ Probability.
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary BinomialDistribution
+instance Binary BinomialDistribution where
+    put (BD x y) = put x >> put y
+    get = BD <$> get <*> get
 
 instance D.Distribution BinomialDistribution where
     cumulative = cumulative
@@ -99,7 +103,7 @@
 {-# INLINE variance #-}
 
 directEntropy :: BinomialDistribution -> Double
-directEntropy d@(BD n _) =   
+directEntropy d@(BD n _) =
   negate . sum $
   takeWhile (< negate m_epsilon) $
   dropWhile (not . (< negate m_epsilon)) $
diff --git a/Statistics/Distribution/CauchyLorentz.hs b/Statistics/Distribution/CauchyLorentz.hs
--- a/Statistics/Distribution/CauchyLorentz.hs
+++ b/Statistics/Distribution/CauchyLorentz.hs
@@ -25,6 +25,8 @@
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
 import qualified Statistics.Distribution as D
+import Data.Binary (put, get)
+import Control.Applicative ((<$>), (<*>))
 
 -- | Cauchy-Lorentz distribution.
 data CauchyDistribution = CD {
@@ -39,7 +41,9 @@
   }
   deriving (Eq, Show, Read, Typeable, Data, Generic)
 
-instance Binary CauchyDistribution
+instance Binary CauchyDistribution where
+    put (CD x y) = put x >> put y
+    get = CD <$> get <*> get
 
 -- | Cauchy distribution
 cauchyDistribution :: Double    -- ^ Central point
@@ -72,6 +76,6 @@
 
 instance D.Entropy CauchyDistribution where
   entropy (CD _ s) = log s + log (4*pi)
-    
+
 instance D.MaybeEntropy CauchyDistribution where
   maybeEntropy = Just . D.entropy
diff --git a/Statistics/Distribution/ChiSquared.hs b/Statistics/Distribution/ChiSquared.hs
--- a/Statistics/Distribution/ChiSquared.hs
+++ b/Statistics/Distribution/ChiSquared.hs
@@ -26,13 +26,16 @@
 
 import qualified Statistics.Distribution         as D
 import qualified System.Random.MWC.Distributions as MWC
+import Data.Binary (put, get)
 
 
 -- | Chi-squared distribution
 newtype ChiSquared = ChiSquared Int
                      deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary ChiSquared
+instance Binary ChiSquared where
+    get = fmap ChiSquared get
+    put (ChiSquared x) = put x
 
 -- | Get number of degrees of freedom
 chiSquaredNDF :: ChiSquared -> Int
@@ -69,12 +72,12 @@
 instance D.MaybeVariance ChiSquared where
     maybeStdDev   = Just . D.stdDev
     maybeVariance = Just . D.variance
-    
+
 instance D.Entropy ChiSquared where
   entropy (ChiSquared ndf) =
     let kHalf = 0.5 * fromIntegral ndf in
-    kHalf 
-    + log 2 
+    kHalf
+    + log 2
     + logGamma kHalf
     + (1-kHalf) * digamma kHalf
 
diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs
--- a/Statistics/Distribution/Exponential.hs
+++ b/Statistics/Distribution/Exponential.hs
@@ -31,13 +31,16 @@
 import qualified Statistics.Sample               as S
 import qualified System.Random.MWC.Distributions as MWC
 import Statistics.Types (Sample)
+import Data.Binary (put, get)
 
 
 newtype ExponentialDistribution = ED {
       edLambda :: Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary ExponentialDistribution
+instance Binary ExponentialDistribution where
+    put = put . edLambda
+    get = fmap ED get
 
 instance D.Distribution ExponentialDistribution where
     cumulative      = cumulative
diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs
--- a/Statistics/Distribution/FDistribution.hs
+++ b/Statistics/Distribution/FDistribution.hs
@@ -23,6 +23,8 @@
 import qualified Statistics.Distribution as D
 import Numeric.SpecFunctions (
   logBeta, incompleteBeta, invIncompleteBeta, digamma)
+import Data.Binary (put, get)
+import Control.Applicative ((<$>), (<*>))
 
 
 
@@ -33,7 +35,9 @@
                        }
                    deriving (Eq, Show, Read, Typeable, Data, Generic)
 
-instance Binary FDistribution
+instance Binary FDistribution where
+    get = F <$> get <*> get <*> get
+    put (F x y z) = put x >> put y >> put z
 
 fDistribution :: Int -> Int -> FDistribution
 fDistribution n m
@@ -89,9 +93,9 @@
   entropy (F n m _) =
     let nHalf = 0.5 * n
         mHalf = 0.5 * m in
-    log (n/m) 
+    log (n/m)
     + logBeta nHalf mHalf
-    + (1 - nHalf) * digamma nHalf 
+    + (1 - nHalf) * digamma nHalf
     - (1 + mHalf) * digamma mHalf
     + (nHalf + mHalf) * digamma (nHalf + mHalf)
 
diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs
--- a/Statistics/Distribution/Gamma.hs
+++ b/Statistics/Distribution/Gamma.hs
@@ -25,14 +25,15 @@
     , gdScale
     ) where
 
+import Control.Applicative ((<$>), (<*>))
 import Data.Binary (Binary)
+import Data.Binary (put, get)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
 import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN, m_neg_inf)
-import Numeric.SpecFunctions (
-  incompleteGamma, invIncompleteGamma, logGamma, digamma)
-import Statistics.Distribution.Poisson.Internal  as Poisson
-import qualified Statistics.Distribution         as D
+import Numeric.SpecFunctions (incompleteGamma, invIncompleteGamma, logGamma, digamma)
+import Statistics.Distribution.Poisson.Internal as Poisson
+import qualified Statistics.Distribution as D
 import qualified System.Random.MWC.Distributions as MWC
 
 -- | The gamma distribution.
@@ -41,7 +42,9 @@
     , gdScale :: {-# UNPACK #-} !Double -- ^ Scale parameter, &#977;.
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary GammaDistribution
+instance Binary GammaDistribution where
+    put (GD x y) = put x >> put y
+    get = GD <$> get <*> get
 
 -- | Create gamma distribution. Both shape and scale parameters must
 -- be positive.
@@ -90,11 +93,11 @@
 
 instance D.MaybeEntropy GammaDistribution where
   maybeEntropy (GD a l)
-    | a > 0 && l > 0 = 
+    | a > 0 && l > 0 =
       Just $
-      a 
-      + log l 
-      + logGamma a 
+      a
+      + log l
+      + logGamma a
       + (1-a) * digamma a
     | otherwise = Nothing
 
diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs
--- a/Statistics/Distribution/Geometric.hs
+++ b/Statistics/Distribution/Geometric.hs
@@ -30,12 +30,14 @@
     , gdSuccess0
     ) where
 
-import Control.Monad  (liftM)
-import Data.Binary    (Binary)
-import Data.Data      (Data, Typeable)
-import GHC.Generics   (Generic)
-import Numeric.MathFunctions.Constants(m_pos_inf,m_neg_inf)
-import qualified Statistics.Distribution         as D
+import Control.Applicative ((<$>))
+import Control.Monad (liftM)
+import Data.Binary (Binary)
+import Data.Binary (put, get)
+import Data.Data (Data, Typeable)
+import GHC.Generics (Generic)
+import Numeric.MathFunctions.Constants (m_pos_inf, m_neg_inf)
+import qualified Statistics.Distribution as D
 import qualified System.Random.MWC.Distributions as MWC
 
 ----------------------------------------------------------------
@@ -45,7 +47,9 @@
       gdSuccess :: Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary GeometricDistribution
+instance Binary GeometricDistribution where
+    get = GD <$> get
+    put (GD x) = put x
 
 instance D.Distribution GeometricDistribution where
     cumulative = cumulative
@@ -115,7 +119,9 @@
       gdSuccess0 :: Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary GeometricDistribution0
+instance Binary GeometricDistribution0 where
+    get = GD0 <$> get
+    put (GD0 x) = put x
 
 instance D.Distribution GeometricDistribution0 where
     cumulative (GD0 s) x = cumulative (GD s) (x + 1)
diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs
--- a/Statistics/Distribution/Hypergeometric.hs
+++ b/Statistics/Distribution/Hypergeometric.hs
@@ -33,6 +33,8 @@
 import Numeric.MathFunctions.Constants (m_epsilon)
 import Numeric.SpecFunctions (choose)
 import qualified Statistics.Distribution as D
+import Data.Binary (put, get)
+import Control.Applicative ((<$>), (<*>))
 
 data HypergeometricDistribution = HD {
       hdM :: {-# UNPACK #-} !Int
@@ -40,7 +42,9 @@
     , hdK :: {-# UNPACK #-} !Int
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary HypergeometricDistribution
+instance Binary HypergeometricDistribution where
+    get = HD <$> get <*> get <*> get
+    put (HD x y z) = put x >> put y >> put z
 
 instance D.Distribution HypergeometricDistribution where
     cumulative = cumulative
@@ -63,7 +67,7 @@
 
 instance D.Entropy HypergeometricDistribution where
   entropy = directEntropy
-  
+
 instance D.MaybeEntropy HypergeometricDistribution where
   maybeEntropy = Just . D.entropy
 
diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs
--- a/Statistics/Distribution/Normal.hs
+++ b/Statistics/Distribution/Normal.hs
@@ -20,17 +20,18 @@
     , standard
     ) where
 
+import Control.Applicative ((<$>), (<*>))
 import Data.Binary (Binary)
+import Data.Binary (put, get)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
 import Numeric.MathFunctions.Constants (m_sqrt_2, m_sqrt_2_pi)
-import Numeric.SpecFunctions           (erfc, invErfc)
+import Numeric.SpecFunctions (erfc, invErfc)
 import qualified Statistics.Distribution as D
-import qualified Statistics.Sample       as S
+import qualified Statistics.Sample as S
 import qualified System.Random.MWC.Distributions as MWC
 
 
-
 -- | The normal distribution.
 data NormalDistribution = ND {
       mean       :: {-# UNPACK #-} !Double
@@ -39,7 +40,9 @@
     , ndCdfDenom :: {-# UNPACK #-} !Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary NormalDistribution
+instance Binary NormalDistribution where
+    put (ND w x y z) = put w >> put x >> put y >> put z
+    get = ND <$> get <*> get <*> get <*> get
 
 instance D.Distribution NormalDistribution where
     cumulative      = cumulative
diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs
--- a/Statistics/Distribution/Poisson.hs
+++ b/Statistics/Distribution/Poisson.hs
@@ -31,13 +31,16 @@
 import qualified Statistics.Distribution.Poisson.Internal as I
 import Numeric.SpecFunctions (incompleteGamma,logFactorial)
 import Numeric.MathFunctions.Constants (m_neg_inf)
+import Data.Binary (put, get)
 
 
 newtype PoissonDistribution = PD {
       poissonLambda :: Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary PoissonDistribution
+instance Binary PoissonDistribution where
+    get = fmap PD get
+    put = put . poissonLambda
 
 instance D.Distribution PoissonDistribution where
     cumulative (PD lambda) x
@@ -78,8 +81,9 @@
 poisson :: Double -> PoissonDistribution
 poisson l
   | l >=  0   = PD l
-  | otherwise = error $ "Statistics.Distribution.Poisson.poisson:\
-                        \ lambda must be non-negative. Got " ++ show l
+  | otherwise = error $
+    "Statistics.Distribution.Poisson.poisson: lambda must be non-negative. Got "
+    ++ show l
 {-# INLINE poisson #-}
 
 -- $references
diff --git a/Statistics/Distribution/Poisson/Internal.hs b/Statistics/Distribution/Poisson/Internal.hs
--- a/Statistics/Distribution/Poisson/Internal.hs
+++ b/Statistics/Distribution/Poisson/Internal.hs
@@ -14,10 +14,10 @@
       probability, poissonEntropy
     ) where
 
-import Data.List(unfoldr)
+import Data.List (unfoldr)
 import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny, m_epsilon)
 import Numeric.SpecFunctions (logGamma, stirlingError, choose, logFactorial)
-import Numeric.SpecFunctions.Extra     (bd0)
+import Numeric.SpecFunctions.Extra (bd0)
 
 -- | An unchecked, non-integer-valued version of Loader's saddle point
 -- algorithm.
@@ -50,7 +50,7 @@
   where parity j
           | even (k-j) = -1
           | otherwise  = 1
-                         
+
 -- | Returns [x, x^2, x^3, x^4, ...]
 powers :: Double -> [Double]
 powers x = unfoldr (\y -> Just (y*x,y*x)) 1
@@ -68,7 +68,7 @@
 alyThm2 lambda upper lower =
   alyThm2Upper lambda upper + 0.5 * (zipCoefficients lambda lower)
 
-zipCoefficients :: Double -> [Double] -> Double 
+zipCoefficients :: Double -> [Double] -> Double
 zipCoefficients lambda coefficients =
   (sum $ map (uncurry (*)) (zip (powers $ recip lambda) coefficients))
 
@@ -76,12 +76,12 @@
 --
 -- poissonMoment[0, s_] := 1
 -- poissonMoment[1, s_] := 0
--- poissonMoment[k_, s_] := 
+-- poissonMoment[k_, s_] :=
 --   Sum[s * Binomial[k - 1, j] * poissonMoment[j, s], {j, 0, k - 2}]
 --
 -- upperSeries[m_]  :=
 --  Distribute[Integrate[
---    Sum[(-1)^(j - 1) * 
+--    Sum[(-1)^(j - 1) *
 --      poissonMoment[j, \[Lambda]] / (j * (j - 1)* \[Lambda]^j),
 --     {j, 3, 2 m - 1}],
 --    \[Lambda]]]
@@ -89,10 +89,10 @@
 -- lowerSeries[m_] :=
 --  Distribute[Integrate[
 --    poissonMoment[
---      2 m + 2, \[Lambda]] / ((2 m + 
+--      2 m + 2, \[Lambda]] / ((2 m +
 --         1)*\[Lambda]^(2 m + 2)), \[Lambda]]]
 --
--- upperBound[m_] := upperSeries[m] + (Log[2*Pi*\[Lambda]] + 1)/2 
+-- upperBound[m_] := upperSeries[m] + (Log[2*Pi*\[Lambda]] + 1)/2
 --
 -- lowerBound[m_] := upperBound[m] + lowerSeries[m]
 
@@ -120,7 +120,7 @@
 lowerCoefficients8 = [0,0,0,0,0,0,0, -2027025/8, -15315300, -105252147,
                       -178127950, -343908565/4, -10929270, -3721149/14,
                       -7709/15, -1/272]
-  
+
 upperCoefficients10 :: [Double]
 upperCoefficients10 = [1/12, 1/24, 19/360, 9,80, 863/2520, 1375/1008,
                        33953/5040, 57281/1440, -2271071617/11880,
@@ -128,13 +128,13 @@
                        -7531072742237/131040, -1405507544003/65520,
                        -21001919627/10080, -1365808297/36720,
                        -26059/544, -1/5814]
-                      
+
 lowerCoefficients10 :: [Double]
 lowerCoefficients10 = [0,0,0,0,0,0,0,0,0,-130945815/2, -7638505875,
                        -438256243425/4, -435477637540, -3552526473925/6,
                        -857611717105/3, -545654955967/12, -5794690528/3,
                        -578334559/42, -699043/133, -1/420]
-                 
+
 upperCoefficients12 :: [Double]
 upperCoefficients12 = [1/12, 1/24, 19/360, 863/2520, 1375/1008,
                        33953/5040, 57281/1440, 3250433/11880,
@@ -153,13 +153,13 @@
                        -347362037754732, -2205885452434521/100,
                        -12237195698286/35, -16926981721/22,
                        -6710881/155, -1/600]
-                      
+
 -- | Compute entropy directly from its definition. This is just as accurate
 -- as 'alyThm1' for lambda <= 1 and is faster, but is slow for large lambda,
 -- and produces some underestimation due to accumulation of floating point
 -- error.
 directEntropy :: Double -> Double
-directEntropy lambda =   
+directEntropy lambda =
   negate . sum $
   takeWhile (< negate m_epsilon * lambda) $
   dropWhile (not . (< negate m_epsilon * lambda)) $
diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs
--- a/Statistics/Distribution/StudentT.hs
+++ b/Statistics/Distribution/StudentT.hs
@@ -23,12 +23,15 @@
 import Statistics.Distribution.Transform (LinearTransform (..))
 import Numeric.SpecFunctions (
   logBeta, incompleteBeta, invIncompleteBeta, digamma)
+import Data.Binary (put, get)
 
 -- | Student-T distribution
 newtype StudentT = StudentT { studentTndf :: Double }
                    deriving (Eq, Show, Read, Typeable, Data, Generic)
 
-instance Binary StudentT
+instance Binary StudentT where
+    put = put . studentTndf
+    get = fmap StudentT get
 
 -- | Create Student-T distribution. Number of parameters must be positive.
 studentT :: Double -> StudentT
@@ -76,7 +79,7 @@
 instance D.Entropy StudentT where
   entropy (StudentT ndf) =
     0.5 * (ndf+1) * (digamma ((1+ndf)/2) - digamma(ndf/2))
-    + log (sqrt ndf) 
+    + log (sqrt ndf)
     + logBeta (ndf/2) 0.5
 
 instance D.MaybeEntropy StudentT where
diff --git a/Statistics/Distribution/Transform.hs b/Statistics/Distribution/Transform.hs
--- a/Statistics/Distribution/Transform.hs
+++ b/Statistics/Distribution/Transform.hs
@@ -10,16 +10,19 @@
 -- Portability : portable
 --
 -- Transformations over distributions
-module Statistics.Distribution.Transform (
-    LinearTransform (..)
-  , linTransFixedPoint
-  , scaleAround
-  ) where
+module Statistics.Distribution.Transform
+    (
+      LinearTransform (..)
+    , linTransFixedPoint
+    , scaleAround
+    ) where
 
+import Control.Applicative ((<*>))
 import Data.Binary (Binary)
+import Data.Binary (put, get)
 import Data.Data (Data, Typeable)
+import Data.Functor ((<$>))
 import GHC.Generics (Generic)
-import Data.Functor          ((<$>))
 import qualified Statistics.Distribution as D
 
 -- | Linear transformation applied to distribution.
@@ -35,7 +38,9 @@
     -- ^ Distribution being transformed.
   } deriving (Eq, Show, Read, Typeable, Data, Generic)
 
-instance (Binary d) => Binary (LinearTransform d)
+instance (Binary d) => Binary (LinearTransform d) where
+    get = LinearTransform <$> get <*> get <*> get
+    put (LinearTransform x y z) = put x >> put y >> put z
 
 -- | Apply linear transformation to distribution.
 scaleAround :: Double           -- ^ Fixed point
@@ -73,11 +78,11 @@
   variance (LinearTransform _ sc dist) = sc * sc * D.variance dist
   stdDev   (LinearTransform _ sc dist) = sc * D.stdDev dist
 
-instance (D.MaybeEntropy d, D.DiscreteDistr d) 
+instance (D.MaybeEntropy d, D.DiscreteDistr d)
          => D.MaybeEntropy (LinearTransform d) where
   maybeEntropy (LinearTransform _ _ dist) = D.maybeEntropy dist
 
-instance (D.Entropy d, D.DiscreteDistr d) 
+instance (D.Entropy d, D.DiscreteDistr d)
          => D.Entropy (LinearTransform d) where
   entropy (LinearTransform _ _ dist) = D.entropy dist
 
diff --git a/Statistics/Distribution/Uniform.hs b/Statistics/Distribution/Uniform.hs
--- a/Statistics/Distribution/Uniform.hs
+++ b/Statistics/Distribution/Uniform.hs
@@ -24,6 +24,8 @@
 import GHC.Generics (Generic)
 import qualified Statistics.Distribution as D
 import qualified System.Random.MWC       as MWC
+import Data.Binary (put, get)
+import Control.Applicative ((<$>), (<*>))
 
 
 -- | Uniform distribution from A to B
@@ -32,7 +34,9 @@
     , uniformB :: {-# UNPACK #-} !Double -- ^ Upper boundary of distribution
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary UniformDistribution
+instance Binary UniformDistribution where
+    put (UniformDistribution x y) = put x >> put y
+    get = UniformDistribution <$> get <*> get
 
 -- | Create uniform distribution.
 uniformDistr :: Double -> Double -> UniformDistribution
diff --git a/Statistics/Function/Comparison.hs b/Statistics/Function/Comparison.hs
--- a/Statistics/Function/Comparison.hs
+++ b/Statistics/Function/Comparison.hs
@@ -16,9 +16,9 @@
       within
     ) where
 
-import Control.Monad.ST         (runST)
+import Control.Monad.ST (runST)
+import Data.Int (Int64)
 import Data.Primitive.ByteArray (newByteArray, readByteArray, writeByteArray)
-import Data.Int                 (Int64)
 
 -- | Compare two 'Double' values for approximate equality, using
 -- Dawson's method.
diff --git a/Statistics/Math.hs b/Statistics/Math.hs
--- a/Statistics/Math.hs
+++ b/Statistics/Math.hs
@@ -15,7 +15,7 @@
 -- 'Numeric.SpecFunctions.Extra' and 'Numeric.Polynomial.Chebyshev'.
 
 module Statistics.Math
-{-# DEPRECATED "Use package math-function" #-} 
+{-# DEPRECATED "Use package math-function" #-}
     ( module Numeric.Polynomial.Chebyshev
     , module Numeric.SpecFunctions
     , module Numeric.SpecFunctions.Extra
@@ -24,4 +24,3 @@
 import Numeric.Polynomial.Chebyshev
 import Numeric.SpecFunctions
 import Numeric.SpecFunctions.Extra
-
diff --git a/Statistics/Math/RootFinding.hs b/Statistics/Math/RootFinding.hs
--- a/Statistics/Math/RootFinding.hs
+++ b/Statistics/Math/RootFinding.hs
@@ -20,13 +20,15 @@
     -- $references
     ) where
 
-import Statistics.Function.Comparison
-
+import Control.Applicative (Alternative(..), Applicative(..))
+import Control.Monad (MonadPlus(..), ap)
 import Data.Binary (Binary)
-import Control.Applicative
-import Control.Monad       (MonadPlus(..), ap)
+import Data.Binary (put, get)
+import Data.Binary.Get (getWord8)
+import Data.Binary.Put (putWord8)
 import Data.Data (Data, Typeable)
 import GHC.Generics (Generic)
+import Statistics.Function.Comparison (within)
 
 
 -- | The result of searching for a root of a mathematical function.
@@ -40,7 +42,18 @@
             -- ^ A root was successfully found.
               deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance (Binary a) => Binary (Root a)
+instance (Binary a) => Binary (Root a) where
+    put NotBracketed = putWord8 0
+    put SearchFailed = putWord8 1
+    put (Root a) = putWord8 2 >> put a
+
+    get = do
+        i <- getWord8
+        case i of
+            0 -> return NotBracketed
+            1 -> return SearchFailed
+            2 -> fmap Root get
+            _ -> fail $ "Root.get: Invalid value: " ++ show i
 
 instance Functor Root where
     fmap _ NotBracketed = NotBracketed
diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs
--- a/Statistics/Quantile.hs
+++ b/Statistics/Quantile.hs
@@ -37,14 +37,14 @@
     -- $references
     ) where
 
-import Data.Vector.Generic             ((!))
+import Data.Vector.Generic ((!))
 import Numeric.MathFunctions.Constants (m_epsilon)
-import Statistics.Function             (partialSort)
+import Statistics.Function (partialSort)
 import qualified Data.Vector.Generic as G
 
 -- | O(/n/ log /n/). Estimate the /k/th /q/-quantile of a sample,
 -- using the weighted average method.
-weightedAvg :: G.Vector v Double => 
+weightedAvg :: G.Vector v Double =>
                Int        -- ^ /k/, the desired quantile.
             -> Int        -- ^ /q/, the number of quantiles.
             -> v Double   -- ^ /x/, the sample data.
diff --git a/Statistics/Resampling.hs b/Statistics/Resampling.hs
--- a/Statistics/Resampling.hs
+++ b/Statistics/Resampling.hs
@@ -15,7 +15,12 @@
     (
       Resample(..)
     , jackknife
+    , jackknifeMean
+    , jackknifeVariance
+    , jackknifeVarianceUnb
+    , jackknifeStdDev
     , resample
+    , estimate
     ) where
 
 import Control.Concurrent (forkIO, newChan, readChan, writeChan)
@@ -29,9 +34,12 @@
 import Data.Word (Word32)
 import GHC.Conc (numCapabilities)
 import GHC.Generics (Generic)
+import Numeric.Sum (Summation(..), kbn)
 import Statistics.Function (indices)
-import Statistics.Types (Estimator, Sample)
+import Statistics.Sample (mean, stdDev, variance, varianceUnbiased)
+import Statistics.Types (Estimator(..), Sample)
 import System.Random.MWC (Gen, initialize, uniform, uniformVector)
+import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
 
@@ -42,7 +50,9 @@
       fromResample :: U.Vector Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary Resample
+instance Binary Resample where
+    put = put . fromResample
+    get = fmap Resample get
 
 -- | /O(e*r*s)/ Resample a data set repeatedly, with replacement,
 -- computing each estimate over the resampled data.
@@ -80,18 +90,84 @@
             forM_ ers $ \(est,arr) ->
                 MU.write arr k . est $ re
             loop (k+1) ers
-      loop start (zip ests results)
+      loop start (zip ests' results)
   replicateM_ numCapabilities $ readChan done
   mapM_ sort results
   mapM (liftM Resample . unsafeFreeze) results
+ where
+  ests' = map estimate ests
 
--- | Compute a statistical estimate repeatedly over a sample, each
--- time omitting a successive element.
+-- | Run an 'Estimator' over a sample.
+estimate :: Estimator -> Sample -> Double
+estimate Mean             = mean
+estimate Variance         = variance
+estimate VarianceUnbiased = varianceUnbiased
+estimate StdDev           = stdDev
+estimate (Function est) = est
+
+-- | /O(n) or O(n^2)/ Compute a statistical estimate repeatedly over a
+-- sample, each time omitting a successive element.
 jackknife :: Estimator -> Sample -> U.Vector Double
-jackknife est sample = U.map f . indices $ sample
-    where f i = est (dropAt i sample)
-{- INLINE jackknife #-}
+jackknife Mean sample             = jackknifeMean sample
+jackknife Variance sample         = jackknifeVariance sample
+jackknife VarianceUnbiased sample = jackknifeVarianceUnb sample
+jackknife StdDev sample = jackknifeStdDev sample
+jackknife (Function est) sample
+  | G.length sample == 1 = singletonErr "jackknife"
+  | otherwise            = U.map f . indices $ sample
+  where f i = est (dropAt i sample)
 
+-- | /O(n)/ Compute the jackknife mean of a sample.
+jackknifeMean :: Sample -> U.Vector Double
+jackknifeMean samp
+  | len == 1  = singletonErr "jackknifeMean"
+  | otherwise = G.map (/l) $ G.zipWith (+) (pfxSumL samp) (pfxSumR samp)
+  where
+    l   = fromIntegral (len - 1)
+    len = G.length samp
+
+-- | /O(n)/ Compute the jackknife variance of a sample with a
+-- correction factor @c@, so we can get either the regular or
+-- \"unbiased\" variance.
+jackknifeVariance_ :: Double -> Sample -> U.Vector Double
+jackknifeVariance_ c samp
+  | len == 1  = singletonErr "jackknifeVariance"
+  | otherwise = G.zipWith4 go als ars bls brs
+  where
+    als = pfxSumL . G.map goa $ samp
+    ars = pfxSumR . G.map goa $ samp
+    goa x = v * v where v = x - m
+    bls = pfxSumL . G.map (subtract m) $ samp
+    brs = pfxSumR . G.map (subtract m) $ samp
+    m = mean samp
+    n = fromIntegral len
+    go al ar bl br = (al + ar - (b * b) / q) / q
+      where b = bl + br
+            q = n - (1 + c)
+    len = G.length samp
+
+-- | /O(n)/ Compute the unbiased jackknife variance of a sample.
+jackknifeVarianceUnb :: Sample -> U.Vector Double
+jackknifeVarianceUnb = jackknifeVariance_ 1
+
+-- | /O(n)/ Compute the jackknife variance of a sample.
+jackknifeVariance :: Sample -> U.Vector Double
+jackknifeVariance = jackknifeVariance_ 0
+
+-- | /O(n)/ Compute the jackknife standard deviation of a sample.
+jackknifeStdDev :: Sample -> U.Vector Double
+jackknifeStdDev = G.map sqrt . jackknifeVarianceUnb
+
+pfxSumL :: U.Vector Double -> U.Vector Double
+pfxSumL = G.map kbn . G.scanl add zero
+
+pfxSumR :: U.Vector Double -> U.Vector Double
+pfxSumR = G.tail . G.map kbn . G.scanr (flip add) zero
+
 -- | Drop the /k/th element of a vector.
 dropAt :: U.Unbox e => Int -> U.Vector e -> U.Vector e
 dropAt n v = U.slice 0 n v U.++ U.slice (n+1) (U.length v - n - 1) v
+
+singletonErr :: String -> a
+singletonErr func = error $
+                    "Statistics.Resampling." ++ func ++ ": singleton input"
diff --git a/Statistics/Resampling/Bootstrap.hs b/Statistics/Resampling/Bootstrap.hs
--- a/Statistics/Resampling/Bootstrap.hs
+++ b/Statistics/Resampling/Bootstrap.hs
@@ -21,20 +21,23 @@
     -- $references
     ) where
 
+import Control.Applicative ((<$>), (<*>))
 import Control.DeepSeq (NFData)
 import Control.Exception (assert)
-import Control.Monad.Par               (parMap,runPar)
+import Control.Monad.Par (parMap, runPar)
 import Data.Binary (Binary)
+import Data.Binary (put, get)
 import Data.Data (Data)
 import Data.Typeable (Typeable)
 import Data.Vector.Unboxed ((!))
-import GHC.Generics
+import GHC.Generics (Generic)
 import Statistics.Distribution (cumulative, quantile)
 import Statistics.Distribution.Normal
 import Statistics.Resampling (Resample(..), jackknife)
 import Statistics.Sample (mean)
 import Statistics.Types (Estimator, Sample)
 import qualified Data.Vector.Unboxed as U
+import qualified Statistics.Resampling as R
 
 -- | A point and interval estimate computed via an 'Estimator'.
 data Estimate = Estimate {
@@ -50,7 +53,9 @@
     -- ^ Confidence level of the confidence intervals.
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary Estimate
+instance Binary Estimate where
+    put (Estimate w x y z) = put w >> put x >> put y >> put z
+    get = Estimate <$> get <*> get <*> get <*> get
 instance NFData Estimate
 
 -- | Multiply the point, lower bound, and upper bound in an 'Estimate'
@@ -93,7 +98,7 @@
       | otherwise =
           estimate pt (resample ! lo) (resample ! hi) confidenceLevel
       where
-        pt    = est sample
+        pt    = R.estimate est sample
         lo    = max (cumn a1) 0
           where a1 = bias + b1 / (1 - accel * b1)
                 b1 = bias + z1
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -21,6 +21,7 @@
 
     -- * Statistics of location
     , mean
+    , welfordMean
     , meanWeighted
     , harmonicMean
     , geometricMean
@@ -54,11 +55,14 @@
     ) where
 
 import Statistics.Function (minMax)
+import Statistics.Sample.Internal (robustSumVar, sum)
 import Statistics.Types (Sample,WeightedSample)
+import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
 
 -- Operator ^ will be overriden
-import Prelude hiding ((^))
+import Prelude hiding ((^), sum)
 
 -- | /O(n)/ Range. The difference between the largest and smallest
 -- elements of a sample.
@@ -67,19 +71,31 @@
     where (lo , hi) = minMax s
 {-# INLINE range #-}
 
+-- | /O(n)/ Arithmetic mean.  This uses Kahan-Babuška-Neumaier
+-- summation, so is more accurate than 'welfordMean' unless the input
+-- values are very large.
+mean :: (G.Vector v Double) => v Double -> Double
+mean xs = sum xs / fromIntegral (G.length xs)
+{-# SPECIALIZE mean :: U.Vector Double -> Double #-}
+{-# SPECIALIZE mean :: V.Vector Double -> Double #-}
+
 -- | /O(n)/ Arithmetic mean.  This uses Welford's algorithm to provide
 -- numerical stability, using a single pass over the sample data.
-mean :: (G.Vector v Double) => v Double -> Double
-mean = fini . G.foldl' go (T 0 0)
+--
+-- Compared to 'mean', this loses a surprising amount of precision
+-- unless the inputs are very large.
+welfordMean :: (G.Vector v Double) => v Double -> Double
+welfordMean = fini . G.foldl' go (T 0 0)
   where
     fini (T a _) = a
     go (T m n) x = T m' n'
         where m' = m + (x - m) / fromIntegral n'
               n' = n + 1
-{-# INLINE mean #-}
+{-# SPECIALIZE welfordMean :: U.Vector Double -> Double #-}
+{-# SPECIALIZE welfordMean :: V.Vector Double -> Double #-}
 
 -- | /O(n)/ Arithmetic mean for weighted sample. It uses a single-pass
--- algorithm analogous to the one used by 'mean'.
+-- algorithm analogous to the one used by 'welfordMean'.
 meanWeighted :: (G.Vector v (Double,Double)) => v (Double,Double) -> Double
 meanWeighted = fini . G.foldl' go (V 0 0)
     where
@@ -117,7 +133,7 @@
     | a < 0  = error "Statistics.Sample.centralMoment: negative input"
     | a == 0 = 1
     | a == 1 = 0
-    | otherwise = G.sum (G.map go xs) / fromIntegral (G.length xs)
+    | otherwise = sum (G.map go xs) / fromIntegral (G.length xs)
   where
     go x = (x-m) ^ a
     m    = mean xs
@@ -202,11 +218,6 @@
 
 data V = V {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
-robustSumVar :: (G.Vector v Double) => Double -> v Double -> Double
-robustSumVar m = G.sum . G.map (square . subtract m)
-  where square x = x * x
-{-# INLINE robustSumVar #-}
-
 -- | Maximum likelihood estimate of a sample's variance.  Also known
 -- as the population variance, where the denominator is /n/.
 variance :: (G.Vector v Double) => v Double -> Double
@@ -242,7 +253,7 @@
 -- | Calculate mean and unbiased estimate of variance. This
 -- function should be used if both mean and variance are required
 -- since it will calculate mean only once.
-meanVarianceUnb ::  (G.Vector v Double) => v Double -> (Double,Double)
+meanVarianceUnb :: (G.Vector v Double) => v Double -> (Double,Double)
 meanVarianceUnb samp
   | n > 1     = (m, robustSumVar m samp / fromIntegral (n-1))
   | otherwise = (m, 0)
diff --git a/Statistics/Sample/Internal.hs b/Statistics/Sample/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Sample/Internal.hs
@@ -0,0 +1,30 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+-- |
+-- Module    : Statistics.Sample.Internal
+-- Copyright : (c) 2013 Bryan O'Sullivan
+-- License   : BSD3
+--
+-- Maintainer  : bos@serpentine.com
+-- Stability   : experimental
+-- Portability : portable
+--
+-- Internal functions for computing over samples.
+module Statistics.Sample.Internal
+    (
+      robustSumVar
+    , sum
+    ) where
+
+import Numeric.Sum (kbn, sumVector)
+import Prelude hiding (sum)
+import qualified Data.Vector.Generic as G
+
+robustSumVar :: (G.Vector v Double) => Double -> v Double -> Double
+robustSumVar m = sum . G.map (square . subtract m)
+  where square x = x * x
+{-# INLINE robustSumVar #-}
+
+sum :: (G.Vector v Double) => v Double -> Double
+sum = sumVector kbn
+{-# INLINE sum #-}
diff --git a/Statistics/Sample/KernelDensity.hs b/Statistics/Sample/KernelDensity.hs
--- a/Statistics/Sample/KernelDensity.hs
+++ b/Statistics/Sample/KernelDensity.hs
@@ -25,17 +25,17 @@
     -- $references
     ) where
 
-import Prelude hiding (const,min,max)
 import Numeric.MathFunctions.Constants (m_sqrt_2_pi)
-import Statistics.Function             (minMax, nextHighestPowerOfTwo)
-import Statistics.Math.RootFinding     (fromRoot, ridders)
-import Statistics.Sample.Histogram     (histogram_)
-import Statistics.Transform            (dct, idct)
+import Prelude hiding (const, min, max, sum)
+import Statistics.Function (minMax, nextHighestPowerOfTwo)
+import Statistics.Math.RootFinding (fromRoot, ridders)
+import Statistics.Sample.Histogram (histogram_)
+import Statistics.Sample.Internal (sum)
+import Statistics.Transform (dct, idct)
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
 
 
-
 -- | Gaussian kernel density estimator for one-dimensional data, using
 -- the method of Botev et al.
 --
@@ -87,13 +87,13 @@
     !n  = fromIntegral ni
     !ni = nextHighestPowerOfTwo n0
     !r  = max - min
-    a   = dct . G.map (/ G.sum h) $ h
+    a   = dct . G.map (/ 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
-        f q t = 2 * pi ** (q*2) * G.sum (G.zipWith g iv a2v)
+        f q t = 2 * pi ** (q*2) * sum (G.zipWith g iv a2v)
           where g i a2 = i ** q * a2 * exp ((-i) * sqr pi * t)
                 a2v = G.map (sqr . (*0.5)) $ G.tail a
                 iv = G.map sqr $ G.enumFromTo 1 (n-1)
diff --git a/Statistics/Sample/KernelDensity/Simple.hs b/Statistics/Sample/KernelDensity/Simple.hs
--- a/Statistics/Sample/KernelDensity/Simple.hs
+++ b/Statistics/Sample/KernelDensity/Simple.hs
@@ -51,17 +51,21 @@
 import Data.Vector.Binary ()
 import GHC.Generics (Generic)
 import Numeric.MathFunctions.Constants (m_1_sqrt_2, m_2_sqrt_pi)
+import Prelude hiding (sum)
 import Statistics.Function (minMax)
-import Statistics.Sample   (stdDev)
-import qualified Data.Vector.Unboxed as U
+import Statistics.Sample (stdDev)
+import Statistics.Sample.Internal (sum)
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
 
 -- | Points from the range of a 'Sample'.
 newtype Points = Points {
       fromPoints :: U.Vector Double
     } deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary Points
+instance Binary Points where
+    get = fmap Points get
+    put = put . fromPoints
 
 -- | Bandwidth estimator for an Epanechnikov kernel.
 epanechnikovBW :: Double -> Bandwidth
@@ -148,7 +152,7 @@
     | n < 2     = errorShort "estimatePDF"
     | otherwise = U.map k . fromPoints
   where
-    k p = G.sum . G.map (kernel f h p) $ sample
+    k p = sum . G.map (kernel f h p) $ sample
     f   = 1 / (h * fromIntegral n)
     n   = G.length sample
 {-# INLINE estimatePDF #-}
diff --git a/Statistics/Sample/Powers.hs b/Statistics/Sample/Powers.hs
--- a/Statistics/Sample/Powers.hs
+++ b/Statistics/Sample/Powers.hs
@@ -50,22 +50,25 @@
 import Data.Binary (Binary(..))
 import Data.Data (Data, Typeable)
 import Data.Vector.Binary ()
-import Data.Vector.Generic   (unsafeFreeze)
-import Data.Vector.Unboxed   ((!))
+import Data.Vector.Generic (unsafeFreeze)
+import Data.Vector.Unboxed ((!))
 import GHC.Generics (Generic)
-import Prelude hiding (sum)
-import Statistics.Function   (indexed)
-import Statistics.Internal   (inlinePerformIO)
 import Numeric.SpecFunctions (choose)
-import System.IO.Unsafe      (unsafePerformIO)
-import qualified Data.Vector.Unboxed as U
+import Prelude hiding (sum)
+import Statistics.Function (indexed)
+import Statistics.Internal (inlinePerformIO)
+import System.IO.Unsafe (unsafePerformIO)
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as MU
+import qualified Statistics.Sample.Internal as S
 
 newtype Powers = Powers (U.Vector Double)
     deriving (Eq, Read, Show, Typeable, Data, Generic)
 
-instance Binary Powers
+instance Binary Powers where
+    put (Powers v) = put v
+    get = fmap Powers get
 
 -- | O(/n/) Collect the /n/ simple powers of a sample.
 --
@@ -112,7 +115,7 @@
                   error ("Statistics.Sample.Powers.centralMoment: "
                          ++ "invalid argument")
     | k == 0    = 1
-    | otherwise = (/n) . U.sum . U.map go . indexed . U.take (k+1) $ pa
+    | otherwise = (/n) . S.sum . U.map go . indexed . U.take (k+1) $ pa
   where
     go (i , e) = (k `choose` i) * ((-m) ^ (k-i)) * e
     n = U.head pa
diff --git a/Statistics/Test/ChiSquared.hs b/Statistics/Test/ChiSquared.hs
--- a/Statistics/Test/ChiSquared.hs
+++ b/Statistics/Test/ChiSquared.hs
@@ -7,11 +7,12 @@
   , TestResult(..)
   ) where
 
-import qualified Data.Vector.Generic as G
-
+import Prelude hiding (sum)
 import Statistics.Distribution
 import Statistics.Distribution.ChiSquared
+import Statistics.Sample.Internal (sum)
 import Statistics.Test.Types
+import qualified Data.Vector.Generic as G
 
 
 -- | Generic form of Pearson chi squared tests for binned data. Data
@@ -33,7 +34,7 @@
   | otherwise      = error $ "Statistics.Test.ChiSquare.chi2test: bad p-value: " ++ show p
   where
     n     = G.length vec - ndf - 1
-    chi2  = G.sum $ G.map (\(o,e) -> sqr (fromIntegral o - e) / e) vec
+    chi2  = sum $ G.map (\(o,e) -> sqr (fromIntegral o - e) / e) vec
     d     = chiSquared n
     sqr x = x * x
 {-# INLINE chi2test #-}
diff --git a/Statistics/Test/KolmogorovSmirnov.hs b/Statistics/Test/KolmogorovSmirnov.hs
--- a/Statistics/Test/KolmogorovSmirnov.hs
+++ b/Statistics/Test/KolmogorovSmirnov.hs
@@ -29,21 +29,19 @@
     -- $references
   ) where
 
-import Control.Monad
-import Control.Monad.ST  (ST)
-
-import qualified Data.Vector.Unboxed         as U
+import Control.Monad (when)
+import Control.Monad.ST (ST)
+import Prelude hiding (sum)
+import Statistics.Distribution (Distribution(..))
+import Statistics.Function (sort)
+import Statistics.Sample.Internal (sum)
+import Statistics.Test.Types (TestResult(..), TestType(..), significant)
+import Statistics.Types (Sample)
+import Text.Printf (printf)
+import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Unboxed.Mutable as M
 
-import Statistics.Distribution        (Distribution(..))
-import Statistics.Types               (Sample)
-import Statistics.Function            (sort)
-import Statistics.Test.Types
 
-import Text.Printf
-
-
-
 ----------------------------------------------------------------
 -- Test
 ----------------------------------------------------------------
@@ -267,7 +265,7 @@
 matrixMultiply (Matrix n xs e1) (Matrix _ ys e2) =
   Matrix n (U.generate (n*n) go) (e1 + e2)
   where
-    go i = U.sum $ U.zipWith (*) row col
+    go i = sum $ U.zipWith (*) row col
       where
         nCol = i `rem` n
         row  = U.slice (i - nCol) n xs
diff --git a/Statistics/Test/MannWhitneyU.hs b/Statistics/Test/MannWhitneyU.hs
--- a/Statistics/Test/MannWhitneyU.hs
+++ b/Statistics/Test/MannWhitneyU.hs
@@ -27,20 +27,18 @@
   ) where
 
 import Control.Applicative ((<$>))
-import Data.List           (findIndex)
-import Data.Ord            (comparing)
-import qualified Data.Vector.Unboxed as U
-
-import Numeric.SpecFunctions          (choose)
-
-import Statistics.Distribution        (quantile)
+import Data.List (findIndex)
+import Data.Ord (comparing)
+import Numeric.SpecFunctions (choose)
+import Prelude hiding (sum)
+import Statistics.Distribution (quantile)
 import Statistics.Distribution.Normal (standard)
-import Statistics.Types               (Sample)
-import Statistics.Function            (sortBy)
-import Statistics.Test.Types
-import Statistics.Test.Internal
-
-
+import Statistics.Function (sortBy)
+import Statistics.Sample.Internal (sum)
+import Statistics.Test.Internal (rank, splitByTags)
+import Statistics.Test.Types (TestResult(..), TestType(..), significant)
+import Statistics.Types (Sample)
+import qualified Data.Vector.Unboxed as U
 
 -- | The Wilcoxon Rank Sums Test.
 --
@@ -54,8 +52,7 @@
 -- and the related functions for testing significance, but this function is exposed
 -- for completeness.
 wilcoxonRankSums :: Sample -> Sample -> (Double, Double)
-wilcoxonRankSums xs1 xs2 = ( U.sum ranks1 , U.sum ranks2
-                           )
+wilcoxonRankSums xs1 xs2 = (sum ranks1, sum ranks2)
   where
     -- Ranks for each sample
     (ranks1,ranks2) = splitByTags $ U.zip tags (rank (==) joinSample)
@@ -138,12 +135,12 @@
     n = bigN - m
 -}
 
--- Memoised version of the original a function, above. 
+-- Memoised version of the original a function, above.
 --
 -- Doubles are stored to avoid integer overflow. 32-bit Ints begin to
 -- overflow for bigN as small as 33 (64-bit one at 66) while Double to
 -- go to infinity till bigN=1029
--- 
+--
 --
 -- outer list is indexed by big N - 2
 -- inner list by (m-1) (we know m < bigN)
diff --git a/Statistics/Test/WilcoxonT.hs b/Statistics/Test/WilcoxonT.hs
--- a/Statistics/Test/WilcoxonT.hs
+++ b/Statistics/Test/WilcoxonT.hs
@@ -25,34 +25,34 @@
   , TestResult(..)
   ) where
 
-import Control.Applicative ((<$>))
-import Data.Function       (on)
-import Data.List           (findIndex)
-import Data.Ord            (comparing)
-import qualified Data.Vector.Unboxed as U
 
-import Statistics.Types               (Sample)
-import Statistics.Function            (sortBy)
-import Statistics.Test.Types
-import Statistics.Test.Internal
 
-
--- | The Wilcoxon matched-pairs signed-rank test.
 --
+--
+--
+-- Note that: wilcoxonMatchedPairSignedRank == (\(x, y) -> (y, x)) . flip wilcoxonMatchedPairSignedRank
+-- The samples are zipped together: if one is longer than the other, both are truncated
 -- The value returned is the pair (T+, T-).  T+ is the sum of positive ranks (the
--- ranks of the differences where the first parameter is higher) whereas T- is
--- the sum of negative ranks (the ranks of the differences where the second parameter is higher).
 -- These values mean little by themselves, and should be combined with the 'wilcoxonSignificant'
 -- function in this module to get a meaningful result.
---
--- The samples are zipped together: if one is longer than the other, both are truncated
+-- ranks of the differences where the first parameter is higher) whereas T- is
+-- the sum of negative ranks (the ranks of the differences where the second parameter is higher).
 -- to the the length of the shorter sample.
---
--- Note that: wilcoxonMatchedPairSignedRank == (\(x, y) -> (y, x)) . flip wilcoxonMatchedPairSignedRank
+
+import Control.Applicative ((<$>))
+import Data.Function (on)
+import Data.List (findIndex)
+import Data.Ord (comparing)
+import Prelude hiding (sum)
+import Statistics.Function (sortBy)
+import Statistics.Sample.Internal (sum)
+import Statistics.Test.Internal (rank, splitByTags)
+import Statistics.Test.Types (TestResult(..), TestType(..), significant)
+import Statistics.Types (Sample)
+import qualified Data.Vector.Unboxed as U
+
 wilcoxonMatchedPairSignedRank :: Sample -> Sample -> (Double, Double)
-wilcoxonMatchedPairSignedRank a b = (          U.sum ranks1
-                                    , negate $ U.sum ranks2
-                                    )
+wilcoxonMatchedPairSignedRank a b = (sum ranks1, negate (sum ranks2))
   where
     (ranks1, ranks2) = splitByTags
                      $ U.zip tags (rank ((==) `on` abs) diffs)
diff --git a/Statistics/Transform.hs b/Statistics/Transform.hs
--- a/Statistics/Transform.hs
+++ b/Statistics/Transform.hs
@@ -29,14 +29,14 @@
     , ifft
     ) where
 
-import Control.Monad         (when)
-import Control.Monad.ST      (ST)
-import Data.Bits             (shiftL, shiftR)
-import Data.Complex          (Complex(..), conjugate, realPart)
+import Control.Monad (when)
+import Control.Monad.ST (ST)
+import Data.Bits (shiftL, shiftR)
+import Data.Complex (Complex(..), conjugate, realPart)
 import Numeric.SpecFunctions (log2)
-import qualified Data.Vector.Generic         as G
+import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Generic.Mutable as M
-import qualified Data.Vector.Unboxed         as U
+import qualified Data.Vector.Unboxed as U
 
 
 type CD = Complex Double
@@ -90,7 +90,7 @@
     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 
+    weights
       = G.cons n
       $ G.generate (len - 1) $ \x -> 2 * n * exp ((0:+1) * fi (x+1) * pi/(2*n))
       where n = fi len
diff --git a/Statistics/Types.hs b/Statistics/Types.hs
--- a/Statistics/Types.hs
+++ b/Statistics/Types.hs
@@ -11,7 +11,7 @@
 
 module Statistics.Types
     (
-      Estimator
+      Estimator(..)
     , Sample
     , WeightedSample
     , Weights
@@ -25,9 +25,16 @@
 -- | Sample with weights. First element of sample is data, second is weight
 type WeightedSample = U.Vector (Double,Double)
 
--- | A function that estimates a property of a sample, such as its
--- 'mean'.
-type Estimator = Sample -> Double
+-- | An estimator of a property of a sample, such as its 'mean'.
+--
+-- The use of an algebraic data type here allows functions such as
+-- 'jackknife' and 'bootstrapBCA' to use more efficient algorithms
+-- when possible.
+data Estimator = Mean
+               | Variance
+               | VarianceUnbiased
+               | StdDev
+               | Function (Sample -> Double)
 
 -- | Weights for affecting the importance of elements of a sample.
 type Weights = U.Vector Double
diff --git a/benchmark/bench.hs b/benchmark/bench.hs
--- a/benchmark/bench.hs
+++ b/benchmark/bench.hs
@@ -1,12 +1,10 @@
 import Control.Monad.ST (runST)
-import Data.Complex
-import qualified Data.Vector.Unboxed as U
-
-import System.Random.MWC
 import Criterion.Main
-
+import Data.Complex
 import Statistics.Sample
 import Statistics.Transform
+import System.Random.MWC
+import qualified Data.Vector.Unboxed as U
 
 
 -- Test sample
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,5 @@
 name:           statistics
-version:        0.10.5.2
+version:        0.11.0.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -87,15 +87,16 @@
     Statistics.Distribution.Poisson.Internal
     Statistics.Function.Comparison
     Statistics.Internal
+    Statistics.Sample.Internal
     Statistics.Test.Internal
   build-depends:
     base < 5,
-    binary >= 0.6.3.0,
+    binary >= 0.5.1.0,
     deepseq >= 1.1.0.2,
     erf,
     monad-par         >= 0.3.4,
     mwc-random        >= 0.13.0.0,
-    math-functions    >= 0.1.2,
+    math-functions    >= 0.1.5.2,
     primitive         >= 0.3,
     vector            >= 0.7.1,
     vector-algorithms >= 0.4,
@@ -122,8 +123,8 @@
     Tests.Distribution
     Tests.Helpers
     Tests.Function
-    Tests.NonparametricTest
-    Tests.NonparametricTest.Table
+    Tests.NonParametric
+    Tests.NonParametric.Table
     Tests.Transform
     Tests.KDE
 
@@ -132,6 +133,7 @@
 
   build-depends:
     base,
+    binary,
     ieee754 >= 0.7.3,
     HUnit,
     QuickCheck >= 2,
diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs
--- a/tests/Tests/Distribution.hs
+++ b/tests/Tests/Distribution.hs
@@ -1,51 +1,41 @@
 {-# 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
+{-# LANGUAGE FlexibleInstances, OverlappingInstances, ScopedTypeVariables,
+    ViewPatterns #-}
+module Tests.Distribution (tests) where
 
-import Data.List     (find)
+import Control.Applicative ((<$>), (<*>))
+import Data.Binary (Binary, decode, encode)
+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.Beta (BetaDistribution, betaDistr)
+import Statistics.Distribution.Binomial (BinomialDistribution, binomial)
 import Statistics.Distribution.CauchyLorentz
-import Statistics.Distribution.Exponential
-import Statistics.Distribution.FDistribution
-import Statistics.Distribution.Gamma
+import Statistics.Distribution.ChiSquared (ChiSquared, chiSquared)
+import Statistics.Distribution.Exponential (ExponentialDistribution, exponential)
+import Statistics.Distribution.FDistribution (FDistribution, fDistribution)
+import Statistics.Distribution.Gamma (GammaDistribution, gammaDistr)
 import Statistics.Distribution.Geometric
 import Statistics.Distribution.Hypergeometric
-import Statistics.Distribution.Normal
-import Statistics.Distribution.Poisson
+import Statistics.Distribution.Normal (NormalDistribution, normalDistr)
+import Statistics.Distribution.Poisson (PoissonDistribution, poisson)
 import Statistics.Distribution.StudentT
-import Statistics.Distribution.Transform
-import Statistics.Distribution.Uniform
-
-import Prelude hiding (catch)
-
-import Tests.Helpers
+import Statistics.Distribution.Transform (LinearTransform, linTransDistr)
+import Statistics.Distribution.Uniform (UniformDistribution, uniformDistr)
+import Test.Framework (Test, testGroup)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.QuickCheck as QC
+import Test.QuickCheck.Monadic as QC
+import Tests.Helpers (T(..), (=~), eq, testAssertion, typeName)
+import Tests.Helpers (monotonicallyIncreasesIEEE)
+import Text.Printf (printf)
+import qualified Control.Exception as E
+import qualified Numeric.IEEE as IEEE
 
 
 -- | Tests for all distributions
-distributionTests :: Test
-distributionTests = testGroup "Tests for all distributions"
+tests :: Test
+tests = testGroup "Tests for all distributions"
   [ contDistrTests (T :: T BetaDistribution        )
   , contDistrTests (T :: T CauchyDistribution      )
   , contDistrTests (T :: T ChiSquared              )
@@ -71,7 +61,7 @@
 ----------------------------------------------------------------
 
 -- Tests for continous distribution
-contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
+contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d, Binary d, Eq d) => T d -> Test
 contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
   cdfTests t ++
   [ testProperty "PDF sanity"              $ pdfSanityCheck     t
@@ -81,7 +71,7 @@
   ]
 
 -- Tests for discrete distribution
-discreteDistrTests :: (Param d, DiscreteDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> Test
+discreteDistrTests :: (Param d, DiscreteDistr d, QC.Arbitrary d, Typeable d, Show d, Binary d, Eq d) => T d -> Test
 discreteDistrTests t = testGroup ("Tests for: " ++ typeName t) $
   cdfTests t ++
   [ testProperty "Prob. sanity"         $ probSanityCheck       t
@@ -91,7 +81,7 @@
   ]
 
 -- Tests for distributions which have CDF
-cdfTests :: (Param d, Distribution d, QC.Arbitrary d, Show d) => T d -> [Test]
+cdfTests :: (Param d, Distribution d, QC.Arbitrary d, Show d, Binary d, Eq d) => T d -> [Test]
 cdfTests t =
   [ testProperty "C.D.F. sanity"        $ cdfSanityCheck         t
   , testProperty "CDF limit at +inf"    $ cdfLimitAtPosInfinity  t
@@ -100,12 +90,15 @@
   , testProperty "CDF at -inf = 1"      $ cdfAtNegInfinity       t
   , testProperty "CDF is nondecreasing" $ cdfIsNondecreasing     t
   , testProperty "1-CDF is correct"     $ cdfComplementIsCorrect t
+  , testProperty "Binary OK"            $ p_binary t
   ]
+
+
 ----------------------------------------------------------------
 
 -- CDF is in [0,1] range
 cdfSanityCheck :: (Distribution d) => T d -> d -> Double -> Bool
-cdfSanityCheck _ d x = c >= 0 && c <= 1 
+cdfSanityCheck _ d x = c >= 0 && c <= 1
   where c = cumulative d x
 
 -- CDF never decreases
@@ -148,7 +141,7 @@
 cdfDiscreteIsCorrect :: (DiscreteDistr d) => T d -> d -> Property
 cdfDiscreteIsCorrect _ d
   = printTestCase (unlines badN)
-  $ null badN  
+  $ null badN
   where
     -- We are checking that:
     --
@@ -202,15 +195,15 @@
 -- 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
+  p < 0 || p > 1 ==> QC.monadicIO $ do r <- QC.run $ E.catch
                                               (do { return $! quantile d p; return False })
-                                              (\(e :: SomeException) -> return True)
+                                              (\(_ :: 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 
+probSanityCheck _ d x = p >= 0 && p <= 1
   where p = probability d x
 
 -- Check that discrete CDF is correct
@@ -245,7 +238,11 @@
     logP = logProbability d x
 
 
-    
+p_binary :: (Eq a, Show a, Binary a) => T a -> a -> Bool
+p_binary _ a = a == (decode . encode) a
+
+
+
 ----------------------------------------------------------------
 -- Arbitrary instances for ditributions
 ----------------------------------------------------------------
@@ -289,7 +286,7 @@
            <*> ((abs <$> arbitrary))
            <*> ((abs <$> arbitrary) `suchThat` (>0))
 instance QC.Arbitrary FDistribution where
-  arbitrary =  fDistribution 
+  arbitrary =  fDistribution
            <$> ((abs <$> arbitrary) `suchThat` (>0))
            <*> ((abs <$> arbitrary) `suchThat` (>0))
 
diff --git a/tests/Tests/Function.hs b/tests/Tests/Function.hs
--- a/tests/Tests/Function.hs
+++ b/tests/Tests/Function.hs
@@ -1,15 +1,11 @@
 module Tests.Function ( tests ) where
 
-import qualified Data.Vector.Unboxed as U
-import           Data.Vector.Unboxed   ((!))
-
-import Test.QuickCheck
+import Statistics.Function
 import Test.Framework
 import Test.Framework.Providers.QuickCheck2
-
+import Test.QuickCheck
 import Tests.Helpers
-import Statistics.Function
-
+import qualified Data.Vector.Unboxed as U
 
 
 tests :: Test
@@ -29,5 +25,5 @@
 p_nextHighestPowerOfTwo
   = all (\(good, is) -> all ((==good) . nextHighestPowerOfTwo) is) lists
   where
-    pows  = [1 .. 17]
+    pows  = [1 .. 17 :: Int]
     lists = [ (2^m, [2^n+1 .. 2^m]) | (n,m) <- pows `zip` tail pows ]
diff --git a/tests/Tests/KDE.hs b/tests/Tests/KDE.hs
--- a/tests/Tests/KDE.hs
+++ b/tests/Tests/KDE.hs
@@ -1,18 +1,16 @@
 -- | Tests for Kernel density estimates.
-module Tests.KDE ( 
-  tests 
+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 Numeric.Sum (kbn, sumVector)
 import Statistics.Sample.KernelDensity
-
+import Test.Framework (Test, testGroup)
+import Test.Framework.Providers.QuickCheck2 (testProperty)
+import Test.QuickCheck (Property, (==>), printTestCase)
+import Text.Printf (printf)
+import qualified Data.Vector.Unboxed as U
 
 
 tests :: Test
@@ -21,7 +19,7 @@
   ]
 
 t_densityIsPDF :: [Double] -> Property
-t_densityIsPDF vec 
+t_densityIsPDF vec
   = not (null vec) ==> test
   where
     (xs,ys)  = kde 4096 (U.fromList vec)
@@ -31,11 +29,11 @@
     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)
+integratePDF step vec
+  = step * sumVector kbn (U.zipWith (*) vec weights)
   where
     n       = U.length vec
     weights = U.generate n go
diff --git a/tests/Tests/NonParametric.hs b/tests/Tests/NonParametric.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/NonParametric.hs
@@ -0,0 +1,226 @@
+-- Tests for Statistics.Test.NonParametric
+module Tests.NonParametric (tests) where
+
+import Statistics.Distribution.Normal (standard)
+import Statistics.Test.KolmogorovSmirnov
+import Statistics.Test.MannWhitneyU
+import Statistics.Test.WilcoxonT
+import Test.Framework (Test, testGroup)
+import Test.Framework.Providers.HUnit
+import Test.HUnit (assertEqual)
+import Tests.Helpers (eq, testAssertion, testEquality)
+import Tests.NonParametric.Table (tableKSD, tableKS2D)
+import qualified Data.Vector.Unboxed as U
+
+
+tests :: Test
+tests = 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/NonParametric/Table.hs b/tests/Tests/NonParametric/Table.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/NonParametric/Table.hs
@@ -0,0 +1,39 @@
+module Tests.NonParametric.Table (
+      tableKSD
+    , tableKS2D
+    ) 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/NonparametricTest.hs b/tests/Tests/NonparametricTest.hs
deleted file mode 100644
--- a/tests/Tests/NonparametricTest.hs
+++ /dev/null
@@ -1,233 +0,0 @@
--- 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
deleted file mode 100644
--- a/tests/Tests/NonparametricTest/Table.hs
+++ /dev/null
@@ -1,36 +0,0 @@
-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
--- a/tests/Tests/Transform.hs
+++ b/tests/Tests/Transform.hs
@@ -6,26 +6,21 @@
       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 Data.Bits ((.&.), shiftL)
+import Data.Complex (Complex((:+)))
+import Data.Functor ((<$>))
+import Numeric.Sum (kbn, sumVector)
+import Statistics.Function (within)
+import Statistics.Transform (CD, dct, fft, idct, ifft)
+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 Test.QuickCheck (Positive(..), Property, Arbitrary(..), Gen, choose, vectorOf, printTestCase)
+import Tests.Helpers (testAssertion)
+import Text.Printf (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
@@ -138,10 +133,10 @@
   vectorNorm :: a -> Double
 
 instance HasNorm (U.Vector Double) where
-  vectorNorm = sqrt . U.sum . U.map (\x -> x*x)
+  vectorNorm = sqrt . sumVector kbn . U.map (\x -> x*x)
 
 instance HasNorm (U.Vector CD) where
-  vectorNorm = sqrt . U.sum . U.map (\(x :+ y) -> x*x + y*y)
+  vectorNorm = sqrt . sumVector kbn . U.map (\(x :+ y) -> x*x + y*y)
 
 -- Approximate equality for vectors
 vecEqual :: Double -> U.Vector Double -> U.Vector Double -> Bool
diff --git a/tests/tests.hs b/tests/tests.hs
--- a/tests/tests.hs
+++ b/tests/tests.hs
@@ -1,15 +1,14 @@
-import Test.Framework       (defaultMain)
-
-import Tests.Distribution
-import Tests.NonparametricTest
-import qualified Tests.Transform
-import qualified Tests.Function
-import qualified Tests.KDE
+import Test.Framework (defaultMain)
+import qualified Tests.Distribution as Distribution
+import qualified Tests.Function as Function
+import qualified Tests.KDE as KDE
+import qualified Tests.NonParametric as NonParametric
+import qualified Tests.Transform as Transform
 
 main :: IO ()
-main = defaultMain [ distributionTests 
-                   , nonparametricTests
-                   , Tests.Transform.tests
-                   , Tests.Function.tests
-                   , Tests.KDE.tests
+main = defaultMain [ Distribution.tests
+                   , Function.tests
+                   , KDE.tests
+                   , NonParametric.tests
+                   , Transform.tests
                    ]