monoid-statistics 1.0.1.0 → 1.1.0

raw patch · 9 files changed

+1105/−389 lines, 9 filesdep +criteriondep +doctestdep +exceptionsdep ~basedep ~math-functionsPVP ok

version bump matches the API change (PVP)

Dependencies added: criterion, doctest, exceptions, mwc-random, tasty-expected-failure, tasty-hunit

Dependency ranges changed: base, math-functions

API changes (from Hackage documentation)

Files

• Changelog.md +32/−0
• Data/Monoid/Statistics.hs +1/−0
• Data/Monoid/Statistics/Class.hs +307/−34
• Data/Monoid/Statistics/Extra.hs +175/−0
• Data/Monoid/Statistics/Numeric.hs +216/−215
• bench/Main.hs +68/−0
• monoid-statistics.cabal +53/−14
• tests/Main.hs +247/−126
• tests/doctests.hs +6/−0

@@ -0,0 +1,32 @@+# Changes in 1.1.0.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.

@@ -12,3 +12,4 @@  import Data.Monoid.Statistics.Class import Data.Monoid.Statistics.Numeric+

@@ -1,14 +1,17 @@-{-# LANGUAGE BangPatterns          #-}-{-# LANGUAGE CPP                   #-}-{-# LANGUAGE DeriveDataTypeable    #-}-{-# LANGUAGE DeriveGeneric         #-}-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE RankNTypes            #-}-{-# LANGUAGE TemplateHaskell       #-}-{-# LANGUAGE TypeFamilies          #-}----{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# 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>@@ -17,31 +20,43 @@ -- Stability  : experimental -- module Data.Monoid.Statistics.Class-  ( -- * Type class and helpers+  ( -- * Monoid Type class and helpers     StatMonoid(..)   , reduceSample   , reduceSampleVec+    -- * Ad-hoc type classes for select statistics+    -- $adhoc+  , CalcCount(..)+  , CalcMean(..)+  , HasMean(..)+  , CalcVariance(..)+  , HasVariance(..)+  , CalcViaHas(..)+    -- * Exception handling+  , Partial(..)+  , partial+  , SampleError(..)     -- * Data types   , Pair(..)   ) where -import           Data.Data    (Typeable,Data)-#if MIN_VERSION_base(4,9,0)-import qualified Data.Semigroup as SG (Semigroup(..))-#endif-import           Data.Monoid    (Monoid(..),(<>),Sum(..),Product(..))+import           Control.Exception+import           Control.Monad.Catch (MonadThrow(..))+import           Data.Data           (Typeable,Data)+import           Data.Monoid 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. By definitions---   /statistic/ is some measure of sample which doesn't depend on---   order of elements (for example: mean, sum, number of elements,---   variance, etc).+--   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@@ -49,7 +64,7 @@ --   fold accumulator to get statistic for union of two samples. -- --   Thus for such algorithm we have value which corresponds to empty---   sample, merge function which which corresponds to merging of two+--   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@@ -57,7 +72,7 @@ -- --   Instance must satisfy following laws. If floating point --   arithmetics is used then equality should be understood as---   approximate. +--   approximate. -- --   > 1. addValue (addValue y mempty) x  == addValue mempty x <> addValue mempty y --   > 2. x <> y == y <> x@@ -73,17 +88,66 @@   {-# MINIMAL addValue | singletonMonoid #-}  -- | Calculate statistic over 'Foldable'. It's implemented in terms of---   foldl'.-reduceSample :: (F.Foldable f, StatMonoid m a) => f a -> m+--   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. It's implemented in terms of---   foldl'.-reduceSampleVec :: (G.Vector v a, StatMonoid m a) => v a -> m+-- | 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 @@ -98,8 +162,189 @@   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) #-}+++++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 ---------------------------------------------------------------- @@ -107,16 +352,15 @@ data Pair a b = Pair !a !b               deriving (Show,Eq,Ord,Typeable,Data,Generic) -#if MIN_VERSION_base(4,9,0)-instance (SG.Semigroup a, SG.Semigroup b) => SG.Semigroup (Pair a b) where-  Pair x y <> Pair x' y' = Pair (x SG.<> x') (y SG.<> y')-#endif+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 (Pair x y) (Pair x' y') = Pair (x <> x') (y <> y')+  mappend = (<>)   {-# INLINABLE mempty  #-}-  {-# INLINE mappend #-}+  {-# 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)@@ -124,7 +368,36 @@   {-# 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

@@ -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

@@ -1,26 +1,39 @@ {-# LANGUAGE BangPatterns          #-}-{-# LANGUAGE CPP                   #-} {-# LANGUAGE DeriveDataTypeable    #-}+{-# LANGUAGE DeriveFoldable        #-} {-# LANGUAGE DeriveGeneric         #-}+{-# LANGUAGE DeriveTraversable     #-} {-# LANGUAGE FlexibleContexts      #-} {-# LANGUAGE FlexibleInstances     #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE RankNTypes            #-}+{-# LANGUAGE ScopedTypeVariables   #-} {-# LANGUAGE TemplateHaskell       #-} {-# LANGUAGE TypeFamilies          #-}+{-# LANGUAGE ViewPatterns          #-}+-- |+-- Monoids for calculating various statistics in constant space module Data.Monoid.Statistics.Numeric (     -- * Mean & Variance     -- ** Number of elements     CountG(..)   , Count   , asCount-    -- ** Mean+    -- ** Mean algorithms+    -- ** Default algorithms+  , Mean+  , asMean+  , WMean+  , asWMean+    -- *** Mean+  , MeanNaive(..)+  , asMeanNaive   , MeanKBN(..)   , asMeanKBN-  , WelfordMean(..)-  , asWelfordMean-  , MeanKahan(..)-  , asMeanKahan+    -- *** Weighted mean+  , WMeanNaive(..)+  , asWMeanNaive+  , WMeanKBN(..)+  , asWMeanKBN     -- ** Variance   , Variance(..)   , asVariance@@ -32,27 +45,22 @@     -- * Binomial trials   , BinomAcc(..)   , asBinomAcc-    -- * Accessors-  , CalcCount(..)-  , CalcMean(..)-  , CalcVariance(..)-  , calcStddev-  , calcStddevML+    -- * Rest+  , Weighted(..)     -- * References     -- $references   ) where -import Data.Monoid                  ((<>))-import Data.Monoid.Statistics.Class-#if MIN_VERSION_base(4,9,0)-import qualified Data.Semigroup as SG (Semigroup(..))-#endif+import Control.Monad.Catch          (MonadThrow(..)) import Data.Data                    (Typeable,Data) import Data.Vector.Unboxed          (Unbox) import Data.Vector.Unboxed.Deriving (derivingUnbox) import Numeric.Sum import GHC.Generics                 (Generic) +import Data.Monoid.Statistics.Class++ ---------------------------------------------------------------- -- Statistical monoids ----------------------------------------------------------------@@ -67,144 +75,188 @@ asCount :: CountG a -> CountG a asCount = id -#if MIN_VERSION_base(4,9,0)-instance Integral a => SG.Semigroup (CountG a) where-  (<>) = mappend-#endif+instance Integral a => Semigroup (CountG a) where+  CountG i <> CountG j = CountG (i + j)  instance Integral a => Monoid (CountG a) where-  mempty                      = CountG 0-  CountG i `mappend` CountG j = CountG (i + j)-  {-# INLINE mempty  #-}-  {-# INLINE mappend #-}+  mempty  = CountG 0+  mappend = (<>)  instance (Integral a) => StatMonoid (CountG a) b where   singletonMonoid _            = CountG 1   addValue        (CountG n) _ = CountG (n + 1)-  {-# INLINE singletonMonoid #-}-  {-# INLINE addValue        #-} + instance CalcCount (CountG Int) where   calcCount = calcCountN-  {-# INLINE calcCount #-}    ---------------------------------------------------------------- --- | Incremental calculation of mean. Sum of elements is calculated---   using compensated Kahan summation.-data MeanKahan = MeanKahan !Int !KahanSum+-- | Type alias for currently recommended algorithms for calculation+--   of mean. It should be default choice+type Mean = MeanKBN++asMean :: Mean -> Mean+asMean = id++-- | 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 (Show,Eq,Typeable,Data,Generic) -asMeanKahan :: MeanKahan -> MeanKahan-asMeanKahan = id+asMeanNaive :: MeanNaive -> MeanNaive+asMeanNaive = id -#if MIN_VERSION_base(4,9,0)-instance SG.Semigroup MeanKahan where-  (<>) = mappend-#endif -instance Monoid MeanKahan where-  mempty = MeanKahan 0 mempty-  MeanKahan 0  _  `mappend` m               = m-  m               `mappend` MeanKahan 0  _  = m-  MeanKahan n1 s1 `mappend` MeanKahan n2 s2 = MeanKahan (n1+n2) (s1<>s2)+instance Semigroup MeanNaive where+  MeanNaive 0  _  <> m               = m+  m               <> MeanNaive 0  _  = m+  MeanNaive n1 s1 <> MeanNaive n2 s2 = MeanNaive (n1+n2) (s1 + s2) -instance Real a => StatMonoid MeanKahan a where-  addValue (MeanKahan n m) x = MeanKahan (n+1) (addValue m x)+instance Monoid MeanNaive where+  mempty  = MeanNaive 0 0+  mappend = (<>) -instance CalcCount MeanKahan where-  calcCount (MeanKahan n _) = n-instance CalcMean MeanKahan where-  calcMean (MeanKahan 0 _) = Nothing-  calcMean (MeanKahan n s) = Just (kahan s / fromIntegral n)+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. Sum of elements is calculated---   using Kahan-Babuška-Neumaier summation.-data MeanKBN = MeanKBN !Int !KBNSum+----------------------------------------------------------------++-- | 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 (Show,Eq,Typeable,Data,Generic)  asMeanKBN :: MeanKBN -> MeanKBN asMeanKBN = id -#if MIN_VERSION_base(4,9,0)-instance SG.Semigroup MeanKBN where-  (<>) = mappend-#endif -instance Monoid MeanKBN where-  mempty = MeanKBN 0 mempty-  MeanKBN 0  _  `mappend` m             = m-  m             `mappend` MeanKBN 0  _  = m-  MeanKBN n1 s1 `mappend` MeanKBN n2 s2 = MeanKBN (n1+n2) (s1<>s2)+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 _) = Nothing-  calcMean (MeanKBN n s) = Just (kbn s / fromIntegral n)+  calcMean (MeanKBN 0 _) = throwM $ EmptySample "Data.Monoid.Statistics.Numeric.MeanKBN: calcMean"+  calcMean (MeanKBN n s) = return (kbn s / fromIntegral n)  +---------------------------------------------------------------- --- | Incremental calculation of mean. One of algorithm's advantage is---   protection against double overflow:------   > λ> calcMean $ asMeanKBN     $ reduceSample (replicate 100 1e308)---   > Just NaN---   > λ> calcMean $ asWelfordMean $ reduceSample (replicate 100 1e308)---   > Just 1.0e308------   Algorithm is due to Welford [Welford1962]-data WelfordMean = WelfordMean !Int    -- Number of entries-                               !Double -- Current mean+-- | Incremental calculation of weighed mean.+data WMeanNaive = WMeanNaive+  !Double  -- Weight+  !Double  -- Weighted sum   deriving (Show,Eq,Typeable,Data,Generic) --- | Type restricted 'id'-asWelfordMean :: WelfordMean -> WelfordMean-asWelfordMean = id+asWMeanNaive :: WMeanNaive -> WMeanNaive+asWMeanNaive = id -#if MIN_VERSION_base(4,9,0)-instance SG.Semigroup WelfordMean where-  (<>) = mappend-#endif -instance Monoid WelfordMean where-  mempty = WelfordMean 0 0-  mappend (WelfordMean 0 _) m = m-  mappend m (WelfordMean 0 _) = m-  mappend (WelfordMean n x) (WelfordMean k y)-    = WelfordMean (n + k) ((x*n' + y*k') / (n' + k'))-    where-      n' = fromIntegral n-      k' = fromIntegral k-  {-# INLINE mempty  #-}-  {-# INLINE mappend #-}+instance Semigroup WMeanNaive where+  WMeanNaive w1 s1 <> WMeanNaive w2 s2 = WMeanNaive (w1 + w2) (s1 + s2) --- | \[ s_n = s_{n-1} + \frac{x_n - s_{n-1}}{n} \]-instance Real a => StatMonoid WelfordMean a where-  addValue (WelfordMean n m) !x-    = WelfordMean n' (m + (realToFrac x - m) / fromIntegral n')+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' = n+1+      w' = realToFrac w+      a' = realToFrac a   {-# INLINE addValue #-} -instance CalcCount WelfordMean where-  calcCount (WelfordMean n _) = n-instance CalcMean WelfordMean where-  calcMean (WelfordMean 0 _) = Nothing-  calcMean (WelfordMean _ m) = Just m+instance CalcMean WMeanNaive where+  calcMean (WMeanNaive w s)+    | w <= 0    = throwM $ EmptySample "Data.Monoid.Statistics.Numeric.WMeanNaive: calcMean"+    | otherwise = return (s / w) +---------------------------------------------------------------- +-- | 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 (Show,Eq,Typeable,Data,Generic) +asWMeanKBN :: WMeanKBN -> WMeanKBN+asWMeanKBN = id+++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)++ ---------------------------------------------------------------- --- | Incremental algorithms for calculation the standard deviation.+-- | 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@@ -213,17 +265,9 @@ -- | Type restricted 'id ' asVariance :: Variance -> Variance asVariance = id-{-# INLINE asVariance #-} -#if MIN_VERSION_base(4,9,0)-instance SG.Semigroup Variance where-  (<>) = mappend-#endif---- | Iterative algorithm for calculation of variance [Chan1979]-instance Monoid Variance where-  mempty = Variance 0 0 0-  mappend (Variance n1 ta sa) (Variance n2 tb sb)+instance Semigroup Variance where+  Variance n1 ta sa <> Variance n2 tb sb     = Variance (n1+n2) (ta+tb) sumsq     where       na = fromIntegral n1@@ -232,10 +276,18 @@       sumsq | n1 == 0   = sb             | n2 == 0   = sa             | otherwise = sa + sb + nom / ((na + nb) * na * nb)-  {-# INLINE mempty #-}-  {-# INLINE mappend #-} +instance Monoid Variance where+  mempty  = Variance 0 0 0+  mappend = (<>)+ instance Real a => StatMonoid Variance a where+  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 #-} @@ -243,17 +295,20 @@   calcCount (Variance n _ _) = n  instance CalcMean Variance where-  calcMean (Variance 0 _ _) = Nothing-  calcMean (Variance n s _) = Just (s / fromIntegral n)+  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)-    | n < 2     = Nothing-    | otherwise = Just $! s / fromIntegral (n - 1)+    | 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     = Nothing-    | otherwise = Just $! s / fromIntegral n-+    | n < 1     = throwM $ InvalidSample+                    "Data.Monoid.Statistics.Numeric.Variance: calcVarianceML"+                    "Need at least 1 element"+    | otherwise = return $! s / fromIntegral n   @@ -263,36 +318,33 @@ newtype Min a = Min { calcMin :: Maybe a }               deriving (Show,Eq,Ord,Typeable,Data,Generic) -#if MIN_VERSION_base(4,9,0)-instance Ord a => SG.Semigroup (Min a) where-  (<>) = mappend-#endif+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-  Min (Just a) `mappend` Min (Just b) = Min (Just $! min a b)-  Min a        `mappend` Min Nothing  = Min a-  Min Nothing  `mappend` Min b        = Min b+  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 (Show,Eq,Ord,Typeable,Data,Generic) -#if MIN_VERSION_base(4,9,0)-instance Ord a => SG.Semigroup (Max a) where-  (<>) = mappend-#endif+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-  Max (Just a) `mappend` Max (Just b) = Max (Just $! min a b)-  Max a        `mappend` Max Nothing  = Max a-  Max Nothing  `mappend` Max b        = Max b+  mempty  = Max Nothing+  mappend = (<>)  instance (Ord a, a ~ a') => StatMonoid (Max a) a' where   singletonMonoid a = Max (Just a)@@ -310,21 +362,17 @@     | isNaN a && isNaN b = True     | otherwise          = a == b -#if MIN_VERSION_base(4,9,0)-instance SG.Semigroup MinD where-  (<>) = mappend-#endif---- N.B. forall (x :: Double) (x <= NaN) == False-instance Monoid MinD where-  mempty = MinD (0/0)-  mappend (MinD x) (MinD y)+instance Semigroup MinD where+  MinD x <> MinD y     | isNaN x   = MinD y     | isNaN y   = MinD x     | otherwise = MinD (min x y)-  {-# INLINE mempty  #-}-  {-# INLINE mappend #-} +-- N.B. forall (x :: Double) (x <= NaN) == False+instance Monoid MinD where+  mempty  = MinD (0/0)+  mappend = (<>)+ instance a ~ Double => StatMonoid MinD a where   singletonMonoid = MinD @@ -340,20 +388,16 @@     | isNaN a && isNaN b = True     | otherwise          = a == b -#if MIN_VERSION_base(4,9,0)-instance SG.Semigroup MaxD where-  (<>) = mappend-#endif--instance Monoid MaxD where-  mempty = MaxD (0/0)-  mappend (MaxD x) (MaxD y)+instance Semigroup MaxD where+  MaxD x <> MaxD y     | isNaN x   = MaxD y     | isNaN y   = MaxD x     | otherwise = MaxD (max x y)-  {-# INLINE mempty  #-}-  {-# INLINE mappend #-} +instance Monoid MaxD where+  mempty  = MaxD (0/0)+  mappend = (<>)+ instance a ~ Double => StatMonoid MaxD a where   singletonMonoid = MaxD @@ -370,64 +414,21 @@ asBinomAcc :: BinomAcc -> BinomAcc asBinomAcc = id -#if MIN_VERSION_base(4,9,0)-instance SG.Semigroup BinomAcc where-  (<>) = mappend-#endif+instance Semigroup BinomAcc where+  BinomAcc n1 m1 <> BinomAcc n2 m2 = BinomAcc (n1+n2) (m1+m2)  instance Monoid BinomAcc where-  mempty = BinomAcc 0 0-  mappend (BinomAcc n1 m1) (BinomAcc n2 m2) = BinomAcc (n1+n2) (m1+m2)+  mempty  = BinomAcc 0 0+  mappend = (<>)  instance StatMonoid BinomAcc Bool where   addValue (BinomAcc nS nT) True  = BinomAcc (nS+1) (nT+1)   addValue (BinomAcc nS nT) False = BinomAcc  nS    (nT+1)  --------------------------------------------------------------------- Ad-hoc type class--------------------------------------------------------------------- | Accumulator could be used to evaluate number of elements in---   sample.-class CalcCount m where-  -- | Number of elements in sample-  calcCount :: m -> Int---- | Monoids which could be used to calculate sample mean:------   \[ \bar{x} = \frac{1}{N}\sum_{i=1}^N{x_i} \]-class CalcMean m where-  -- | Returns @Nothing@ if there isn't enough data to make estimate.-  calcMean :: m -> Maybe Double---- | Monoids which could be used to calculate sample variance. Both---   methods return @Nothing@ if there isn't enough data to make---   estimate.-class CalcVariance m where-  -- | Calculate unbiased estimate of variance:-  ---  --   \[ \sigma^2 = \frac{1}{N-1}\sum_{i=1}^N(x_i - \bar{x})^2 \]-  calcVariance   :: m -> Maybe Double-  -- | Calculate maximum likelihood estimate of variance:-  ---  --   \[ \sigma^2 = \frac{1}{N}\sum_{i=1}^N(x_i - \bar{x})^2 \]-  calcVarianceML :: m -> Maybe Double---- | Calculate sample standard deviation from unbiased estimation of---   variance:------   \[ \sigma = \sqrt{\frac{1}{N-1}\sum_{i=1}^N(x_i - \bar{x})^2 } \]-calcStddev :: CalcVariance m => m -> Maybe Double-calcStddev = fmap sqrt . calcVariance---- | Calculate sample standard deviation from maximum likelihood---   estimation of variance:------   \[ \sigma = \sqrt{\frac{1}{N}\sum_{i=1}^N(x_i - \bar{x})^2 } \]-calcStddevML :: CalcVariance m => m -> Maybe Double-calcStddevML = fmap sqrt . calcVarianceML+-- | Value @a@ weighted by weight @w@+data Weighted w a = Weighted w a+              deriving (Show,Eq,Ord,Typeable,Data,Generic,Functor,Foldable,Traversable)   @@ -454,11 +455,6 @@   [| \(MeanKBN a (KBNSum b c)) -> (a,b,c)   |]   [| \(a,b,c) -> MeanKBN a (KBNSum b c) |] -derivingUnbox "WelfordMean"-  [t| WelfordMean -> (Int,Double) |]-  [| \(WelfordMean a b) -> (a,b)  |]-  [| \(a,b) -> WelfordMean a b    |]- derivingUnbox "Variance"   [t| Variance -> (Int,Double,Double) |]   [| \(Variance a b c) -> (a,b,c)  |]@@ -473,6 +469,11 @@   [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   |]  -- $references --

@@ -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+      ]+  ]

@@ -1,5 +1,5 @@ Name:           monoid-statistics-Version:        1.0.1.0+Version:        1.1.0 Cabal-Version:  >= 1.10 License:        BSD3 License-File:   LICENSE@@ -17,12 +17,15 @@   possibility to parallelize calculations. However not all statistics    could be calculated this way. +Extra-Source-Files:+  Changelog.md+ tested-with:-    GHC ==7.10.3-     || ==8.0.2-     || ==8.2.2-     || ==8.4.4-     || ==8.6.5+    GHC ==8.6.5+     || ==8.8.4+     || ==8.10.7+     || ==9.0.1+     || ==9.2.1  extra-source-files:   README.md@@ -34,15 +37,19 @@ Library   default-language: Haskell2010   ghc-options:      -Wall -O2-  Build-Depends:    base            >=4.8  && <5+  --+  Build-Depends:    base            >=4.12  && <5+                  , exceptions      >=0.10                   , vector          >=0.11 && <1                   , vector-th-unbox >=0.2.1.6-                  , math-functions  >=0.2.1.0+                  , math-functions  >=0.3+  --   Exposed-modules: Data.Monoid.Statistics                    Data.Monoid.Statistics.Class                    Data.Monoid.Statistics.Numeric+                   Data.Monoid.Statistics.Extra -test-suite tests+test-suite monoid-statistics-tests   default-language: Haskell2010   type:             exitcode-stdio-1.0   ghc-options:      -Wall -threaded@@ -51,11 +58,43 @@   if arch(i386)     ghc-options:  -msse2   hs-source-dirs: tests-  main-is:        Main.hs-  other-modules:-  build-depends: monoid-statistics-               , base             >=4.8 && <5-               , math-functions   >=0.2.1+  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.21+      , 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

@@ -1,31 +1,72 @@-{-# LANGUAGE LambdaCase          #-}+{-# LANGUAGE AllowAmbiguousTypes #-}+{-# LANGUAGE FlexibleContexts    #-}+{-# LANGUAGE FlexibleInstances   #-} {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications    #-} -- {-# OPTIONS_GHC -fno-warn-orphans #-}-import Data.Monoid+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  -data T a = T+----------------------------------------------------------------+-- 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, Arbitrary m, Show m, Eq m)-  => T m -> TestTree-p_memptyIsNeutral _-  = testProperty "mempty is neutral" $ \(m :: m) ->-       (m <> mempty) == m-    && (mempty <> m) == m+  :: 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. (Monoid m, Arbitrary m, Show m, Eq m)-  => T m -> TestTree-p_associativity _-  = testProperty "associativity" $ \(a :: m) b c ->+  :: 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)@@ -33,29 +74,35 @@      $ val1 == val2  p_commutativity-  :: forall m. (Monoid m, Arbitrary m, Show m, Eq m)-  => T m -> TestTree-p_commutativity _-  = testProperty "commutativity" $ \(a :: m) b ->-    (a <> b) == (b <> a)+  :: 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 a m. ( StatMonoid m a-                 , Arbitrary m, Show m, Eq m-                 , Arbitrary a, Show a, Eq a)-  => T a -> T m -> TestTree-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 a m. ( StatMonoid m a-                 , Arbitrary m, Show m, Eq m-                 , Arbitrary a, Show a, Eq a)-  => T a -> T m -> TestTree-p_addValue2 _ _-  = testProperty "addValue law" $ \(x :: a) (y :: a) ->+  :: 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)@@ -66,85 +113,132 @@  ---------------------------------------------------------------- -testType :: forall m. Typeable m => T m -> [T m -> TestTree] -> TestTree-testType t props = testGroup (show (typeRep (Proxy :: Proxy m)))-                             (fmap ($ t) props)+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"-  [ testType (T :: T (CountG Int))-      [ p_memptyIsNeutral-      , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Int)-      , p_addValue2 (T :: T Int)-      ]-  , testType (T :: T (Min Int))-      [ p_memptyIsNeutral-      , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Int)-      , p_addValue2 (T :: T Int)-      ]-  , testType (T :: T (Max Int))-      [ p_memptyIsNeutral-      , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Int)-      , p_addValue2 (T :: T Int)-      ]-  , testType (T :: T MinD)-      [ p_memptyIsNeutral-      , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Double)-      , p_addValue2 (T :: T Double)-      ]-  , testType (T :: T MaxD)-      [ p_memptyIsNeutral-      , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Double)-      , p_addValue2 (T :: T Double)-      ]-  , testType (T :: T BinomAcc)-      [ p_memptyIsNeutral-      , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Bool)-      , p_addValue2 (T :: T Bool)-      ]-  , testType (T :: T WelfordMean)-      [ p_memptyIsNeutral-      -- , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Double)-      -- , p_addValue2 (T :: T Double)-      ]-  , testType (T :: T MeanKBN)-      [ p_memptyIsNeutral-      -- , p_associativity-      -- , p_commutativity-      , p_addValue1 (T :: T Double)-      , p_addValue2 (T :: T Double)-      ]-  , testType (T :: T MeanKahan)-      [ p_memptyIsNeutral-      -- , p_associativity-      -- , p_commutativity-      , p_addValue1 (T :: T Double)-      -- , p_addValue2 (T :: T Double)-      ]-  , testType (T :: T Variance)-      [ p_memptyIsNeutral-      -- , p_associativity-      , p_commutativity-      , p_addValue1 (T :: T Double)-      , p_addValue2 (T :: T Double)-      ]+  [ 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@@ -153,19 +247,19 @@     return (CountG n)  instance (Arbitrary a) => Arbitrary (Max a) where-  arbitrary = Max <$> arbitrary+  arbitrary = fmap Max arbitrary  instance (Arbitrary a) => Arbitrary (Min a) where-  arbitrary = Min <$> arbitrary+  arbitrary = fmap Min arbitrary  instance Arbitrary MinD where-  arbitrary = frequency [ (1, pure mempty)-                        , (4, MinD <$> arbitrary)+  arbitrary = frequency [ (1, return mempty)+                        , (4, fmap MinD arbitrary)                         ]  instance Arbitrary MaxD where-  arbitrary = frequency [ (1, pure mempty)-                        , (4, MaxD <$> arbitrary)+  arbitrary = frequency [ (1, return mempty)+                        , (4, fmap MaxD arbitrary)                         ]  instance Arbitrary BinomAcc where@@ -174,34 +268,61 @@     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 >>= \case+  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 >>= \case+  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 MeanKBN where-  arbitrary = arbitrary >>= \case-    NonNegative 0 -> return mempty-    NonNegative n -> do-      x1 <- arbitrary-      x2 <- arbitrary-      x3 <- arbitrary-      return $ MeanKBN n (((zero `add` x1) `add` x2) `add` x3)+instance (Arbitrary a, Arbitrary w) => Arbitrary (Weighted w a) where+  arbitrary = Weighted <$> arbitrary <*> arbitrary -instance Arbitrary MeanKahan where-  arbitrary = arbitrary >>= \case-    NonNegative 0 -> return mempty-    NonNegative n -> do-      x1 <- arbitrary-      x2 <- arbitrary-      x3 <- arbitrary-      return $ MeanKahan n (((zero `add` x1) `add` x2) `add` x3)+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

@@ -0,0 +1,6 @@+module Main where++import Test.DocTest (doctest)++main :: IO ()+main = doctest ["Data"]