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,11 +27,18 @@
 
 -- | 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 #-}
 
+-- | 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
@@ -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)
             )
@@ -63,6 +70,27 @@
   where
     (x, y) = G.unzip xy
 {-# INLINE spearman #-}
+
+-- | 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
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/Regression.hs b/Statistics/Regression.hs
--- a/Statistics/Regression.hs
+++ b/Statistics/Regression.hs
@@ -41,8 +41,15 @@
 --   element than the list of predictors; the last element is the
 --   /y/-intercept value.
 --
--- * /R²/, 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
@@ -65,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
@@ -92,7 +122,7 @@
   where n = rows r
         l = U.length b
 
--- | Compute /R²/, 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
@@ -101,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/Sample.hs b/Statistics/Sample.hs
--- a/Statistics/Sample.hs
+++ b/Statistics/Sample.hs
@@ -20,6 +20,7 @@
     , range
 
     -- * Statistics of location
+    , expectation
     , mean
     , welfordMean
     , meanWeighted
@@ -54,17 +55,20 @@
     -- * 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 +80,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 +134,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 +150,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
@@ -354,40 +366,77 @@
 
 -- | 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
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/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)
diff --git a/Statistics/Types.hs b/Statistics/Types.hs
--- a/Statistics/Types.hs
+++ b/Statistics/Types.hs
@@ -94,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@.
@@ -134,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")
@@ -144,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)
@@ -155,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
 
@@ -163,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
@@ -299,6 +299,7 @@
 -- >                       , confIntUDX = 6
 -- >                       , confIntCL  = cl95
 -- >                       }
+-- >          }
 --
 -- Prior to statistics 0.14 @Estimate@ data type used following definition:
 --
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/bench.hs b/benchmark/bench.hs
--- a/benchmark/bench.hs
+++ b/benchmark/bench.hs
@@ -1,24 +1,26 @@
-import Control.Monad.ST (runST)
-import Criterion.Main
 import Data.Complex
 import Statistics.Sample
 import Statistics.Transform
-import Statistics.Correlation.Pearson
+import Statistics.Correlation
 import System.Random.MWC
-import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed as VU
+import qualified Data.Vector.Unboxed.Mutable as MVU
 
+import Bench
 
+
 -- Test sample
-sample :: U.Vector Double
-sample = runST $ flip uniformVector 10000 =<< create
+sample :: VU.Vector Double
+sample = VU.create $ do g <- create
+                        MVU.replicateM 10000 (uniform g)
 
 -- Weighted test sample
-sampleW :: U.Vector (Double,Double)
-sampleW = U.zip sample (U.reverse sample)
+sampleW :: VU.Vector (Double,Double)
+sampleW = VU.zip sample (VU.reverse sample)
 
 -- Complex vector for FFT tests
-sampleC :: U.Vector (Complex Double)
-sampleC = U.zipWith (:+) sample (U.reverse sample)
+sampleC :: VU.Vector (Complex Double)
+sampleC = VU.zipWith (:+) sample (VU.reverse sample)
 
 
 -- Simple benchmark for functions from Statistics.Sample
@@ -37,9 +39,11 @@
     , 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
+    , 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
@@ -52,17 +56,17 @@
     ]
   , bgroup "FFT"
     [ bgroup "fft"
-      [ bench  (show n) $ whnf fft   (U.take n sampleC) | n <- fftSizes ]
+      [ bench  (show n) $ whnf fft   (VU.take n sampleC) | n <- fftSizes ]
     , bgroup "ifft"
-      [ bench  (show n) $ whnf ifft  (U.take n sampleC) | n <- fftSizes ]
+      [ bench  (show n) $ whnf ifft  (VU.take n sampleC) | n <- fftSizes ]
     , bgroup "dct"
-      [ bench  (show n) $ whnf dct   (U.take n sample)  | n <- fftSizes ]
+      [ bench  (show n) $ whnf dct   (VU.take n sample)  | n <- fftSizes ]
     , bgroup "dct_"
-      [ bench  (show n) $ whnf dct_  (U.take n sampleC) | n <- fftSizes ]
+      [ bench  (show n) $ whnf dct_  (VU.take n sampleC) | n <- fftSizes ]
     , bgroup "idct"
-      [ bench  (show n) $ whnf idct  (U.take n sample)  | n <- fftSizes ]
+      [ bench  (show n) $ whnf idct  (VU.take n sample)  | n <- fftSizes ]
     , bgroup "idct_"
-      [ bench  (show n) $ whnf idct_ (U.take n sampleC) | n <- fftSizes ]
+      [ bench  (show n) $ whnf idct_ (VU.take n sampleC) | n <- fftSizes ]
     ]
   ]
 
diff --git a/changelog.md b/changelog.md
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,25 @@
+## 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`.
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.16.2.1
+version:        0.16.3.0
 synopsis:       A library of statistical types, data, and functions
 description:
   This library provides a number of common functions and types useful
@@ -22,7 +25,7 @@
   * Common statistical tests for significant differences between
     samples.
 
-license:        BSD2
+license:        BSD-2-Clause
 license-file:   LICENSE
 homepage:       https://github.com/haskell/statistics
 bug-reports:    https://github.com/haskell/statistics/issues
@@ -30,32 +33,39 @@
 maintainer:     Alexey Khudaykov <alexey.skladnoy@gmail.com>
 copyright:      2009-2014 Bryan O'Sullivan
 category:       Math, Statistics
-build-type:     Simple
-cabal-version:  >= 1.10
+
 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
 
 tested-with:
-    GHC ==8.4.4
-    GHC ==8.6.5
-    GHC ==8.8.4
-    GHC ==8.10.7
-    GHC ==9.0.2
-    GHC ==9.2.8
-    GHC ==9.4.6
-    GHC ==9.6.2
+  GHC ==8.4.4
+   || ==8.6.5
+   || ==8.8.4
+   || ==8.10.7
+   || ==9.0.2
+   || ==9.2.8
+   || ==9.4.8
+   || ==9.6.6
+   || ==9.8.4
+   || ==9.10.1
 
+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:
@@ -176,6 +186,49 @@
                , vector
                , vector-algorithms
 
-source-repository head
-  type:     git
-  location: https://github.com/haskell/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.24
+
+-- 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:        bench.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:        bench.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
@@ -102,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
@@ -117,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/ExactDistribution.hs b/tests/Tests/ExactDistribution.hs
--- a/tests/Tests/ExactDistribution.hs
+++ b/tests/Tests/ExactDistribution.hs
@@ -1,8 +1,9 @@
-{-# LANGUAGE BangPatterns,
-             FlexibleInstances,
-             FlexibleContexts,
-             ScopedTypeVariables
-  #-}
+{-# LANGUAGE BangPatterns        #-}
+{-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE FlexibleInstances   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications    #-}
+{-# LANGUAGE TypeFamilies        #-}
 -- |
 -- Module    : Tests.ExactDistribution
 -- Copyright : (c) 2022 Lorenz Minder
@@ -64,8 +65,6 @@
     , ExactHypergeomDistr(..)
 
     -- * Linking to production distributions
-    , ProductionProbFuncs(..)
-    , productionProbFuncs
     , ProductionLinkage
 
     -- * Individual test routines
@@ -88,6 +87,7 @@
 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
@@ -295,92 +295,84 @@
 --
 ----------------------------------------------------------------
 
--- | Distribution evaluation functions.
---
--- This is used to store a
-data ProductionProbFuncs = ProductionProbFuncs {
-        prodProb            :: Int -> Double
-    ,   prodCumulative      :: Double -> Double
-    ,   prodComplCumulative :: Double -> Double
-    }
-
-productionProbFuncs :: (DiscreteDistr a) => a -> ProductionProbFuncs
-productionProbFuncs d = ProductionProbFuncs {
-        prodProb = probability d
-    ,   prodCumulative = cumulative d
-    ,   prodComplCumulative = complCumulative d
-    }
-
-class (ExactDiscreteDistr a) => ProductionLinkage a where
-    productionLinkage :: a -> ProductionProbFuncs
+class (ExactDiscreteDistr a, DiscreteDistr (ProdDistrib a)
+      ) => ProductionLinkage a where
+  type ProdDistrib a
+  toProd :: a -> ProdDistrib a
 
 instance ProductionLinkage ExactBinomialDistr where
-    productionLinkage (ExactBD n p) =
-        let d = binomial (fromIntegral n) (fromRational p)
-        in  productionProbFuncs d
+  type ProdDistrib ExactBinomialDistr = BinomialDistribution
+  toProd (ExactBD n p) = binomial (fromIntegral n) (fromRational p)
 
 instance ProductionLinkage ExactDiscreteUniformDistr where
-    productionLinkage (ExactDU lower upper) =
-        let d = discreteUniformAB (fromIntegral lower) (fromIntegral upper)
-        in  productionProbFuncs d
+  type ProdDistrib ExactDiscreteUniformDistr = DiscreteUniform
+  toProd (ExactDU lower upper) = discreteUniformAB (fromIntegral lower) (fromIntegral upper)
 
 instance ProductionLinkage ExactGeometricDistr where
-    productionLinkage (ExactGeom p) =
-        let d = geometric $ fromRational p
-        in  productionProbFuncs d
+  type ProdDistrib ExactGeometricDistr = GeometricDistribution
+  toProd (ExactGeom p) = geometric $ fromRational p
 
 instance ProductionLinkage ExactHypergeomDistr where
-    productionLinkage (ExactHG nK nN n) =
-        let d = hypergeometric (fromIntegral nK) (fromIntegral nN) (fromIntegral n)
-        in  productionProbFuncs d
+  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 -> Bool
-pmfMatch tol (TestCase dExact k) =
-    let dProd = productionLinkage dExact
-        pe = fromRational $ exactProb dExact k
-        pa = prodProb dProd k'
-        k' = fromIntegral k
-    in  relativeError pe pa < tol
+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) =
-    let dProd = productionLinkage dExact
-        pe = fromRational $ exactCumulative dExact k
-        pa = prodCumulative dProd k'
-        k' = fromIntegral k
-    in  relativeError pe pa < tol
+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) =
-    let dProd = productionLinkage dExact
-        pe = fromRational $ 1 - exactCumulative dExact k
-        pa = prodComplCumulative dProd k'
-        k' = fromIntegral k
-    in  relativeError pe pa < tol
+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 :: (Show a, ProductionLinkage a, Arbitrary (TestCase a)) =>
+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 tol) :: TestCase a -> Bool)
-    ,   testProperty "CDF match" $ ((cdfMatch tol) :: TestCase a -> Bool)
-    ,   testProperty "1 - CDF match" $ ((complCdfMatch tol) :: TestCase a -> Bool)
+  testGroup ("Exact tests for " ++ name)
+    [ testProperty "PMF match"     $ pmfMatch      @a tol
+    , testProperty "CDF match"     $ cdfMatch      @a tol
+    , testProperty "1 - CDF match" $ complCdfMatch @a tol
     ]
 
 
@@ -388,9 +380,8 @@
 
 exactDistributionTests :: TestTree
 exactDistributionTests = testGroup "Test distributions against exact"
-  [
-    distTests (Tag :: Tag ExactBinomialDistr)       "Binomial"          1.0e-12
-  , distTests (Tag :: Tag ExactDiscreteUniformDistr) "DiscreteUniform"  1.0e-12
-  , distTests (Tag :: Tag ExactGeometricDistr)      "Geometric"         1.0e-13
-  , distTests (Tag :: Tag ExactHypergeomDistr)      "Hypergeometric"    1.0e-12
+  [ 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/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/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"]
+