diff --git a/Changelog.md b/Changelog.md new file mode 100644 --- /dev/null +++ b/Changelog.md @@ -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. diff --git a/Data/Monoid/Statistics.hs b/Data/Monoid/Statistics.hs --- a/Data/Monoid/Statistics.hs +++ b/Data/Monoid/Statistics.hs @@ -12,3 +12,4 @@ import Data.Monoid.Statistics.Class import Data.Monoid.Statistics.Numeric + diff --git a/Data/Monoid/Statistics/Class.hs b/Data/Monoid/Statistics/Class.hs --- a/Data/Monoid/Statistics/Class.hs +++ b/Data/Monoid/Statistics/Class.hs @@ -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 @@ -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 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. + +-- $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,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 -- 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,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 diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -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 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"]