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