diff --git a/README.markdown b/README.markdown
--- a/README.markdown
+++ b/README.markdown
@@ -18,18 +18,11 @@
 # Get involved!
 
 Please report bugs via the
-[github issue tracker](https://github.com/bos/statistics/issues).
-
-Master [git mirror](https://github.com/bos/statistics):
-
-* `git clone git://github.com/bos/statistics.git`
-
-There's also a [Mercurial mirror](https://bitbucket.org/bos/statistics):
-
-* `hg clone https://bitbucket.org/bos/statistics`
+[github issue tracker](https://github.com/haskell/statistics/issues).
 
-(You can create and contribute changes using either Mercurial or git.)
+Master [git mirror](https://github.com/haskell/statistics):
 
+* `git clone git://github.com/haskell/statistics.git`
 
 # Authors
 
diff --git a/Setup.lhs b/Setup.lhs
deleted file mode 100644
--- a/Setup.lhs
+++ /dev/null
@@ -1,3 +0,0 @@
-#!/usr/bin/env runhaskell
-> import Distribution.Simple
-> main = defaultMain
diff --git a/Statistics/Correlation.hs b/Statistics/Correlation.hs
--- a/Statistics/Correlation.hs
+++ b/Statistics/Correlation.hs
@@ -6,9 +6,11 @@
 module Statistics.Correlation
     ( -- * Pearson correlation
       pearson
+    , pearson2
     , pearsonMatByRow
       -- * Spearman correlation
     , spearman
+    , spearman2
     , spearmanMatByRow
     ) where
 
@@ -25,12 +27,19 @@
 
 -- | Pearson correlation for sample of pairs. Exactly same as
 -- 'Statistics.Sample.correlation'
-pearson :: (G.Vector v (Double, Double), G.Vector v Double)
+pearson :: (G.Vector v (Double, Double))
         => v (Double, Double) -> Double
 pearson = correlation
 {-# INLINE pearson #-}
 
--- | Compute pairwise pearson correlation between rows of a matrix
+-- | Pearson correlation for sample of pairs. Exactly same as
+-- 'Statistics.Sample.correlation'
+pearson2 :: (G.Vector v Double)
+         => v Double -> v Double -> Double
+pearson2 = correlation2
+{-# INLINE pearson2 #-}
+
+-- | Compute pairwise Pearson correlation between rows of a matrix
 pearsonMatByRow :: Matrix -> Matrix
 pearsonMatByRow m
   = generateSym (rows m)
@@ -43,15 +52,13 @@
 -- Spearman
 ----------------------------------------------------------------
 
--- | compute spearman correlation between two samples
+-- | Compute Spearman correlation between two samples
 spearman :: ( Ord a
             , Ord b
             , G.Vector v a
             , G.Vector v b
             , G.Vector v (a, b)
             , G.Vector v Int
-            , G.Vector v Double
-            , G.Vector v (Double, Double)
             , G.Vector v (Int, a)
             , G.Vector v (Int, b)
             )
@@ -64,7 +71,28 @@
     (x, y) = G.unzip xy
 {-# INLINE spearman #-}
 
--- | compute pairwise spearman correlation between rows of a matrix
+-- | Compute Spearman correlation between two samples. Samples must
+--   have same length.
+spearman2 :: ( Ord a
+            , Ord b
+            , G.Vector v a
+            , G.Vector v b
+            , G.Vector v Int
+            , G.Vector v (Int, a)
+            , G.Vector v (Int, b)
+            )
+         => v a
+         -> v b
+         -> Double
+spearman2 xs ys
+  | nx /= ny  = error "Statistics.Correlation.spearman2: samples must have same length"
+  | otherwise = pearson $ G.zip (rankUnsorted xs) (rankUnsorted ys)
+  where
+    nx = G.length xs
+    ny = G.length ys
+{-# INLINE spearman2 #-}
+
+-- | compute pairwise Spearman correlation between rows of a matrix
 spearmanMatByRow :: Matrix -> Matrix
 spearmanMatByRow
   = pearsonMatByRow . fromRows . fmap rankUnsorted . toRows
diff --git a/Statistics/Correlation/Kendall.hs b/Statistics/Correlation/Kendall.hs
--- a/Statistics/Correlation/Kendall.hs
+++ b/Statistics/Correlation/Kendall.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE BangPatterns, CPP, FlexibleContexts #-}
+{-# LANGUAGE BangPatterns, FlexibleContexts #-}
 -- |
 -- Module      : Statistics.Correlation.Kendall
 --
@@ -130,11 +130,6 @@
         _  -> do GM.unsafeWrite src iIns eLow
                  wroteLow low (iLow+1) high iHigh eHigh (iIns+1)
 {-# INLINE merge #-}
-
-#if !MIN_VERSION_base(4,6,0)
-modifySTRef' :: STRef s a -> (a -> a) -> ST s ()
-modifySTRef' = modifySTRef
-#endif
 
 -- $references
 --
diff --git a/Statistics/Distribution.hs b/Statistics/Distribution.hs
--- a/Statistics/Distribution.hs
+++ b/Statistics/Distribution.hs
@@ -29,18 +29,15 @@
     , ContGen(..)
     , DiscreteGen(..)
     , genContinuous
-    , genContinous
       -- * Helper functions
     , findRoot
     , sumProbabilities
     ) where
 
-import Control.Applicative ((<$>), Applicative(..))
-import Control.Monad.Primitive (PrimMonad,PrimState)
 import Prelude hiding (sum)
-import Statistics.Function (square)
+import Statistics.Function        (square)
 import Statistics.Sample.Internal (sum)
-import System.Random.MWC (Gen, uniform)
+import System.Random.Stateful     (StatefulGen, uniformDouble01M)
 import qualified Data.Vector.Unboxed as U
 import qualified Data.Vector.Generic as G
 
@@ -50,13 +47,13 @@
 class Distribution d where
     -- | Cumulative distribution function.  The probability that a
     -- random variable /X/ is less or equal than /x/,
-    -- i.e. P(/X/&#8804;/x/). Cumulative should be defined for
+    -- i.e. P(/X/≤/x/). Cumulative should be defined for
     -- infinities as well:
     --
     -- > cumulative d +∞ = 1
     -- > cumulative d -∞ = 0
     cumulative :: d -> Double -> Double
-
+    cumulative d x = 1 - complCumulative d x
     -- | One's complement of cumulative distribution:
     --
     -- > complCumulative d x = 1 - cumulative d x
@@ -67,17 +64,18 @@
     -- encouraged to provide more precise implementation.
     complCumulative :: d -> Double -> Double
     complCumulative d x = 1 - cumulative d x
+    {-# MINIMAL (cumulative | complCumulative) #-}
 
+
 -- | Discrete probability distribution.
 class Distribution  d => DiscreteDistr d where
     -- | Probability of n-th outcome.
     probability :: d -> Int -> Double
     probability d = exp . logProbability d
-
     -- | Logarithm of probability of n-th outcome
     logProbability :: d -> Int -> Double
     logProbability d = log . probability d
-
+    {-# MINIMAL (probability | logProbability) #-}
 
 -- | Continuous probability distribution.
 --
@@ -86,25 +84,23 @@
 class Distribution d => ContDistr d where
     -- | Probability density function. Probability that random
     -- variable /X/ lies in the infinitesimal interval
-    -- [/x/,/x+/&#948;/x/) equal to /density(x)/&#8901;&#948;/x/
+    -- [/x/,/x+/δ/x/) equal to /density(x)/⋅δ/x/
     density :: d -> Double -> Double
     density d = exp . logDensity d
-
+    -- | Natural logarithm of density.
+    logDensity :: d -> Double -> Double
+    logDensity d = log . density d
     -- | Inverse of the cumulative distribution function. The value
-    -- /x/ for which P(/X/&#8804;/x/) = /p/. If probability is outside
+    -- /x/ for which P(/X/≤/x/) = /p/. If probability is outside
     -- of [0,1] range function should call 'error'
     quantile :: d -> Double -> Double
-
+    quantile d x = complQuantile d (1 - x)
     -- | 1-complement of @quantile@:
     --
     -- > complQuantile x ≡ quantile (1 - x)
     complQuantile :: d -> Double -> Double
     complQuantile d x = quantile d (1 - x)
-
-    -- | Natural logarithm of density.
-    logDensity :: d -> Double -> Double
-    logDensity d = log . density d
-
+    {-# MINIMAL (density | logDensity), (quantile | complQuantile) #-}
 
 -- | Type class for distributions with mean. 'maybeMean' should return
 --   'Nothing' if it's undefined for current value of data
@@ -126,9 +122,10 @@
 --   Minimal complete definition is 'maybeVariance' or 'maybeStdDev'
 class MaybeMean d => MaybeVariance d where
     maybeVariance :: d -> Maybe Double
-    maybeVariance d = (*) <$> x <*> x where x = maybeStdDev d
+    maybeVariance = fmap square . maybeStdDev
     maybeStdDev   :: d -> Maybe Double
-    maybeStdDev = fmap sqrt . maybeVariance
+    maybeStdDev   = fmap sqrt . maybeVariance
+    {-# MINIMAL (maybeVariance | maybeStdDev) #-}
 
 -- | Type class for distributions with variance. If distribution have
 --   finite variance for all valid parameter values it should be
@@ -140,7 +137,9 @@
     variance d = square (stdDev d)
     stdDev   :: d -> Double
     stdDev = sqrt . variance
+    {-# MINIMAL (variance | stdDev) #-}
 
+
 -- | 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.  'maybeEntropy' should return 'Nothing' if
@@ -161,35 +160,30 @@
 -- | Generate discrete random variates which have given
 --   distribution.
 class Distribution d => ContGen d where
-  genContVar :: PrimMonad m => d -> Gen (PrimState m) -> m Double
+  genContVar :: (StatefulGen g m) => d -> g -> m Double
 
 -- | Generate discrete random variates which have given
 --   distribution. 'ContGen' is superclass because it's always possible
 --   to generate real-valued variates from integer values
 class (DiscreteDistr d, ContGen d) => DiscreteGen d where
-  genDiscreteVar :: PrimMonad m => d -> Gen (PrimState m) -> m Int
+  genDiscreteVar :: (StatefulGen g m) => d -> g -> m Int
 
 -- | Estimate distribution from sample. First parameter in sample is
 --   distribution type and second is element type.
 class FromSample d a where
-  -- | Estimate distribution from sample. Returns nothing is there's
-  --   not enough data to estimate or sample clearly doesn't come from
-  --   distribution in question. For example if there's negative
-  --   samples in exponential distribution.
+  -- | Estimate distribution from sample. Returns 'Nothing' if there is
+  --   not enough data, or if no usable fit results from the method
+  --   used, e.g., the estimated distribution parameters would be
+  --   invalid or inaccurate.
   fromSample :: G.Vector v a => v a -> Maybe d
 
 
 -- | Generate variates from continuous distribution using inverse
 --   transform rule.
-genContinuous :: (ContDistr d, PrimMonad m) => d -> Gen (PrimState m) -> m Double
+genContinuous :: (ContDistr d, StatefulGen g m) => d -> g -> m Double
 genContinuous d gen = do
-  x <- uniform gen
+  x <- uniformDouble01M gen
   return $! quantile d x
-
--- | Backwards compatibility with genContinuous.
-genContinous :: (ContDistr d, PrimMonad m) => d -> Gen (PrimState m) -> m Double
-genContinous = genContinuous
-{-# DEPRECATED genContinous "Use genContinuous" #-}
 
 data P = P {-# UNPACK #-} !Double {-# UNPACK #-} !Double
 
diff --git a/Statistics/Distribution/Binomial.hs b/Statistics/Distribution/Binomial.hs
--- a/Statistics/Distribution/Binomial.hs
+++ b/Statistics/Distribution/Binomial.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE PatternGuards     #-}
 {-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
 -- |
 -- Module    : Statistics.Distribution.Binomial
@@ -31,7 +32,7 @@
 import Data.Data             (Data, Typeable)
 import GHC.Generics          (Generic)
 import Numeric.SpecFunctions           (choose,logChoose,incompleteBeta,log1p)
-import Numeric.MathFunctions.Constants (m_epsilon)
+import Numeric.MathFunctions.Constants (m_epsilon,m_tiny)
 
 import qualified Statistics.Distribution as D
 import qualified Statistics.Distribution.Poisson.Internal as I
@@ -70,6 +71,7 @@
 
 instance D.Distribution BinomialDistribution where
     cumulative = cumulative
+    complCumulative = complCumulative
 
 instance D.DiscreteDistr BinomialDistribution where
     probability    = probability
@@ -104,9 +106,16 @@
   | n == 0         = 1
     -- choose could overflow Double for n >= 1030 so we switch to
     -- log-domain to calculate probability
-  | n < 1000       = choose n k * p^k * (1-p)^(n-k)
-  | otherwise      = exp $ logChoose n k + log p * k' + log1p (-p) * nk'
+    --
+    -- We also want to avoid underflow when computing p^k &
+    -- (1-p)^(n-k).
+  | n < 1000
+  , pK  >= m_tiny
+  , pNK >= m_tiny = choose n k * pK * pNK
+  | otherwise     = exp $ logChoose n k + log p * k' + log1p (-p) * nk'
   where
+    pK  = p^k
+    pNK = (1-p)^(n-k)
     k'  = fromIntegral k
     nk' = fromIntegral $ n - k
 
@@ -119,7 +128,6 @@
     k'  = fromIntegral   k
     nk' = fromIntegral $ n - k
 
--- Summation from different sides required to reduce roundoff errors
 cumulative :: BinomialDistribution -> Double -> Double
 cumulative (BD n p) x
   | isNaN x      = error "Statistics.Distribution.Binomial.cumulative: NaN input"
@@ -130,6 +138,16 @@
   where
     k = floor x
 
+complCumulative :: BinomialDistribution -> Double -> Double
+complCumulative (BD n p) x
+  | isNaN x      = error "Statistics.Distribution.Binomial.complCumulative: NaN input"
+  | isInfinite x = if x > 0 then 0 else 1
+  | k <  0       = 1
+  | k >= n       = 0
+  | otherwise    = incompleteBeta (fromIntegral (k+1)) (fromIntegral (n-k)) p
+  where
+    k = floor x
+
 mean :: BinomialDistribution -> Double
 mean (BD n p) = fromIntegral n * p
 
@@ -157,7 +175,7 @@
           -> Maybe BinomialDistribution
 binomialE n p
   | n < 0            = Nothing
-  | p >= 0 || p <= 1 = Just (BD n p)
+  | p >= 0 && p <= 1 = Just (BD n p)
   | otherwise        = Nothing
 
 errMsg :: Int -> Double -> String
diff --git a/Statistics/Distribution/CauchyLorentz.hs b/Statistics/Distribution/CauchyLorentz.hs
--- a/Statistics/Distribution/CauchyLorentz.hs
+++ b/Statistics/Distribution/CauchyLorentz.hs
@@ -94,23 +94,43 @@
 
 
 instance D.Distribution CauchyDistribution where
-  cumulative (CD m s) x = 0.5 + atan( (x - m) / s ) / pi
+  cumulative (CD m s) x
+    | y < -1    = atan (-1/y) / pi
+    | otherwise = 0.5 + atan y / pi
+    where
+       y = (x - m) / s
+  complCumulative (CD m s) x
+    | y > 1     = atan (1/y) / pi
+    | otherwise = 0.5 - atan y / pi
+    where
+       y = (x - m) / s
 
 instance D.ContDistr CauchyDistribution where
   density (CD m s) x = (1 / pi) / (s * (1 + y*y))
     where y = (x - m) / s
   quantile (CD m s) p
-    | p > 0 && p < 1 = m + s * tan( pi * (p - 0.5) )
-    | p == 0         = -1 / 0
-    | p == 1         =  1 / 0
-    | otherwise      =
-      error $ "Statistics.Distribution.CauchyLorentz.quantile: p must be in [0,1] range. Got: "++show p
+    | p == 0    = -1 / 0
+    | p == 1    =  1 / 0
+    | p == 0.5  = m
+    | p < 0     = err
+    | p < 0.5   = m - s / tan( pi * p )
+    | p < 1     = m + s / tan( pi * (1 - p) )
+    | otherwise = err
+    where
+      err = error
+          $ "Statistics.Distribution.CauchyLorentz.quantile: p must be in [0,1] range. Got: "++show p
   complQuantile (CD m s) p
-    | p > 0 && p < 1 = m + s * tan( pi * (0.5 - p) )
-    | p == 0         =  1 / 0
-    | p == 1         = -1 / 0
-    | otherwise      =
-      error $ "Statistics.Distribution.CauchyLorentz.complQuantile: p must be in [0,1] range. Got: "++show p
+    | p == 0    =  1 / 0
+    | p == 1    = -1 / 0
+    | p == 0.5  = m
+    | p < 0     = err
+    | p < 0.5   = m + s / tan( pi * p )
+    | p < 1     = m - s / tan( pi * (1 - p) )
+    | otherwise = err
+    where
+      err = error
+          $ "Statistics.Distribution.CauchyLorentz.quantile: p must be in [0,1] range. Got: "++show p
+
 
 instance D.ContGen CauchyDistribution where
   genContVar = D.genContinuous
diff --git a/Statistics/Distribution/DiscreteUniform.hs b/Statistics/Distribution/DiscreteUniform.hs
--- a/Statistics/Distribution/DiscreteUniform.hs
+++ b/Statistics/Distribution/DiscreteUniform.hs
@@ -28,10 +28,11 @@
     , rangeTo
     ) where
 
-import Control.Applicative ((<$>), (<*>), empty)
+import Control.Applicative (empty)
 import Data.Aeson   (FromJSON(..), ToJSON, Value(..), (.:))
 import Data.Binary  (Binary(..))
 import Data.Data    (Data, Typeable)
+import System.Random.Stateful (uniformRM)
 import GHC.Generics (Generic)
 
 import qualified Statistics.Distribution as D
@@ -93,6 +94,12 @@
 
 instance D.MaybeEntropy DiscreteUniform where
   maybeEntropy = Just . D.entropy
+
+instance D.ContGen DiscreteUniform where
+  genContVar d = fmap fromIntegral . D.genDiscreteVar d
+
+instance D.DiscreteGen DiscreteUniform where
+  genDiscreteVar (U a b) = uniformRM (a,b)
 
 -- | Construct discrete uniform distribution on support {1, ..., n}.
 --   Range /n/ must be >0.
diff --git a/Statistics/Distribution/Exponential.hs b/Statistics/Distribution/Exponential.hs
--- a/Statistics/Distribution/Exponential.hs
+++ b/Statistics/Distribution/Exponential.hs
@@ -11,7 +11,7 @@
 -- Portability : portable
 --
 -- The exponential distribution.  This is the continuous probability
--- distribution of the times between events in a poisson process, in
+-- distribution of the times between events in a Poisson process, in
 -- which events occur continuously and independently at a constant
 -- average rate.
 
@@ -33,7 +33,6 @@
 import Numeric.SpecFunctions           (log1p,expm1)
 import Numeric.MathFunctions.Constants (m_neg_inf)
 import qualified System.Random.MWC.Distributions as MWC
-import qualified Data.Vector.Generic as G
 
 import qualified Statistics.Distribution         as D
 import qualified Statistics.Sample               as S
@@ -136,11 +135,9 @@
 errMsg :: Double -> String
 errMsg l = "Statistics.Distribution.Exponential.exponential: scale parameter must be positive. Got " ++ show l
 
--- | Create exponential distribution from sample. Returns @Nothing@ if
---   sample is empty or contains negative elements. No other tests are
---   made to check whether it truly is exponential.
+-- | Create exponential distribution from sample.  Estimates the rate
+--   with the maximum likelihood estimator, which is biased. Returns
+--   @Nothing@ if the sample mean does not exist or is not positive.
 instance D.FromSample ExponentialDistribution Double where
-  fromSample xs
-    | G.null xs       = Nothing
-    | G.all (>= 0) xs = Nothing
-    | otherwise       = Just $! ED (S.mean xs)
+  fromSample xs = let m = S.mean xs
+                  in  if m > 0 then Just (ED (1/m)) else Nothing
diff --git a/Statistics/Distribution/FDistribution.hs b/Statistics/Distribution/FDistribution.hs
--- a/Statistics/Distribution/FDistribution.hs
+++ b/Statistics/Distribution/FDistribution.hs
@@ -109,15 +109,36 @@
 cumulative :: FDistribution -> Double -> Double
 cumulative (F n m _) x
   | x <= 0       = 0
-  | isInfinite x = 1            -- Only matches +∞
-  | otherwise    = let y = n*x in incompleteBeta (0.5 * n) (0.5 * m) (y / (m + y))
+  -- Only matches +∞
+  | isInfinite x = 1
+  -- NOTE: Here we rely on implementation detail of incompleteBeta. It
+  --       computes using series expansion for sufficiently small x
+  --       and uses following identity otherwise:
+  --
+  --           I(x; a, b) = 1 - I(1-x; b, a)
+  --
+  --       Point is we can compute 1-x as m/(m+y) without loss of
+  --       precision for large x. Sadly this switchover point is
+  --       implementation detail.
+  | n >= (n+m)*bx = incompleteBeta (0.5 * n) (0.5 * m) bx
+  | otherwise     = 1 - incompleteBeta (0.5 * m) (0.5 * n) bx1
+  where
+    y   = n * x
+    bx  = y / (m + y)
+    bx1 = m / (m + y)
 
 complCumulative :: FDistribution -> Double -> Double
 complCumulative (F n m _) x
-  | x <= 0       = 1
-  | isInfinite x = 0            -- Only matches +∞
-  | otherwise    = let y = n*x
-                   in incompleteBeta (0.5 * m) (0.5 * n) (m / (m + y))
+  | x <= 0        = 1
+  -- Only matches +∞
+  | isInfinite x  = 0
+  -- See NOTE at cumulative
+  | m >= (n+m)*bx = incompleteBeta (0.5 * m) (0.5 * n) bx
+  | otherwise     = 1 - incompleteBeta (0.5 * n) (0.5 * m) bx1
+  where
+    y   = n*x
+    bx  = m / (m + y)
+    bx1 = y / (m + y)
 
 logDensity :: FDistribution -> Double -> Double
 logDensity (F n m fac) x
diff --git a/Statistics/Distribution/Gamma.hs b/Statistics/Distribution/Gamma.hs
--- a/Statistics/Distribution/Gamma.hs
+++ b/Statistics/Distribution/Gamma.hs
@@ -36,6 +36,7 @@
 import Numeric.MathFunctions.Constants (m_pos_inf, m_NaN, m_neg_inf)
 import Numeric.SpecFunctions (incompleteGamma, invIncompleteGamma, logGamma, digamma)
 import qualified System.Random.MWC.Distributions as MWC
+import qualified Numeric.Sum as Sum
 
 import Statistics.Distribution.Poisson.Internal as Poisson
 import qualified Statistics.Distribution as D
@@ -126,7 +127,11 @@
     density    = density
     logDensity (GD k theta) x
       | x <= 0    = m_neg_inf
-      | otherwise = log x * (k - 1) - (x / theta) - logGamma k - log theta * k
+      | otherwise = Sum.sum Sum.kbn [ log x * (k - 1)
+                                    , - (x / theta)
+                                    , - logGamma k
+                                    , - log theta * k
+                                    ]
     quantile   = quantile
 
 instance D.Variance GammaDistribution where
diff --git a/Statistics/Distribution/Geometric.hs b/Statistics/Distribution/Geometric.hs
--- a/Statistics/Distribution/Geometric.hs
+++ b/Statistics/Distribution/Geometric.hs
@@ -39,7 +39,7 @@
 import Data.Binary         (Binary(..))
 import Data.Data           (Data, Typeable)
 import GHC.Generics        (Generic)
-import Numeric.MathFunctions.Constants (m_pos_inf, m_neg_inf)
+import Numeric.MathFunctions.Constants (m_neg_inf)
 import Numeric.SpecFunctions           (log1p,expm1)
 import qualified System.Random.MWC.Distributions as MWC
 
@@ -81,10 +81,11 @@
 instance D.DiscreteDistr GeometricDistribution where
     probability (GD s) n
       | n < 1     = 0
-      | otherwise = s * (1-s) ** (fromIntegral n - 1)
+      | s >= 0.5  = s * (1 - s)^(n - 1)
+      | otherwise = s * (exp $ log1p (-s) * (fromIntegral n - 1))
     logProbability (GD s) n
        | n < 1     = m_neg_inf
-       | otherwise = log s + log (1-s) * (fromIntegral n - 1)
+       | otherwise = log s + log1p (-s) * (fromIntegral n - 1)
 
 
 instance D.Mean GeometricDistribution where
@@ -102,9 +103,8 @@
 
 instance D.Entropy GeometricDistribution where
   entropy (GD s)
-    | s == 0 = m_pos_inf
     | s == 1 = 0
-    | otherwise = negate $ (s * log s + (1-s) * log (1-s)) / s
+    | otherwise = -(s * log s + (1-s) * log1p (-s)) / s
 
 instance D.MaybeEntropy GeometricDistribution where
   maybeEntropy = Just . D.entropy
@@ -120,14 +120,18 @@
   | x < 1        = 0
   | isInfinite x = 1
   | isNaN      x = error "Statistics.Distribution.Geometric.cumulative: NaN input"
-  | otherwise    = negate $ expm1 $ fromIntegral (floor x :: Int) * log1p (-s)
+  | s >= 0.5     = 1 - (1 - s)^k
+  | otherwise    = negate $ expm1 $ fromIntegral k * log1p (-s)
+    where k = floor x :: Int
 
 complCumulative :: GeometricDistribution -> Double -> Double
 complCumulative (GD s) x
   | x < 1        = 1
   | isInfinite x = 0
-  | isNaN      x = error "Statistics.Distribution.Geometric.cumulative: NaN input"
-  | otherwise    = exp $ fromIntegral (floor x :: Int) * log1p (-s)
+  | isNaN      x = error "Statistics.Distribution.Geometric.complCumulative: NaN input"
+  | s >= 0.5     = (1 - s)^k
+  | otherwise    = exp $ fromIntegral k * log1p (-s)
+    where k = floor x :: Int
 
 
 -- | Create geometric distribution.
@@ -139,11 +143,11 @@
 geometricE :: Double                -- ^ Success rate
            -> Maybe GeometricDistribution
 geometricE x
-  | x >= 0 && x <= 1 = Just (GD x)
+  | x > 0 && x <= 1  = Just (GD x)
   | otherwise        = Nothing
 
 errMsg :: Double -> String
-errMsg x = "Statistics.Distribution.Geometric.geometric: probability must be in [0,1] range. Got " ++ show x
+errMsg x = "Statistics.Distribution.Geometric.geometric: probability must be in (0,1] range. Got " ++ show x
 
 
 ----------------------------------------------------------------
@@ -215,8 +219,8 @@
 geometric0E :: Double                -- ^ Success rate
             -> Maybe GeometricDistribution0
 geometric0E x
-  | x >= 0 && x <= 1 = Just (GD0 x)
+  | x > 0 && x <= 1  = Just (GD0 x)
   | otherwise        = Nothing
 
 errMsg0 :: Double -> String
-errMsg0 x = "Statistics.Distribution.Geometric.geometric0: probability must be in [0,1] range. Got " ++ show x
+errMsg0 x = "Statistics.Distribution.Geometric.geometric0: probability must be in (0,1] range. Got " ++ show x
diff --git a/Statistics/Distribution/Hypergeometric.hs b/Statistics/Distribution/Hypergeometric.hs
--- a/Statistics/Distribution/Hypergeometric.hs
+++ b/Statistics/Distribution/Hypergeometric.hs
@@ -71,6 +71,7 @@
 
 instance D.Distribution HypergeometricDistribution where
     cumulative = cumulative
+    complCumulative = complCumulative
 
 instance D.DiscreteDistr HypergeometricDistribution where
     probability    = probability
@@ -133,10 +134,10 @@
 
 errMsg :: Int -> Int -> Int -> String
 errMsg m l k
-  =  "Statistics.Distribution.Hypergeometric.hypergeometric: "
-  ++ "m=" ++ show m
-  ++ "l=" ++ show l
-  ++ "k=" ++ show k
+  =  "Statistics.Distribution.Hypergeometric.hypergeometric:"
+  ++ " m=" ++ show m
+  ++ " l=" ++ show l
+  ++ " k=" ++ show k
   ++ " should hold: l>0 & m in [0,l] & k in (0,l]"
 
 -- Naive implementation
@@ -164,6 +165,18 @@
   | n <  minN    = 0
   | n >= maxN    = 1
   | otherwise    = D.sumProbabilities d minN n
+  where
+    n    = floor x
+    minN = max 0 (mi+ki-li)
+    maxN = min mi ki
+
+complCumulative :: HypergeometricDistribution -> Double -> Double
+complCumulative d@(HD mi li ki) x
+  | isNaN x      = error "Statistics.Distribution.Hypergeometric.complCumulative: NaN argument"
+  | isInfinite x = if x > 0 then 0 else 1
+  | n <  minN    = 1
+  | n >= maxN    = 0
+  | otherwise    = D.sumProbabilities d (n + 1) maxN
   where
     n    = floor x
     minN = max 0 (mi+ki-li)
diff --git a/Statistics/Distribution/Laplace.hs b/Statistics/Distribution/Laplace.hs
--- a/Statistics/Distribution/Laplace.hs
+++ b/Statistics/Distribution/Laplace.hs
@@ -151,9 +151,9 @@
 errMsg _ s = "Statistics.Distribution.Laplace.laplace: scale parameter must be positive. Got " ++ show s
 
 
--- | Create Laplace distribution from sample. No tests are made to
---   check whether it truly is Laplace. Location of distribution
---   estimated as median of sample.
+-- | Create Laplace distribution from sample.  The location is estimated
+--   as the median of the sample, and the scale as the mean absolute
+--   deviation of the median.
 instance D.FromSample LaplaceDistribution Double where
   fromSample xs
     | G.null xs = Nothing
diff --git a/Statistics/Distribution/Lognormal.hs b/Statistics/Distribution/Lognormal.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Distribution/Lognormal.hs
@@ -0,0 +1,172 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
+-- |
+-- Module    : Statistics.Distribution.Lognormal
+-- Copyright : (c) 2020 Ximin Luo
+-- License   : BSD3
+--
+-- Maintainer  : infinity0@pwned.gg
+-- Stability   : experimental
+-- Portability : portable
+--
+-- The log normal distribution.  This is a continuous probability
+-- distribution that describes data whose log is clustered around a
+-- mean. For example, the multiplicative product of many independent
+-- positive random variables.
+
+module Statistics.Distribution.Lognormal
+    (
+      LognormalDistribution
+      -- * Constructors
+    , lognormalDistr
+    , lognormalDistrErr
+    , lognormalDistrMeanStddevErr
+    , lognormalStandard
+    ) where
+
+import Data.Aeson            (FromJSON, ToJSON)
+import Data.Binary           (Binary (..))
+import Data.Data             (Data, Typeable)
+import GHC.Generics          (Generic)
+import Numeric.MathFunctions.Constants (m_huge, m_sqrt_2_pi)
+import Numeric.SpecFunctions (expm1, log1p)
+import qualified Data.Vector.Generic as G
+
+import qualified Statistics.Distribution as D
+import qualified Statistics.Distribution.Normal as N
+import Statistics.Internal
+
+
+-- | The lognormal distribution.
+newtype LognormalDistribution = LND N.NormalDistribution
+    deriving (Eq, Typeable, Data, Generic)
+
+instance Show LognormalDistribution where
+  showsPrec i (LND d) = defaultShow2 "lognormalDistr" m s i
+   where
+    m = D.mean d
+    s = D.stdDev d
+instance Read LognormalDistribution where
+  readPrec = defaultReadPrecM2 "lognormalDistr" $
+    (either (const Nothing) Just .) . lognormalDistrErr
+
+instance ToJSON LognormalDistribution
+instance FromJSON LognormalDistribution
+
+instance Binary LognormalDistribution where
+  put (LND d) = put m >> put s
+   where
+    m = D.mean d
+    s = D.stdDev d
+  get = do
+    m  <- get
+    sd <- get
+    either fail return $ lognormalDistrErr m sd
+
+instance D.Distribution LognormalDistribution where
+  cumulative      = cumulative
+  complCumulative = complCumulative
+
+instance D.ContDistr LognormalDistribution where
+  logDensity    = logDensity
+  quantile      = quantile
+  complQuantile = complQuantile
+
+instance D.MaybeMean LognormalDistribution where
+  maybeMean = Just . D.mean
+
+instance D.Mean LognormalDistribution where
+  mean (LND d) = exp (m + v / 2)
+   where
+    m = D.mean d
+    v = D.variance d
+
+instance D.MaybeVariance LognormalDistribution where
+  maybeStdDev   = Just . D.stdDev
+  maybeVariance = Just . D.variance
+
+instance D.Variance LognormalDistribution where
+  variance (LND d) = expm1 v * exp (2 * m + v)
+   where
+    m = D.mean d
+    v = D.variance d
+
+instance D.Entropy LognormalDistribution where
+  entropy (LND d) = logBase 2 (s * exp (m + 0.5) * m_sqrt_2_pi)
+   where
+    m = D.mean d
+    s = D.stdDev d
+
+instance D.MaybeEntropy LognormalDistribution where
+  maybeEntropy = Just . D.entropy
+
+instance D.ContGen LognormalDistribution where
+  genContVar d = D.genContinuous d
+
+-- | Standard log normal distribution with mu 0 and sigma 1.
+--
+-- Mean is @sqrt e@ and variance is @(e - 1) * e@.
+lognormalStandard :: LognormalDistribution
+lognormalStandard = LND N.standard
+
+-- | Create log normal distribution from parameters.
+lognormalDistr
+  :: Double            -- ^ Mu
+  -> Double            -- ^ Sigma
+  -> LognormalDistribution
+lognormalDistr mu sig = either error id $ lognormalDistrErr mu sig
+
+-- | Create log normal distribution from parameters.
+lognormalDistrErr
+  :: Double            -- ^ Mu
+  -> Double            -- ^ Sigma
+  -> Either String LognormalDistribution
+lognormalDistrErr mu sig
+  | sig >= sqrt (log m_huge - 2 * mu) = Left $ errMsg mu sig
+  | otherwise = LND <$> N.normalDistrErr mu sig
+
+errMsg :: Double -> Double -> String
+errMsg mu sig =
+  "Statistics.Distribution.Lognormal.lognormalDistr: sigma must be > 0 && < "
+    ++ show lim ++ ". Got " ++ show sig
+  where lim = sqrt (log m_huge - 2 * mu)
+
+-- | Create log normal distribution from mean and standard deviation.
+lognormalDistrMeanStddevErr
+  :: Double            -- ^ Mu
+  -> Double            -- ^ Sigma
+  -> Either String LognormalDistribution
+lognormalDistrMeanStddevErr m sd = LND <$> N.normalDistrErr mu sig
+  where r = sd / m
+        sig2 = log1p (r * r)
+        sig = sqrt sig2
+        mu = log m - sig2 / 2
+
+-- | Variance is estimated using maximum likelihood method
+--   (biased estimation) over the log of the data.
+--
+--   Returns @Nothing@ if sample contains less than one element or
+--   variance is zero (all elements are equal)
+instance D.FromSample LognormalDistribution Double where
+  fromSample = fmap LND . D.fromSample . G.map log
+
+logDensity :: LognormalDistribution -> Double -> Double
+logDensity (LND d) x
+  | x > 0 = let lx = log x in D.logDensity d lx - lx
+  | otherwise = 0
+
+cumulative :: LognormalDistribution -> Double -> Double
+cumulative (LND d) x
+  | x > 0 = D.cumulative d $ log x
+  | otherwise = 0
+
+complCumulative :: LognormalDistribution -> Double -> Double
+complCumulative (LND d) x
+  | x > 0 = D.complCumulative d $ log x
+  | otherwise = 1
+
+quantile :: LognormalDistribution -> Double -> Double
+quantile (LND d) = exp . D.quantile d
+
+complQuantile :: LognormalDistribution -> Double -> Double
+complQuantile (LND d) = exp . D.complQuantile d
diff --git a/Statistics/Distribution/NegativeBinomial.hs b/Statistics/Distribution/NegativeBinomial.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Distribution/NegativeBinomial.hs
@@ -0,0 +1,188 @@
+{-# LANGUAGE OverloadedStrings, PatternGuards,
+             DeriveDataTypeable, DeriveGeneric #-}
+-- |
+-- Module    : Statistics.Distribution.NegativeBinomial
+-- Copyright : (c) 2022 Lorenz Minder
+-- License   : BSD3
+--
+-- Maintainer  : lminder@gmx.net
+-- Stability   : experimental
+-- Portability : portable
+--
+-- The negative binomial distribution.  This is the discrete probability
+-- distribution of the number of failures in a sequence of independent
+-- yes\/no experiments before a specified number of successes /r/.  Each
+-- Bernoulli trial has success probability /p/ in the range (0, 1].  The
+-- parameter /r/ must be positive, but does not have to be integer.
+
+module Statistics.Distribution.NegativeBinomial (
+      NegativeBinomialDistribution
+    -- * Constructors
+    , negativeBinomial
+    , negativeBinomialE
+    -- * Accessors
+    , nbdSuccesses
+    , nbdProbability
+) where
+
+import Control.Applicative
+import Data.Aeson                       (FromJSON(..), ToJSON, Value(..), (.:))
+import Data.Binary                      (Binary(..))
+import Data.Data                        (Data, Typeable)
+import Data.Foldable                    (foldl')
+import GHC.Generics                     (Generic)
+import Numeric.SpecFunctions            (incompleteBeta, log1p)
+import Numeric.SpecFunctions.Extra      (logChooseFast)
+import Numeric.MathFunctions.Constants  (m_epsilon, m_tiny)
+
+import qualified Statistics.Distribution as D
+import Statistics.Internal
+
+-- Math helper functions
+
+-- | Generalized binomial coefficients.
+--
+--   These computes binomial coefficients with the small generalization
+--   that the /n/ need not be integer, but can be real.
+gChoose :: Double -> Int -> Double
+gChoose n k
+    | k < 0             = 0
+    | k' >= 50          = exp $ logChooseFast n k'
+    | otherwise         = foldl' (*) 1 factors
+    where   factors = [ (n - k' + j) / j | j <- [1..k'] ]
+            k' = fromIntegral k
+
+
+-- Implementation of Negative Binomial
+
+-- | The negative binomial distribution.
+data NegativeBinomialDistribution = NBD {
+      nbdSuccesses   :: {-# UNPACK #-} !Double
+    -- ^ Number of successes until stop
+    , nbdProbability :: {-# UNPACK #-} !Double
+    -- ^ Success probability.
+    } deriving (Eq, Typeable, Data, Generic)
+
+instance Show NegativeBinomialDistribution where
+  showsPrec i (NBD r p) = defaultShow2 "negativeBinomial" r p i
+instance Read NegativeBinomialDistribution where
+  readPrec = defaultReadPrecM2 "negativeBinomial" negativeBinomialE
+
+instance ToJSON NegativeBinomialDistribution
+instance FromJSON NegativeBinomialDistribution where
+  parseJSON (Object v) = do
+    r <- v .: "nbdSuccesses"
+    p <- v .: "nbdProbability"
+    maybe (fail $ errMsg r p) return $ negativeBinomialE r p
+  parseJSON _ = empty
+
+instance Binary NegativeBinomialDistribution where
+  put (NBD r p) = put r >> put p
+  get = do
+    r <- get
+    p <- get
+    maybe (fail $ errMsg r p) return $ negativeBinomialE r p
+
+instance D.Distribution NegativeBinomialDistribution where
+    cumulative = cumulative
+    complCumulative = complCumulative
+
+instance D.DiscreteDistr NegativeBinomialDistribution where
+    probability    = probability
+    logProbability = logProbability
+
+instance D.Mean NegativeBinomialDistribution where
+    mean = mean
+
+instance D.Variance NegativeBinomialDistribution where
+    variance = variance
+
+instance D.MaybeMean NegativeBinomialDistribution where
+    maybeMean = Just . D.mean
+
+instance D.MaybeVariance NegativeBinomialDistribution where
+    maybeStdDev   = Just . D.stdDev
+    maybeVariance = Just . D.variance
+
+instance D.Entropy NegativeBinomialDistribution where
+   entropy = directEntropy
+
+instance D.MaybeEntropy NegativeBinomialDistribution where
+   maybeEntropy = Just . D.entropy
+
+-- This could be slow for big n
+probability :: NegativeBinomialDistribution -> Int -> Double
+probability d@(NBD r p) k
+  | k < 0          = 0
+    -- Switch to log domain for large k + r to avoid overflows.
+    --
+    -- We also want to avoid underflow when computing (1-p)^k &
+    -- p^r.
+  | k' + r < 1000
+  , pK >= m_tiny
+  , pR >= m_tiny  = gChoose (k' + r - 1) k * pK * pR
+  | otherwise     = exp $ logProbability d k
+  where
+    pK  = exp $ log1p (-p) * k'
+    pR  = p**r
+    k'  = fromIntegral k
+
+logProbability :: NegativeBinomialDistribution -> Int -> Double
+logProbability (NBD r p) k
+  | k < 0                   = (-1)/0
+  | otherwise               = logChooseFast (k' + r - 1) k'
+                              + log1p (-p) * k'
+                              + log p * r
+  where k' = fromIntegral k
+
+cumulative :: NegativeBinomialDistribution -> Double -> Double
+cumulative (NBD r p) x
+  | isNaN x      = error "Statistics.Distribution.NegativeBinomial.cumulative: NaN input"
+  | isInfinite x = if x > 0 then 1 else 0
+  | k < 0        = 0
+  | otherwise    = incompleteBeta r (fromIntegral (k+1)) p
+  where
+    k = floor x :: Integer
+
+complCumulative :: NegativeBinomialDistribution -> Double -> Double
+complCumulative (NBD r p) x
+  | isNaN x      = error "Statistics.Distribution.NegativeBinomial.complCumulative: NaN input"
+  | isInfinite x = if x > 0 then 0 else 1
+  | k < 0        = 1
+  | otherwise    = incompleteBeta (fromIntegral (k+1)) r (1 - p)
+  where
+    k = floor x :: Integer
+
+mean :: NegativeBinomialDistribution -> Double
+mean (NBD r p) = r * (1 - p)/p
+
+variance :: NegativeBinomialDistribution -> Double
+variance (NBD r p) = r * (1 - p)/(p * p)
+
+directEntropy :: NegativeBinomialDistribution -> Double
+directEntropy d =
+  negate . sum $
+  takeWhile (< -m_epsilon) $
+  dropWhile (>= -m_epsilon) $
+  [ let x = probability d k in x * log x | k <- [0..]]
+
+-- | Construct negative binomial distribution. Number of successes /r/
+--   must be positive and probability must be in (0,1] range
+negativeBinomial :: Double              -- ^ Number of successes.
+                 -> Double              -- ^ Success probability.
+                 -> NegativeBinomialDistribution
+negativeBinomial r p = maybe (error $ errMsg r p) id $ negativeBinomialE r p
+
+-- | Construct negative binomial distribution. Number of successes /r/
+--   must be positive and probability must be in (0,1] range
+negativeBinomialE :: Double              -- ^ Number of successes.
+                  -> Double              -- ^ Success probability.
+                  -> Maybe NegativeBinomialDistribution
+negativeBinomialE r p
+  | r > 0 && 0 < p && p <= 1            = Just (NBD r p)
+  | otherwise                           = Nothing
+
+errMsg :: Double -> Double -> String
+errMsg r p
+  = "Statistics.Distribution.NegativeBinomial.negativeBinomial: r=" ++ show r
+  ++ " p=" ++ show p ++ ", but need r>0 and p in (0,1]"
diff --git a/Statistics/Distribution/Normal.hs b/Statistics/Distribution/Normal.hs
--- a/Statistics/Distribution/Normal.hs
+++ b/Statistics/Distribution/Normal.hs
@@ -19,6 +19,7 @@
     -- * Constructors
     , normalDistr
     , normalDistrE
+    , normalDistrErr
     , standard
     ) where
 
@@ -55,7 +56,7 @@
   parseJSON (Object v) = do
     m  <- v .: "mean"
     sd <- v .: "stdDev"
-    maybe (fail $ errMsg m sd) return $ normalDistrE m sd
+    either fail return $ normalDistrErr m sd
   parseJSON _ = empty
 
 instance Binary NormalDistribution where
@@ -63,7 +64,7 @@
     get = do
       m  <- get
       sd <- get
-      maybe (fail $ errMsg m sd) return $ normalDistrE m sd
+      either fail return $ normalDistrErr m sd
 
 instance D.Distribution NormalDistribution where
     cumulative      = cumulative
@@ -111,7 +112,7 @@
 normalDistr :: Double            -- ^ Mean of distribution
             -> Double            -- ^ Standard deviation of distribution
             -> NormalDistribution
-normalDistr m sd = maybe (error $ errMsg m sd) id $ normalDistrE m sd
+normalDistr m sd = either error id $ normalDistrErr m sd
 
 -- | Create normal distribution from parameters.
 --
@@ -120,13 +121,20 @@
 normalDistrE :: Double            -- ^ Mean of distribution
              -> Double            -- ^ Standard deviation of distribution
              -> Maybe NormalDistribution
-normalDistrE m sd
-  | sd > 0    = Just ND { mean       = m
-                        , stdDev     = sd
-                        , ndPdfDenom = log $ m_sqrt_2_pi * sd
-                        , ndCdfDenom = m_sqrt_2 * sd
-                        }
-  | otherwise = Nothing
+normalDistrE m sd = either (const Nothing) Just $ normalDistrErr m sd
+
+-- | Create normal distribution from parameters.
+--
+normalDistrErr :: Double            -- ^ Mean of distribution
+               -> Double            -- ^ Standard deviation of distribution
+               -> Either String NormalDistribution
+normalDistrErr m sd
+  | sd > 0    = Right $ ND { mean       = m
+                           , stdDev     = sd
+                           , ndPdfDenom = log $ m_sqrt_2_pi * sd
+                           , ndCdfDenom = m_sqrt_2 * sd
+                           }
+  | otherwise = Left $ errMsg m sd
 
 errMsg :: Double -> Double -> String
 errMsg _ sd = "Statistics.Distribution.Normal.normalDistr: standard deviation must be positive. Got " ++ show sd
diff --git a/Statistics/Distribution/Poisson.hs b/Statistics/Distribution/Poisson.hs
--- a/Statistics/Distribution/Poisson.hs
+++ b/Statistics/Distribution/Poisson.hs
@@ -31,15 +31,18 @@
 import Data.Binary          (Binary(..))
 import Data.Data            (Data, Typeable)
 import GHC.Generics         (Generic)
+
+import qualified System.Random.MWC.Distributions as MWC
+
 import Numeric.SpecFunctions (incompleteGamma,logFactorial)
 import Numeric.MathFunctions.Constants (m_neg_inf)
 
+
 import qualified Statistics.Distribution as D
 import qualified Statistics.Distribution.Poisson.Internal as I
 import Statistics.Internal
 
 
-
 newtype PoissonDistribution = PD {
       poissonLambda :: Double
     } deriving (Eq, Typeable, Data, Generic)
@@ -92,6 +95,14 @@
 
 instance D.MaybeEntropy PoissonDistribution where
   maybeEntropy = Just . D.entropy
+
+-- | @since 0.16.5.0
+instance D.DiscreteGen PoissonDistribution where
+  genDiscreteVar (PD lambda) = MWC.poisson lambda
+
+-- | @since 0.16.5.0
+instance D.ContGen PoissonDistribution where
+  genContVar (PD lambda) gen = fromIntegral <$> MWC.poisson lambda gen
 
 -- | Create Poisson distribution.
 poisson :: Double -> PoissonDistribution
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
@@ -33,7 +33,7 @@
                            (m_sqrt_2_pi * sqrt x)
 
 -- -- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy
--- -- of the Poisson Law".  This function is unused because 'directEntorpy'
+-- -- of the Poisson Law".  This function is unused because 'directEntropy'
 -- -- is just as accurate and is faster by about a factor of 4.
 -- alyThm1 :: Double -> Double
 -- alyThm1 lambda =
@@ -61,7 +61,7 @@
   1.4189385332046727 + 0.5 * log lambda +
   zipCoefficients lambda coefficients
 
--- | Returns the average of the upper and lower bounds accounding to
+-- | Returns the average of the upper and lower bounds according to
 -- theorem 2.
 alyThm2 :: Double -> [Double] -> [Double] -> Double
 alyThm2 lambda upper lower =
@@ -164,7 +164,7 @@
   dropWhile (not . (< negate m_epsilon * lambda)) $
   [ let x = probability lambda k in x * log x | k <- [0..]]
 
--- | Compute the entropy of a poisson distribution using the best available
+-- | Compute the entropy of a Poisson distribution using the best available
 -- method.
 poissonEntropy :: Double -> Double
 poissonEntropy lambda
diff --git a/Statistics/Distribution/StudentT.hs b/Statistics/Distribution/StudentT.hs
--- a/Statistics/Distribution/StudentT.hs
+++ b/Statistics/Distribution/StudentT.hs
@@ -26,7 +26,7 @@
 import Data.Data           (Data, Typeable)
 import GHC.Generics        (Generic)
 import Numeric.SpecFunctions (
-  logBeta, incompleteBeta, invIncompleteBeta, digamma)
+  logBeta, incompleteBeta, invIncompleteBeta, digamma, log1p)
 
 import qualified Statistics.Distribution as D
 import Statistics.Distribution.Transform (LinearTransform (..))
@@ -94,8 +94,9 @@
 
 
 logDensityUnscaled :: StudentT -> Double -> Double
-logDensityUnscaled (StudentT ndf) x =
-    log (ndf / (ndf + x*x)) * (0.5 * (1 + ndf)) - logBeta 0.5 (0.5 * ndf)
+logDensityUnscaled (StudentT ndf) x
+  = log1p (x*x/ndf) * (-(0.5 * (1 + ndf)))
+  - logBeta 0.5 (0.5 * ndf)
 
 quantile :: StudentT -> Double -> Double
 quantile (StudentT ndf) p
diff --git a/Statistics/Distribution/Transform.hs b/Statistics/Distribution/Transform.hs
--- a/Statistics/Distribution/Transform.hs
+++ b/Statistics/Distribution/Transform.hs
@@ -18,11 +18,9 @@
     ) where
 
 import Data.Aeson (FromJSON, ToJSON)
-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 qualified Statistics.Distribution as D
 
diff --git a/Statistics/Distribution/Uniform.hs b/Statistics/Distribution/Uniform.hs
--- a/Statistics/Distribution/Uniform.hs
+++ b/Statistics/Distribution/Uniform.hs
@@ -22,11 +22,11 @@
     ) where
 
 import Control.Applicative
-import Data.Aeson          (FromJSON(..), ToJSON, Value(..), (.:))
-import Data.Binary         (Binary(..))
-import Data.Data           (Data, Typeable)
-import GHC.Generics        (Generic)
-import qualified System.Random.MWC       as MWC
+import Data.Aeson             (FromJSON(..), ToJSON, Value(..), (.:))
+import Data.Binary            (Binary(..))
+import Data.Data              (Data, Typeable)
+import System.Random.Stateful (uniformRM)
+import GHC.Generics           (Generic)
 
 import qualified Statistics.Distribution as D
 import Statistics.Internal
@@ -117,4 +117,4 @@
   maybeEntropy = Just . D.entropy
 
 instance D.ContGen UniformDistribution where
-    genContVar (UniformDistribution a b) gen = MWC.uniformR (a,b) gen
+    genContVar (UniformDistribution a b) = uniformRM (a,b)
diff --git a/Statistics/Distribution/Weibull.hs b/Statistics/Distribution/Weibull.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Distribution/Weibull.hs
@@ -0,0 +1,224 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
+-- |
+-- Module    : Statistics.Distribution.Lognormal
+-- Copyright : (c) 2020 Ximin Luo
+-- License   : BSD3
+--
+-- Maintainer  : infinity0@pwned.gg
+-- Stability   : experimental
+-- Portability : portable
+--
+-- The Weibull distribution.  This is a continuous probability
+-- distribution that describes the occurrence of a single event whose
+-- probability changes over time, controlled by the shape parameter.
+
+module Statistics.Distribution.Weibull
+    (
+      WeibullDistribution
+      -- * Constructors
+    , weibullDistr
+    , weibullDistrErr
+    , weibullStandard
+    , weibullDistrApproxMeanStddevErr
+    ) where
+
+import Control.Applicative
+import Data.Aeson            (FromJSON(..), ToJSON, Value(..), (.:))
+import Data.Binary           (Binary(..))
+import Data.Data             (Data, Typeable)
+import GHC.Generics          (Generic)
+import Numeric.MathFunctions.Constants (m_eulerMascheroni)
+import Numeric.SpecFunctions (expm1, log1p, logGamma)
+import qualified Data.Vector.Generic as G
+
+import qualified Statistics.Distribution as D
+import qualified Statistics.Sample as S
+import Statistics.Internal
+
+
+-- | The Weibull distribution.
+data WeibullDistribution = WD {
+      wdShape  :: {-# UNPACK #-} !Double
+    , wdLambda :: {-# UNPACK #-} !Double
+    } deriving (Eq, Typeable, Data, Generic)
+
+instance Show WeibullDistribution where
+  showsPrec i (WD k l) = defaultShow2 "weibullDistr" k l i
+instance Read WeibullDistribution where
+  readPrec = defaultReadPrecM2 "weibullDistr" $
+    (either (const Nothing) Just .) . weibullDistrErr
+
+instance ToJSON WeibullDistribution
+instance FromJSON WeibullDistribution where
+  parseJSON (Object v) = do
+    k <- v .: "wdShape"
+    l <- v .: "wdLambda"
+    either fail return $ weibullDistrErr k l
+  parseJSON _ = empty
+
+instance Binary WeibullDistribution where
+  put (WD k l) = put k >> put l
+  get = do
+    k <- get
+    l <- get
+    either fail return $ weibullDistrErr k l
+
+instance D.Distribution WeibullDistribution where
+  cumulative      = cumulative
+  complCumulative = complCumulative
+
+instance D.ContDistr WeibullDistribution where
+  logDensity    = logDensity
+  quantile      = quantile
+  complQuantile = complQuantile
+
+instance D.MaybeMean WeibullDistribution where
+  maybeMean = Just . D.mean
+
+instance D.Mean WeibullDistribution where
+  mean (WD k l) = l * exp (logGamma (1 + 1 / k))
+
+instance D.MaybeVariance WeibullDistribution where
+  maybeStdDev   = Just . D.stdDev
+  maybeVariance = Just . D.variance
+
+instance D.Variance WeibullDistribution where
+  variance (WD k l) = l * l * (exp (logGamma (1 + 2 * invk)) - q * q)
+   where
+    invk = 1 / k
+    q    = exp (logGamma (1 + invk))
+
+instance D.Entropy WeibullDistribution where
+  entropy (WD k l) = m_eulerMascheroni * (1 - 1 / k) + log (l / k) + 1
+
+instance D.MaybeEntropy WeibullDistribution where
+  maybeEntropy = Just . D.entropy
+
+instance D.ContGen WeibullDistribution where
+  genContVar d = D.genContinuous d
+
+-- | Standard Weibull distribution with scale factor (lambda) 1.
+weibullStandard :: Double -> WeibullDistribution
+weibullStandard k = weibullDistr k 1.0
+
+-- | Create Weibull distribution from parameters.
+--
+-- If the shape (first) parameter is @1.0@, the distribution is equivalent to a
+-- 'Statistics.Distribution.Exponential.ExponentialDistribution' with parameter
+-- @1 / lambda@ the scale (second) parameter.
+weibullDistr
+  :: Double            -- ^ Shape
+  -> Double            -- ^ Lambda (scale)
+  -> WeibullDistribution
+weibullDistr k l = either error id $ weibullDistrErr k l
+
+-- | Create Weibull distribution from parameters.
+--
+-- If the shape (first) parameter is @1.0@, the distribution is equivalent to a
+-- 'Statistics.Distribution.Exponential.ExponentialDistribution' with parameter
+-- @1 / lambda@ the scale (second) parameter.
+weibullDistrErr
+  :: Double            -- ^ Shape
+  -> Double            -- ^ Lambda (scale)
+  -> Either String WeibullDistribution
+weibullDistrErr k l | k <= 0     = Left $ errMsg k l
+                    | l <= 0     = Left $ errMsg k l
+                    | otherwise = Right $ WD k l
+
+errMsg :: Double -> Double -> String
+errMsg k l =
+  "Statistics.Distribution.Weibull.weibullDistr: both shape and lambda must be positive. Got shape "
+    ++ show k
+    ++ " and lambda "
+    ++ show l
+
+-- | Create Weibull distribution from mean and standard deviation.
+--
+-- The algorithm is from "Methods for Estimating Wind Speed Frequency
+-- Distributions", C. G. Justus, W. R. Hargreaves, A. Mikhail, D. Graber, 1977.
+-- Given the identity:
+--
+-- \[
+-- (\frac{\sigma}{\mu})^2 = \frac{\Gamma(1+2/k)}{\Gamma(1+1/k)^2} - 1
+-- \]
+--
+-- \(k\) can be approximated by
+--
+-- \[
+-- k \approx (\frac{\sigma}{\mu})^{-1.086}
+-- \]
+--
+-- \(\lambda\) is then calculated straightforwardly via the identity
+--
+-- \[
+-- \lambda = \frac{\mu}{\Gamma(1+1/k)}
+-- \]
+--
+-- Numerically speaking, the approximation for \(k\) is accurate only within a
+-- certain range. We arbitrarily pick the range \(0.033 \le \frac{\sigma}{\mu} \le 1.45\)
+-- where it is good to ~6%, and will refuse to create a distribution outside of
+-- this range. The paper does not cover these details but it is straightforward
+-- to check them numerically.
+weibullDistrApproxMeanStddevErr
+  :: Double            -- ^ Mean
+  -> Double            -- ^ Stddev
+  -> Either String WeibullDistribution
+weibullDistrApproxMeanStddevErr m s = if r > 1.45 || r < 0.033
+    then Left msg
+    else weibullDistrErr k l
+  where r = s / m
+        k = (s / m) ** (-1.086)
+        l = m / exp (logGamma (1 + 1/k))
+        msg = "Statistics.Distribution.Weibull.weibullDistr: stddev-mean ratio "
+          ++ "outside approximation accuracy range [0.033, 1.45]. Got "
+          ++ "stddev " ++ show s ++ " and mean " ++ show m
+
+-- | Uses an approximation based on the mean and standard deviation in
+--   'weibullDistrEstMeanStddevErr', with standard deviation estimated
+--   using maximum likelihood method (unbiased estimation).
+--
+--   Returns @Nothing@ if sample contains less than one element or
+--   variance is zero (all elements are equal), or if the estimated mean
+--   and standard-deviation lies outside the range for which the
+--   approximation is accurate.
+instance D.FromSample WeibullDistribution Double where
+  fromSample xs
+    | G.length xs <= 1 = Nothing
+    | v == 0           = Nothing
+    | otherwise        = either (const Nothing) Just $
+      weibullDistrApproxMeanStddevErr m (sqrt v)
+    where
+      (m,v) = S.meanVarianceUnb xs
+
+logDensity :: WeibullDistribution -> Double -> Double
+logDensity (WD k l) x
+  | x < 0     = 0
+  | otherwise = log k + (k - 1) * log x - k * log l - (x / l) ** k
+
+cumulative :: WeibullDistribution -> Double -> Double
+cumulative (WD k l) x | x < 0     = 0
+                      | otherwise = -expm1 (-(x / l) ** k)
+
+complCumulative :: WeibullDistribution -> Double -> Double
+complCumulative (WD k l) x | x < 0     = 1
+                           | otherwise = exp (-(x / l) ** k)
+
+quantile :: WeibullDistribution -> Double -> Double
+quantile (WD k l) p
+  | p == 0         = 0
+  | p == 1         = inf
+  | p > 0 && p < 1 = l * (-log1p (-p)) ** (1 / k)
+  | otherwise      =
+    error $ "Statistics.Distribution.Weibull.quantile: p must be in [0,1] range. Got: " ++ show p
+  where inf = 1 / 0
+
+complQuantile :: WeibullDistribution -> Double -> Double
+complQuantile (WD k l) q
+  | q == 0         = inf
+  | q == 1         = 0
+  | q > 0 && q < 1 = l * (-log q) ** (1 / k)
+  | otherwise      =
+    error $ "Statistics.Distribution.Weibull.complQuantile: q must be in [0,1] range. Got: " ++ show q
+  where inf = 1 / 0
diff --git a/Statistics/Function.hs b/Statistics/Function.hs
--- a/Statistics/Function.hs
+++ b/Statistics/Function.hs
@@ -76,8 +76,8 @@
 {-# INLINE indices #-}
 
 -- | Zip a vector with its indices.
-indexed :: (G.Vector v e, G.Vector v Int, G.Vector v (Int,e)) => v e -> v (Int,e)
-indexed a = G.zip (indices a) a
+indexed :: (G.Vector v e, G.Vector v (Int,e)) => v e -> v (Int,e)
+indexed xs = G.imap (,) xs
 {-# INLINE indexed #-}
 
 data MM = MM {-# UNPACK #-} !Double {-# UNPACK #-} !Double
diff --git a/Statistics/Internal.hs b/Statistics/Internal.hs
--- a/Statistics/Internal.hs
+++ b/Statistics/Internal.hs
@@ -25,8 +25,6 @@
 import Control.Applicative
 import Control.Monad
 import Text.Read
-import Data.Orphans ()
-
 
 
 ----------------------------------------------------------------
diff --git a/Statistics/Quantile.hs b/Statistics/Quantile.hs
--- a/Statistics/Quantile.hs
+++ b/Statistics/Quantile.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE CPP                #-}
 {-# LANGUAGE DeriveDataTypeable #-}
 {-# LANGUAGE DeriveFoldable     #-}
 {-# LANGUAGE DeriveFunctor      #-}
@@ -55,7 +54,6 @@
 import           Data.Aeson             (ToJSON,FromJSON)
 import           Data.Data              (Data,Typeable)
 import           Data.Default.Class
-import           Data.Functor
 import qualified Data.Foldable        as F
 import           Data.Vector.Generic ((!))
 import qualified Data.Vector          as V
@@ -196,7 +194,7 @@
   | F.any (badQ nQ) qs = modErr "quantiles" "Wrong quantile number"
   | G.any isNaN xs     = modErr "quantiles" "Sample contains NaNs"
   -- Doesn't matter what we put into empty container
-  | fnull qs           = 0 <$ qs
+  | null qs            = 0 <$ qs
   | otherwise          = fmap (estimateQuantile sortedXs) ks'
   where
     ks'      = fmap (\q -> toPk param n q nQ) qs
@@ -209,14 +207,6 @@
       :: (Functor f, F.Foldable f) => ContParam -> f Int -> Int -> U.Vector Double -> f Double #-}
 {-# SPECIALIZE quantiles
       :: (Functor f, F.Foldable f) => ContParam -> f Int -> Int -> S.Vector Double -> f Double #-}
-
--- COMPAT
-fnull :: F.Foldable f => f a -> Bool
-#if !MIN_VERSION_base(4,8,0)
-fnull = F.foldr (\_ _ -> False) True
-#else
-fnull = null
-#endif
 
 -- | O(/k·n/·log /n/). Same as quantiles but uses 'G.Vector' container
 --   instead of 'Foldable' one.
diff --git a/Statistics/Regression.hs b/Statistics/Regression.hs
--- a/Statistics/Regression.hs
+++ b/Statistics/Regression.hs
@@ -13,7 +13,6 @@
     , bootstrapRegress
     ) where
 
-import Control.Applicative ((<$>))
 import Control.Concurrent.Async (forConcurrently)
 import Control.DeepSeq (rnf)
 import Control.Monad (when)
@@ -42,8 +41,15 @@
 --   element than the list of predictors; the last element is the
 --   /y/-intercept value.
 --
--- * /R&#0178;/, the coefficient of determination (see 'rSquare' for
+-- * /R²/, the coefficient of determination (see 'rSquare' for
 --   details).
+--
+-- >>> import qualified Data.Vector.Unboxed as VU
+-- >>> :{
+--  olsRegress [ VU.fromList [0,1,2,3]
+--             ] (VU.fromList [1000, 1001, 1002, 1003])
+-- :}
+-- ([1.0000000000000218,999.9999999999999],1.0)
 olsRegress :: [Vector]
               -- ^ Non-empty list of predictor vectors.  Must all have
               -- the same length.  These will become the columns of
@@ -66,7 +72,30 @@
     lss@(n:ls) = map G.length preds
 olsRegress _ _ = error "no predictors given"
 
--- | Compute the ordinary least-squares solution to /A x = b/.
+-- | Compute the ordinary least-squares solution to overdetermined
+--   linear system \(Ax = b\). In other words it finds
+--
+--   \[ \operatorname{argmin}|Ax-b|^2 \].
+--
+--   All columns of \(A\) must be linearly independent. It's not
+--   checked function will return nonsensical result if resulting
+--   linear system is poorly conditioned.
+--
+-- >>> import qualified Data.Vector.Unboxed as VU
+-- >>> :{
+--  ols (fromColumns [ VU.fromList [0,1,2,3]
+--                   , VU.fromList [1,1,1,1]
+--                   ]) (VU.fromList [1000, 1001, 1002, 1003])
+-- :}
+-- [1.0000000000000218,999.9999999999999]
+--
+-- >>> :{
+--  ols (fromColumns [ VU.fromList [0,1,2,3]
+--                   , VU.fromList [4,2,1,1]
+--                   , VU.fromList [1,1,1,1]
+--                   ]) (VU.fromList [1000, 1001, 1002, 1003])
+-- :}
+-- [1.0000000000005393,4.2290644612446807e-13,999.9999999999983]
 ols :: Matrix     -- ^ /A/ has at least as many rows as columns.
     -> Vector     -- ^ /b/ has the same length as columns in /A/.
     -> Vector
@@ -93,7 +122,7 @@
   where n = rows r
         l = U.length b
 
--- | Compute /R&#0178;/, the coefficient of determination that
+-- | Compute /R²/, the coefficient of determination that
 -- indicates goodness-of-fit of a regression.
 --
 -- This value will be 1 if the predictors fit perfectly, dropping to 0
@@ -102,11 +131,19 @@
         -> Vector               -- ^ Responders.
         -> Vector               -- ^ Regression coefficients.
         -> Double
-rSquare pred resp coeff = 1 - r / t
+rSquare pred resp coeff
+  -- Data has zero variance. If fit is perfect we set R² to 1 else to
+  -- 0. This is not perfect heuristic. Fit residuals may be nonzero
+  -- due to rounding.
+  | t == 0             = if r == 0 then 1 else 0
+  -- If fit residuals are worse than average we simply set R² to 0
+  | r2 >= 0 && r2 <= 1 = r2
+  | otherwise          = 0
   where
-    r   = sum $ flip U.imap resp $ \i x -> square (x - p i)
-    t   = sum $ flip U.map resp $ \x -> square (x - mean resp)
-    p i = sum . flip U.imap coeff $ \j -> (* unsafeIndex pred i j)
+    r2  = 1 - r / t
+    r   = sum $ flip U.imap resp  $ \i x -> square (x - p i)
+    t   = sum $ flip U.map  resp  $ \x   -> square (x - mean resp)
+    p i = sum $ flip U.imap coeff $ \j x -> x * unsafeIndex pred i j
 
 -- | Bootstrap a regression function.  Returns both the results of the
 -- regression and the requested confidence interval values.
diff --git a/Statistics/Resampling.hs b/Statistics/Resampling.hs
--- a/Statistics/Resampling.hs
+++ b/Statistics/Resampling.hs
@@ -1,5 +1,4 @@
 {-# LANGUAGE BangPatterns       #-}
-{-# LANGUAGE CPP                #-}
 {-# LANGUAGE DeriveDataTypeable #-}
 {-# LANGUAGE DeriveFoldable     #-}
 {-# LANGUAGE DeriveFunctor      #-}
@@ -40,7 +39,6 @@
     ) where
 
 import Data.Aeson (FromJSON, ToJSON)
-import Control.Applicative
 import Control.Concurrent.Async (forConcurrently_)
 import Control.Monad (forM_, forM, replicateM, liftM2)
 import Control.Monad.Primitive (PrimMonad(..))
@@ -88,9 +86,7 @@
   , resamples  :: v a
   }
   deriving (Eq, Read, Show , Generic, Functor, T.Foldable, T.Traversable
-#if __GLASGOW_HASKELL__ >= 708
            , Typeable, Data
-#endif
            )
 
 instance (Binary a,   Binary   (v a)) => Binary   (Bootstrap v a) where
diff --git a/Statistics/Resampling/Bootstrap.hs b/Statistics/Resampling/Bootstrap.hs
--- a/Statistics/Resampling/Bootstrap.hs
+++ b/Statistics/Resampling/Bootstrap.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE CPP #-}
 -- |
 -- Module    : Statistics.Resampling.Bootstrap
 -- Copyright : (c) 2009, 2011 Bryan O'Sullivan
@@ -31,9 +30,7 @@
 
 import qualified Statistics.Resampling as R
 
-#if !defined(__GHCJS__)
-import Control.Monad.Par (parMap, runPar)
-#endif
+import Control.Parallel.Strategies (parMap, rdeepseq)
 
 data T = {-# UNPACK #-} !Double :< {-# UNPACK #-} !Double
 infixl 2 :<
@@ -51,15 +48,7 @@
   --   this.
   -> [Estimate ConfInt Double]
 bootstrapBCA confidenceLevel sample resampledData
-#if defined(__GHCJS__)
-  -- monad-par causes seems to cause "thread blocked indefinitely on MVar"
-  -- on GHCJS still
-  --
-  -- I (phadej) would change the interface to return IO, and use mapConcurrently from async
-  = map e resampledData
-#else
-  = runPar $ parMap e resampledData
-#endif
+  = parMap rdeepseq e resampledData
   where
     e (est, Bootstrap pt resample)
       | U.length sample == 1 || isInfinite bias =
diff --git a/Statistics/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE BangPatterns #-}
 -- |
 -- Module    : Statistics.Sample
 -- Copyright : (c) 2008 Don Stewart, 2009 Bryan O'Sullivan
@@ -20,6 +21,7 @@
     , range
 
     -- * Statistics of location
+    , expectation
     , mean
     , welfordMean
     , meanWeighted
@@ -51,20 +53,23 @@
     , fastVarianceUnbiased
     , fastStdDev
 
-    -- * Joint distirbutions
+    -- * Joint distributions
     , covariance
     , correlation
+    , covariance2
+    , correlation2
     , pair
     -- * References
     -- $references
     ) where
 
-import Statistics.Function (minMax)
+import Statistics.Function (minMax,square)
 import Statistics.Sample.Internal (robustSumVar, sum)
 import Statistics.Types.Internal  (Sample,WeightedSample)
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
+import Numeric.Sum (kbn, Summation(zero,add))
 
 -- Operator ^ will be overridden
 import Prelude hiding ((^), sum)
@@ -76,9 +81,17 @@
     where (lo , hi) = minMax s
 {-# INLINE range #-}
 
+-- | /O(n)/ Compute expectation of function over for sample. This is
+--   simply @mean . map f@ but won't create intermediate vector.
+expectation :: (G.Vector v a) => (a -> Double) -> v a -> Double
+expectation f xs = kbn (G.foldl' (\s -> add s . f) zero xs)
+                 / fromIntegral (G.length xs)
+{-# INLINE expectation #-}
+
 -- | /O(n)/ Arithmetic mean.  This uses Kahan-Babuška-Neumaier
 -- summation, so is more accurate than 'welfordMean' unless the input
--- values are very large.
+-- values are very large. This function is not subject to stream
+-- fusion.
 mean :: (G.Vector v Double) => v Double -> Double
 mean xs = sum xs / fromIntegral (G.length xs)
 {-# SPECIALIZE mean :: U.Vector Double -> Double #-}
@@ -122,7 +135,7 @@
 
 -- | /O(n)/ Geometric mean of a sample containing no negative values.
 geometricMean :: (G.Vector v Double) => v Double -> Double
-geometricMean = exp . mean . G.map log
+geometricMean = exp . expectation log
 {-# INLINE geometricMean #-}
 
 -- | Compute the /k/th central moment of a sample.  The central moment
@@ -138,7 +151,7 @@
     | a < 0  = error "Statistics.Sample.centralMoment: negative input"
     | a == 0 = 1
     | a == 1 = 0
-    | otherwise = sum (G.map go xs) / fromIntegral (G.length xs)
+    | otherwise = expectation go xs
   where
     go x = (x-m) ^ a
     m    = mean xs
@@ -215,7 +228,7 @@
 
 -- $variance
 --
--- The variance&#8212;and hence the standard deviation&#8212;of a
+-- The variance — and hence the standard deviation — of a
 -- sample of fewer than two elements are both defined to be zero.
 
 -- $robust
@@ -354,42 +367,79 @@
 
 -- | Covariance of sample of pairs. For empty sample it's set to
 --   zero
-covariance :: (G.Vector v (Double,Double), G.Vector v Double)
+covariance :: (G.Vector v (Double,Double))
            => v (Double,Double)
            -> Double
 covariance xy
   | n == 0    = 0
-  | otherwise = mean $ G.zipWith (*)
-                         (G.map (\x -> x - muX) xs)
-                         (G.map (\y -> y - muY) ys)
+  | otherwise = expectation (\(x,y) -> (x - muX)*(y - muY)) xy
   where
-    n       = G.length xy
-    (xs,ys) = G.unzip xy
-    muX     = mean xs
-    muY     = mean ys
+    n   = G.length xy
+    muX = expectation fst xy
+    muY = expectation snd xy
 {-# SPECIALIZE covariance :: U.Vector (Double,Double) -> Double #-}
 {-# SPECIALIZE covariance :: V.Vector (Double,Double) -> Double #-}
 
 -- | Correlation coefficient for sample of pairs. Also known as
 --   Pearson's correlation. For empty sample it's set to zero.
-correlation :: (G.Vector v (Double,Double), G.Vector v Double)
+correlation :: (G.Vector v (Double,Double))
            => v (Double,Double)
            -> Double
 correlation xy
   | n == 0    = 0
   | otherwise = cov / sqrt (varX * varY)
   where
-    n       = G.length xy
-    (xs,ys) = G.unzip xy
-    (muX,varX) = meanVariance xs
-    (muY,varY) = meanVariance ys
-    cov = mean $ G.zipWith (*)
-            (G.map (\x -> x - muX) xs)
-            (G.map (\y -> y - muY) ys)
+    n    = G.length xy
+    muX  = expectation (\(x,_) -> x) xy
+    muY  = expectation (\(_,y) -> y) xy
+    varX = expectation (\(x,_) -> square (x - muX))    xy
+    varY = expectation (\(_,y) -> square (y - muY))    xy
+    cov  = expectation (\(x,y) -> (x - muX)*(y - muY)) xy
 {-# SPECIALIZE correlation :: U.Vector (Double,Double) -> Double #-}
 {-# SPECIALIZE correlation :: V.Vector (Double,Double) -> Double #-}
 
 
+-- | Covariance of two samples. Both vectors must be of the same
+--   length. If both are empty it's set to zero
+covariance2 :: (G.Vector v Double)
+           => v Double
+           -> v Double
+           -> Double
+covariance2 xs ys
+  | nx /= ny  = error $ "Statistics.Sample.covariance2: both samples must have same length"
+  | nx == 0   = 0
+  | otherwise = sum (G.zipWith (\x y -> (x - muX)*(y - muY)) xs ys)
+              / fromIntegral nx
+  where
+    nx  = G.length xs
+    ny  = G.length ys
+    muX = mean xs
+    muY = mean ys
+{-# SPECIALIZE covariance2 :: U.Vector Double -> U.Vector Double -> Double #-}
+{-# SPECIALIZE covariance2 :: V.Vector Double -> V.Vector Double -> Double #-}
+
+-- | Correlation coefficient for two samples. Both vector must have
+--   same length Also known as Pearson's correlation. For empty sample
+--   it's set to zero.
+correlation2 :: (G.Vector v Double)
+             => v Double
+             -> v Double
+             -> Double
+correlation2 xs ys
+  | nx /= ny  = error $ "Statistics.Sample.correlation2: both samples must have same length"
+  | nx == 0   = 0
+  | otherwise = cov / sqrt (varX * varY)
+  where
+    nx         = G.length xs
+    ny         = G.length ys
+    (muX,varX) = meanVariance xs
+    (muY,varY) = meanVariance ys
+    cov = sum (G.zipWith (\x y -> (x - muX)*(y - muY)) xs ys)
+        / fromIntegral nx
+{-# SPECIALIZE correlation2 :: U.Vector Double -> U.Vector Double -> Double #-}
+{-# SPECIALIZE correlation2 :: V.Vector Double -> V.Vector Double -> Double #-}
+
+
 -- | Pair two samples. It's like 'G.zip' but requires that both
 --   samples have equal size.
 pair :: (G.Vector v a, G.Vector v b, G.Vector v (a,b)) => v a -> v b -> v (a,b)
@@ -403,8 +453,9 @@
 
 -- (^) operator from Prelude is just slow.
 (^) :: Double -> Int -> Double
-x ^ 1 = x
-x ^ n = x * (x ^ (n-1))
+x0 ^ n0 = go (n0-1) x0 where
+    go 0 !acc = acc
+    go n  acc = go (n-1) (acc*x0)
 {-# INLINE (^) #-}
 
 -- don't support polymorphism, as we can't get unboxed returns if we use it.
diff --git a/Statistics/Sample/Histogram.hs b/Statistics/Sample/Histogram.hs
--- a/Statistics/Sample/Histogram.hs
+++ b/Statistics/Sample/Histogram.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE FlexibleContexts, BangPatterns #-}
+{-# LANGUAGE FlexibleContexts, BangPatterns, ScopedTypeVariables #-}
 
 -- |
 -- Module    : Statistics.Sample.Histogram
@@ -19,6 +19,7 @@
     , range
     ) where
 
+import Control.Monad.ST
 import Numeric.MathFunctions.Constants (m_epsilon,m_tiny)
 import Statistics.Function (minMax)
 import qualified Data.Vector.Generic as G
@@ -49,7 +50,7 @@
 --
 -- Interval (bin) sizes are uniform, based on the supplied upper
 -- and lower bounds.
-histogram_ :: (Num b, RealFrac a, G.Vector v0 a, G.Vector v1 b) =>
+histogram_ :: forall b a v0 v1. (Num b, RealFrac a, G.Vector v0 a, G.Vector v1 b) =>
               Int
            -- ^ Number of bins.  This value must be positive.  A zero
            -- or negative value will cause an error.
@@ -65,6 +66,7 @@
            -> v1 b
 histogram_ numBins lo hi xs0 = G.create (GM.replicate numBins 0 >>= bin xs0)
   where
+    bin :: forall s. v0 a -> G.Mutable v1 s b -> ST s (G.Mutable v1 s b)
     bin xs bins = go 0
      where
        go i | i >= len = return bins
@@ -75,7 +77,7 @@
          go (i+1)
        write' bins' b !e = GM.write bins' b e
        len = G.length xs
-       d = ((hi - lo) * (1 + realToFrac m_epsilon)) / fromIntegral numBins
+       d = ((hi - lo) / fromIntegral numBins) * (1 + realToFrac m_epsilon)
 {-# INLINE histogram_ #-}
 
 -- | /O(n)/ Compute decent defaults for the lower and upper bounds of
diff --git a/Statistics/Sample/Internal.hs b/Statistics/Sample/Internal.hs
--- a/Statistics/Sample/Internal.hs
+++ b/Statistics/Sample/Internal.hs
@@ -14,9 +14,10 @@
     (
       robustSumVar
     , sum
+    , sumF
     ) where
 
-import Numeric.Sum (kbn, sumVector)
+import qualified Numeric.Sum as Sum
 import Prelude hiding (sum)
 import Statistics.Function (square)
 import qualified Data.Vector.Generic as G
@@ -26,5 +27,9 @@
 {-# INLINE robustSumVar #-}
 
 sum :: (G.Vector v Double) => v Double -> Double
-sum = sumVector kbn
+sum = Sum.sumVector Sum.kbn
 {-# INLINE sum #-}
+
+sumF :: Foldable f => f Double -> Double
+sumF = Sum.sum Sum.kbn
+{-# INLINE sumF #-}
diff --git a/Statistics/Test/Bartlett.hs b/Statistics/Test/Bartlett.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Test/Bartlett.hs
@@ -0,0 +1,99 @@
+{-# LANGUAGE CPP              #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-|
+Module      : Statistics.Test.Bartlett
+Description : Bartlett's test for homogeneity of variances.
+Copyright   : (c) Praneya Kumar, Alexey Khudyakov, 2025
+License     : BSD-3-Clause
+
+Bartlett's test is used to check that multiple groups of observations
+come from distributions with equal variances. This test assumes that
+samples come from normal distribution. If this is not the case it may
+simple test for non-normality and Levene's ("Statistics.Test.Levene")
+is preferred
+
+>>> import qualified Data.Vector.Unboxed as VU
+>>> import Statistics.Test.Bartlett
+>>> :{
+let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
+    c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
+in bartlettTest [a,b,c]
+:}
+Right (Test {testSignificance = mkPValue 1.1254782518843598e-5, testStatistics = 22.789434813726768, testDistribution = chiSquared 2})
+
+-}
+module Statistics.Test.Bartlett (
+    bartlettTest,
+    module Statistics.Distribution.ChiSquared
+) where
+
+import qualified Data.Vector           as V
+import qualified Data.Vector.Unboxed   as VU
+import qualified Data.Vector.Generic   as VG
+import qualified Data.Vector.Storable  as VS
+import qualified Data.Vector.Primitive as VP
+#if MIN_VERSION_vector(0,13,2)
+import qualified Data.Vector.Strict    as VV
+#endif
+
+import Statistics.Distribution (complCumulative)
+import Statistics.Distribution.ChiSquared (chiSquared, ChiSquared(..))
+import Statistics.Sample (varianceUnbiased)
+import Statistics.Types (mkPValue)
+import Statistics.Test.Types (Test(..))
+
+-- | Perform Bartlett's test for equal variances. The input is a list
+--   of vectors, where each vector represents a group of observations.
+bartlettTest :: VG.Vector v Double => [v Double] -> Either String (Test ChiSquared)
+bartlettTest groups
+  | length groups < 2                 = Left "At least two groups are required for Bartlett's test."
+  | any ((< 2) . VG.length) groups    = Left "Each group must have at least two observations."
+  | any ((<= 0) . var) groupVariances = Left "All groups must have positive variance."
+  | otherwise = Right Test
+      { testSignificance = pValue
+      , testStatistics   = tStatistic
+      , testDistribution = chiDist
+      }
+  where
+    -- Number of groups
+    k = length groups
+    -- Sample sizes for each group
+    ni  = map (fromIntegral . VG.length) groups
+    -- Total number of observations across all groups
+    n_tot = sum $ fromIntegral . VG.length <$> groups
+    -- Variance estimates
+    groupVariances = toVar <$> groups
+    sumWeightedVars = sum [ (n - 1) * v | Var{sampleN=n, var=v} <- groupVariances ]
+    pooledVariance  = sumWeightedVars / fromIntegral (n_tot - k)
+    -- Numerator of Bartlett's statistic
+    numerator =
+      fromIntegral (n_tot - k) * log pooledVariance -
+      sum [ (n - 1) * log v | Var{sampleN=n, var=v} <- groupVariances ]
+    -- Denominator correction term
+    sumReciprocals = sum [1 / (n - 1) | n <- ni]
+    denomCorrection =
+      1 + (sumReciprocals - 1 / fromIntegral (n_tot - k)) / (3 * (fromIntegral k - 1))
+
+    -- Test statistic and test distrubution
+    tStatistic = max 0 $ numerator / denomCorrection
+    chiDist    = chiSquared (k - 1)
+    pValue     = mkPValue $ complCumulative chiDist tStatistic
+{-# SPECIALIZE bartlettTest :: [V.Vector  Double] -> Either String (Test ChiSquared) #-}
+{-# SPECIALIZE bartlettTest :: [VU.Vector Double] -> Either String (Test ChiSquared) #-}
+{-# SPECIALIZE bartlettTest :: [VS.Vector Double] -> Either String (Test ChiSquared) #-}
+{-# SPECIALIZE bartlettTest :: [VP.Vector Double] -> Either String (Test ChiSquared) #-}
+#if MIN_VERSION_vector(0,13,2)
+{-# SPECIALIZE bartlettTest :: [VV.Vector Double] -> Either String (Test ChiSquared) #-}
+#endif
+
+-- Estimate of variance
+data Var = Var
+  { sampleN :: !Double -- ^ N of elements
+  , var     :: !Double -- ^ Sample variance
+  }
+
+toVar :: VG.Vector v Double => v Double -> Var
+toVar xs = Var { sampleN = fromIntegral $ VG.length xs
+               , var     = varianceUnbiased xs
+               }
diff --git a/Statistics/Test/ChiSquared.hs b/Statistics/Test/ChiSquared.hs
--- a/Statistics/Test/ChiSquared.hs
+++ b/Statistics/Test/ChiSquared.hs
@@ -17,8 +17,8 @@
 import qualified Data.Vector as V
 import qualified Data.Vector.Generic as G
 import qualified Data.Vector.Unboxed as U
-
-
+import qualified Data.Vector.Fusion.Bundle as F
+import qualified Numeric.Sum as Sum
 
 -- | Generic form of Pearson chi squared tests for binned data. Data
 --   sample is supplied in form of tuples (observed quantity,
@@ -26,7 +26,7 @@
 --
 --   This test should be used only if all bins have expected values of
 --   at least 5.
-chi2test :: (G.Vector v (Int,Double), G.Vector v Double)
+chi2test :: (G.Vector v (Int,Double))
          => Int                 -- ^ Number of additional degrees of
                                 --   freedom. One degree of freedom
                                 --   is due to the fact that the are
@@ -39,12 +39,15 @@
   | n   > 0   = Just Test
               { testSignificance = mkPValue $ complCumulative d chi2
               , testStatistics   = chi2
-              , testDistribution = chiSquared ndf
+              , testDistribution = chiSquared n
               }
   | otherwise = Nothing
   where
     n     = G.length vec - ndf - 1
-    chi2  = sum $ G.map (\(o,e) -> square (fromIntegral o - e) / e) vec
+    chi2  = Sum.kbn
+          $ F.foldl' Sum.add Sum.zero
+          $ F.map (\(o,e) -> square (fromIntegral o - e) / e)
+          $ G.stream vec
     d     = chiSquared n
 {-# INLINABLE  chi2test #-}
 {-# SPECIALIZE
@@ -56,7 +59,7 @@
 -- | Chi squared test for data with normal errors. Data is supplied in
 --   form of pair (observation with error, and expectation).
 chi2testCont
-  :: (G.Vector v (Estimate NormalErr Double, Double), G.Vector v Double)
+  :: (G.Vector v (Estimate NormalErr Double, Double))
   => Int                                   -- ^ Number of additional
                                            --   degrees of freedom.
   -> v (Estimate NormalErr Double, Double) -- ^ Observation and expectation.
@@ -66,10 +69,13 @@
   | n   > 0   = Just Test
               { testSignificance = mkPValue $ complCumulative d chi2
               , testStatistics   = chi2
-              , testDistribution = chiSquared ndf
+              , testDistribution = chiSquared n
               }
   | otherwise = Nothing
   where
     n     = G.length vec - ndf - 1
-    chi2  = sum $ G.map (\(Estimate o (NormalErr s),e) -> square (o - e) / s) vec
+    chi2  = Sum.kbn
+          $ F.foldl' Sum.add Sum.zero
+          $ F.map (\(Estimate o (NormalErr s),e) -> square (o - e) / s)
+          $ G.stream vec
     d     = chiSquared n
diff --git a/Statistics/Test/Internal.hs b/Statistics/Test/Internal.hs
--- a/Statistics/Test/Internal.hs
+++ b/Statistics/Test/Internal.hs
@@ -8,6 +8,7 @@
 import Data.Ord
 import           Data.Vector.Generic           ((!))
 import qualified Data.Vector.Generic         as G
+import qualified Data.Vector.Unboxed         as U
 import qualified Data.Vector.Generic.Mutable as M
 import Statistics.Function
 
@@ -27,15 +28,16 @@
 --   In case of ties average of ranks of equal elements is assigned
 --   to each
 --
--- >>> rank (==) (fromList [10,20,30::Int])
--- > fromList [1.0,2.0,3.0]
+-- >>> import qualified Data.Vector.Unboxed as VU
+-- >>> rank (==) (VU.fromList [10,20,30::Int])
+-- [1.0,2.0,3.0]
 --
--- >>> rank (==) (fromList [10,10,10,30::Int])
--- > fromList [2.0,2.0,2.0,4.0]
-rank :: (G.Vector v a, G.Vector v Double)
+-- >>> rank (==) (VU.fromList [10,10,10,30::Int])
+-- [2.0,2.0,2.0,4.0]
+rank :: (G.Vector v a)
      => (a -> a -> Bool)        -- ^ Equivalence relation
      -> v a                     -- ^ Vector to rank
-     -> v Double
+     -> U.Vector Double
 rank eq vec = G.unfoldr go (Rank 0 (-1) 1 vec)
   where
     go (Rank 0 _ r v)
@@ -58,11 +60,10 @@
 rankUnsorted :: ( Ord a
                 , G.Vector v a
                 , G.Vector v Int
-                , G.Vector v Double
                 , G.Vector v (Int, a)
                 )
              => v a
-             -> v Double
+             -> U.Vector Double
 rankUnsorted xs = G.create $ do
     -- Put ranks into their original positions
     -- NOTE: backpermute will do wrong thing
diff --git a/Statistics/Test/KolmogorovSmirnov.hs b/Statistics/Test/KolmogorovSmirnov.hs
--- a/Statistics/Test/KolmogorovSmirnov.hs
+++ b/Statistics/Test/KolmogorovSmirnov.hs
@@ -21,7 +21,7 @@
   , kolmogorovSmirnovCdfD
   , kolmogorovSmirnovD
   , kolmogorovSmirnov2D
-    -- * Probablities
+    -- * Probabilities
   , kolmogorovSmirnovProbability
     -- * References
     -- $references
diff --git a/Statistics/Test/KruskalWallis.hs b/Statistics/Test/KruskalWallis.hs
--- a/Statistics/Test/KruskalWallis.hs
+++ b/Statistics/Test/KruskalWallis.hs
@@ -17,7 +17,6 @@
   ) where
 
 import Data.Ord (comparing)
-import Data.Foldable (foldMap)
 import qualified Data.Vector.Unboxed as U
 import Statistics.Function (sort, sortBy, square)
 import Statistics.Distribution (complCumulative)
diff --git a/Statistics/Test/Levene.hs b/Statistics/Test/Levene.hs
new file mode 100644
--- /dev/null
+++ b/Statistics/Test/Levene.hs
@@ -0,0 +1,153 @@
+{-# LANGUAGE CPP              #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-|
+Module      : Statistics.Test.Levene
+Description : Levene's test for homogeneity of variances.
+Copyright   : (c) Praneya Kumar, Alexey Khudyakov, 2025
+License     : BSD-3-Clause
+
+Levene's test used to check whether samples have equal variance. Null
+hypothesis is all samples are from distributions with same variance
+(homoscedacity). Test is robust to non-normality, and versatile with
+mean or median centering.
+
+>>> import qualified Data.Vector.Unboxed as VU
+>>> import Statistics.Test.Levene
+>>> :{
+let a = VU.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = VU.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
+    c = VU.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
+in levenesTest Median [a, b, c]
+:}
+Right (Test {testSignificance = mkPValue 2.4315059672496814e-3, testStatistics = 7.584952754501659, testDistribution = fDistributionReal 2.0 27.0})
+-}
+module Statistics.Test.Levene (
+    Center(..),
+    levenesTest
+) where
+
+import Control.Monad
+import qualified Data.Vector           as V
+import qualified Data.Vector.Unboxed   as VU
+import qualified Data.Vector.Generic   as VG
+import qualified Data.Vector.Storable  as VS
+import qualified Data.Vector.Primitive as VP
+#if MIN_VERSION_vector(0,13,2)
+import qualified Data.Vector.Strict    as VV
+#endif
+import Statistics.Distribution (complCumulative)
+import Statistics.Distribution.FDistribution (fDistribution, FDistribution)
+import Statistics.Types      (mkPValue)
+import Statistics.Test.Types (Test(..))
+import Statistics.Function   (gsort)
+import Statistics.Sample     (mean)
+
+import qualified Statistics.Sample.Internal as IS
+import Statistics.Quantile
+
+
+-- | Center calculation method
+data Center
+  = Mean             -- ^ Use arithmetic mean
+  | Median           -- ^ Use median
+  | Trimmed !Double  -- ^ Trimmed mean with given proportion to cut from each end
+  deriving (Eq, Show)
+
+-- | Main Levene's test function with full error handling
+levenesTest
+  :: (VG.Vector v Double)
+  => Center      -- ^ Centering method
+  -> [v Double]  -- ^ Input samples
+  -> Either String (Test FDistribution)
+{-# INLINABLE levenesTest #-}
+levenesTest center samples
+  | length samples < 2 = Left "At least two samples required"
+  -- NOTE: We don't have nice way of computing mean of a list!
+  | otherwise = do
+      let residuals = computeResiduals center <$> samples
+      -- Average of all Z
+      let n_tot = sum $ VG.length . vecZ <$> residuals -- Total number of samples
+      let zbar = IS.sumF [ meanZ z * sampleN z
+                         | z <- residuals]
+               / fromIntegral n_tot
+      -- Numerator: Sum over (ni * (Z[i] - Z)^2)
+      let numerator = IS.sumF [ sampleN z * sqr (meanZ z - zbar)
+                              | z <- residuals]
+      -- Denominator: Sum over Σ((dev_ij - zbari)^2)
+      let denominator = IS.sumF
+            [ IS.sum $ VU.map (sqr . subtract (meanZ z)) (vecZ z)
+            | z <- residuals
+            ]
+      -- Handle division by zero and invalid values
+      when (denominator <= 0 || isNaN denominator || isInfinite denominator)
+        $ Left "Invalid denominator in W-statistic calculation"
+      let wStat = (fromIntegral (n_tot - k) / fromIntegral (k - 1)) * (numerator / denominator)
+          fDist = fDistribution (k - 1) (n_tot - k)
+      Right Test { testStatistics   = wStat
+                 , testSignificance = mkPValue $ complCumulative fDist wStat
+                 , testDistribution = fDist
+                 }
+  where
+    k = length samples -- Number of groups
+{-# SPECIALIZE levenesTest :: Center -> [V.Vector  Double] -> Either String (Test FDistribution) #-}
+{-# SPECIALIZE levenesTest :: Center -> [VU.Vector Double] -> Either String (Test FDistribution) #-}
+{-# SPECIALIZE levenesTest :: Center -> [VS.Vector Double] -> Either String (Test FDistribution) #-}
+{-# SPECIALIZE levenesTest :: Center -> [VP.Vector Double] -> Either String (Test FDistribution) #-}
+#if MIN_VERSION_vector(0,13,2)
+{-# SPECIALIZE levenesTest :: Center -> [VV.Vector Double] -> Either String (Test FDistribution) #-}
+#endif
+
+----------------------------------------------------------------
+-- Implementation
+----------------------------------------------------------------
+
+-- | Trim data from both ends with error handling and performance optimization
+trimboth :: (Ord a, Fractional a, VG.Vector v a)
+         => v a
+         -> Double
+         -> v a
+{-# INLINE trimboth #-}
+trimboth vec p
+  | p < 0 || p >= 0.5 = error "Statistics.Test.Levene: trimming: proportion must be between 0 and 0.5"
+  | VG.null vec       = vec
+  | otherwise         = VG.slice lowerCut (upperCut - lowerCut) sorted
+  where
+    n        = VG.length vec
+    sorted   = gsort vec
+    lowerCut = ceiling $ p * fromIntegral n
+    upperCut = n - lowerCut
+
+data Residuals = Residuals
+  { sampleN :: !Double
+  , meanZ   :: !Double
+  , vecZ    :: !(VU.Vector Double)
+  }
+
+computeResiduals
+  :: VG.Vector v Double
+  => Center
+  -> v Double
+  -> Residuals
+{-# INLINE computeResiduals #-}
+computeResiduals method xs = case method of
+  Mean   ->
+    let c  = mean xs
+        zs = VU.map (\x -> abs (x - c)) $ VU.convert xs
+    in makeR zs
+  Median ->
+    let c  = median medianUnbiased xs
+        zs = VU.map (\x -> abs (x - c)) $ VU.convert xs
+    in makeR zs
+  Trimmed p ->
+    let trimmed = trimboth xs p
+        c       = mean trimmed
+        zs      = VU.map (\x -> abs (x - c)) $ VU.convert trimmed
+    in makeR zs
+  where
+    makeR zs = Residuals { sampleN = fromIntegral $ VU.length zs
+                         , meanZ   = mean zs
+                         , vecZ    = zs
+                         }
+
+sqr :: Double -> Double
+sqr x = x * x
diff --git a/Statistics/Test/MannWhitneyU.hs b/Statistics/Test/MannWhitneyU.hs
--- a/Statistics/Test/MannWhitneyU.hs
+++ b/Statistics/Test/MannWhitneyU.hs
@@ -24,7 +24,6 @@
     -- $references
   ) where
 
-import Control.Applicative ((<$>))
 import Data.List (findIndex)
 import Data.Ord (comparing)
 import Numeric.SpecFunctions (choose)
diff --git a/Statistics/Test/StudentT.hs b/Statistics/Test/StudentT.hs
--- a/Statistics/Test/StudentT.hs
+++ b/Statistics/Test/StudentT.hs
@@ -71,7 +71,7 @@
 -- | Paired two-sample t-test. Two samples are paired in a
 -- within-subject design. Returns @Nothing@ if sample size is not
 -- sufficient.
-pairedTTest :: forall v. (G.Vector v (Double, Double), G.Vector v Double)
+pairedTTest :: forall v. (G.Vector v (Double, Double))
             => PositionTest          -- ^ one- or two-tailed test
             -> v (Double, Double)    -- ^ paired samples
             -> Maybe (Test StudentT)
@@ -134,16 +134,16 @@
 
 
 -- Calculate T-statistics for paired sample
-tStatisticsPaired :: (G.Vector v (Double, Double), G.Vector v Double)
+tStatisticsPaired :: (G.Vector v (Double, Double))
                   => v (Double, Double)
                   -> (Double, Double)
 {-# INLINE tStatisticsPaired #-}
 tStatisticsPaired sample = (t, ndf)
   where
     -- t-statistics
-    t = let d    = G.map (uncurry (-)) sample
-            sumd = G.sum d
-        in sumd / sqrt ((n * G.sum (G.map square d) - square sumd) / ndf)
+    t = let d    = U.map (uncurry (-)) $ G.convert sample
+            sumd = U.sum d
+        in sumd / sqrt ((n * U.sum (U.map square d) - square sumd) / ndf)
     -- degree of freedom
     ndf = n - 1
     n   = fromIntegral $ G.length sample
diff --git a/Statistics/Test/WilcoxonT.hs b/Statistics/Test/WilcoxonT.hs
--- a/Statistics/Test/WilcoxonT.hs
+++ b/Statistics/Test/WilcoxonT.hs
@@ -39,7 +39,6 @@
 -- the sum of negative ranks (the ranks of the differences where the second parameter is higher).
 -- to the length of the shorter sample.
 
-import Control.Applicative ((<$>))
 import Data.Function (on)
 import Data.List (findIndex)
 import Data.Ord (comparing)
diff --git a/Statistics/Types.hs b/Statistics/Types.hs
--- a/Statistics/Types.hs
+++ b/Statistics/Types.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE CPP #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE TypeFamilies #-}
@@ -72,12 +71,6 @@
 import Data.Vector.Unboxed          (Unbox)
 import Data.Vector.Unboxed.Deriving (derivingUnbox)
 import GHC.Generics                 (Generic)
-
-#if __GLASGOW_HASKELL__ == 704
-import qualified Data.Vector.Generic
-import qualified Data.Vector.Generic.Mutable
-#endif
-
 import Statistics.Internal
 import Statistics.Types.Internal
 import Statistics.Distribution
@@ -101,7 +94,7 @@
 -- second from @1 - CL@ or significance level.
 --
 -- >>> cl95
--- mkCLFromSignificance 0.05
+-- mkCLFromSignificance 5.0e-2
 --
 -- Prior to 0.14 confidence levels were passed to function as plain
 -- @Doubles@. Use 'mkCL' to convert them to @CL@.
@@ -141,7 +134,7 @@
 --   exception if parameter is out of [0,1] range
 --
 -- >>> mkCL 0.95    -- same as cl95
--- mkCLFromSignificance 0.05
+-- mkCLFromSignificance 5.0000000000000044e-2
 mkCL :: (Ord a, Num a) => a -> CL a
 mkCL
   = fromMaybe (error "Statistics.Types.mkCL: probability is out if [0,1] range")
@@ -151,7 +144,7 @@
 --   parameter is out of [0,1] range
 --
 -- >>> mkCLE 0.95    -- same as cl95
--- Just (mkCLFromSignificance 0.05)
+-- Just (mkCLFromSignificance 5.0000000000000044e-2)
 mkCLE :: (Ord a, Num a) => a -> Maybe (CL a)
 mkCLE p
   | p >= 0 && p <= 1 = Just $ CL (1 - p)
@@ -162,7 +155,7 @@
 --   throw exception if parameter is out of [0,1] range
 --
 -- >>> mkCLFromSignificance 0.05    -- same as cl95
--- mkCLFromSignificance 0.05
+-- mkCLFromSignificance 5.0e-2
 mkCLFromSignificance :: (Ord a, Num a) => a -> CL a
 mkCLFromSignificance = fromMaybe (error errMkCL) . mkCLFromSignificanceE
 
@@ -170,7 +163,7 @@
 --   parameter is out of [0,1] range
 --
 -- >>> mkCLFromSignificanceE 0.05    -- same as cl95
--- Just (mkCLFromSignificance 0.05)
+-- Just (mkCLFromSignificance 5.0e-2)
 mkCLFromSignificanceE :: (Ord a, Num a) => a -> Maybe (CL a)
 mkCLFromSignificanceE p
   | p >= 0 && p <= 1 = Just $ CL p
@@ -306,6 +299,7 @@
 -- >                       , confIntUDX = 6
 -- >                       , confIntCL  = cl95
 -- >                       }
+-- >          }
 --
 -- Prior to statistics 0.14 @Estimate@ data type used following definition:
 --
@@ -324,9 +318,7 @@
     , estError           :: !(e a)
       -- ^ Confidence interval for estimate.
     } deriving (Eq, Read, Show, Generic
-#if __GLASGOW_HASKELL__ >= 708
                , Typeable, Data
-#endif
                )
 
 instance (Binary   (e a), Binary   a) => Binary   (Estimate e a) where
diff --git a/bench-papi/Bench.hs b/bench-papi/Bench.hs
new file mode 100644
--- /dev/null
+++ b/bench-papi/Bench.hs
@@ -0,0 +1,14 @@
+-- |
+-- Here we reexport definitions of tasty-bench
+module Bench
+  ( whnf
+  , nf
+  , nfIO
+  , whnfIO
+  , bench
+  , bgroup
+  , defaultMain
+  , benchIngredients
+  ) where
+
+import Test.Tasty.PAPI
diff --git a/bench-time/Bench.hs b/bench-time/Bench.hs
new file mode 100644
--- /dev/null
+++ b/bench-time/Bench.hs
@@ -0,0 +1,14 @@
+-- |
+-- Here we reexport definitions of tasty-bench
+module Bench
+  ( whnf
+  , nf
+  , nfIO
+  , whnfIO
+  , bench
+  , bgroup
+  , defaultMain
+  , benchIngredients
+  ) where
+
+import Test.Tasty.Bench
diff --git a/benchmark/Main.hs b/benchmark/Main.hs
new file mode 100644
--- /dev/null
+++ b/benchmark/Main.hs
@@ -0,0 +1,77 @@
+module Main where
+
+import Data.Complex
+import Statistics.Sample
+import Statistics.Transform
+import Statistics.Correlation
+import System.Random.MWC
+import qualified Data.Vector.Unboxed as VU
+import qualified Data.Vector.Unboxed.Mutable as MVU
+
+import Bench
+
+
+-- Test sample
+sample :: VU.Vector Double
+sample = VU.create $ do g <- create
+                        MVU.replicateM 10000 (uniform g)
+
+-- Weighted test sample
+sampleW :: VU.Vector (Double,Double)
+sampleW = VU.zip sample (VU.reverse sample)
+
+-- Complex vector for FFT tests
+sampleC :: VU.Vector (Complex Double)
+sampleC = VU.zipWith (:+) sample (VU.reverse sample)
+
+
+-- Simple benchmark for functions from Statistics.Sample
+main :: IO ()
+main =
+  defaultMain
+  [ bgroup "sample"
+    [ bench "range"            $ nf (\x -> range x)            sample
+      -- Mean
+    , bench "mean"             $ nf (\x -> mean x)             sample
+    , bench "meanWeighted"     $ nf (\x -> meanWeighted x)     sampleW
+    , bench "harmonicMean"     $ nf (\x -> harmonicMean x)     sample
+    , bench "geometricMean"    $ nf (\x -> geometricMean x)    sample
+      -- Variance
+    , bench "variance"         $ nf (\x -> variance x)         sample
+    , bench "varianceUnbiased" $ nf (\x -> varianceUnbiased x) sample
+    , bench "varianceWeighted" $ nf (\x -> varianceWeighted x) sampleW
+      -- Correlation
+    , bench "pearson"          $ nf pearson     sampleW
+    , bench "covariance"       $ nf covariance  sampleW
+    , bench "correlation"      $ nf correlation sampleW
+    , bench "covariance2"      $ nf (covariance2  sample) sample
+    , bench "correlation2"     $ nf (correlation2 sample) sample
+      -- Other
+    , bench "stdDev"           $ nf (\x -> stdDev x)           sample
+    , bench "skewness"         $ nf (\x -> skewness x)         sample
+    , bench "kurtosis"         $ nf (\x -> kurtosis x)         sample
+      -- Central moments
+    , bench "C.M. 2"           $ nf (\x -> centralMoment 2 x)  sample
+    , bench "C.M. 3"           $ nf (\x -> centralMoment 3 x)  sample
+    , bench "C.M. 4"           $ nf (\x -> centralMoment 4 x)  sample
+    , bench "C.M. 5"           $ nf (\x -> centralMoment 5 x)  sample
+    ]
+  , bgroup "FFT"
+    [ bgroup "fft"
+      [ bench  (show n) $ whnf fft   (VU.take n sampleC) | n <- fftSizes ]
+    , bgroup "ifft"
+      [ bench  (show n) $ whnf ifft  (VU.take n sampleC) | n <- fftSizes ]
+    , bgroup "dct"
+      [ bench  (show n) $ whnf dct   (VU.take n sample)  | n <- fftSizes ]
+    , bgroup "dct_"
+      [ bench  (show n) $ whnf dct_  (VU.take n sampleC) | n <- fftSizes ]
+    , bgroup "idct"
+      [ bench  (show n) $ whnf idct  (VU.take n sample)  | n <- fftSizes ]
+    , bgroup "idct_"
+      [ bench  (show n) $ whnf idct_ (VU.take n sampleC) | n <- fftSizes ]
+    ]
+  ]
+
+
+fftSizes :: [Int]
+fftSizes = [32,128,512,2048]
diff --git a/benchmark/bench.hs b/benchmark/bench.hs
deleted file mode 100644
--- a/benchmark/bench.hs
+++ /dev/null
@@ -1,71 +0,0 @@
-import Control.Monad.ST (runST)
-import Criterion.Main
-import Data.Complex
-import Statistics.Sample
-import Statistics.Transform
-import Statistics.Correlation.Pearson
-import System.Random.MWC
-import qualified Data.Vector.Unboxed as U
-
-
--- Test sample
-sample :: U.Vector Double
-sample = runST $ flip uniformVector 10000 =<< create
-
--- Weighted test sample
-sampleW :: U.Vector (Double,Double)
-sampleW = U.zip sample (U.reverse sample)
-
--- Comlex vector for FFT tests
-sampleC :: U.Vector (Complex Double)
-sampleC = U.zipWith (:+) sample (U.reverse sample)
-
-
--- Simple benchmark for functions from Statistics.Sample
-main :: IO ()
-main =
-  defaultMain
-  [ bgroup "sample"
-    [ bench "range"            $ nf (\x -> range x)            sample
-      -- Mean
-    , bench "mean"             $ nf (\x -> mean x)             sample
-    , bench "meanWeighted"     $ nf (\x -> meanWeighted x)     sampleW
-    , bench "harmonicMean"     $ nf (\x -> harmonicMean x)     sample
-    , bench "geometricMean"    $ nf (\x -> geometricMean x)    sample
-      -- Variance
-    , bench "variance"         $ nf (\x -> variance x)         sample
-    , bench "varianceUnbiased" $ nf (\x -> varianceUnbiased x) sample
-    , bench "varianceWeighted" $ nf (\x -> varianceWeighted x) sampleW
-      -- Correlation
-    , bench "pearson"          $ nf (\x -> pearson (U.reverse sample) x) sample
-    , bench "pearson'"          $ nf (\x -> pearson' (U.reverse sample) x) sample
-    , bench "pearsonFast"      $ nf (\x -> pearsonFast (U.reverse sample) x) sample
-      -- Other
-    , bench "stdDev"           $ nf (\x -> stdDev x)           sample
-    , bench "skewness"         $ nf (\x -> skewness x)         sample
-    , bench "kurtosis"         $ nf (\x -> kurtosis x)         sample
-      -- Central moments
-    , bench "C.M. 2"           $ nf (\x -> centralMoment 2 x)  sample
-    , bench "C.M. 3"           $ nf (\x -> centralMoment 3 x)  sample
-    , bench "C.M. 4"           $ nf (\x -> centralMoment 4 x)  sample
-    , bench "C.M. 5"           $ nf (\x -> centralMoment 5 x)  sample
-    ]
-  , bgroup "FFT"
-    [ bgroup "fft"
-      [ bench  (show n) $ whnf fft   (U.take n sampleC) | n <- fftSizes ]
-    , bgroup "ifft"
-      [ bench  (show n) $ whnf ifft  (U.take n sampleC) | n <- fftSizes ]
-    , bgroup "dct"
-      [ bench  (show n) $ whnf dct   (U.take n sample)  | n <- fftSizes ]
-    , bgroup "dct_"
-      [ bench  (show n) $ whnf dct_  (U.take n sampleC) | n <- fftSizes ]
-    , bgroup "idct"
-      [ bench  (show n) $ whnf idct  (U.take n sample)  | n <- fftSizes ]
-    , bgroup "idct_"
-      [ bench  (show n) $ whnf idct_ (U.take n sampleC) | n <- fftSizes ]
-    ]
-  ]
-
-
-fftSizes :: [Int]
-fftSizes = [32,128,512,2048]
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,11 +1,97 @@
+## Changes in 0.16.5.0 [2026.01.09]
+
+ * `ContGen` and `DiscreteGen` instances for `Poisson` distributions are added.
+
+
+## Changes in 0.16.4.0 [2025.10.23]
+
+ * Bartlett's test (`Statistics.Test.Bartlett`) and Levene's test
+   (`Statistics.Test.Levene`) for homogeneity of variances is added.
+
+ * Improved performance in calculation of moments.
+
+ * Improved precision in calculation of `logDensity` of Student T distribution.
+
+
+## Changes in 0.16.3.0
+
+ * `S.Sample.correlation`, `S.Sample.covariance`,
+   `S.Correlation.pearson` do not allocate temporary arrays.
+
+ * Variants of correlation which take two vectors as input are added:
+   `S.Sample.correlation2`, `S.Sample.covariance2`, `S.Correlation.pearson2`,
+   `S.Correlation.spearman2`.
+
+ * Contexts for `S.Function.indexed`, `S.Correlation.spearman`, `S.pairedTTest`,
+   `S.Sample.correlation`, `S.Sample.covariance`, reduced.
+
+ * Computation of `rSquare` in linear regression has special case for case when
+   data variation is 0.
+
+ * Doctests added.
+
+ * Benchmarks using `tasty-bench` and `tasty-papi` added.
+
+ * Spurious test failures fixed.
+
+
+## Changes in 0.16.2.1
+
+ * Unnecessary constraint dropped from `tStatisticsPaired`.
+
+ * Compatibility with QuickCheck-2.14. Test suite doesn't fail every time.
+
+
+## Changes in 0.16.2.0
+
+ * Improved precision for `complCumulative` for hypergeometric and binomial
+   distributions. Precision improvements of geometric distribution
+
+ * Negative binomial distribution added.
+
+
+## Changes in 0.16.1.2
+
+ * Fixed bug in `fromSample` for exponential distribudion (#190)
+
+
+## Changes in 0.16.1.0
+
+ * Dependency on monad-par is dropped. `parMap` from `parallel` is used instead.
+
+
+## Changes in 0.16.0.2
+
+ * Bug in constructor of binomial distribution is fixed (#181). It accepted
+   out-of range probability before.
+
+
+## Changes in 0.16.0.0
+
+ * Random number generation switched to API introduced in random-1.2
+
+ * Support of GHC<7.10 is dropped
+
+ * Fix for chi-squared test (#167) which was completely wrong
+
+ * Computation of CDF and quantiles of Cauchy distribution is now numerically
+   stable.
+
+ * Fix loss of precision in computing of CDF of gamma distribution
+
+ * Log-normal and Weibull distributions added.
+
+ * `DiscreteGen` instance added for `DiscreteUniform`
+
+
 ## Changes in 0.15.2.0
 
  * Test suite is finally fixed (#42, #123). It took very-very-very long
    time but finally happened.
 
- * Avoid loss of precision when computing CDF for exponential districution.
+ * Avoid loss of precision when computing CDF for exponential distribution.
 
- * Avoid loss of precision when computing CDF for geometric districution. Add
+ * Avoid loss of precision when computing CDF for geometric distribution. Add
    complement of CDF.
 
  * Correctly handle case of n=0 in poissonCI
@@ -59,7 +145,7 @@
 ## Changes in 0.14.0.0
 
 Breaking update. It seriously changes parts of API. It adds new data types for
-dealing with with estimates, confidence intervals, confidence levels and
+dealing with estimates, confidence intervals, confidence levels and
 p-value. Also API for statistical tests is changed.
 
  * Module `Statistis.Types` now contains new data types for estimates,
@@ -231,19 +317,19 @@
 
   * Bugs in DCT and IDCT are fixed.
 
-  * Accesors for uniform distribution are added.
+  * Accessors for uniform distribution are added.
 
-  * ContGen instances for all continuous distribtuions are added.
+  * ContGen instances for all continuous distributions are added.
 
   * Beta distribution is added.
 
-  * Constructor for improper gamma distribtuion is added.
+  * Constructor for improper gamma distribution is added.
 
   * Binomial distribution allows zero trials.
 
   * Poisson distribution now accept zero parameter.
 
-  * Integer overflow in caculation of Wilcoxon-T test is fixed.
+  * Integer overflow in calculation of Wilcoxon-T test is fixed.
 
   * Bug in 'ContGen' instance for normal distribution is fixed.
 
@@ -284,7 +370,7 @@
   * Root finding is added, in S.Math.RootFinding.
 
   * The complCumulative function is added to the Distribution
-    class in order to accurately assess probalities P(X>x) which are
+    class in order to accurately assess probabilities P(X>x) which are
     used in one-tailed tests.
 
   * A stdDev function is added to the Variance class for
@@ -299,7 +385,7 @@
   * Bugs in quantile estimations for chi-square and gamma distribution
     are fixed.
 
-  * Integer overlow in mannWhitneyUCriticalValue is fixed. It
+  * Integer overflow in mannWhitneyUCriticalValue is fixed. It
     produced incorrect critical values for moderately large
     samples. Something around 20 for 32-bit machines and 40 for 64-bit
     ones.
@@ -321,18 +407,18 @@
 
   * Mean and variance for gamma distribution are fixed.
 
-  * Much faster cumulative probablity functions for Poisson and
+  * Much faster cumulative probability functions for Poisson and
     hypergeometric distributions.
 
   * Better density functions for gamma and Poisson distributions.
 
   * Student-T, Fisher-Snedecor F-distributions and Cauchy-Lorentz
-    distrbution are added.
+    distribution are added.
 
   * The function S.Function.create is removed. Use generateM from
     the vector package instead.
 
-  * Function to perform approximate comparion of doubles is added to
+  * Function to perform approximate comparison of doubles is added to
     S.Function.Comparison
 
   * Regularized incomplete beta function and its inverse are added to
diff --git a/statistics.cabal b/statistics.cabal
--- a/statistics.cabal
+++ b/statistics.cabal
@@ -1,5 +1,8 @@
+cabal-version:  3.0
+build-type:     Simple
+
 name:           statistics
-version:        0.15.2.0
+version:        0.16.5.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -22,43 +25,52 @@
   * Common statistical tests for significant differences between
     samples.
 
-license:        BSD2
+license:        BSD-2-Clause
 license-file:   LICENSE
-homepage:       https://github.com/bos/statistics
-bug-reports:    https://github.com/bos/statistics/issues
+homepage:       https://github.com/haskell/statistics
+bug-reports:    https://github.com/haskell/statistics/issues
 author:         Bryan O'Sullivan <bos@serpentine.com>, Alexey Khudaykov <alexey.skladnoy@gmail.com>
-maintainer:     Bryan O'Sullivan <bos@serpentine.com>, Alexey Khudaykov <alexey.skladnoy@gmail.com>
+maintainer:     Alexey Khudaykov <alexey.skladnoy@gmail.com>
 copyright:      2009-2014 Bryan O'Sullivan
 category:       Math, Statistics
-build-type:     Simple
-cabal-version:  >= 1.8
+
 extra-source-files:
   README.markdown
-  benchmark/bench.hs
-  changelog.md
   examples/kde/KDE.hs
   examples/kde/data/faithful.csv
   examples/kde/kde.html
   examples/kde/kde.tpl
-  tests/Tests/Math/Tables.hs
-  tests/Tests/Math/gen.py
   tests/utils/Makefile
   tests/utils/fftw.c
 
+extra-doc-files:
+  changelog.md
+
 tested-with:
-    GHC ==7.4.2
-     || ==7.6.3
-     || ==7.8.4
-     || ==7.10.3
-     || ==8.0.2
-     || ==8.2.2
-     || ==8.4.4
-     || ==8.6.5
-     || ==8.8.1
-  , GHCJS ==8.4
+  GHC ==8.4.4
+   || ==8.6.5
+   || ==8.8.4
+   || ==8.10.7
+   || ==9.0.2
+   || ==9.2.8
+   || ==9.4.8
+   || ==9.6.7
+   || ==9.8.4
+   || ==9.10.2
+   || ==9.12.2
 
+source-repository head
+  type:     git
+  location: https://github.com/haskell/statistics
 
+flag BenchPAPI
+  Description: Enable building of benchmarks which use instruction counters.
+               It requires libpapi and only works on Linux so it's protected by flag
+  Default: False
+  Manual:  True
+
 library
+  default-language: Haskell2010
   exposed-modules:
     Statistics.Autocorrelation
     Statistics.ConfidenceInt
@@ -76,11 +88,14 @@
     Statistics.Distribution.Geometric
     Statistics.Distribution.Hypergeometric
     Statistics.Distribution.Laplace
+    Statistics.Distribution.Lognormal
+    Statistics.Distribution.NegativeBinomial
     Statistics.Distribution.Normal
     Statistics.Distribution.Poisson
     Statistics.Distribution.StudentT
     Statistics.Distribution.Transform
     Statistics.Distribution.Uniform
+    Statistics.Distribution.Weibull
     Statistics.Function
     Statistics.Quantile
     Statistics.Regression
@@ -93,6 +108,8 @@
     Statistics.Sample.KernelDensity.Simple
     Statistics.Sample.Normalize
     Statistics.Sample.Powers
+    Statistics.Test.Bartlett
+    Statistics.Test.Levene
     Statistics.Test.ChiSquared
     Statistics.Test.KolmogorovSmirnov
     Statistics.Test.KruskalWallis
@@ -108,11 +125,11 @@
     Statistics.Internal
     Statistics.Test.Internal
     Statistics.Types.Internal
-  build-depends: base                    >= 4.5 && < 5
-               , base-orphans            >= 0.6 && <0.9
+  build-depends: base                    >= 4.9 && < 5
                  --
-               , math-functions          >= 0.3
-               , mwc-random              >= 0.13.0.0
+               , math-functions          >= 0.3.4.1
+               , mwc-random              >= 0.15.3.0
+               , random                  >= 1.2
                  --
                , aeson                   >= 0.6.0.0
                , async                   >= 2.2.2 && <2.3
@@ -120,21 +137,21 @@
                , binary                  >= 0.5.1.0
                , primitive               >= 0.3
                , dense-linear-algebra    >= 0.1 && <0.2
+               , parallel                >= 3.2.2.0 && <3.4
                , vector                  >= 0.10
                , vector-algorithms       >= 0.4
                , vector-th-unbox
                , vector-binary-instances >= 0.2.1
                , data-default-class      >= 0.1.2
-  if !impl(ghcjs)
-    build-depends:
-      monad-par               >= 0.3.4
+
   -- Older GHC
   if impl(ghc < 7.6)
     build-depends:
       ghc-prim
   ghc-options: -O2 -Wall -fwarn-tabs -funbox-strict-fields
 
-test-suite tests
+test-suite statistics-tests
+  default-language: Haskell2010
   type:           exitcode-stdio-1.0
   hs-source-dirs: tests
   main-is:        tests.hs
@@ -142,6 +159,7 @@
     Tests.ApproxEq
     Tests.Correlation
     Tests.Distribution
+    Tests.ExactDistribution
     Tests.Function
     Tests.Helpers
     Tests.KDE
@@ -156,6 +174,8 @@
     Tests.Quantile
   ghc-options:
     -Wall -threaded -rtsopts -fsimpl-tick-factor=500
+  if impl(ghc >= 9.8)
+    ghc-options: -Wno-x-partial
   build-depends: base
                , statistics
                , dense-linear-algebra
@@ -165,7 +185,6 @@
                , aeson
                , ieee754 >= 0.7.3
                , math-functions
-               , mwc-random
                , primitive
                , tasty
                , tasty-hunit
@@ -174,6 +193,49 @@
                , vector
                , vector-algorithms
 
-source-repository head
-  type:     git
-  location: https://github.com/bos/statistics
+test-suite statistics-doctests
+  default-language: Haskell2010
+  type:             exitcode-stdio-1.0
+  hs-source-dirs:   tests
+  main-is:          doctest.hs
+  if impl(ghcjs) || impl(ghc < 8.0)
+    Buildable: False
+  -- Linker on macos prints warnings to console which confuses doctests.
+  -- We simply disable doctests on ma for older GHC
+  -- > warning: -single_module is obsolete
+  if os(darwin) && impl(ghc < 9.6)
+    buildable: False
+  build-depends:
+            base       -any
+          , statistics -any
+          , doctest    >=0.15 && <0.25
+
+-- We want to be able to build benchmarks using both tasty-bench and tasty-papi.
+-- They have similar API so we just create two shim modules which reexport
+-- definitions from corresponding library and pick one in cabal file.
+common bench-stanza
+  ghc-options:      -Wall
+  default-language: Haskell2010
+  build-depends: base < 5
+               , vector          >= 0.12.3
+               , statistics
+               , mwc-random
+               , tasty           >=1.3.1
+
+benchmark statistics-bench
+  import:         bench-stanza
+  type:           exitcode-stdio-1.0
+  hs-source-dirs: benchmark bench-time
+  main-is:        Main.hs
+  Other-modules:  Bench
+  build-depends:  tasty-bench >= 0.3
+
+benchmark statistics-bench-papi
+  import:         bench-stanza
+  type:           exitcode-stdio-1.0
+  if impl(ghcjs) || !flag(BenchPAPI)
+     buildable: False
+  hs-source-dirs: benchmark bench-papi
+  main-is:        Main.hs
+  Other-modules:  Bench
+  build-depends:  tasty-papi >= 0.1.2
diff --git a/tests/Tests/Correlation.hs b/tests/Tests/Correlation.hs
--- a/tests/Tests/Correlation.hs
+++ b/tests/Tests/Correlation.hs
@@ -5,11 +5,11 @@
 
 import Control.Arrow (Arrow(..))
 import qualified Data.Vector as V
+import Data.Maybe
 import Statistics.Correlation
 import Statistics.Correlation.Kendall
-import Test.QuickCheck ((==>),Property,counterexample)
 import Test.Tasty
-import Test.Tasty.QuickCheck
+import Test.Tasty.QuickCheck hiding (sample)
 import Test.Tasty.HUnit
 
 import Tests.ApproxEq
@@ -34,15 +34,19 @@
 
 testPearson :: [(Double,Double)] -> Property
 testPearson sample
-  = (length sample > 1) ==> (exact ~= fast)
+  = (length sample > 1 && isJust exact) ==> (case exact of
+                                               Just e  -> e ~= fast
+                                               Nothing -> property False
+                                            )
   where
     (~=) = eql 1e-12
     exact = exactPearson $ map (realToFrac *** realToFrac) sample
     fast  = pearson $ V.fromList sample
 
-exactPearson :: [(Rational,Rational)] -> Double
+exactPearson :: [(Rational,Rational)] -> Maybe Double
 exactPearson sample
-  = realToFrac cov / sqrt (realToFrac (varX * varY))
+  | varX == 0 || varY == 0 = Nothing
+  | otherwise              = Just $ realToFrac cov / sqrt (realToFrac (varX * varY))
   where
     (xs,ys) = unzip sample
     n       = fromIntegral $ length sample
@@ -98,11 +102,11 @@
         , not (isNaN c3)
         , not (isNaN c4)
         ]
-  ==> ( counterexample (show sample0)
-      $ counterexample (show sample1)
-      $ counterexample (show sample2)
-      $ counterexample (show sample3)
-      $ counterexample (show sample4)
+  ==> ( counterexample ("S0 = " ++ show sample0)
+      $ counterexample ("S1 = " ++ show sample1)
+      $ counterexample ("S2 = " ++ show sample2)
+      $ counterexample ("S3 = " ++ show sample3)
+      $ counterexample ("S4 = " ++ show sample4)
       $ counterexample (show (c1,c2,c3,c4))
       $ and [ c1 == c2
             , c1 == c3
@@ -113,8 +117,8 @@
     -- We need to stretch sample into [-10 .. 10] range to avoid
     -- problems with under/overflows etc.
     stretch xs
-      | a == b = xs
-      | otherwise = [ (x - a - 10) * 20 / (a - b) | x <- xs ]
+      | a == b    = xs
+      | otherwise = [ ((x - a)/(b - a) - 0.5) * 20 | x <- xs ]
       where
         a = minimum xs
         b = maximum xs
diff --git a/tests/Tests/Distribution.hs b/tests/Tests/Distribution.hs
--- a/tests/Tests/Distribution.hs
+++ b/tests/Tests/Distribution.hs
@@ -2,10 +2,10 @@
     ViewPatterns #-}
 module Tests.Distribution (tests) where
 
-import Control.Applicative ((<$), (<$>), (<*>))
 import qualified Control.Exception as E
 import Data.List (find)
 import Data.Typeable (Typeable)
+import Data.Word
 import Numeric.MathFunctions.Constants (m_tiny,m_huge,m_epsilon)
 import Numeric.MathFunctions.Comparison
 import Statistics.Distribution
@@ -19,11 +19,14 @@
 import Statistics.Distribution.Geometric
 import Statistics.Distribution.Hypergeometric
 import Statistics.Distribution.Laplace        (LaplaceDistribution)
+import Statistics.Distribution.Lognormal      (LognormalDistribution)
+import Statistics.Distribution.NegativeBinomial (NegativeBinomialDistribution)
 import Statistics.Distribution.Normal         (NormalDistribution)
 import Statistics.Distribution.Poisson        (PoissonDistribution)
 import Statistics.Distribution.StudentT
 import Statistics.Distribution.Transform      (LinearTransform)
 import Statistics.Distribution.Uniform        (UniformDistribution)
+import Statistics.Distribution.Weibull        (WeibullDistribution)
 import Statistics.Distribution.DiscreteUniform (DiscreteUniform)
 import Test.Tasty                 (TestTree, testGroup)
 import Test.Tasty.QuickCheck      (testProperty)
@@ -33,6 +36,7 @@
 import Text.Printf (printf)
 
 import Tests.ApproxEq  (ApproxEq(..))
+import Tests.ExactDistribution (exactDistributionTests)
 import Tests.Helpers   (T(..), Double01(..), testAssertion, typeName)
 import Tests.Helpers   (monotonicallyIncreasesIEEE,isDenorm)
 import Tests.Orphanage ()
@@ -46,8 +50,10 @@
   , contDistrTests (T :: T ExponentialDistribution )
   , contDistrTests (T :: T GammaDistribution       )
   , contDistrTests (T :: T LaplaceDistribution     )
+  , contDistrTests (T :: T LognormalDistribution   )
   , contDistrTests (T :: T NormalDistribution      )
   , contDistrTests (T :: T UniformDistribution     )
+  , contDistrTests (T :: T WeibullDistribution     )
   , contDistrTests (T :: T StudentT                )
   , contDistrTests (T :: T (LinearTransform NormalDistribution))
   , contDistrTests (T :: T FDistribution           )
@@ -56,9 +62,11 @@
   , discreteDistrTests (T :: T GeometricDistribution      )
   , discreteDistrTests (T :: T GeometricDistribution0     )
   , discreteDistrTests (T :: T HypergeometricDistribution )
+  , discreteDistrTests (T :: T NegativeBinomialDistribution )
   , discreteDistrTests (T :: T PoissonDistribution        )
   , discreteDistrTests (T :: T DiscreteUniform            )
 
+  , exactDistributionTests
   , unitTests
   ]
 
@@ -71,8 +79,7 @@
 contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
   cdfTests t ++
   [ testProperty "PDF sanity"              $ pdfSanityCheck     t
-  ] ++
-  [ (if quantileIsInvCDF_enabled t then id else ignoreTest)
+  , (if quantileIsInvCDF_enabled t then id else ignoreTest)
   $ testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t
   , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t
   , testProperty "log density check"       $ logDensityCheck    t
@@ -86,7 +93,7 @@
   [ testProperty "Prob. sanity"         $ probSanityCheck       t
   , testProperty "CDF is sum of prob."  $ discreteCDFcorrect    t
   , testProperty "Discrete CDF is OK"   $ cdfDiscreteIsCorrect  t
-  , testProperty "log probabilty check" $ logProbabilityCheck   t
+  , testProperty "log probability check" $ logProbabilityCheck   t
   ]
 
 -- Tests for distributions which have CDF
@@ -94,7 +101,8 @@
 cdfTests t =
   [ testProperty "C.D.F. sanity"        $ cdfSanityCheck         t
   , testProperty "CDF limit at +inf"    $ cdfLimitAtPosInfinity  t
-  , testProperty "CDF limit at -inf"    $ cdfLimitAtNegInfinity  t
+  , (if cdfLimitAtNegInfinity_enabled t then id else ignoreTest)
+  $ testProperty "CDF limit at -inf"    $ cdfLimitAtNegInfinity  t
   , testProperty "CDF at +inf = 1"      $ cdfAtPosInfinity       t
   , testProperty "CDF at -inf = 1"      $ cdfAtNegInfinity       t
   , testProperty "CDF is nondecreasing" $ cdfIsNondecreasing     t
@@ -143,11 +151,14 @@
 
 
 -- CDF's complement is implemented correctly
-cdfComplementIsCorrect :: (Distribution d, Param d) => T d -> d -> Double -> Bool
+cdfComplementIsCorrect :: (Distribution d, Param d) => T d -> d -> Double -> Property
 cdfComplementIsCorrect _ d x
-  = 1 - (cumulative d x + complCumulative d x) <= tol
+  = counterexample ("err. tolerance = " ++ show tol)
+  $ counterexample ("difference     = " ++ show delta)
+  $ delta <= tol
   where
-    tol = prec_complementCDF d
+    tol   = prec_complementCDF d
+    delta = 1 - (cumulative d x + complCumulative d x)
 
 -- CDF for discrete distribution uses <= for comparison
 cdfDiscreteIsCorrect :: (Param d, DiscreteDistr d) => T d -> d -> Property
@@ -159,38 +170,43 @@
     --
     -- > CDF(i) - CDF(i-e) = P(i)
     --
-    -- Apporixmate equality is tricky here. Scale is set by maximum
-    -- value of CDF and probability. Case when all proabilities are
-    -- zero should be trated specially.
+    -- Approximate equality is tricky here. Scale is set by maximum
+    -- value of CDF and probability. Case when all probabilities are
+    -- zero should be treated specially.
     badN = [ printf "N=%3i    p[i]=%g\tp[i+1]=%g\tdP=%g\trelerr=%g" i p p1 dp ((p1-p-dp) / max p1 dp)
            | i <- [0 .. 100]
            , let p      = cumulative d $ fromIntegral i - 1e-6
                  p1     = cumulative d $ fromIntegral i
                  dp     = probability d i
                  relerr = ((p1 - p) - dp) / max p1 dp
-           ,  not (p == 0 && p1 == 0 && dp == 0)
-           && relerr > tol
+           , p  > m_tiny || p == 0
+           , p1 > m_tiny
+           , dp > m_tiny
+           , relerr > tol
            ]
     tol = prec_discreteCDF d
 
-logDensityCheck :: (ContDistr d) => T d -> d -> Double -> Property
+logDensityCheck :: (Param d, ContDistr d) => T d -> d -> Double -> Property
 logDensityCheck _ d x
   = not (isDenorm x)
   ==> ( counterexample (printf "density    = %g" p)
       $ counterexample (printf "logDensity = %g" logP)
       $ counterexample (printf "log p      = %g" (log p))
-      $ counterexample (printf "eps        = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
+      $ counterexample (printf "ulps[log]  = %i" ulpsLog)
+      $ counterexample (printf "ulps[lin]  = %i" ulpsLin)
       $ or [ p == 0      && logP == (-1/0)
            , p <= m_tiny && logP < log m_tiny
              -- To avoid problems with roundtripping error in case
              -- when density is computed as exponent of logDensity we
              -- accept either inequality
-           ,  (ulpDistance (log p) logP <= 32)
-           || (ulpDistance p (exp logP) <= 32)
+           ,  (ulpsLog <= n) || (ulpsLin <= n)
            ])
   where
-    p    = density d x
-    logP = logDensity d x
+    p       = density d x
+    logP    = logDensity d x
+    n       = prec_logDensity d
+    ulpsLog = ulpDistance (log p) logP
+    ulpsLin = ulpDistance p       (exp logP)
 
 -- PDF is positive
 pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool
@@ -201,6 +217,8 @@
 complQuantileCheck _ d (Double01 p)
   = counterexample (printf "x0 = %g" x0)
   $ counterexample (printf "x1 = %g" x1)
+  $ counterexample (printf "abs err = %g" $ abs (x1 - x0))
+  $ counterexample (printf "rel err = %g" $ relativeError x1 x0)
   -- We avoid extreme tails of distributions
   --
   -- FIXME: all parameters are arbitrary at the moment
@@ -208,7 +226,7 @@
         , p < 0.99
         , not $ isInfinite x0
         , not $ isInfinite x1
-        ] ==> (abs (x1 - x0) < 1e-6)
+        ] ==> (if x0 < 1e6 then abs (x1 - x0) < 1e-6 else relativeError x1 x0 < 1e-12)
   where
     x0 = quantile      d (1 - p)
     x1 = complQuantile d p
@@ -263,9 +281,9 @@
   $ counterexample (printf "Sum   = %g" p2)
   $ counterexample (printf "Delta = %g" (abs (p1 - p2)))
   $ abs (p1 - p2) < 3e-10
-  -- Avoid too large differeneces. Otherwise there is to much to sum
+  -- Avoid too large differences. Otherwise there is to much to sum
   --
-  -- Absolute difference is used guard againist precision loss when
+  -- Absolute difference is used guard against precision loss when
   -- close values of CDF are subtracted
   where
     n  = min a b
@@ -273,23 +291,26 @@
     p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)
     p2 = sum $ map (probability d) [n .. m]
 
-logProbabilityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Property
+logProbabilityCheck :: (Param d, DiscreteDistr d) => T d -> d -> Int -> Property
 logProbabilityCheck _ d x
   = counterexample (printf "probability    = %g" p)
   $ counterexample (printf "logProbability = %g" logP)
   $ counterexample (printf "log p          = %g" (log p))
-  $ counterexample (printf "eps            = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
+  $ counterexample (printf "ulps[log]      = %i" ulpsLog)
+  $ counterexample (printf "ulps[lin]      = %i" ulpsLin)
   $ or [ p == 0     && logP == (-1/0)
        , p < 1e-308 && logP < 609
          -- To avoid problems with roundtripping error in case
          -- when density is computed as exponent of logDensity we
          -- accept either inequality
-       ,  (ulpDistance (log p) logP <= 32)
-       || (ulpDistance p (exp logP) <= 32)
+       ,  (ulpsLog <= n) || (ulpsLin <= n)
        ]
   where
     p    = probability d x
     logP = logProbability d x
+    n    = prec_logDensity d
+    ulpsLog = ulpDistance (log p) logP
+    ulpsLin = ulpDistance p       (exp logP)
 
 
 -- | Parameters for distribution testing. Some distribution require
@@ -298,6 +319,9 @@
   -- | Whether quantileIsInvCDF is enabled
   quantileIsInvCDF_enabled :: T a -> Bool
   quantileIsInvCDF_enabled _ = True
+  -- | Whether cdfLimitAtNegInfinity is enabled
+  cdfLimitAtNegInfinity_enabled :: T a -> Bool
+  cdfLimitAtNegInfinity_enabled _ = True
   -- | Precision for 'quantileIsInvCDF' test
   prec_quantile_CDF :: a -> (Double,Double)
   prec_quantile_CDF _ = (16,16)
@@ -307,36 +331,57 @@
   -- | Precision of CDF's complement
   prec_complementCDF :: a -> Double
   prec_complementCDF _ = 1e-14
+  -- | Precision for logDensity check
+  prec_logDensity :: a -> Word64
+  prec_logDensity _ = 32
 
 instance Param StudentT where
   -- FIXME: disabled unless incompleteBeta troubles are sorted out
   quantileIsInvCDF_enabled _ = False
+
 instance Param BetaDistribution where
-  -- FIXME: See https://github.com/bos/statistics/issues/161 for details
+  -- FIXME: See https://github.com/haskell/statistics/issues/161 for details
   quantileIsInvCDF_enabled _ = False
+
 instance Param FDistribution where
   -- FIXME: disabled unless incompleteBeta troubles are sorted out
   quantileIsInvCDF_enabled _ = False
+  -- We compute CDF and complement using same method so precision
+  -- should be very good here.
+  prec_complementCDF _ = 64 * m_epsilon
 
 instance Param ChiSquared where
   prec_quantile_CDF _ = (32,32)
 
 instance Param BinomialDistribution where
-  prec_discreteCDF _ = 1e-13
-instance Param CauchyDistribution
+  prec_discreteCDF _ = 1e-12
+  prec_logDensity  _ = 48
+instance Param CauchyDistribution where
+  -- Distribution is long-tailed enough that we may never get to zero
+  cdfLimitAtNegInfinity_enabled _ = False
+
 instance Param DiscreteUniform
 instance Param ExponentialDistribution
-instance Param GammaDistribution
+instance Param GammaDistribution where
+  -- We lose precision near `incompleteGamma 10` because of error
+  -- introduced by exp . logGamma.  This could only be fixed in
+  -- math-function by implementing gamma
+  prec_quantile_CDF _ = (24,24)
+  prec_logDensity   _ = 512
 instance Param GeometricDistribution
 instance Param GeometricDistribution0
 instance Param HypergeometricDistribution
 instance Param LaplaceDistribution
+instance Param LognormalDistribution where
+  prec_quantile_CDF _ = (64,64)
+instance Param NegativeBinomialDistribution where
+  prec_discreteCDF  _ = 1e-12
+  prec_logDensity   _ = 48
 instance Param NormalDistribution
 instance Param PoissonDistribution
 instance Param UniformDistribution
+instance Param WeibullDistribution
 instance Param a => Param (LinearTransform a)
-
-
 
 ----------------------------------------------------------------
 -- Unit tests
diff --git a/tests/Tests/ExactDistribution.hs b/tests/Tests/ExactDistribution.hs
new file mode 100644
--- /dev/null
+++ b/tests/Tests/ExactDistribution.hs
@@ -0,0 +1,387 @@
+{-# LANGUAGE BangPatterns        #-}
+{-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE FlexibleInstances   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications    #-}
+{-# LANGUAGE TypeFamilies        #-}
+-- |
+-- Module    : Tests.ExactDistribution
+-- Copyright : (c) 2022 Lorenz Minder
+-- License   : BSD3
+--
+-- Maintainer  : lminder@gmx.net
+-- Stability   : experimental
+-- Portability : portable
+--
+-- Tests comparing distributions to exact versions.
+--
+-- This module provides exact versions of some distributions, and tests
+-- to compare them to the production implementations in
+-- Statistics.Distribution.*.  It also contains the functionality to
+-- test the production distributions against the exact versions.  Errors
+-- are flagged if data points are discovered where the probability mass
+-- function, the cumulative probability function, or its complement
+-- deviates too far (more than a prescribed tolerance) from the exact
+-- calculation.
+--
+-- The distributions here are implemented with rational integer
+-- arithmetic, using pretty much the textbook definitions formulas.
+-- Numerical problems like overflow or rounding errors cannot occur with
+-- this approach, making them are easy to write, read and verify.  They
+-- are, of course, substantially slower than the production
+-- distributions in Statistics.Distribution.*.  This makes them
+-- unsuitable for most uses other than testing and debugging.  (Also,
+-- only a handful of distributions can be implemented exactly with
+-- rational arithmetic.)
+--
+-- This module has the following sub-components:
+-- 
+-- * Exact (rational) definitions of some distribution functions,
+--   including both the probability mass as well as the CDF.
+--
+-- * QC.Arbitrary implementations to sample test cases (i.e.,
+--   distribution parameters and evaluation points).
+--
+-- * "Linkage": a mechanism to construct a production distribution
+--   corresponding to a test case for an exact distribution.
+--
+-- * A set of tests for the distributions derived using all of the above
+--   components.
+--
+-- This module exports a number symbols which can be useful for
+-- debugging and experimentation.  For use in a test suite, only the
+-- `exactDistributionTests` function is needed.
+
+module Tests.ExactDistribution (
+    -- * Exact math functions
+      exactChoose
+
+    -- * Exact distributions
+    , ExactDiscreteDistr(..)
+
+    , ExactBinomialDistr(..)
+    , ExactDiscreteUniformDistr(..)
+    , ExactGeometricDistr(..)
+    , ExactHypergeomDistr(..)
+
+    -- * Linking to production distributions
+    , ProductionLinkage
+
+    -- * Individual test routines
+    , pmfMatch
+    , cdfMatch
+    , complCdfMatch
+
+    -- * Test groups
+    , Tag(..)
+    , distTests
+    , exactDistributionTests
+) where
+
+----------------------------------------------------------------
+
+import Data.Foldable
+import Data.Ratio
+
+import Test.Tasty                       (TestTree, testGroup)
+import Test.Tasty.QuickCheck            (testProperty)
+import Test.QuickCheck as QC
+import Numeric.MathFunctions.Comparison (relativeError)
+import Numeric.MathFunctions.Constants  (m_tiny)
+
+import Statistics.Distribution
+import Statistics.Distribution.Binomial
+import Statistics.Distribution.DiscreteUniform
+import Statistics.Distribution.Geometric
+import Statistics.Distribution.Hypergeometric
+
+----------------------------------------------------------------
+--
+-- Math functions.
+--
+-- Used for implementing the distributions below.
+--
+----------------------------------------------------------------
+
+-- | Exactly compute binomial coefficient.
+--
+-- /n/ need not be an integer, can be fractional.
+exactChoose :: Ratio Integer -> Integer -> Ratio Integer
+exactChoose n k
+    | k < 0     = 0
+    | otherwise = foldl' (*) 1 factors
+    where   factors = [ (n - k' + j) / j | j <- [1..k'] ]
+            k' = fromInteger k :: Ratio Integer
+
+----------------------------------------------------------------
+--
+-- Exact distributions.
+--
+----------------------------------------------------------------
+
+-- | Exact discrete distribution.
+class ExactDiscreteDistr a where
+    -- | Probability mass function.
+    exactProb :: a -> Integer -> Ratio Integer
+    exactProb d x = exactCumulative d x - exactCumulative d (x - 1)
+
+    -- | Cumulative distribution function.
+    exactCumulative :: a -> Integer -> Ratio Integer
+
+-- | Exact Binomial distribution.
+data ExactBinomialDistr = ExactBD Integer (Ratio Integer)
+    deriving(Show)
+
+instance ExactDiscreteDistr ExactBinomialDistr where
+    -- Probability mass, computed with textbook formula.
+    exactProb (ExactBD n p) k
+        | k < 0 || k > n    = 0
+        | otherwise         = exactChoose n' k * p^k * (1-p)^(n-k)
+        where n' = fromIntegral n
+    -- CDF 
+    --
+    -- Computed iteratively by summing up all the probabilities
+    -- <= /k/.  Rather than computing everything from scratch for each
+    -- probability, we reuse previous results.  The meanings of the
+    -- variables in the "update" function are:
+    -- 
+    -- bc   is the binomial coefficient (n choose j),
+    -- pj   is the term p^j,
+    -- pnj  is the term (1 - p)^(n - j)
+    -- r    is the (partial) sum of the probabilities 
+    --
+    exactCumulative (ExactBD n p) k
+        | k < 0             = 0
+        | k >= n            = 1
+        -- Special case for p = 1, since in the below fold we
+        -- divide by (1 - p).
+        | p == 1            = if k == n then 1 else 0
+        | otherwise
+          = result $ foldl' update (1, 1, (1 - p)^n, (1 - p)^n) [1..k]
+          where update (!bc, !pj, !pnj, !r) !j =
+                    let bc' = bc * (n - j + 1) `div` j 
+                        pj' = pj * p
+                        pnj' = pnj / (1 - p)
+                        r' = r + (fromIntegral bc') * pj' * pnj'
+                    in  (bc', pj', pnj', r')
+                result (_, _, _, r) = r
+
+-- | Exact Discrete Uniform distribution.
+data ExactDiscreteUniformDistr = ExactDU Integer Integer
+    deriving(Show)
+
+instance ExactDiscreteDistr ExactDiscreteUniformDistr  where
+    exactProb (ExactDU lower upper) k
+        | k < lower || k > upper    = 0
+        | otherwise                 = 1 % (upper - lower + 1)
+    exactCumulative (ExactDU lower upper) k
+        | k < lower                 = 0
+        | k > upper                 = 1
+        | otherwise                 =
+            let d = (k - lower + 1)
+            in  d % (upper - lower + 1)
+
+-- | Geometric distribution.
+data ExactGeometricDistr = ExactGeom (Ratio Integer)
+    deriving(Show)
+
+instance ExactDiscreteDistr ExactGeometricDistr where
+    exactProb (ExactGeom p) k
+        | k < 1                     = 0
+        | otherwise                 = (1 - p)^(k - 1) * p
+
+    exactCumulative (ExactGeom p) k = 1 - (1 - p)^k
+
+-- | Hypergeometric distribution.
+--
+--   Parameters are /K/, /N/ and /n/, where:
+--   - /N/ is the total sample space size.
+--   - /K/ is number of "good" objects among /N/.
+--   - /n/ is the number of draws without replacement.
+data ExactHypergeomDistr = ExactHG Integer Integer Integer
+    deriving(Show)
+
+instance ExactDiscreteDistr ExactHypergeomDistr where
+    exactProb (ExactHG nK nN n) k
+        | k < 0                     = 0
+        | k > n || k > nN           = 0
+        | otherwise                 =
+            exactChoose nK' k * exactChoose (nN' - nK') (n - k)
+                / exactChoose nN' n
+            where nN' = fromIntegral nN
+                  nK' = fromIntegral nK
+
+    exactCumulative d k = sum [ exactProb d i | i <- [0..k] ]
+
+----------------------------------------------------------------
+--
+-- TestCase construction.
+--
+-- Contains the TestCase data type which encapsulates an instance of an
+-- exact distribution together with an evaluation point.
+--
+-- Then in contains the QC.Arbitrary implementations for TestCases of
+-- the different exact distributions.  As a general rule, we try the
+-- sampling to be relatively efficient, i.e., we only want to sample
+-- valid distribution parameters.  The evaluation points are sampled
+-- such that most points are within the support of the distribution.
+--
+----------------------------------------------------------------
+
+-- Divisor to compute a rational number from an integer.
+--
+-- We want input parameters to be exactly representable as
+-- Double values.  This is so that the production distribution does not
+-- mismatch the exact one simply because the input values don't exactly
+-- match.  (This can happen if the derivative of the distribution
+-- function is large.)   For this reason, the gd value needs to be a
+-- power of 2, and <= 2^53, since the mantissa of a Double is 53 bits.
+--
+-- A value of 2^53 gives the most accurate and diverse tests, but the
+-- cost is increased running times, as the computed numerators and
+-- denominators will become quite large.
+gd :: Integer
+gd = 2^(16 :: Int)
+
+-- TestCase
+--
+-- Combination of an exact distribution together with an evaluation point.
+data TestCase a = TestCase a Integer deriving (Show)
+
+instance QC.Arbitrary (TestCase ExactBinomialDistr) where
+    arbitrary = do
+        -- This somewhat odd sampling of /n/ is done so that lower
+        -- values (<1000) are more often represented as the larger ones.
+        n <- (*) <$> chooseInteger (1,1000) <*> chooseInteger(1,2)
+        p <- (% gd) <$> chooseInteger (0, gd)
+        k <- chooseInteger (-1, n + 1)
+        return $ TestCase (ExactBD n p) k
+    shrink _ = []
+
+instance QC.Arbitrary (TestCase ExactDiscreteUniformDistr) where
+    arbitrary = do
+        a <- chooseInteger (-1000, 1000)
+        sz <- chooseInteger (1, 1000)
+        let b = a + sz
+        k <- chooseInteger (a - 10, b + 10)
+        return $ TestCase (ExactDU a b) k
+    shrink _ = []
+
+instance QC.Arbitrary (TestCase ExactGeometricDistr) where
+    arbitrary = do
+        p <- (% gd) <$> chooseInteger (1, gd)
+        let lim = (floor $ 100 / p) :: Integer
+        k <- chooseInteger (0, lim)
+        return $ TestCase (ExactGeom p) k
+    shrink _ = []
+
+instance QC.Arbitrary (TestCase ExactHypergeomDistr) where
+    arbitrary = do
+        nN <- chooseInteger (1, 100)        -- XXX lower bound should be 0
+        nK <- chooseInteger (0, nN)
+        n  <- chooseInteger (1, nN)         -- XXX lower bound should be 0
+        k  <- chooseInteger (0, min n nK)
+        return $ TestCase (ExactHG nK nN n) k
+    shrink _ = []
+
+----------------------------------------------------------------
+--
+-- Linking to the production distributions
+--
+-- This section contains the ProductionLinkage typeclass and
+-- implementation, that allows to obtain a functions for evaluating
+-- the production distribution functions for a corresponding exact
+-- distribution.
+--
+----------------------------------------------------------------
+
+class (ExactDiscreteDistr a, DiscreteDistr (ProdDistrib a)
+      ) => ProductionLinkage a where
+  type ProdDistrib a
+  toProd :: a -> ProdDistrib a
+
+instance ProductionLinkage ExactBinomialDistr where
+  type ProdDistrib ExactBinomialDistr = BinomialDistribution
+  toProd (ExactBD n p) = binomial (fromIntegral n) (fromRational p)
+
+instance ProductionLinkage ExactDiscreteUniformDistr where
+  type ProdDistrib ExactDiscreteUniformDistr = DiscreteUniform
+  toProd (ExactDU lower upper) = discreteUniformAB (fromIntegral lower) (fromIntegral upper)
+
+instance ProductionLinkage ExactGeometricDistr where
+  type ProdDistrib ExactGeometricDistr = GeometricDistribution
+  toProd (ExactGeom p) = geometric $ fromRational p
+
+instance ProductionLinkage ExactHypergeomDistr where
+  type ProdDistrib ExactHypergeomDistr = HypergeometricDistribution
+  toProd (ExactHG nK nN n) =
+    hypergeometric (fromIntegral nK) (fromIntegral nN) (fromIntegral n)
+
+
+----------------------------------------------------------------
+-- Tests
+----------------------------------------------------------------
+
+-- Compare that probabilities agree. If they are denormalized just
+-- return True. You can't say much about precision
+probabilityAgree :: Double -> Double -> Double -> Bool
+probabilityAgree tol pe pa
+  | pa < 0      = False
+  | pe < 0      = False
+  | pe < m_tiny = True
+  | otherwise   = relativeError pe pa < tol
+
+-- Check production probability mass function accuracy.
+--
+-- Inputs: tolerance (max relative error) and test case
+pmfMatch :: (Show a, ProductionLinkage a) => Double -> TestCase a -> Property
+pmfMatch tol (TestCase dExact k)
+  = counterexample ("Exact  = " ++ show pe)
+  $ counterexample ("Approx = " ++ show pa)
+  $ probabilityAgree tol pe pa
+  where
+    pe = fromRational $ exactProb dExact k
+    pa = probability (toProd dExact) (fromIntegral k)
+
+-- Check production cumulative probability function accuracy.
+--
+-- Inputs:  tolerance (max relative error) and test case.
+cdfMatch :: (Show a, ProductionLinkage a) => Double -> TestCase a -> Bool
+cdfMatch tol (TestCase dExact k)
+  = probabilityAgree tol pe pa
+  where
+    pe = fromRational $ exactCumulative dExact k
+    pa = cumulative (toProd dExact) (fromIntegral k)
+
+-- Check production complement cumulative function accuracy.
+--
+-- Inputs:  tolerance (max relative error) and test case.
+complCdfMatch :: (Show a, ProductionLinkage a) => Double -> TestCase a -> Bool
+complCdfMatch tol (TestCase dExact k)
+  = probabilityAgree tol pe pa
+  where
+    pe = fromRational $ 1 - exactCumulative dExact k
+    pa = complCumulative (toProd dExact) (fromIntegral k)
+
+-- Phantom type to encode an exact distribution.
+data Tag a = Tag
+
+distTests :: forall a. (Show a, ProductionLinkage a, Arbitrary (TestCase a)) =>
+    Tag a -> String -> Double -> TestTree
+distTests (Tag :: Tag a) name tol =
+  testGroup ("Exact tests for " ++ name)
+    [ testProperty "PMF match"     $ pmfMatch      @a tol
+    , testProperty "CDF match"     $ cdfMatch      @a tol
+    , testProperty "1 - CDF match" $ complCdfMatch @a tol
+    ]
+
+
+-- Test driver -------------------------------------------------
+
+exactDistributionTests :: TestTree
+exactDistributionTests = testGroup "Test distributions against exact"
+  [ distTests (Tag @ExactBinomialDistr)        "Binomial"         1.0e-12
+  , distTests (Tag @ExactDiscreteUniformDistr) "DiscreteUniform"  1.0e-12
+  , distTests (Tag @ExactGeometricDistr)       "Geometric"        1.0e-13
+  , distTests (Tag @ExactHypergeomDistr)       "Hypergeometric"   1.0e-12
+  ]
diff --git a/tests/Tests/Helpers.hs b/tests/Tests/Helpers.hs
--- a/tests/Tests/Helpers.hs
+++ b/tests/Tests/Helpers.hs
@@ -60,8 +60,8 @@
 -- Check that function is nondecreasing taking rounding errors into
 -- account.
 --
--- In fact funstion is allowed to decrease less than one ulp in order
--- to guard againist problems with excess precision. On x86 FPU works
+-- In fact function is allowed to decrease less than one ulp in order
+-- to guard against problems with excess precision. On x86 FPU works
 -- with 80-bit numbers but doubles are 64-bit so rounding happens
 -- whenever values are moved from registers to memory
 monotonicallyIncreasesIEEE :: (Ord a, IEEE.IEEE b)  => (a -> b) -> a -> a -> Bool
diff --git a/tests/Tests/Math/Tables.hs b/tests/Tests/Math/Tables.hs
deleted file mode 100644
--- a/tests/Tests/Math/Tables.hs
+++ /dev/null
@@ -1,47 +0,0 @@
-module Tests.Math.Tables where
-
-tableLogGamma :: [(Double,Double)]
-tableLogGamma =
-  [(0.000001250000000, 13.592366285131769033)
-  , (0.000068200000000, 9.5930266308318756785)
-  , (0.000246000000000, 8.3100370767447966358)
-  , (0.000880000000000, 7.03508133735248542)
-  , (0.003120000000000, 5.768129358365567505)
-  , (0.026700000000000, 3.6082588918892977148)
-  , (0.077700000000000, 2.5148371858768232556)
-  , (0.234000000000000, 1.3579557559432759994)
-  , (0.860000000000000, 0.098146578027685615897)
-  , (1.340000000000000, -0.11404757557207759189)
-  , (1.890000000000000, -0.0425116422978701336)
-  , (2.450000000000000, 0.25014296569217625565)
-  , (3.650000000000000, 1.3701041997380685178)
-  , (4.560000000000000, 2.5375143317949580002)
-  , (6.660000000000000, 5.9515377269550207018)
-  , (8.250000000000000, 9.0331869196051233217)
-  , (11.300000000000001, 15.814180681373947834)
-  , (25.600000000000001, 56.711261598328121636)
-  , (50.399999999999999, 146.12815158702164808)
-  , (123.299999999999997, 468.85500075897556371)
-  , (487.399999999999977, 2526.9846647543727158)
-  , (853.399999999999977, 4903.9359135978220365)
-  , (2923.300000000000182, 20402.93198938705973)
-  , (8764.299999999999272, 70798.268343590112636)
-  , (12630.000000000000000, 106641.77264982508495)
-  , (34500.000000000000000, 325976.34838781820145)
-  , (82340.000000000000000, 849629.79603036714252)
-  , (234800.000000000000000, 2668846.4390507959761)
-  , (834300.000000000000000, 10540830.912557534873)
-  , (1230000.000000000000000, 16017699.322315014899)
-  ]
-tableIncompleteBeta :: [(Double,Double,Double,Double)]
-tableIncompleteBeta =
-  [(2.000000000000000, 3.000000000000000, 0.030000000000000, 0.0051864299999999996862)
-  , (2.000000000000000, 3.000000000000000, 0.230000000000000, 0.22845923000000001313)
-  , (2.000000000000000, 3.000000000000000, 0.760000000000000, 0.95465728000000005249)
-  , (4.000000000000000, 2.300000000000000, 0.890000000000000, 0.93829812158347802864)
-  , (1.000000000000000, 1.000000000000000, 0.550000000000000, 0.55000000000000004441)
-  , (0.300000000000000, 12.199999999999999, 0.110000000000000, 0.95063000053947077639)
-  , (13.100000000000000, 9.800000000000001, 0.120000000000000, 1.3483109941962659385e-07)
-  , (13.100000000000000, 9.800000000000001, 0.420000000000000, 0.071321857831804780226)
-  , (13.100000000000000, 9.800000000000001, 0.920000000000000, 0.99999578339197081611)
-  ]
diff --git a/tests/Tests/Math/gen.py b/tests/Tests/Math/gen.py
deleted file mode 100644
--- a/tests/Tests/Math/gen.py
+++ /dev/null
@@ -1,51 +0,0 @@
-#!/usr/bin/python
-"""
-"""
-
-from mpmath import *
-
-def printListLiteral(lines) :
-    print "  [" + "\n  , ".join(lines) + "\n  ]"
-
-################################################################
-# Generate header
-print "module Tests.Math.Tables where"
-print
-
-################################################################
-## Generate table for logGamma
-print "tableLogGamma :: [(Double,Double)]"
-print "tableLogGamma ="
-
-gammaArg = [ 1.25e-6, 6.82e-5, 2.46e-4, 8.8e-4,  3.12e-3, 2.67e-2,
-             7.77e-2, 0.234,   0.86,    1.34,    1.89,    2.45,
-             3.65,    4.56,    6.66,    8.25,    11.3,    25.6,
-             50.4,    123.3,   487.4,   853.4,   2923.3,  8764.3,
-             1.263e4, 3.45e4,  8.234e4, 2.348e5, 8.343e5, 1.23e6,
-             ]
-printListLiteral(
-    [ '(%.15f, %.20g)' % (x, log(gamma(x))) for x in gammaArg ]
-    )
-
-
-################################################################
-## Generate table for incompleteBeta
-
-print "tableIncompleteBeta :: [(Double,Double,Double,Double)]"
-print "tableIncompleteBeta ="
-
-incompleteBetaArg = [
-    (2,    3,    0.03),
-    (2,    3,    0.23),
-    (2,    3,    0.76),
-    (4,    2.3,  0.89),
-    (1,    1,    0.55),
-    (0.3,  12.2, 0.11),
-    (13.1, 9.8,  0.12),
-    (13.1, 9.8,  0.42),
-    (13.1, 9.8,  0.92),
-    ]
-printListLiteral(
-    [ '(%.15f, %.15f, %.15f, %.20g)' % (p,q,x, betainc(p,q,0,x, regularized=True))
-      for (p,q,x) in incompleteBetaArg
-      ])
diff --git a/tests/Tests/Matrix.hs b/tests/Tests/Matrix.hs
--- a/tests/Tests/Matrix.hs
+++ b/tests/Tests/Matrix.hs
@@ -5,7 +5,6 @@
 import Test.Tasty (TestTree, testGroup)
 import Test.Tasty.QuickCheck (testProperty)
 import Test.QuickCheck
-import Tests.ApproxEq (ApproxEq(..))
 import Tests.Matrix.Types
 import qualified Data.Vector.Unboxed as U
 
@@ -27,9 +26,20 @@
 t_transpose m = U.concat (map (column n) [0..rows m-1]) === toVector m
   where n = transpose m
 
-t_qr :: Matrix -> Property
-t_qr a = hasNaN p .||. eql 1e-10 a p
-  where p = uncurry multiply (qr a)
+t_qr :: Property
+t_qr = property $ do
+  a <- do (r,c) <- arbitrary
+          fromMat <$> arbMatWith r c (fromIntegral <$> choose (-10, 10::Int))
+  let (q,r) = qr a
+      a'    = multiply q r
+  pure $ counterexample ("A  = \n"++show a)
+       $ counterexample ("A' = \n"++show a')
+       $ counterexample ("Q  = \n"++show q)
+       $ counterexample ("R  = \n"++show r)
+       $ dimension a == dimension a'
+      && ( hasNaN a'
+        || and (zipWith (\x y -> abs (x - y) < 1e-12) (toList a) (toList a'))
+         )
 
 tests :: TestTree
 tests = testGroup "Matrix"
diff --git a/tests/Tests/Matrix/Types.hs b/tests/Tests/Matrix/Types.hs
--- a/tests/Tests/Matrix/Types.hs
+++ b/tests/Tests/Matrix/Types.hs
@@ -6,6 +6,8 @@
       Mat(..)
     , fromMat
     , toMat
+    , arbMat
+    , arbMatWith
     ) where
 
 import Control.Monad (join)
@@ -32,10 +34,21 @@
     arbitrary = small $ join (arbMat <$> arbitrary <*> arbitrary)
     shrink (Mat r c xs) = Mat r c <$> shrinkFixedList (shrinkFixedList shrink) xs
 
-arbMat :: (Arbitrary a) => Positive (Small Int) -> Positive (Small Int)
-       -> Gen (Mat a)
-arbMat (Positive (Small r)) (Positive (Small c)) =
-    Mat r c <$> vectorOf r (vector c)
+arbMat
+  :: (Arbitrary a)
+  => Positive (Small Int)
+  -> Positive (Small Int)
+  -> Gen (Mat a)
+arbMat r c = arbMatWith r c arbitrary
+
+arbMatWith
+  :: (Arbitrary a)
+  => Positive (Small Int)
+  -> Positive (Small Int)
+  -> Gen a
+  -> Gen (Mat a)
+arbMatWith (Positive (Small r)) (Positive (Small c)) genA =
+    Mat r c <$> vectorOf r (vectorOf c genA)
 
 instance Arbitrary Matrix where
     arbitrary = fromMat <$> arbitrary
diff --git a/tests/Tests/Orphanage.hs b/tests/Tests/Orphanage.hs
--- a/tests/Tests/Orphanage.hs
+++ b/tests/Tests/Orphanage.hs
@@ -16,11 +16,14 @@
 import Statistics.Distribution.Geometric
 import Statistics.Distribution.Hypergeometric
 import Statistics.Distribution.Laplace         (LaplaceDistribution, laplace)
+import Statistics.Distribution.Lognormal       (LognormalDistribution, lognormalDistr)
+import Statistics.Distribution.NegativeBinomial (NegativeBinomialDistribution, negativeBinomial)
 import Statistics.Distribution.Normal          (NormalDistribution, normalDistr)
 import Statistics.Distribution.Poisson         (PoissonDistribution, poisson)
 import Statistics.Distribution.StudentT
 import Statistics.Distribution.Transform       (LinearTransform, scaleAround)
 import Statistics.Distribution.Uniform         (UniformDistribution, uniformDistr)
+import Statistics.Distribution.Weibull         (WeibullDistribution, weibullDistr)
 import Statistics.Distribution.DiscreteUniform (DiscreteUniform, discreteUniformAB)
 import Statistics.Types
 
@@ -28,7 +31,7 @@
 
 
 ----------------------------------------------------------------
--- Arbitrary instances for ditributions
+-- Arbitrary instances for distributions
 ----------------------------------------------------------------
 
 instance QC.Arbitrary BinomialDistribution where
@@ -38,18 +41,23 @@
 instance QC.Arbitrary LaplaceDistribution where
   arbitrary = laplace <$> QC.choose (-10,10) <*> QC.choose (0, 2)
 instance QC.Arbitrary GammaDistribution where
-  arbitrary = gammaDistr <$> QC.choose (0.1,10) <*> QC.choose (0.1,10)
+  arbitrary = gammaDistr <$> QC.choose (0.1,100) <*> QC.choose (0.1,100)
 instance QC.Arbitrary BetaDistribution where
   arbitrary = betaDistr <$> QC.choose (1e-3,10) <*> QC.choose (1e-3,10)
 instance QC.Arbitrary GeometricDistribution where
-  arbitrary = geometric <$> QC.choose (0,1)
+  arbitrary = geometric <$> QC.choose (1e-10,1)
 instance QC.Arbitrary GeometricDistribution0 where
-  arbitrary = geometric0 <$> QC.choose (0,1)
+  arbitrary = geometric0 <$> QC.choose (1e-10,1)
 instance QC.Arbitrary HypergeometricDistribution where
   arbitrary = do l <- QC.choose (1,20)
                  m <- QC.choose (0,l)
                  k <- QC.choose (1,l)
                  return $ hypergeometric m l k
+instance QC.Arbitrary LognormalDistribution where
+  -- can't choose sigma too big, otherwise goes outside of double-float limit
+  arbitrary = lognormalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-10, 20)
+instance QC.Arbitrary NegativeBinomialDistribution where
+  arbitrary = negativeBinomial <$> QC.choose (1,100) <*> QC.choose (1e-10,1)
 instance QC.Arbitrary NormalDistribution where
   arbitrary = normalDistr <$> QC.choose (-100,100) <*> QC.choose (1e-3, 1e3)
 instance QC.Arbitrary PoissonDistribution where
@@ -60,6 +68,8 @@
   arbitrary = do a <- QC.arbitrary
                  b <- QC.arbitrary `suchThat` (/= a)
                  return $ uniformDistr a b
+instance QC.Arbitrary WeibullDistribution where
+  arbitrary = weibullDistr <$> QC.choose (1e-3,1e3) <*> QC.choose (1e-3, 1e3)
 instance QC.Arbitrary CauchyDistribution where
   arbitrary = cauchyDistribution
                 <$> arbitrary
diff --git a/tests/Tests/Parametric.hs b/tests/Tests/Parametric.hs
--- a/tests/Tests/Parametric.hs
+++ b/tests/Tests/Parametric.hs
@@ -4,12 +4,18 @@
 import Statistics.Test.StudentT
 import Statistics.Types
 import qualified Data.Vector.Unboxed as U
-import Test.Tasty (testGroup)
-import Tests.Helpers  (testEquality)
+import qualified Data.Vector as V
+import Test.Tasty (testGroup, TestTree)
+import Test.Tasty.HUnit (testCase, assertBool)
+import Tests.Helpers (testEquality)
 import qualified Test.Tasty as Tst
 
+import Statistics.Test.Levene
+import Statistics.Test.Bartlett
+
+
 tests :: Tst.TestTree
-tests = testGroup "Parametric tests" studentTTests
+tests = testGroup "Parametric tests" [studentTTests, bartlettTests, leveneTests]
 
 -- 2 samples x 20 obs data
 --
@@ -77,9 +83,9 @@
   , testEquality name (isSignificant (mkPValue $ pValue pVal + 1e-5) test)
     Significant
   ]
-  
-studentTTests :: [Tst.TestTree]
-studentTTests = concat
+
+studentTTests :: Tst.TestTree
+studentTTests = testGroup "StudentT test" $ concat
   [ -- R: t.test(sample1, sample2, alt="two.sided", var.equal=T)
     testTTest "two-sample t-test SamplesDiffer Student"
       (mkPValue 0.03410) (fromJust $ studentTTest SamplesDiffer sample1 sample2)
@@ -100,3 +106,119 @@
       (mkPValue 0.01705) (fromJust $ pairedTTest BGreater sample12)
   ]
   where sample12 = U.zip sample1 sample2
+
+
+------------------------------------------------------------
+-- Bartlett's Test
+------------------------------------------------------------
+
+bartlettTests :: TestTree
+bartlettTests = testGroup "Bartlett's test"
+  [ testCase "a,b,c" $ testBartlettTest [a,b,c] 1.8027132567760222   0.40601846976301237
+  , testCase "a,b"   $ testBartlettTest [a,b]   0.005221063776321886 0.9423974408021293
+  , testCase "a,c"   $ testBartlettTest [a,c]   1.1531619271845452   0.2828882244527482
+  , testCase "a,a"   $ testBartlettTest [a,a]   0.0                  1.0
+  ]
+  where
+    a = U.fromList [9.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = U.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 9.36, 9.18, 8.67, 9.05]
+    c = U.fromList [8.95, 8.12, 8.95, 8.85, 8.03, 8.84, 8.07, 8.98, 8.86, 8.98]
+
+testBartlettTest
+  :: [U.Vector Double]
+  -> Double
+  -> Double
+  -> IO ()
+testBartlettTest samples w p = do
+  r <- case bartlettTest samples of
+    Left  _ -> error "Bartlett's test failed"
+    Right r -> pure r
+  approxEqual "W" 1e-9 (testStatistics r)            w
+  approxEqual "p" 1e-9 (pValue $ testSignificance r) p
+
+------------------------------------------------------------
+-- Levene's Test (Trimmed Mean)
+------------------------------------------------------------
+
+leveneTests :: TestTree
+leveneTests = testGroup "Levene test"
+  -- Statistics' value and p-values are computed using 
+  [ testCase "a,b,c Mean"    $ testLeveneTest [a,b,c] Mean   7.905194483442054 0.001983795817472731
+  , testCase "a,b   Mean"    $ testLeveneTest [a,b]   Mean   8.83873787256358  0.008149720958328811
+  , testCase "a,a   Mean"    $ testLeveneTest [a,a]   Mean   0.0               1.0
+  , testCase "a,b,c Median"  $ testLeveneTest [a,b,c] Median 7.584952754501659 0.002431505967249681
+  , testCase "a,b   Median"  $ testLeveneTest [a,b]   Median 8.461374333228711 0.009364737715584399
+  , testCase "aL,bL Mean"    $ testLeveneTest [aL,bL] Mean   5.84424549939465  0.01653410652558999
+  , testCase "aL,bL Trimmed" $ testLeveneTest [aL,bL] (Trimmed 0.05) 8.368311226366314 0.004294953946529551
+  ]
+  where
+    a = V.fromList [8.88, 9.12, 9.04, 8.98, 9.00, 9.08, 9.01, 8.85, 9.06, 8.99]
+    b = V.fromList [8.88, 8.95, 9.29, 9.44, 9.15, 9.58, 8.36, 9.18, 8.67, 9.05]
+    c = V.fromList [8.95, 9.12, 8.95, 8.85, 9.03, 8.84, 9.07, 8.98, 8.86, 8.98]
+    -- Large samples for testing trimmed
+    aL = V.fromList [
+      -0.18919252, -1.62837673,  5.21332355, -0.00962043, -0.28417847,
+      -0.88128233,  1.49698436,  6.1780359 , -1.22301348,  3.34598245,
+       5.33227264, -0.88732069,  0.14487346,  2.61060215,  4.22033907,
+       2.53139215, -0.72131061,  0.53063607, -0.60510374, -0.73230842,
+       1.54037043, -2.81103963,  3.40763063,  0.49005324,  2.13085513,
+       5.68650547,  4.16397279, -0.17325097,  1.12664972,  4.23297516,
+       4.15943436, -1.01452078,  2.40391646,  0.83019962,  0.29665879,
+      -3.83031046, -1.98576933,  1.5356527 ,  1.30773365,  0.292818  ,
+       2.45877828,  1.06482289, -0.63241873,  1.58465379,  1.96577614,
+       2.25791943,  4.13769848, -2.38595767, -0.65801423, -2.54007791,
+       3.17428087,  4.32096964,  0.92240335, -2.38101319,  1.35692587,
+       1.48279101, -0.04438309,  0.50296642,  2.08261495,  1.33181215,
+      -1.95427198,  4.95406809,  1.51294898, -2.68536129, -0.2441218 ,
+       2.41142613,  4.71051493,  2.66618697,  1.12668301, -0.25732583,
+       1.25021838, -1.27523641,  5.01638744,  3.38864442,  0.17979744,
+      -0.88481645,  3.89346357, -0.51512217, -1.60542888,  0.88378679,
+      -2.12962732, -1.35989539,  5.09215112, -1.37442481,  0.83578405,
+       0.13829571,  1.25171481,  3.60552158, -3.24051591, -0.44301834,
+       0.78253445,  1.76098254,  1.79677434, -0.19010505,  3.07640466,
+       3.02853882,  1.24849063,  4.84505382,  6.82274999,  2.24063474]
+    bL = V.fromList [
+        2.15584101, -2.74876744, -0.82231894,  1.97518087,  2.59280595,
+        1.28703417,  2.40450278,  1.9761031 ,  2.35186598,  1.15611047,
+        2.26709318,  1.2832138 , -2.1486074 ,  0.27563011, -0.51816861,
+        0.89658424,  3.27069545,  1.72846646,  3.84454277,  5.58301459,
+       -0.40878188,  3.41602853,  1.1281526 ,  0.9665913 ,  0.76567084,
+        1.69522855,  1.69133014,  0.70529264,  2.65243202, -1.0088019 ,
+       -0.62431026,  3.76667396,  3.66225181,  0.73217579,  0.04478736,
+        0.4169833 ,  0.77065631, -1.31484093,  1.23858618, -0.08339456,
+        3.14154286,  1.84358218, -0.53511423, -3.4919477 ,  0.24076997,
+        3.59381684,  1.99497806,  2.95499775,  1.67157731,  0.0214764 ,
+        3.32161612, -2.64762427,  0.06486472,  0.19653897,  1.34954235,
+        1.18568747, -0.54434597, -3.35544223,  1.41933109,  0.95100195,
+        2.7182116 ,  1.1334068 , -0.95297806, -0.05421818,  1.42248799,
+       -3.96201277, -3.21309254, -0.21209211,  0.9689551 ,  0.13526401,
+       -0.88656198,  0.41331783, -3.18766064,  4.34948246,  1.35656384,
+        0.41920101, -0.46578994,  1.55181583,  2.43937014,  2.49040644,
+        4.10505494,  1.68856296,  1.31503895,  0.41123368,  0.73242999,
+        0.2804349 , -1.83494592, -0.31073195,  2.61185513,  2.91645094,
+        1.26097638,  2.64197134,  3.88931972,  0.03783002,  2.55209729,
+        3.46869549,  0.96348003,  2.27658242,  2.7613171 , -0.1372434 ]
+
+    
+testLeveneTest
+  :: [V.Vector Double]
+  -> Center
+  -> Double
+  -> Double
+  -> IO ()
+testLeveneTest samples center w p = do
+  r <- case levenesTest center samples of
+    Left  _ -> error "Levene's test failed"
+    Right r -> pure r
+  approxEqual "W" 1e-9 (testStatistics r)            w
+  approxEqual "p" 1e-9 (pValue $ testSignificance r) p
+
+
+----------------------------------------------------------------
+
+approxEqual :: String -> Double -> Double -> Double -> IO ()
+approxEqual name epsilon actual expected =
+  assertBool (name ++ ": expected ≈ " ++ show expected ++ ", got " ++ show actual)
+             (diff < epsilon)
+  where
+    diff = abs (actual - expected)
diff --git a/tests/Tests/Serialization.hs b/tests/Tests/Serialization.hs
--- a/tests/Tests/Serialization.hs
+++ b/tests/Tests/Serialization.hs
@@ -16,11 +16,14 @@
 import Statistics.Distribution.Geometric
 import Statistics.Distribution.Hypergeometric
 import Statistics.Distribution.Laplace        (LaplaceDistribution)
+import Statistics.Distribution.Lognormal      (LognormalDistribution)
+import Statistics.Distribution.NegativeBinomial (NegativeBinomialDistribution)
 import Statistics.Distribution.Normal         (NormalDistribution)
 import Statistics.Distribution.Poisson        (PoissonDistribution)
 import Statistics.Distribution.StudentT
 import Statistics.Distribution.Transform      (LinearTransform)
 import Statistics.Distribution.Uniform        (UniformDistribution)
+import Statistics.Distribution.Weibull        (WeibullDistribution)
 import Statistics.Types
 
 import Test.Tasty            (TestTree, testGroup)
@@ -50,8 +53,11 @@
   , serializationTests (T :: T ExponentialDistribution )
   , serializationTests (T :: T GammaDistribution       )
   , serializationTests (T :: T LaplaceDistribution     )
+  , serializationTests (T :: T LognormalDistribution   )
+  , serializationTests (T :: T NegativeBinomialDistribution         )
   , serializationTests (T :: T NormalDistribution      )
   , serializationTests (T :: T UniformDistribution     )
+  , serializationTests (T :: T WeibullDistribution     )
   , serializationTests (T :: T StudentT                )
   , serializationTests (T :: T (LinearTransform NormalDistribution))
   , serializationTests (T :: T FDistribution           )
diff --git a/tests/Tests/Transform.hs b/tests/Tests/Transform.hs
--- a/tests/Tests/Transform.hs
+++ b/tests/Tests/Transform.hs
@@ -8,7 +8,6 @@
 
 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)
diff --git a/tests/doctest.hs b/tests/doctest.hs
new file mode 100644
--- /dev/null
+++ b/tests/doctest.hs
@@ -0,0 +1,5 @@
+import Test.DocTest (doctest)
+
+main :: IO ()
+main = doctest ["-XHaskell2010", "Statistics"]
+