diff --git a/Changelog.md b/Changelog.md
new file mode 100644
--- /dev/null
+++ b/Changelog.md
@@ -0,0 +1,59 @@
+# Changes in 1.1.5
+
+- Bifunctor instance is added to `Weighted`
+
+
+# Changes in 1.1.4
+
+- Actually export `CountW`
+
+
+# Changes in 1.1.3
+
+- `Data` and `Storable` instances for `CountG`.
+
+- `CalcNEvt` type class added and `CountW` accumulator for counting weighted
+  events.
+
+
+# Changes in 1.1.2
+
+- `Unbox` instances for `MeanNaive`, `WMeanNaive`, `WMeanKBN`.
+
+# Changes in 1.1.1
+
+- `Unbox` instance for `BinomAcc` is added.
+
+
+# Changes in 1.1.0
+
+- Type classes `CalcMean` and `CalcVar` are generalized to use `MonadThrow` to
+  signal failure instead of using `Maybe` only
+
+- Functions for computing standard deviation are placed into type
+  classes. Sometimes we have standard deviation at hand, if distribution is
+  parameterized by it for example.
+
+- `Mean` now type synonym for `MeanKBN`.
+
+- `WelfordMean` and `KahanMean` are moved to `D.M.S.Extra` module.
+
+- Support for calculating weighted mean.
+
+- `StatMonoid` instances for up to 4-tuples.
+
+- `Max` now works correctly (#2).
+
+- `PPair` for use in parallel computation is added.
+
+
+# Changes in 1.0.0.0
+
+- Type class definition changed: now it has both `addValue :: m → a → m` and
+  `singletonMonoid :: a → m`
+
+- `Mean` renamed as `WelfordMean`
+
+- `Unbox` instances added for all data types.
+
+- `BinomAcc` added.
diff --git a/Data/Monoid/Statistics.hs b/Data/Monoid/Statistics.hs
--- a/Data/Monoid/Statistics.hs
+++ b/Data/Monoid/Statistics.hs
@@ -1,7 +1,3 @@
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE BangPatterns          #-}
-{-# LANGUAGE DeriveDataTypeable    #-}
 -- |
 -- Module     : Data.Monoid.Statistics
 -- Copyright  : Copyright (c) 2010, Alexey Khudyakov <alexey.skladnoy@gmail.com>
@@ -9,136 +5,11 @@
 -- Maintainer : Alexey Khudyakov <alexey.skladnoy@gmail.com>
 -- Stability  : experimental
 -- 
-module Data.Monoid.Statistics ( 
-    -- * Type class
-    StatMonoid(..)
-  , evalStatistic
-    -- ** Examples
-    -- $examples
-    -- * Generic monoid
-  , TwoStats(..)
-    -- * Additional information
-    -- $info
+module Data.Monoid.Statistics (
+    module Data.Monoid.Statistics.Class
+  , module Data.Monoid.Statistics.Numeric
   ) where
 
-
-import Data.Monoid
-import Data.Typeable (Typeable)
-import qualified Data.Foldable as F
-
-
-
--- | Monoid which corresponds to some stattics. In order to do so it
---   must be commutative. In many cases it's not practical to
---   construct monoids for each element so 'papennd' was added.
---   First parameter of type class is monoidal accumulator. Second is
---   type of element over which statistic is calculated. 
---
---   Statistic could be calculated with fold over sample. Since
---   accumulator is 'Monoid' such fold could be easily parralelized.
---   Check examples section for more information.
---
---   Instance must satisfy following law:
---
---   > pappend x (pappend y mempty) == pappend x mempty `mappend` pappend y mempty
---   > mappend x y == mappend y x
---
---   It is very similar to Reducer type class from monoids package but
---   require commutative monoids
-class Monoid m => StatMonoid m a where
-  -- | Add one element to monoid accumulator. P stands for point in
-  --   analogy for Pointed.
-  pappend :: a -> m -> m
-
--- | Calculate statistic over 'Foldable'. It's implemented in terms of
---   foldl'.
-evalStatistic :: (F.Foldable d, StatMonoid m a) => d a -> m
-evalStatistic = F.foldl' (flip pappend) mempty
-  
-
-----------------------------------------------------------------
--- Generic monoids
-----------------------------------------------------------------
-
--- | Monoid which allows to calculate two statistics in parralel
-data TwoStats a b = TwoStats { calcStat1 :: !a
-                             , calcStat2 :: !b
-                             }
-                    deriving (Show,Eq,Typeable)
-
-instance (Monoid a, Monoid b) => Monoid (TwoStats a b) where
-  mempty = TwoStats mempty mempty
-  mappend !(TwoStats x y) !(TwoStats x' y') = 
-    TwoStats (mappend x x') (mappend y y')
-  {-# INLINE mempty  #-}
-  {-# INLINE mappend #-}
-
-instance (StatMonoid a x, StatMonoid b x) => StatMonoid (TwoStats a b) x where
-  pappend !x !(TwoStats a b) = TwoStats (pappend x a) (pappend x b)
-  {-# INLINE pappend #-}
-
-            
--- $info
---
--- Statistic is function of a sample which does not depend on order of
--- elements in a sample. For each statistics corresponding monoid
--- could be constructed:
---
--- > f :: [A] -> B
--- >
--- > data F = F [A]
--- >
--- > evalF (F xs) = f xs
--- >
--- > instance Monoid F here
--- >   mempty = F []
--- >   (F a) `mappend` (F b) = F (a ++ b)
---
--- This indeed proves that monoid could be constructed. Monoid above
--- is completely impractical. It runs in O(n) space. However for some
--- statistics monoids which runs in O(1) space could be
--- implemented. Simple examples of such statistics are number of
--- elements in sample or mean of a sample.
---
--- On the other hand some statistics could not be implemented in such
--- way. For example calculation of median require O(n) space. Variance
--- could be implemented in O(1) but such implementation will have
--- problems with numberical stability.
-
-
+import Data.Monoid.Statistics.Class
+import Data.Monoid.Statistics.Numeric
 
--- $examples
---
--- These examples show how to find maximum and minimum of a sample in
--- one pass over data.
--- 
--- This is test data. It's not limited to list but could be anything
--- what could be folded.
---
--- > > let xs = [1..100] :: [Double]
--- 
--- Now let calculate maximum of test sample using two methods. First
--- one is to use generic function 'evalStatistic' and another one is
--- fold.
---
--- > > evalStatistic xs :: Max
--- > Max {calcMax = 100.0}
--- > > foldl (flip pappend) mempty xs :: Max
--- > Max {calcMax = 100.0}
---
--- More complicated example allows to combine several monoids
--- together. It allows to calculate two statistics in one pass:
---
--- > > evalStatistic xs :: TwoStats Min Max
--- > TwoStats {calcStat1 = Min {calcMin = 1.0}, calcStat2 = Max {calcMax = 100.0}}
---
--- Last example shows how to calculate nuber of elements, mean and
--- variance at once:
---
--- > > let v = evalStatistic xs :: Variance
--- > > calcCount v
--- > 100
--- > > calcMean v
--- > 50.5
--- > > calcStddev v
--- > 28.86607004772212
diff --git a/Data/Monoid/Statistics/Class.hs b/Data/Monoid/Statistics/Class.hs
new file mode 100644
--- /dev/null
+++ b/Data/Monoid/Statistics/Class.hs
@@ -0,0 +1,449 @@
+{-# LANGUAGE BangPatterns               #-}
+{-# LANGUAGE DefaultSignatures          #-}
+{-# LANGUAGE DeriveDataTypeable         #-}
+{-# LANGUAGE DeriveFoldable             #-}
+{-# LANGUAGE DeriveFunctor              #-}
+{-# LANGUAGE DeriveGeneric              #-}
+{-# LANGUAGE DeriveTraversable          #-}
+{-# LANGUAGE DerivingStrategies         #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses      #-}
+{-# LANGUAGE RankNTypes                 #-}
+{-# LANGUAGE TemplateHaskell            #-}
+{-# LANGUAGE TypeFamilies               #-}
+-- |
+-- Module     : Data.Monoid.Statistics
+-- Copyright  : Copyright (c) 2010,2017, Alexey Khudyakov <alexey.skladnoy@gmail.com>
+-- License    : BSD3
+-- Maintainer : Alexey Khudyakov <alexey.skladnoy@gmail.com>
+-- Stability  : experimental
+--
+module Data.Monoid.Statistics.Class
+  ( -- * Monoid Type class and helpers
+    StatMonoid(..)
+  , reduceSample
+  , reduceSampleVec
+    -- * Ad-hoc type classes for select statistics
+    -- $adhoc
+  , CalcCount(..)
+  , CalcNEvt(..)
+  , CalcMean(..)
+  , HasMean(..)
+  , CalcVariance(..)
+  , HasVariance(..)
+    -- ** Deriving via
+  , CalcViaHas(..)
+    -- * Exception handling
+  , Partial(..)
+  , partial
+  , SampleError(..)
+    -- * Data types
+  , Pair(..)
+  ) where
+
+import           Control.Exception
+import           Control.Monad.Catch (MonadThrow(..))
+import           Data.Data           (Typeable,Data)
+import           Data.Monoid
+import           Data.Int
+import           Data.Word
+import           Data.Vector.Unboxed          (Unbox)
+import           Data.Vector.Unboxed.Deriving (derivingUnbox)
+import qualified Data.Foldable       as F
+import qualified Data.Vector.Generic as G
+import           Numeric.Sum
+import GHC.Stack    (HasCallStack)
+import GHC.Generics (Generic)
+
+
+-- | This type class is used to express parallelizable constant space
+--   algorithms for calculation of statistics. /Statistic/ is function
+--   of type @[a]→b@ which does not depend on order of elements. (for
+--   example: mean, sum, number of elements, variance, etc).
+--
+--   For many statistics it's possible to possible to construct
+--   constant space algorithm which is expressed as fold. Additionally
+--   it's usually possible to write function which combine state of
+--   fold accumulator to get statistic for union of two samples.
+--
+--   Thus for such algorithm we have value which corresponds to empty
+--   sample, function which which corresponds to merging of two
+--   samples, and single step of fold. Last one allows to evaluate
+--   statistic given data sample and first two form a monoid and allow
+--   parallelization: split data into parts, build estimate for each
+--   by folding and then merge them using mappend.
+--
+--   Instance must satisfy following laws. If floating point
+--   arithmetics is used then equality should be understood as
+--   approximate.
+--
+--   > 1. addValue (addValue y mempty) x  == addValue mempty x <> addValue mempty y
+--   > 2. x <> y == y <> x
+class Monoid m => StatMonoid m a where
+  -- | Add one element to monoid accumulator. It's step of fold.
+  addValue :: m -> a -> m
+  addValue m a = m <> singletonMonoid a
+  {-# INLINE addValue #-}
+  -- | State of accumulator corresponding to 1-element sample.
+  singletonMonoid :: a -> m
+  singletonMonoid = addValue mempty
+  {-# INLINE singletonMonoid #-}
+  {-# MINIMAL addValue | singletonMonoid #-}
+
+-- | Calculate statistic over 'Foldable'. It's implemented in terms of
+--   foldl'. Note that in cases when accumulator is immediately
+--   consumed by polymorphic function such as 'callMeam' its type
+--   becomes ambiguous. @TypeApplication@ then could be used to
+--   disambiguate.
+--
+-- >>> reduceSample @Mean [1,2,3,4]
+-- MeanKBN 4 (KBNSum 10.0 0.0)
+-- >>> calcMean $ reduceSample @Mean [1,2,3,4] :: Maybe Double
+-- Just 2.5
+reduceSample :: forall m a f. (StatMonoid m a, F.Foldable f) => f a -> m
+reduceSample = F.foldl' addValue mempty
+
+-- | Calculate statistic over vector. Works in same was as
+-- 'reduceSample' but works for vectors.
+reduceSampleVec :: forall m a v. (StatMonoid m a, G.Vector v a) => v a -> m
+reduceSampleVec = G.foldl' addValue mempty
+{-# INLINE reduceSampleVec #-}
+
+instance ( StatMonoid m1 a
+         , StatMonoid m2 a
+         ) => StatMonoid (m1,m2) a where
+  addValue (!m1, !m2) a =
+    let !m1' = addValue m1 a
+        !m2' = addValue m2 a
+    in (m1', m2')
+  singletonMonoid a = ( singletonMonoid a
+                      , singletonMonoid a
+                      )
+
+instance ( StatMonoid m1 a
+         , StatMonoid m2 a
+         , StatMonoid m3 a
+         ) => StatMonoid (m1,m2,m3) a where
+  addValue (!m1, !m2, !m3) a =
+    let !m1' = addValue m1 a
+        !m2' = addValue m2 a
+        !m3' = addValue m3 a
+    in (m1', m2', m3')
+  singletonMonoid a = ( singletonMonoid a
+                      , singletonMonoid a
+                      , singletonMonoid a
+                      )
+
+instance ( StatMonoid m1 a
+         , StatMonoid m2 a
+         , StatMonoid m3 a
+         , StatMonoid m4 a
+         ) => StatMonoid (m1,m2,m3,m4) a where
+  addValue (!m1, !m2, !m3, !m4) a =
+    let !m1' = addValue m1 a
+        !m2' = addValue m2 a
+        !m3' = addValue m3 a
+        !m4' = addValue m4 a
+    in (m1', m2', m3', m4')
+  singletonMonoid a = ( singletonMonoid a
+                      , singletonMonoid a
+                      , singletonMonoid a
+                      , singletonMonoid a
+                      )
+
+instance (Num a, a ~ a') => StatMonoid (Sum a) a' where
+  singletonMonoid = Sum
+
+instance (Num a, a ~ a') => StatMonoid (Product a) a' where
+  singletonMonoid = Product
+
+instance Real a => StatMonoid KahanSum a where
+  addValue m x = add m (realToFrac x)
+  {-# INLINE addValue #-}
+
+instance Real a => StatMonoid KBNSum a where
+  addValue m x = add m (realToFrac x)
+  {-# INLINE addValue #-}
+
+instance Real a => StatMonoid KB2Sum a where
+  addValue m x = add m (realToFrac x)
+  {-# INLINE addValue #-}
+
+
+----------------------------------------------------------------
+-- Ad-hoc type class
+----------------------------------------------------------------
+
+-- $adhoc
+--
+-- Type classes defined here allows to extract common statistics from
+-- estimators. it's assumed that quantities in question are already
+-- computed so extraction is cheap.
+--
+--
+-- ==== Error handling
+--
+-- Computation of statistics may fail. For example mean is not defined
+-- for an empty sample. @Maybe@ could be seen as easy way to handle
+-- this situation. But in many cases most convenient way to handle
+-- failure is to throw an exception. So failure is encoded by using
+-- polymorphic function of type @MonadThrow m ⇒ a → m X@.
+--
+-- Maybe types has instance, such as 'Maybe', 'Either'
+-- 'Control.Exception.SomeException', 'IO' and most transformers
+-- wrapping it. Notably this library defines 'Partial' monad which
+-- allows to convert failures to exception in pure setting.
+--
+-- >>> calcMean $ reduceSample @Mean []
+-- *** Exception: EmptySample "Data.Monoid.Statistics.Numeric.MeanKBN: calcMean"
+--
+-- >>> calcMean $ reduceSample @Mean [] :: Maybe Double
+-- Nothing
+--
+-- >>> import Control.Exception
+-- >>> calcMean $ reduceSample @Mean [] :: Either SomeException Double
+-- Left (EmptySample "Data.Monoid.Statistics.Numeric.MeanKBN: calcMean")
+--
+-- Last example uses IO
+--
+-- >>> calcMean $ reduceSample @Mean []
+-- *** Exception: EmptySample "Data.Monoid.Statistics.Numeric.MeanKBN: calcMean"
+--
+--
+-- ==== Deriving instances
+--
+-- Type classes come in two variants, one that allow failure and one
+-- for use in cases when quantity is always defined. This is not the
+-- case for estimators, but true for distributions and intended for
+-- such use cases. In that case 'CalcViaHas' could be used to derive
+-- necessary instances.
+--
+-- >>> :{
+-- data NormalDist = NormalDist !Double !Double
+--   deriving (CalcMean,CalcVariance) via CalcViaHas NormalDist
+-- instance HasMean NormalDist where
+--   getMean (NormalDist mu _) = mu
+-- instance HasVariance NormalDist where
+--   getVariance   (NormalDist _ s) = s
+--   getVarianceML (NormalDist _ s) = s
+-- :}
+
+
+-- | Value from which we can efficiently extract number of elements in
+--   sample it represents.
+class CalcCount a where
+  -- | /Assumed O(1)/. Number of elements in sample.
+  calcCount :: a -> Int
+
+-- | Value from which we can efficiently calculate mean of sample or
+--   distribution.
+class CalcMean a where
+  -- | /Assumed O(1)/ Returns @Nothing@ if there isn't enough data to
+  --   make estimate or distribution doesn't have defined mean.
+  --
+  --   \[ \bar{x} = \frac{1}{N}\sum_{i=1}^N{x_i} \]
+  calcMean :: MonadThrow m => a -> m Double
+
+-- | Same as 'CalcMean' but should never fail
+class CalcMean a => HasMean a where
+  getMean :: a -> Double
+
+
+-- | Values from which we can efficiently compute estimate of sample
+--   variance or distribution variance. It has two methods: one which
+--   applies bias correction to estimate and another that returns
+--   maximul likelyhood estimate. For distribution they should return
+--   same value.
+class CalcVariance a where
+  -- | /Assumed O(1)/ Calculate unbiased estimate of variance:
+  --
+  --   \[ \sigma^2 = \frac{1}{N-1}\sum_{i=1}^N(x_i - \bar{x})^2 \]
+  calcVariance :: MonadThrow m => a -> m Double
+  calcVariance = fmap (\x->x*x) . calcStddev
+  -- | /Assumed O(1)/ Calculate maximum likelihood estimate of variance:
+  --
+  --   \[ \sigma^2 = \frac{1}{N}\sum_{i=1}^N(x_i - \bar{x})^2 \]
+  calcVarianceML :: MonadThrow m => a -> m Double
+  calcVarianceML = fmap (\x->x*x) . calcStddevML
+  -- | Calculate sample standard deviation from unbiased estimation of
+  --   variance.
+  calcStddev :: MonadThrow m => a -> m Double
+  calcStddev = fmap sqrt . calcVariance
+  -- | Calculate sample standard deviation from maximum likelihood
+  --   estimation of variance.
+  calcStddevML :: (MonadThrow m) => a -> m Double
+  calcStddevML = fmap sqrt . calcVarianceML
+  {-# MINIMAL (calcVariance,calcVarianceML) | (calcStddev,calcStddevML) #-}
+
+-- | Same as 'CalcVariance' but never fails
+class CalcVariance a => HasVariance a where
+  getVariance   :: a -> Double
+  getVariance   = (\x -> x*x) . getStddev
+  getVarianceML :: a -> Double
+  getVarianceML = (\x -> x*x) . getStddevML
+  getStddev     :: a -> Double
+  getStddev     = sqrt . getVariance
+  getStddevML   :: a -> Double
+  getStddevML   = sqrt . getVarianceML
+  {-# MINIMAL (getVariance,getVarianceML) | (getStddev,getStddevML) #-}
+
+
+-- | Type class for accumulators that are used for event counting with
+--   possibly weighted events. Those are mostly used as accumulators
+--   in histograms.
+class CalcNEvt a where
+  -- | Calculate sum of events weights.
+  calcEvtsW :: a -> Double
+  -- | Calculate error estimate (1σ or 68% CL). All instances defined
+  --   in library use normal approximation which breaks down for small
+  --   number of events.
+  calcEvtsWErr :: a -> Double
+  -- | Calculate effective number of events which is defined as
+  --   \[N=E(w)^2/\operatorname{Var}(w)\] or as number of events
+  --   which will yield same estimate for mean variance is they all
+  --   have same weight.
+  calcEffNEvt :: a -> Double
+  calcEffNEvt = calcEvtsW
+
+instance CalcNEvt Int where
+  calcEvtsW    = fromIntegral
+  calcEvtsWErr = sqrt . calcEvtsW
+
+instance CalcNEvt Int32 where
+  calcEvtsW    = fromIntegral
+  calcEvtsWErr = sqrt . calcEvtsW
+
+instance CalcNEvt Int64 where
+  calcEvtsW    = fromIntegral
+  calcEvtsWErr = sqrt . calcEvtsW
+
+instance CalcNEvt Word where
+  calcEvtsW    = fromIntegral
+  calcEvtsWErr = sqrt . calcEvtsW
+
+instance CalcNEvt Word32 where
+  calcEvtsW    = fromIntegral
+  calcEvtsWErr = sqrt . calcEvtsW
+
+instance CalcNEvt Word64 where
+  calcEvtsW    = fromIntegral
+  calcEvtsWErr = sqrt . calcEvtsW
+
+
+-- | Derive instances for 'CalcMean' and 'CalcVariance' from 'HasMean'
+--   and 'HasVariance' instances.
+newtype CalcViaHas a = CalcViaHas a
+  deriving newtype (HasMean, HasVariance)
+
+instance HasMean a => CalcMean (CalcViaHas a) where
+  calcMean = pure . getMean
+
+instance HasVariance a => CalcVariance (CalcViaHas a) where
+  calcVariance   = pure . getVariance
+  calcVarianceML = pure . getVarianceML
+
+----------------------------------------------------------------
+-- Exceptions
+----------------------------------------------------------------
+
+-- | Identity monad which is used to encode partial functions for
+--   'MonadThrow' based error handling. Its @MonadThrow@ instance
+--   just throws normal exception.
+newtype Partial a = Partial a
+  deriving (Show, Read, Eq, Ord, Typeable, Data, Generic)
+
+-- | Convert error to IO exception. This way one could for example
+--   convert case when some statistics is not defined to an exception:
+--
+-- >>> calcMean $ reduceSample @Mean []
+-- *** Exception: EmptySample "Data.Monoid.Statistics.Numeric.MeanKBN: calcMean"
+partial :: HasCallStack => Partial a -> a
+partial (Partial x) = x
+
+instance Functor Partial where
+  fmap f (Partial a) = Partial (f a)
+
+instance Applicative Partial where
+  pure = Partial
+  Partial f <*> Partial a = Partial (f a)
+  (!_) *> a   = a
+  a   <* (!_) = a
+instance Monad Partial where
+  return = pure
+  Partial a >>= f = f a
+  (>>) = (*>)
+
+instance MonadThrow Partial where
+  throwM = throw
+
+-- | Exception which is thrown when we can't compute some value
+data SampleError
+  = EmptySample String
+  -- ^ @EmptySample function@: We're trying to compute quantity that
+  --   is undefined for empty sample.
+  | InvalidSample String String
+  -- ^ @InvalidSample function descripton@ quantity in question could
+  --   not be computed for some other reason
+  deriving Show
+
+instance Exception SampleError
+
+
+----------------------------------------------------------------
+-- Generic monoids
+----------------------------------------------------------------
+
+-- | Strict pair. It allows to calculate two statistics in parallel
+data Pair a b = Pair !a !b
+              deriving (Show,Eq,Ord,Typeable,Data,Generic)
+
+instance (Semigroup a, Semigroup b) => Semigroup (Pair a b) where
+  Pair x y <> Pair x' y' = Pair (x <> x') (y <> y')
+  {-# INLINABLE (<>) #-}
+
+instance (Monoid a, Monoid b) => Monoid (Pair a b) where
+  mempty  = Pair mempty mempty
+  mappend = (<>)
+  {-# INLINABLE mempty  #-}
+  {-# INLINABLE mappend #-}
+
+instance (StatMonoid a x, StatMonoid b x) => StatMonoid (Pair a b) x where
+  addValue (Pair a b) !x = Pair (addValue a x) (addValue b x)
+  singletonMonoid x = Pair (singletonMonoid x) (singletonMonoid x)
+  {-# INLINE addValue        #-}
+  {-# INLINE singletonMonoid #-}
+
+
+-- | Strict pair for parallel accumulation
+data PPair a b = PPair !a !b
+
+instance (Semigroup a, Semigroup b) => Semigroup (PPair a b) where
+  PPair x y <> PPair x' y' = PPair (x <> x') (y <> y')
+  {-# INLINABLE (<>) #-}
+
+instance (Monoid a, Monoid b) => Monoid (PPair a b) where
+  mempty  = PPair mempty mempty
+  mappend = (<>)
+  {-# INLINABLE mempty  #-}
+  {-# INLINABLE mappend #-}
+
+instance (StatMonoid a x, StatMonoid b y) => StatMonoid (PPair a b) (x,y) where
+  addValue (PPair a b) (!x,!y) = PPair (addValue a x) (addValue b y)
+  singletonMonoid (!x,!y) = PPair (singletonMonoid x) (singletonMonoid y)
+  {-# INLINE addValue        #-}
+  {-# INLINE singletonMonoid #-}
+
+
+
+
+
+derivingUnbox "Pair"
+  [t| forall a b. (Unbox a, Unbox b) => Pair a b -> (a,b) |]
+  [| \(Pair a b) -> (a,b) |]
+  [| \(a,b) -> Pair a b   |]
+
+-- $setup
+--
+-- >>> :set -XDerivingVia
+-- >>> import Data.Monoid.Statistics.Numeric
diff --git a/Data/Monoid/Statistics/Extra.hs b/Data/Monoid/Statistics/Extra.hs
new file mode 100644
--- /dev/null
+++ b/Data/Monoid/Statistics/Extra.hs
@@ -0,0 +1,175 @@
+{-# LANGUAGE BangPatterns          #-}
+{-# LANGUAGE DeriveDataTypeable    #-}
+{-# LANGUAGE DeriveGeneric         #-}
+{-# LANGUAGE FlexibleContexts      #-}
+{-# LANGUAGE FlexibleInstances     #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TemplateHaskell       #-}
+{-# LANGUAGE TypeFamilies          #-}
+-- |
+-- Monoids for calculating various statistics in constant space. This
+-- module contains algorithms that should be generally avoided unless
+-- there's specific reason to use them.
+module Data.Monoid.Statistics.Extra (
+    -- * Mean
+    WelfordMean(..)
+  , asWelfordMean
+  , MeanKahan(..)
+  , asMeanKahan
+  , MeanKB2(..)
+  , asMeanKB2
+    -- $references
+  ) where
+
+import Control.Monad.Catch          (MonadThrow(..))
+import Data.Data                    (Typeable,Data)
+import Data.Vector.Unboxed.Deriving (derivingUnbox)
+import Numeric.Sum
+import GHC.Generics                 (Generic)
+
+import Data.Monoid.Statistics.Class
+
+
+
+----------------------------------------------------------------
+-- Mean
+----------------------------------------------------------------
+
+
+-- | Incremental calculation of mean which uses second-order
+--   compensated Kahan-Babuška summation. In most cases
+--   'Data.Monoid.Statistics.Numeric.KBNSum' should provide enough
+--   precision.
+data MeanKB2 = MeanKB2 !Int {-# UNPACK #-} !KB2Sum
+             deriving (Show,Eq)
+
+asMeanKB2 :: MeanKB2 -> MeanKB2
+asMeanKB2 = id
+
+instance Semigroup MeanKB2 where
+  MeanKB2 0  _  <> m             = m
+  m             <> MeanKB2 0  _  = m
+  MeanKB2 n1 s1 <> MeanKB2 n2 s2 = MeanKB2 (n1+n2) (s1 <> s2)
+
+instance Monoid MeanKB2 where
+  mempty  = MeanKB2 0 mempty
+  mappend = (<>)
+
+instance Real a => StatMonoid MeanKB2 a where
+  addValue (MeanKB2 n m) x = MeanKB2 (n+1) (addValue m x)
+
+instance CalcMean MeanKB2 where
+  calcMean (MeanKB2 0 _) = throwM $ EmptySample "Data.Monoid.Statistics.Extra.MeanKB2"
+  calcMean (MeanKB2 n s) = return $! kb2 s / fromIntegral n
+
+
+
+-- | Incremental calculation of mean. Sum of elements is calculated
+--   using compensated Kahan summation. It's provided only for sake of
+--   completeness. 'Data.Monoid.Statistics.Numeric.KBNSum' should be used
+--   instead.
+data MeanKahan = MeanKahan !Int !KahanSum
+             deriving (Show,Eq,Typeable,Data,Generic)
+
+asMeanKahan :: MeanKahan -> MeanKahan
+asMeanKahan = id
+
+
+instance Semigroup MeanKahan where
+  MeanKahan 0  _  <> m               = m
+  m               <> MeanKahan 0  _  = m
+  MeanKahan n1 s1 <> MeanKahan n2 s2 = MeanKahan (n1+n2) (s1 <> s2)
+  {-# INLINE (<>) #-}
+
+instance Monoid MeanKahan where
+  mempty  = MeanKahan 0 mempty
+  mappend = (<>)
+
+instance Real a => StatMonoid MeanKahan a where
+  addValue (MeanKahan n m) x = MeanKahan (n+1) (addValue m x)
+
+instance CalcCount MeanKahan where
+  calcCount (MeanKahan n _) = n
+instance CalcMean MeanKahan where
+  calcMean (MeanKahan 0 _) = throwM $ EmptySample "Data.Monoid.Statistics.Extra.WelfordMean"
+  calcMean (MeanKahan n s) = return (kahan s / fromIntegral n)
+
+
+-- | Incremental calculation of mean. Note that this algorithm doesn't
+--   offer better numeric precision than plain summation. Its only
+--   advantage is protection against double overflow:
+--
+-- >>> calcMean $ reduceSample @MeanKBN (replicate 100 1e308) :: Maybe Double
+-- Just NaN
+-- >>> calcMean $ reduceSample @WelfordMean (replicate 100 1e308) :: Maybe Double
+-- Just 1.0e308
+--
+--   Unless this feature is needed 'Data.Monoid.Statistics.Numeric.KBNSum'
+--   should be used. Algorithm is due to Welford [Welford1962]
+--
+-- \[ s_n = s_{n-1} + \frac{x_n - s_{n-1}}{n} \]
+data WelfordMean = WelfordMean !Int    -- Number of entries
+                               !Double -- Current mean
+  deriving (Show,Eq,Typeable,Data,Generic)
+
+-- | Type restricted 'id'
+asWelfordMean :: WelfordMean -> WelfordMean
+asWelfordMean = id
+
+instance Semigroup WelfordMean where
+  WelfordMean 0 _ <> m = m
+  m <> WelfordMean 0 _ = m
+  WelfordMean n x <> WelfordMean k y
+    = WelfordMean (n + k) ((x*n' + y*k') / (n' + k'))
+    where
+      n' = fromIntegral n
+      k' = fromIntegral k
+  {-# INLINE (<>) #-}
+
+instance Monoid WelfordMean where
+  mempty  = WelfordMean 0 0
+  mappend = (<>)
+  {-# INLINE mempty  #-}
+  {-# INLINE mappend #-}
+
+instance Real a => StatMonoid WelfordMean a where
+  addValue (WelfordMean n m) !x
+    = WelfordMean n' (m + (realToFrac x - m) / fromIntegral n')
+    where
+      n' = n+1
+  {-# INLINE addValue #-}
+
+instance CalcCount WelfordMean where
+  calcCount (WelfordMean n _) = n
+instance CalcMean WelfordMean where
+  calcMean (WelfordMean 0 _) = throwM $ EmptySample "Data.Monoid.Statistics.Extra.WelfordMean"
+  calcMean (WelfordMean _ m) = return m
+
+
+
+----------------------------------------------------------------
+-- Unboxed instances
+----------------------------------------------------------------
+
+derivingUnbox "MeanKahan"
+  [t| MeanKahan -> (Int,Double,Double) |]
+  [| \(MeanKahan a (KahanSum b c)) -> (a,b,c)   |]
+  [| \(a,b,c) -> MeanKahan a (KahanSum b c) |]
+
+derivingUnbox "WelfordMean"
+  [t| WelfordMean -> (Int,Double) |]
+  [| \(WelfordMean a b) -> (a,b)  |]
+  [| \(a,b) -> WelfordMean a b    |]
+
+
+-- $references
+--
+-- * [Welford1962] Welford, B.P. (1962) Note on a method for
+--   calculating corrected sums of squares and
+--   products. /Technometrics/
+--   4(3):419-420. <http://www.jstor.org/stable/1266577>
+
+-- $setup
+--
+-- >>> :set -XTypeApplications
+-- >>> import Data.Monoid.Statistics.Numeric
diff --git a/Data/Monoid/Statistics/Numeric.hs b/Data/Monoid/Statistics/Numeric.hs
--- a/Data/Monoid/Statistics/Numeric.hs
+++ b/Data/Monoid/Statistics/Numeric.hs
@@ -1,256 +1,574 @@
-{-# LANGUAGE BangPatterns          #-}
-{-# LANGUAGE FlexibleContexts      #-}
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE DeriveDataTypeable    #-}
-module Data.Monoid.Statistics.Numeric ( 
-    -- * Mean and variance
-    Count(..)
+{-# LANGUAGE BangPatterns               #-}
+{-# LANGUAGE DeriveAnyClass             #-}
+{-# LANGUAGE DeriveDataTypeable         #-}
+{-# LANGUAGE DeriveFoldable             #-}
+{-# LANGUAGE DeriveGeneric              #-}
+{-# LANGUAGE DeriveTraversable          #-}
+{-# LANGUAGE DerivingStrategies         #-}
+{-# LANGUAGE DerivingVia                #-}
+{-# LANGUAGE FlexibleContexts           #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses      #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE StandaloneDeriving         #-}
+{-# LANGUAGE TemplateHaskell            #-}
+{-# LANGUAGE TypeFamilies               #-}
+{-# LANGUAGE TypeOperators              #-}
+{-# LANGUAGE ViewPatterns               #-}
+-- |
+-- Monoids for calculating various statistics in constant space
+module Data.Monoid.Statistics.Numeric (
+    -- * Mean & Variance
+    -- ** Number of elements
+    CountG(..)
+  , Count
   , asCount
-  , Mean(..)
+  , CountW(..)
+    -- ** Mean algorithms
+    -- ** Default algorithms
+  , Mean
   , asMean
+  , WMean
+  , asWMean
+    -- *** Mean
+  , MeanNaive(..)
+  , asMeanNaive
+  , MeanKBN(..)
+  , asMeanKBN
+    -- *** Weighted mean
+  , WMeanNaive(..)
+  , asWMeanNaive
+  , WMeanKBN(..)
+  , asWMeanKBN
+    -- ** Variance
   , Variance(..)
   , asVariance
-    -- ** Ad-hoc accessors
-    -- $accessors
-  , CalcCount(..)
-  , CalcMean(..)
-  , CalcVariance(..)
-  , calcStddev
-  , calcStddevUnbiased
     -- * Maximum and minimum
   , Max(..)
   , Min(..)
+  , MaxD(..)
+  , MinD(..)
+    -- * Binomial trials
+  , BinomAcc(..)
+  , asBinomAcc
+    -- * Rest
+  , Weighted(..)
+    -- * References
+    -- $references
   ) where
 
-import Data.Monoid
-import Data.Monoid.Statistics
-import Data.Typeable (Typeable)
+import Control.Monad.Catch          (MonadThrow(..))
+import Data.Bifunctor
+import Data.Data                    (Typeable,Data)
+import Data.Vector.Unboxed          (Unbox)
+import Data.Vector.Unboxed.Deriving (derivingUnbox)
+import qualified Data.Vector.Unboxed          as VU
+import qualified Data.Vector.Generic          as VG
+import qualified Data.Vector.Generic.Mutable  as VGM
+import Foreign.Storable             (Storable)
+import Numeric.Sum
+import GHC.Generics                 (Generic)
 
+import Data.Monoid.Statistics.Class
+
+
 ----------------------------------------------------------------
 -- Statistical monoids
 ----------------------------------------------------------------
 
--- | Simplest statistics. Number of elements in the sample
-newtype Count a = Count { calcCountI :: a }
-                  deriving (Show,Eq,Ord,Typeable)
+-- | Calculate number of elements in the sample.
+newtype CountG a = CountG { calcCountN :: a }
+  deriving stock   (Show,Eq,Ord,Data)
+  deriving newtype (Storable)
 
--- | Fix type of monoid
-asCount :: Count a -> Count a
+type Count = CountG Int
+
+-- | Type restricted 'id'
+asCount :: CountG a -> CountG a
 asCount = id
-{-# INLINE asCount #-}
 
-instance Integral a => Monoid (Count a) where
-  mempty = Count 0
-  (Count i) `mappend` (Count j) = Count (i + j)
-  {-# INLINE mempty  #-}
-  {-# INLINE mappend #-}
-  
-instance (Integral a) => StatMonoid (Count a) b where
-  pappend _ !(Count n) = Count (n + 1)
-  {-# INLINE pappend #-}
+instance Integral a => Semigroup (CountG a) where
+  CountG i <> CountG j = CountG (i + j)
 
-instance CalcCount (Count Int) where
-  calcCount = calcCountI
-  {-# INLINE calcCount #-}
+instance Integral a => Monoid (CountG a) where
+  mempty  = CountG 0
+  mappend = (<>)
 
+instance (Integral a) => StatMonoid (CountG a) b where
+  singletonMonoid _            = CountG 1
+  addValue        (CountG n) _ = CountG (n + 1)
 
 
+instance CalcCount (CountG Int) where
+  calcCount = calcCountN
 
--- | Mean of sample. Samples of Double,Float and bui;t-in integral
---   types are supported
+instance Real a => CalcNEvt (CountG a) where
+  calcEvtsW    = realToFrac . calcCountN
+  calcEvtsWErr = sqrt . calcEvtsW
+  {-# INLINE calcEvtsW    #-}
+  {-# INLINE calcEvtsWErr #-}
+
+----------------------------------------------------------------
+
+-- | Accumulator type for counting weighted events. Weights are
+--   presumed to be independent and follow same distribution \[W\].
+--   In this case sum of weights follows compound Poisson
+--   distribution. Its expectation could be then estimated as
+--   \[\sum_iw_i\] and variance as \[\sum_iw_i^2\].
 --
--- Numeric stability of 'mappend' is not proven.
-data Mean = Mean {-# UNPACK #-} !Int    -- Number of entries
-                 {-# UNPACK #-} !Double -- Current mean
-            deriving (Show,Eq,Typeable)
+--   Main use of this data type is as accumulator in histograms which
+--   count weighted events.
+data CountW = CountW
+  !Double -- Sum of weights
+  !Double -- Sum of weight squares
+  deriving stock (Show,Eq,Generic)
 
--- | Fix type of monoid
+instance Semigroup CountW where
+  CountW wA w2A <> CountW wB w2B = CountW (wA+wB) (w2A+w2B)
+  {-# INLINE (<>) #-}
+instance Monoid CountW where
+  mempty = CountW 0 0
+
+instance Real a => StatMonoid CountW a where
+  addValue (CountW w w2) a = CountW (w + x) (w2 + x*x)
+    where
+      x = realToFrac a
+
+instance CalcNEvt CountW where
+  calcEvtsW    (CountW w _ ) = w
+  calcEvtsWErr (CountW _ w2) = sqrt w2
+  calcEffNEvt  (CountW w w2) = w * w / w2
+
+----------------------------------------------------------------
+
+-- | Type alias for currently recommended algorithms for calculation
+--   of mean. It should be default choice
+type Mean = MeanKBN
+
 asMean :: Mean -> Mean
 asMean = id
-{-# INLINE asMean #-}
 
-instance Monoid Mean where
-  mempty = Mean 0 0
-  mappend !(Mean n x) !(Mean k y) = Mean (n + k) ((x*n' + y*k') / (n' + k')) 
+-- | Type alias for currently recommended algorithms for calculation
+--   of weighted mean. It should be default choice
+type WMean = WMeanKBN
+
+asWMean :: WMean -> WMean
+asWMean = id
+
+
+----------------------------------------------------------------
+
+-- | Incremental calculation of mean. It tracks separately number of
+--   elements and running sum. Note that summation of floating point
+--   numbers loses precision and genrally use 'MeanKBN' is
+--   recommended.
+data MeanNaive = MeanNaive !Int !Double
+  deriving stock (Show,Eq,Data,Generic)
+
+asMeanNaive :: MeanNaive -> MeanNaive
+asMeanNaive = id
+
+
+instance Semigroup MeanNaive where
+  MeanNaive 0  _  <> m               = m
+  m               <> MeanNaive 0  _  = m
+  MeanNaive n1 s1 <> MeanNaive n2 s2 = MeanNaive (n1+n2) (s1 + s2)
+
+instance Monoid MeanNaive where
+  mempty  = MeanNaive 0 0
+  mappend = (<>)
+
+instance Real a => StatMonoid MeanNaive a where
+  addValue (MeanNaive n m) x = MeanNaive (n+1) (m + realToFrac x)
+  {-# INLINE addValue #-}
+
+instance CalcCount MeanNaive where
+  calcCount (MeanNaive n _) = n
+instance CalcMean MeanNaive where
+  calcMean (MeanNaive 0 _) = throwM $ EmptySample "Data.Monoid.Statistics.Numeric.MeanNaive: calcMean"
+  calcMean (MeanNaive n s) = return (s / fromIntegral n)
+
+
+----------------------------------------------------------------
+
+-- | Incremental calculation of mean. It tracks separately number of
+--   elements and running sum. It uses algorithm for compensated
+--   summation which works with mantissa of double size at cost of
+--   doing more operations. This means that it's usually possible to
+--   compute sum (and therefore mean) within 1 ulp.
+data MeanKBN = MeanKBN !Int {-# UNPACK #-} !KBNSum
+  deriving stock (Show,Eq,Data,Generic)
+
+asMeanKBN :: MeanKBN -> MeanKBN
+asMeanKBN = id
+
+
+instance Semigroup MeanKBN where
+  MeanKBN 0  _  <> m             = m
+  m             <> MeanKBN 0  _  = m
+  MeanKBN n1 s1 <> MeanKBN n2 s2 = MeanKBN (n1+n2) (s1 <> s2)
+
+instance Monoid MeanKBN where
+  mempty  = MeanKBN 0 mempty
+  mappend = (<>)
+  
+instance Real a => StatMonoid MeanKBN a where
+  addValue (MeanKBN n m) x = MeanKBN (n+1) (addValue m x)
+  {-# INLINE addValue #-}
+
+instance CalcCount MeanKBN where
+  calcCount (MeanKBN n _) = n
+instance CalcMean MeanKBN where
+  calcMean (MeanKBN 0 _) = throwM $ EmptySample "Data.Monoid.Statistics.Numeric.MeanKBN: calcMean"
+  calcMean (MeanKBN n s) = return (kbn s / fromIntegral n)
+
+
+----------------------------------------------------------------
+
+-- | Incremental calculation of weighed mean.
+data WMeanNaive = WMeanNaive
+  !Double  -- Weight
+  !Double  -- Weighted sum
+  deriving stock (Show,Eq,Data,Generic)
+
+asWMeanNaive :: WMeanNaive -> WMeanNaive
+asWMeanNaive = id
+
+
+instance Semigroup WMeanNaive where
+  WMeanNaive w1 s1 <> WMeanNaive w2 s2 = WMeanNaive (w1 + w2) (s1 + s2)
+
+instance Monoid WMeanNaive where
+  mempty  = WMeanNaive 0 0
+  mappend = (<>)
+
+instance (Real w, Real a) => StatMonoid WMeanNaive (Weighted w a) where
+  addValue (WMeanNaive n s) (Weighted w a)
+    = WMeanNaive (n + w') (s + (w' * a'))
     where
-      n' = fromIntegral n
-      k' = fromIntegral k
-  {-# INLINE mempty  #-}
-  {-# INLINE mappend #-}
+      w' = realToFrac w
+      a' = realToFrac a
+  {-# INLINE addValue #-}
 
-instance Real a => StatMonoid Mean a where
-  pappend !x !(Mean n m) = Mean n' (m + (realToFrac x - m) / fromIntegral n') where n' = n+1
-  {-# INLINE pappend #-}
+instance CalcMean WMeanNaive where
+  calcMean (WMeanNaive w s)
+    | w <= 0    = throwM $ EmptySample "Data.Monoid.Statistics.Numeric.WMeanNaive: calcMean"
+    | otherwise = return (s / w)
 
-instance CalcCount Mean where
-  calcCount (Mean n _) = n
-  {-# INLINE calcCount #-}
-instance CalcMean Mean where
-  calcMean (Mean _ m) = m
-  {-# INLINE calcMean #-}
+----------------------------------------------------------------
 
+-- | Incremental calculation of weighed mean. Sum of both weights and
+--   elements is calculated using Kahan-Babuška-Neumaier summation.
+data WMeanKBN = WMeanKBN
+  {-# UNPACK #-} !KBNSum  -- Weight
+  {-# UNPACK #-} !KBNSum  -- Weighted sum
+  deriving stock (Show,Eq,Data,Generic)
 
+asWMeanKBN :: WMeanKBN -> WMeanKBN
+asWMeanKBN = id
 
 
--- | Intermediate quantities to calculate the standard deviation.
+instance Semigroup WMeanKBN where
+  WMeanKBN n1 s1 <> WMeanKBN n2 s2 = WMeanKBN (n1 <> n2) (s1 <> s2)
+
+instance Monoid WMeanKBN where
+  mempty  = WMeanKBN mempty mempty
+  mappend = (<>)
+
+instance (Real w, Real a) => StatMonoid WMeanKBN (Weighted w a) where
+  addValue (WMeanKBN n m) (Weighted w a)
+    = WMeanKBN (add n w') (add m (w' * a'))
+    where
+      w' = realToFrac w :: Double
+      a' = realToFrac a :: Double
+  {-# INLINE addValue #-}
+
+instance CalcMean WMeanKBN where
+  calcMean (WMeanKBN (kbn -> w) (kbn -> s))
+    | w <= 0    = throwM $ EmptySample "Data.Monoid.Statistics.Numeric.WMeanKBN: calcMean"
+    | otherwise = return (s / w)
+
+
+----------------------------------------------------------------
+
+-- | This is algorithm for estimation of mean and variance of sample
+--   which uses modified Welford algorithm. It uses KBN summation and
+--   provides approximately 2 additional decimal digits
+data VarWelfordKBN = VarWelfordKBN
+  {-# UNPACK #-} !Int    --  Number of elements in the sample
+  {-# UNPACK #-} !KBNSum -- Current sum of elements of sample
+  {-# UNPACK #-} !KBNSum -- Current sum of squares of deviations from current mean
+
+asVarWelfordKBN :: VarWelfordKBN -> VarWelfordKBN
+asVarWelfordKBN = id
+
+
+-- | Incremental algorithms for calculation the standard deviation [Chan1979].
 data Variance = Variance {-# UNPACK #-} !Int    --  Number of elements in the sample
                          {-# UNPACK #-} !Double -- Current sum of elements of sample
                          {-# UNPACK #-} !Double -- Current sum of squares of deviations from current mean
-                deriving (Show,Eq,Typeable)
+  deriving stock (Show,Eq,Typeable)
 
--- | Fix type of monoid
+-- | Type restricted 'id '
 asVariance :: Variance -> Variance
 asVariance = id
-{-# INLINE asVariance #-}
 
--- | Using parallel algorithm from:
--- 
--- Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1979),
--- Updating Formulae and a Pairwise Algorithm for Computing Sample
--- Variances., Technical Report STAN-CS-79-773, Department of
--- Computer Science, Stanford University. Page 4.
--- 
--- <ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf>
---
-instance Monoid Variance where
-  mempty = Variance 0 0 0
-  mappend !(Variance n1 ta sa) !(Variance n2 tb sb) = Variance (n1+n2) (ta+tb) sumsq
+instance Semigroup Variance where
+  Variance n1 ta sa <> Variance n2 tb sb
+    = Variance (n1+n2) (ta+tb) sumsq
     where
       na = fromIntegral n1
       nb = fromIntegral n2
       nom = sqr (ta * nb - tb * na)
-      sumsq
-        | n1 == 0 || n2 == 0 = sa + sb  -- because either sa or sb should be 0
-        | otherwise          = sa + sb + nom / ((na + nb) * na * nb)
-  {-# INLINE mempty #-}
-  {-# INLINE mappend #-}
+      sumsq | n1 == 0   = sb
+            | n2 == 0   = sa
+            | otherwise = sa + sb + nom / ((na + nb) * na * nb)
 
+instance Monoid Variance where
+  mempty  = Variance 0 0 0
+  mappend = (<>)
+
 instance Real a => StatMonoid Variance a where
-  -- Can be implemented directly as in Welford-Knuth algorithm.
-  pappend !x !s = s `mappend` (Variance 1 (realToFrac x) 0)
-  {-# INLINE pappend #-}
+  addValue (Variance 0 _ _) x = singletonMonoid x
+  addValue (Variance n t s) (realToFrac -> x)
+    = Variance (n + 1) (t + x) (s + sqr (t  - n' * x) / (n' * (n'+1)))
+    where
+      n' = fromIntegral n
+  {-# INLINE addValue #-}
+  singletonMonoid x = Variance 1 (realToFrac x) 0
+  {-# INLINE singletonMonoid #-}
 
 instance CalcCount Variance where
   calcCount (Variance n _ _) = n
-  {-# INLINE calcCount #-}
+
 instance CalcMean Variance where
-  calcMean (Variance n t _) = t / fromIntegral n
-  {-# INLINE calcMean #-}
+  calcMean (Variance 0 _ _) = throwM $ EmptySample "Data.Monoid.Statistics.Numeric.Variance: calcMean"
+  calcMean (Variance n s _) = return (s / fromIntegral n)
+
 instance CalcVariance Variance where
-  calcVariance (Variance n _ s) = s / fromIntegral n
-  calcVarianceUnbiased (Variance n _ s) = s / fromIntegral (n-1)
-  {-# INLINE calcVariance         #-}
-  {-# INLINE calcVarianceUnbiased #-}
+  calcVariance (Variance n _ s)
+    | n < 2     = throwM $ InvalidSample
+                    "Data.Monoid.Statistics.Numeric.Variance: calcVariance"
+                    "Need at least 2 elements"
+    | otherwise = return $! s / fromIntegral (n - 1)
+  calcVarianceML (Variance n _ s)
+    | n < 1     = throwM $ InvalidSample
+                    "Data.Monoid.Statistics.Numeric.Variance: calcVarianceML"
+                    "Need at least 1 element"
+    | otherwise = return $! s / fromIntegral n
 
 
 
+----------------------------------------------------------------
 
+-- | Calculate minimum of sample
+newtype Min a = Min { calcMin :: Maybe a }
+  deriving stock (Show,Eq,Ord,Data,Generic)
 
--- | Calculate minimum of sample. For empty sample returns NaN. Any
--- NaN encountedred will be ignored. 
-newtype Min = Min { calcMin :: Double }
-              deriving (Show,Eq,Ord,Typeable)
+instance Ord a => Semigroup (Min a) where
+  Min (Just a) <> Min (Just b) = Min (Just $! min a b)
+  Min a        <> Min Nothing  = Min a
+  Min Nothing  <> Min b        = Min b
 
+instance Ord a => Monoid (Min a) where
+  mempty  = Min Nothing
+  mappend = (<>)
+
+instance (Ord a, a ~ a') => StatMonoid (Min a) a' where
+  singletonMonoid a = Min (Just a)
+
+
+----------------------------------------------------------------
+
+-- | Calculate maximum of sample
+newtype Max a = Max { calcMax :: Maybe a }
+  deriving stock (Show,Eq,Ord,Data,Generic)
+
+instance Ord a => Semigroup (Max a) where
+  Max (Just a) <> Max (Just b) = Max (Just $! max a b)
+  Max a        <> Max Nothing  = Max a
+  Max Nothing  <> Max b        = Max b
+
+instance Ord a => Monoid (Max a) where
+  mempty  = Max Nothing
+  mappend = (<>)
+
+instance (Ord a, a ~ a') => StatMonoid (Max a) a' where
+  singletonMonoid a = Max (Just a)
+
+
+----------------------------------------------------------------
+
+-- | Calculate minimum of sample of Doubles. For empty sample returns NaN. Any
+--   NaN encountered will be ignored.
+newtype MinD = MinD { calcMinD :: Double }
+  deriving stock (Show,Data,Generic)
+
+instance Eq MinD where
+  MinD a == MinD b
+    | isNaN a && isNaN b = True
+    | otherwise          = a == b
+
+instance Semigroup MinD where
+  MinD x <> MinD y
+    | isNaN x   = MinD y
+    | isNaN y   = MinD x
+    | otherwise = MinD (min x y)
+
 -- N.B. forall (x :: Double) (x <= NaN) == False
-instance Monoid Min where
-  mempty = Min (0/0)
-  mappend !(Min x) !(Min y) 
-    | isNaN x   = Min y
-    | isNaN y   = Min x
-    | otherwise = Min (min x y)
-  {-# INLINE mempty  #-}
-  {-# INLINE mappend #-}  
+instance Monoid MinD where
+  mempty  = MinD (0/0)
+  mappend = (<>)
 
-instance StatMonoid Min Double where
-  pappend !x m = mappend (Min x) m
-  {-# INLINE pappend #-}
+instance a ~ Double => StatMonoid MinD a where
+  singletonMonoid = MinD
 
+
+
 -- | Calculate maximum of sample. For empty sample returns NaN. Any
--- NaN encountedred will be ignored. 
-newtype Max = Max { calcMax :: Double }
-              deriving (Show,Eq,Ord,Typeable)
+--   NaN encountered will be ignored.
+newtype MaxD = MaxD { calcMaxD :: Double }
+  deriving stock (Show,Data,Generic)
 
-instance Monoid Max where
-  mempty = Max (0/0)
-  mappend !(Max x) !(Max y) 
-    | isNaN x   = Max y
-    | isNaN y   = Max x
-    | otherwise = Max (max x y)
-  {-# INLINE mempty  #-}
-  {-# INLINE mappend #-}  
+instance Eq MaxD where
+  MaxD a == MaxD b
+    | isNaN a && isNaN b = True
+    | otherwise          = a == b
 
-instance StatMonoid Max Double where
-  pappend !x m = mappend (Max x) m
-  {-# INLINE pappend #-}
+instance Semigroup MaxD where
+  MaxD x <> MaxD y
+    | isNaN x   = MaxD y
+    | isNaN y   = MaxD x
+    | otherwise = MaxD (max x y)
 
+instance Monoid MaxD where
+  mempty  = MaxD (0/0)
+  mappend = (<>)
 
+instance a ~ Double => StatMonoid MaxD a where
+  singletonMonoid = MaxD
 
 
 ----------------------------------------------------------------
--- Ad-hoc type class
-----------------------------------------------------------------
-  
--- $accessors
---
--- Monoids 'Count', 'Mean' and 'Variance' form some kind of tower.
--- Every successive monoid can calculate every statistics previous
--- monoids can. So to avoid replicating accessors for each statistics
--- a set of ad-hoc type classes was added. 
---
--- This approach have deficiency. It becomes to infer type of monoidal
--- accumulator from accessor function so following expression will be
--- rejected:
--- 
--- > calcCount $ evalStatistics xs
---
--- Indeed type of accumulator is:
---
--- > forall a . (StatMonoid a, CalcMean a) => a
---
--- Therefore it must be fixed by adding explicit type annotation. For
--- example:
---
--- > calcMean (evalStatistics xs :: Mean)
 
-  
+-- | Accumulator for binomial trials.
+data BinomAcc = BinomAcc { binomAccSuccess :: !Int
+                         , binomAccTotal   :: !Int
+                         }
+  deriving stock (Show,Eq,Ord,Data,Generic)
 
--- | Statistics which could count number of elements in the sample
-class CalcCount m where
-  -- | Number of elements in sample
-  calcCount :: m -> Int
+-- | Type restricted 'id'
+asBinomAcc :: BinomAcc -> BinomAcc
+asBinomAcc = id
 
--- | Statistics which could estimate mean of sample
-class CalcMean m where
-  -- | Calculate esimate of mean of a sample
-  calcMean :: m -> Double
-  
--- | Statistics which could estimate variance of sample
-class CalcVariance m where
-  -- | Calculate biased estimate of variance
-  calcVariance         :: m -> Double
-  -- | Calculate unbiased estimate of the variance, where the
-  --   denominator is $n-1$.
-  calcVarianceUnbiased :: m -> Double
+instance Semigroup BinomAcc where
+  BinomAcc n1 m1 <> BinomAcc n2 m2 = BinomAcc (n1+n2) (m1+m2)
 
--- | Calculate sample standard deviation (biased estimator, $s$, where
---   the denominator is $n-1$).
-calcStddev :: CalcVariance m => m -> Double
-calcStddev = sqrt . calcVariance
-{-# INLINE calcStddev #-}
+instance Monoid BinomAcc where
+  mempty  = BinomAcc 0 0
+  mappend = (<>)
 
--- | Calculate standard deviation of the sample
--- (unbiased estimator, $\sigma$, where the denominator is $n$).
-calcStddevUnbiased :: CalcVariance m => m -> Double
-calcStddevUnbiased = sqrt . calcVarianceUnbiased
-{-# INLINE calcStddevUnbiased #-}
+instance StatMonoid BinomAcc Bool where
+  addValue (BinomAcc nS nT) True  = BinomAcc (nS+1) (nT+1)
+  addValue (BinomAcc nS nT) False = BinomAcc  nS    (nT+1)
 
 
+-- | Value @a@ weighted by weight @w@
+data Weighted w a = Weighted w a
+  deriving stock (Show,Eq,Ord,Data,Generic,Functor,Foldable,Traversable)
 
+instance Bifunctor Weighted where
+  first  f   (Weighted w a) = Weighted (f w) a
+  second f   (Weighted w a) = Weighted w (f a)
+  bimap  f g (Weighted w a) =Weighted (f w) (g a)
+  {-# INLINE first  #-}
+  {-# INLINE second #-}
+  {-# INLINE bimap  #-}
+
+
 ----------------------------------------------------------------
 -- Helpers
 ----------------------------------------------------------------
- 
+
 sqr :: Double -> Double
 sqr x = x * x
 {-# INLINE sqr #-}
+
+
+----------------------------------------------------------------
+-- Unboxed instances
+----------------------------------------------------------------
+
+derivingUnbox "CountG"
+  [t| forall a. Unbox a => CountG a -> a |]
+  [| calcCountN |]
+  [| CountG     |]
+
+derivingUnbox "MeanNaive"
+  [t| MeanNaive -> (Int,Double) |]
+  [| \(MeanNaive a b) -> (a,b)  |]
+  [| \(a,b) -> MeanNaive a b    |]
+
+derivingUnbox "MeanKBN"
+  [t| MeanKBN -> (Int,Double,Double)      |]
+  [| \(MeanKBN a (KBNSum b c)) -> (a,b,c) |]
+  [| \(a,b,c) -> MeanKBN a (KBNSum b c)   |]
+
+derivingUnbox "WMeanNaive"
+  [t| WMeanNaive -> (Double,Double) |]
+  [| \(WMeanNaive a b) -> (a,b)     |]
+  [| \(a,b) -> WMeanNaive a b       |]
+
+derivingUnbox "WMeanKBN"
+  [t| WMeanKBN -> (Double,Double,Double,Double)         |]
+  [| \(WMeanKBN (KBNSum a b) (KBNSum c d)) -> (a,b,c,d) |]
+  [| \(a,b,c,d) -> WMeanKBN (KBNSum a b) (KBNSum c d)   |]
+
+derivingUnbox "Variance"
+  [t| Variance -> (Int,Double,Double) |]
+  [| \(Variance a b c) -> (a,b,c)     |]
+  [| \(a,b,c) -> Variance a b c       |]
+
+derivingUnbox "MinD"
+  [t| MinD -> Double |]
+  [| calcMinD |]
+  [| MinD     |]
+
+derivingUnbox "MaxD"
+  [t| MaxD -> Double |]
+  [| calcMaxD |]
+  [| MaxD     |]
+
+derivingUnbox "Weighted"
+  [t| forall w a. (Unbox w, Unbox a) => Weighted w a -> (w,a) |]
+  [| \(Weighted w a) -> (w,a) |]
+  [| \(w,a) -> Weighted w a   |]
+
+derivingUnbox "BinomAcc"
+  [t| BinomAcc -> (Int,Int)   |]
+  [| \(BinomAcc k n) -> (k,n) |]
+  [| \(k,n) -> BinomAcc k n   |]
+
+instance VU.IsoUnbox CountW (Double,Double) where
+  toURepr (CountW w w2) = (w,w2)
+  fromURepr (w,w2) = CountW w w2
+  {-# INLINE toURepr   #-}
+  {-# INLINE fromURepr #-}
+newtype instance VU.MVector s CountW = MV_CountW (VU.MVector s (Double,Double))
+newtype instance VU.Vector    CountW = V_CountW  (VU.Vector    (Double,Double))
+deriving via (CountW `VU.As` (Double,Double)) instance VGM.MVector VU.MVector CountW
+deriving via (CountW `VU.As` (Double,Double)) instance VG.Vector   VU.Vector  CountW
+instance VU.Unbox CountW
+
+-- $references
+--
+-- * [Welford1962] Welford, B.P. (1962) Note on a method for
+--   calculating corrected sums of squares and
+--   products. /Technometrics/
+--   4(3):419-420. <http://www.jstor.org/stable/1266577>
+--
+-- * [Chan1979] Chan, Tony F.; Golub, Gene H.; LeVeque, Randall
+--   J. (1979), Updating Formulae and a Pairwise Algorithm for
+--   Computing Sample Variances., Technical Report STAN-CS-79-773,
+--   Department of Computer Science, Stanford University. Page 4.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,7 @@
+# monoid-statistics parallelizable constant space estimators
+
+[![Build Status](https://travis-ci.org/Shimuuar/monoid-statistics.png?branch=master)](https://travis-ci.org/Shimuuar/monoid-statistics)
+
+Monoids for calculation of statistics of sample. This approach allows to
+calculate many statistics in one pass over data and possibility to parallelize
+calculations. However not all statistics could be calculated this way.
diff --git a/bench/Main.hs b/bench/Main.hs
new file mode 100644
--- /dev/null
+++ b/bench/Main.hs
@@ -0,0 +1,68 @@
+-- |
+module Main where
+
+import Control.Monad
+import Control.Monad.ST (runST)
+
+import qualified Data.Vector.Unboxed as U
+
+import Criterion.Main
+import System.Random.MWC
+--
+import Numeric.Sum
+import Data.Monoid.Statistics
+import Data.Monoid.Statistics.Extra
+
+
+sampleD10,sampleD100,sampleD1000 :: U.Vector Double
+sampleD10   = runST $ U.replicateM 10   . uniform =<< create
+sampleD100  = runST $ U.replicateM 100  . uniform =<< create
+sampleD1000 = runST $ U.replicateM 1000 . uniform =<< create
+
+sampleI10,sampleI100,sampleI1000 :: U.Vector Int
+sampleI10   = runST $ U.replicateM 10   . uniform =<< create
+sampleI100  = runST $ U.replicateM 100  . uniform =<< create
+sampleI1000 = runST $ U.replicateM 1000 . uniform =<< create
+
+main :: IO ()
+main = defaultMain
+  [ bgroup "Count"
+      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> Count) sampleD10
+      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> Count) sampleD100
+      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> Count) sampleD1000
+      ]
+  , bgroup "KBN"
+      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> KBNSum) sampleD10
+      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> KBNSum) sampleD100
+      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> KBNSum) sampleD1000
+      ]
+  , bgroup "Kahan"
+      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> KahanSum) sampleD10
+      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> KahanSum) sampleD100
+      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> KahanSum) sampleD1000
+      ]
+  , bgroup "MeanKBN"
+      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> MeanKBN) sampleD10
+      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> MeanKBN) sampleD100
+      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> MeanKBN) sampleD1000
+      , bench "10   I" $ whnf (reduceSampleVec :: U.Vector Int    -> MeanKBN) sampleI10
+      , bench "100  I" $ whnf (reduceSampleVec :: U.Vector Int    -> MeanKBN) sampleI100
+      , bench "1000 I" $ whnf (reduceSampleVec :: U.Vector Int    -> MeanKBN) sampleI1000
+      ]
+  , bgroup "WelfordMean"
+      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> WelfordMean) sampleD10
+      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> WelfordMean) sampleD100
+      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> WelfordMean) sampleD1000
+      , bench "10   I" $ whnf (reduceSampleVec :: U.Vector Int    -> WelfordMean) sampleI10
+      , bench "100  I" $ whnf (reduceSampleVec :: U.Vector Int    -> WelfordMean) sampleI100
+      , bench "1000 I" $ whnf (reduceSampleVec :: U.Vector Int    -> WelfordMean) sampleI1000
+      ]
+  , bgroup "Variance"
+      [ bench "10   D" $ whnf (reduceSampleVec :: U.Vector Double -> Variance) sampleD10
+      , bench "100  D" $ whnf (reduceSampleVec :: U.Vector Double -> Variance) sampleD100
+      , bench "1000 D" $ whnf (reduceSampleVec :: U.Vector Double -> Variance) sampleD1000
+      , bench "10   I" $ whnf (reduceSampleVec :: U.Vector Int    -> Variance) sampleI10
+      , bench "100  I" $ whnf (reduceSampleVec :: U.Vector Int    -> Variance) sampleI100
+      , bench "1000 I" $ whnf (reduceSampleVec :: U.Vector Int    -> Variance) sampleI1000
+      ]
+  ]
diff --git a/monoid-statistics.cabal b/monoid-statistics.cabal
--- a/monoid-statistics.cabal
+++ b/monoid-statistics.cabal
@@ -1,13 +1,12 @@
-
-
 Name:           monoid-statistics
-Version:        0.3.1
-Cabal-Version:  >= 1.6
+Version:        1.1.5
+Cabal-Version:  >= 1.10
 License:        BSD3
 License-File:   LICENSE
 Author:         Alexey Khudyakov <alexey.skladnoy@gmail.com>
 Maintainer:     Alexey Khudyakov <alexey.skladnoy@gmail.com>
-Homepage:       https://bitbucket.org/Shimuuar/monoid-statistics
+Homepage:       https://github.com/Shimuuar/monoid-statistics
+Bug-reports:    https://github.com/Shimuuar/monoid-statistics/issues
 Category:       Statistics
 Build-Type:     Simple
 Synopsis:       
@@ -17,20 +16,85 @@
   allows to calculate many statistics in one pass over data and
   possibility to parallelize calculations. However not all statistics 
   could be calculated this way.
-  .
-  This packages is quite similar to monoids package but limited to
-  calculation on statistics. In particular it makes use of
-  commutatitvity of statistical monoids.
-  .
-  Changes:
-  .
-  * 0.3.1 Better documentation; Fix in Min/Max monoids
 
+Extra-Source-Files:
+  Changelog.md
+
+tested-with:
+    GHC ==8.6.5
+     || ==8.8.4
+     || ==8.10.7
+     || ==9.0.1
+     || ==9.2.1
+
+extra-source-files:
+  README.md
+
 source-repository head
-  type:     hg
-  location: http://bitbucket.org/Shimuuar/monoid-statistics
+  type:     git
+  location: https://github.com/Shimuuar/monoid-statistics
 
 Library
-  Build-Depends:   base >=3 && <5
+  default-language: Haskell2010
+  ghc-options:      -Wall -O2
+  --
+  Build-Depends:    base            >=4.12  && <5
+                  , exceptions      >=0.10
+                  , vector          >=0.13 && <1
+                  , vector-th-unbox >=0.2.1.6
+                  , math-functions  >=0.3
+  --
   Exposed-modules: Data.Monoid.Statistics
+                   Data.Monoid.Statistics.Class
                    Data.Monoid.Statistics.Numeric
+                   Data.Monoid.Statistics.Extra
+
+test-suite monoid-statistics-tests
+  default-language: Haskell2010
+  type:             exitcode-stdio-1.0
+  ghc-options:      -Wall -threaded
+  -- Tests for math-functions' Sum require SSE2 on i686 to pass
+  -- (because of excess precision)
+  if arch(i386)
+    ghc-options:  -msse2
+  hs-source-dirs: tests
+  Main-is:        Main.hs
+  Other-Modules:
+  Build-Depends: monoid-statistics
+               , base             >=4.9 && <5
+               , math-functions   >=0.3
+               , tasty            >=0.11
+               , tasty-quickcheck >=0.9
+               , tasty-hunit
+               , tasty-expected-failure
+               , QuickCheck
+
+test-suite monoid-statistics-doctest
+  if impl(ghcjs)
+    buildable: False
+  -- It seems GHC 9.0 & 9.2 chokes to death on examples with deriving via
+  if impl(ghc >= 9.0) && impl(ghc < 9.1)
+    buildable: False
+  if impl(ghc >= 9.2) && impl(ghc < 9.3)
+    buildable: False
+  type:             exitcode-stdio-1.0
+  main-is:          doctests.hs
+  hs-source-dirs:   tests
+  default-language: Haskell2010
+  build-depends:
+        base                >=4.9  && <5
+      , doctest             >=0.15 && <0.23
+      , monoid-statistics   -any
+
+benchmark monoid-stat-bench
+  default-language: Haskell2010
+  type:             exitcode-stdio-1.0
+  ghc-options:      -Wall -O2
+  hs-source-dirs:   bench
+  Main-is:          Main.hs
+  Build-Depends: monoid-statistics
+               , base           >=4.9  && <5
+               , vector         >=0.11 && <1
+               , math-functions >=0.3
+               , mwc-random     >=0.13
+               , criterion
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -0,0 +1,328 @@
+{-# LANGUAGE AllowAmbiguousTypes #-}
+{-# LANGUAGE FlexibleContexts    #-}
+{-# LANGUAGE FlexibleInstances   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications    #-}
+--
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+module Main (main) where
+import Data.Typeable
+import Numeric.Sum
+import Test.Tasty
+import Test.Tasty.QuickCheck
+import Test.Tasty.HUnit
+import Test.Tasty.ExpectedFailure (ignoreTest)
+
+import Data.Monoid.Statistics
+import Data.Monoid.Statistics.Extra
+
+
+----------------------------------------------------------------
+-- Properties
+----------------------------------------------------------------
+
+class MonoidProperty m where
+  isAssociative, isCommutative, isMemptyDistribute, isMemptyNeutral :: Bool
+  isMemptyNeutral    = True
+  isMemptyDistribute = True
+  isAssociative      = True
+  isCommutative      = True
+
+instance {-# OVERLAPPABLE #-} MonoidProperty m
+instance MonoidProperty MeanNaive where
+  isAssociative = False
+instance MonoidProperty WelfordMean where
+  isAssociative      = False
+  isMemptyDistribute = False
+instance MonoidProperty MeanKahan where
+  isAssociative      = False
+  isCommutative      = False
+  isMemptyDistribute = False
+instance MonoidProperty MeanKBN where
+  isAssociative = False
+  isCommutative = False
+instance MonoidProperty Variance where
+  isAssociative = False
+instance MonoidProperty WMeanNaive where
+  isAssociative   = False
+instance MonoidProperty WMeanKBN where
+  isMemptyNeutral = False
+  isAssociative   = False
+  isCommutative   = False
+
+p_memptyIsNeutral
+  :: forall m. (Monoid m, MonoidProperty m, Arbitrary m, Show m, Eq m)
+  => TestTree
+p_memptyIsNeutral
+  = (if isMemptyNeutral @m then id else ignoreTest)
+  $ testProperty "mempty is neutral" $ \(m :: m) ->
+    counterexample ("m <> mempty = " ++ show (m <> mempty))
+  $ counterexample ("mempty <> m = " ++ show (mempty <> m))
+  $  (m <> mempty) == m
+  && (mempty <> m) == m
+
+p_associativity
+  :: forall m. (MonoidProperty m, Monoid m, Arbitrary m, Show m, Eq m)
+  => TestTree
+p_associativity
+  = (if isAssociative @m then id else ignoreTest)
+  $ testProperty "associativity" $ \(a :: m) b c ->
+    let val1 = (a <> b) <> c
+        val2 = a <> (b <> c)
+    in counterexample ("left : " ++ show val1)
+     $ counterexample ("right: " ++ show val2)
+     $ val1 == val2
+
+p_commutativity
+  :: forall m. (Monoid m, MonoidProperty m, Arbitrary m, Show m, Eq m)
+  => TestTree
+p_commutativity
+  = (if isCommutative @m then id else ignoreTest)
+  $ testProperty "commutativity" $ \(a :: m) b ->
+    let val1 = a <> b
+        val2 = b <> a
+    in counterexample ("a <> b = " ++ show val1)
+     $ counterexample ("b <> a = " ++ show val2)
+     $ val1 == val2
+
+p_addValue1
+  :: forall m a. ( StatMonoid m a
+                 , Eq m
+                 , Arbitrary a, Show a)
+  => TestTree
+p_addValue1
+  = testProperty "addValue x mempty == singletonMonoid" $ \(a :: a) ->
+    singletonMonoid a == addValue (mempty :: m) a
+
+
+p_addValue2
+  :: forall m a. ( MonoidProperty m, StatMonoid m a
+                 , Show m, Eq m
+                 , Arbitrary a, Show a)
+  => TestTree
+p_addValue2
+  = (if isMemptyDistribute @m then id else ignoreTest)
+  $ testProperty "addValue (addValue m x) y = addValue 0 x <> addValue 0 y" $ \(x :: a) (y :: a) ->
+    let val1 = addValue (addValue mempty y) x
+        val2 = (addValue mempty x <> addValue (mempty :: m) y)
+    in counterexample ("left : " ++ show val1)
+     $ counterexample ("right: " ++ show val2)
+     $ val1 == val2
+
+
+
+----------------------------------------------------------------
+
+testMonoid
+  :: forall m a.
+     ( StatMonoid m a, MonoidProperty m
+     , Typeable a, Typeable m, Arbitrary a, Arbitrary m, Show a, Show m, Eq m)
+  => [TestTree] -> TestTree
+testMonoid tests
+  = testGroup (show (typeOf (undefined :: m)) ++ " <= " ++ show (typeOf (undefined :: a)))
+  $ [ p_memptyIsNeutral @m
+    , p_associativity   @m
+    , p_commutativity   @m
+    , p_addValue1       @m @a
+    , p_addValue2       @m @a
+    ]
+ ++ tests
+
+testMeanMonoid
+  :: forall m.
+     ( StatMonoid m Double, CalcMean m, CalcCount m, MonoidProperty m
+     , Typeable m, Arbitrary m, Show m, Eq m)
+  => [TestTree] -> TestTree
+testMeanMonoid tests
+  = testMonoid @m @Double
+  $ [ testCase "Count" $ do
+        let m = reduceSample @m testSample
+        testSampleCount     @=? calcCount m
+    , testCase "Mean" $ do
+        let m = reduceSample @m testSample
+        Just testSampleMean @=? calcMean  m
+    , testCase "Mean (empty sample)" $ do
+        let m = reduceSample @m @Double []
+        Nothing @=? calcMean m
+    ] ++ tests
+
+testVarianceMonoid
+  :: forall m.
+     ( StatMonoid m Double, CalcVariance m, CalcMean m, CalcCount m, MonoidProperty m
+     , Typeable m, Arbitrary m, Show m, Eq m)
+  => [TestTree] -> TestTree
+testVarianceMonoid tests
+  = testMeanMonoid @m
+  $ [ testCase "Variance (unbiased)" $ do
+        let m = reduceSample @m testSample
+        Just testSampleVariance   @=? calcVariance m
+    , testCase "Variance (ML)" $ do
+        let m = reduceSample @m testSample
+        Just testSampleVarianceML @=? calcVarianceML m
+    ] ++ tests
+
+testWMeanMonoid
+  :: forall m.
+     ( StatMonoid m (Weighted Double Double), CalcMean m, MonoidProperty m
+     , Typeable m, Arbitrary m, Show m, Eq m)
+  => [TestTree] -> TestTree
+testWMeanMonoid tests
+  = testMonoid @m @(Weighted Double Double)
+  $ [ testCase "Mean" $ do
+        let m = reduceSample @m testWSample
+        Just testWSampleMean @=? calcMean  m
+    ] ++ tests
+
+
+main :: IO ()
+main = defaultMain $ testGroup "monoid-statistics"
+  [ testMonoid @(CountG Int) @Int
+    [ testCase "CountG"  $ let xs = "acbdef"
+                               n  = reduceSample xs :: Count
+                           in length xs @=? calcCount n
+    ]
+  , testMonoid @(Min Int)    @Int
+    [ testCase "Min []"  $ let xs = []
+                               n  = reduceSample xs :: Min Int
+                           in Nothing @=? calcMin n
+    , testCase "Min"     $ let xs = [1..10]
+                               n  = reduceSample xs :: Min Int
+                           in Just (minimum xs) @=? calcMin n
+    ]
+  , testMonoid @(Max Int)    @Int
+    [ testCase "Max []"  $ let xs = []
+                               n  = reduceSample xs :: Max Int
+                           in Nothing @=? calcMax n
+    , testCase "Max"     $ let xs = [1..10]
+                               n  = reduceSample xs :: Max Int
+                           in Just (maximum xs) @=? calcMax n
+    ]
+  , testMonoid @MinD         @Double
+    [ testCase "MinD"    $ let xs = [1..10]
+                               n  = reduceSample xs :: MinD
+                           in minimum xs @=? calcMinD n
+    ]
+  , testMonoid @MaxD         @Double
+    [ testCase "MaxD"    $ let xs = [1..10]
+                               n  = reduceSample xs :: MaxD
+                           in maximum xs @=? calcMaxD n
+
+    ]
+  , testMonoid @BinomAcc     @Bool   []
+    -- Numeric accumulators
+  , testMeanMonoid     @MeanNaive   []
+  , testMeanMonoid     @WelfordMean []
+  , testMeanMonoid     @MeanKahan   []
+  , testMeanMonoid     @MeanKBN     []
+  , testWMeanMonoid    @WMeanNaive  []
+  , testWMeanMonoid    @WMeanKBN    []
+  , testVarianceMonoid @Variance    []
+  ]
+
+-- | Test sample for which we could compute statistics exactly, and
+--   any reasonable algorithm should be able to return exact answer as
+--   well
+testSample :: [Double]
+testSample = [1..10]
+
+testWSample :: [Weighted Double Double]
+testWSample = [Weighted x x | x <- [1..10]]
+
+testSampleCount :: Int
+testSampleCount = length testSample
+
+testSampleMean,testWSampleMean :: Double
+testSampleMean  = 5.5
+testWSampleMean = 7.0
+
+testSampleVariance,testSampleVarianceML :: Double
+testSampleVariance   = 9.166666666666666
+testSampleVarianceML = 8.25
+
+----------------------------------------------------------------
+
+instance (Arbitrary a, Num a, Ord a) => Arbitrary (CountG a) where
+  arbitrary = do
+    NonNegative n <- arbitrary
+    return (CountG n)
+
+instance (Arbitrary a) => Arbitrary (Max a) where
+  arbitrary = fmap Max arbitrary
+
+instance (Arbitrary a) => Arbitrary (Min a) where
+  arbitrary = fmap Min arbitrary
+
+instance Arbitrary MinD where
+  arbitrary = frequency [ (1, return mempty)
+                        , (4, fmap MinD arbitrary)
+                        ]
+
+instance Arbitrary MaxD where
+  arbitrary = frequency [ (1, return mempty)
+                        , (4, fmap MaxD arbitrary)
+                        ]
+
+instance Arbitrary BinomAcc where
+  arbitrary = do
+    NonNegative nSucc <- arbitrary
+    NonNegative nFail <- arbitrary
+    return $ BinomAcc nSucc (nFail + nSucc)
+
+instance Arbitrary MeanNaive where
+  arbitrary = arbitrary >>= \x -> case x of
+    NonNegative 0 -> return mempty
+    NonNegative n -> do m <- arbitrary
+                        return (MeanNaive n m)
+
+instance Arbitrary WelfordMean where
+  arbitrary = arbitrary >>= \x -> case x of
+    NonNegative 0 -> return mempty
+    NonNegative n -> do m <- arbitrary
+                        return (WelfordMean n m)
+
+instance Arbitrary MeanKahan where
+  arbitrary = do
+    n <- arbitrary
+    s <- arbitraryKBN n
+    return $! MeanKahan (getNonNegative n) s
+
+instance Arbitrary MeanKBN where
+  arbitrary = do
+    n <- arbitrary
+    s <- arbitraryKBN n
+    return $! MeanKBN (getNonNegative n) s
+
+instance Arbitrary WMeanKBN where
+  arbitrary = do
+    n            <- arbitrary
+    KBNSum w1 w2 <- arbitraryKBN n
+    s            <- arbitraryKBN n
+    return $! WMeanKBN (KBNSum (abs w1) w2) s
+
+instance Arbitrary WMeanNaive where
+  arbitrary = do
+    NonNegative w <- arbitrary
+    s             <- arbitrary
+    return $! WMeanNaive w s
+
+instance Arbitrary Variance where
+  arbitrary = arbitrary >>= \x -> case x of
+    NonNegative 0 -> return mempty
+    NonNegative n -> do
+      m             <- arbitrary
+      NonNegative s <- arbitrary
+      return $ Variance n m s
+
+instance (Arbitrary a, Arbitrary w) => Arbitrary (Weighted w a) where
+  arbitrary = Weighted <$> arbitrary <*> arbitrary
+
+arbitraryKBN :: Summation a => NonNegative Int -> Gen a
+arbitraryKBN (NonNegative 0) = return zero
+arbitraryKBN (NonNegative 1) = do
+  x1 <- arbitrary
+  return $! zero `add` x1
+arbitraryKBN _ = do
+  x1 <- arbitrary
+  x2 <- arbitrary
+  x3 <- arbitrary
+  return $! ((zero `add` x1) `add` x2) `add` x3
diff --git a/tests/doctests.hs b/tests/doctests.hs
new file mode 100644
--- /dev/null
+++ b/tests/doctests.hs
@@ -0,0 +1,6 @@
+module Main where
+
+import Test.DocTest (doctest)
+
+main :: IO ()
+main = doctest ["Data"]