monoid-statistics 0.3.1 → 1.0.0
raw patch · 6 files changed
+707/−300 lines, 6 filesdep +QuickCheckdep +math-functionsdep +monoid-statisticsdep ~base
Dependencies added: QuickCheck, math-functions, monoid-statistics, tasty, tasty-quickcheck, vector, vector-th-unbox
Dependency ranges changed: base
Files
- Data/Monoid/Statistics.hs +5/−135
- Data/Monoid/Statistics/Class.hs +128/−0
- Data/Monoid/Statistics/Numeric.hs +326/−149
- README.md +7/−0
- monoid-statistics.cabal +34/−16
- tests/Main.hs +207/−0
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,10 @@ -- 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.------ $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+import Data.Monoid.Statistics.Class+import Data.Monoid.Statistics.Numeric
+ Data/Monoid/Statistics/Class.hs view
@@ -0,0 +1,128 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+--+{-# OPTIONS_GHC -fno-warn-orphans #-}+-- |+-- 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+ ( -- * Type class and helpers+ StatMonoid(..)+ , reduceSample+ , reduceSampleVec+ -- * Data types+ , Pair(..)+ ) where++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.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).+--+-- 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, merge 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'.+reduceSample :: (F.Foldable f, StatMonoid m a) => 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+reduceSampleVec = G.foldl' addValue mempty+{-# INLINE reduceSampleVec #-}+++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 Monoid KahanSum where+ mempty = zero+ mappend s1 s2 = add s1 (kahan s2)+instance Real a => StatMonoid KahanSum a where+ addValue m x = add m (realToFrac x)+ {-# INLINE addValue #-}++instance Monoid KBNSum where+ mempty = zero+ mappend s1 s2 = add s1 (kbn s2)+instance Real a => StatMonoid KBNSum a where+ addValue m x = add m (realToFrac x)+ {-# INLINE addValue #-}+++----------------------------------------------------------------+-- 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 (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')+ {-# 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 #-}++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 |]
Data/Monoid/Statistics/Numeric.hs view
@@ -1,256 +1,433 @@ {-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE DeriveDataTypeable #-}-module Data.Monoid.Statistics.Numeric ( - -- * Mean and variance- Count(..)+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+module Data.Monoid.Statistics.Numeric (+ -- * Mean & Variance+ -- ** Number of elements+ CountG(..)+ , Count , asCount- , Mean(..)- , asMean+ -- ** Mean+ , MeanKBN(..)+ , asMeanKBN+ , WelfordMean(..)+ , asWelfordMean+ , MeanKahan(..)+ , asMeanKahan+ -- ** Variance , Variance(..) , asVariance- -- ** Ad-hoc accessors- -- $accessors+ -- * Maximum and minimum+ , Max(..)+ , Min(..)+ , MaxD(..)+ , MinD(..)+ -- * Binomial trials+ , BinomAcc(..)+ , asBinomAcc+ -- * Accessors , CalcCount(..) , CalcMean(..) , CalcVariance(..) , calcStddev- , calcStddevUnbiased- -- * Maximum and minimum- , Max(..)- , Min(..)+ , calcStddevML+ -- * References+ -- $references ) where -import Data.Monoid-import Data.Monoid.Statistics-import Data.Typeable (Typeable)+import Data.Monoid ((<>))+import Data.Monoid.Statistics.Class+import Data.Data (Typeable,Data)+import Data.Vector.Unboxed (Unbox)+import Data.Vector.Unboxed.Deriving (derivingUnbox)+import Numeric.Sum+import GHC.Generics (Generic) ---------------------------------------------------------------- -- Statistical monoids ---------------------------------------------------------------- --- | Simplest statistics. Number of elements in the sample-newtype Count a = Count { calcCountI :: a }+-- | Calculate number of elements in the sample.+newtype CountG a = CountG { calcCountN :: a } deriving (Show,Eq,Ord,Typeable) --- | 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)+instance Integral a => Monoid (CountG a) where+ mempty = CountG 0+ CountG i `mappend` CountG j = CountG (i + j) {-# INLINE mempty #-} {-# INLINE mappend #-}- -instance (Integral a) => StatMonoid (Count a) b where- pappend _ !(Count n) = Count (n + 1)- {-# INLINE pappend #-} -instance CalcCount (Count Int) where- calcCount = calcCountI+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 #-} +---------------------------------------------------------------- --- | Mean of sample. Samples of Double,Float and bui;t-in integral--- types are supported+-- | Incremental calculation of mean. Sum of elements is calculated+-- using compensated Kahan summation.+data MeanKahan = MeanKahan !Int !KahanSum+ deriving (Show,Eq,Typeable,Data,Generic)++asMeanKahan :: MeanKahan -> MeanKahan+asMeanKahan = id++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 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 _) = Nothing+ calcMean (MeanKahan n s) = Just (kahan s / fromIntegral n)++++-- | Incremental calculation of mean. Sum of elements is calculated+-- using Kahan-Babuška-Neumaier summation.+data MeanKBN = MeanKBN !Int !KBNSum+ deriving (Show,Eq,Typeable,Data,Generic)++asMeanKBN :: MeanKBN -> MeanKBN+asMeanKBN = id++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 Real a => StatMonoid MeanKBN a where+ addValue (MeanKBN n m) x = MeanKBN (n+1) (addValue m x)++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)++++-- | Incremental calculation of mean. One of algorithm's advantage is+-- protection against double overflow: ----- Numeric stability of 'mappend' is not proven.-data Mean = Mean {-# UNPACK #-} !Int -- Number of entries- {-# UNPACK #-} !Double -- Current mean- deriving (Show,Eq,Typeable)+-- > λ> 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+ deriving (Show,Eq,Typeable,Data,Generic) --- | Fix type of monoid-asMean :: Mean -> Mean-asMean = id-{-# INLINE asMean #-}+-- | Type restricted 'id'+asWelfordMean :: WelfordMean -> WelfordMean+asWelfordMean = id -instance Monoid Mean where- mempty = Mean 0 0- mappend !(Mean n x) !(Mean k y) = Mean (n + k) ((x*n' + y*k') / (n' + k')) +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 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 #-}+-- | \[ 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')+ where+ n' = n+1+ {-# INLINE addValue #-} -instance CalcCount Mean where- calcCount (Mean n _) = n- {-# INLINE calcCount #-}-instance CalcMean Mean where- calcMean (Mean _ m) = m- {-# INLINE calcMean #-}+instance CalcCount WelfordMean where+ calcCount (WelfordMean n _) = n+instance CalcMean WelfordMean where+ calcMean (WelfordMean 0 _) = Nothing+ calcMean (WelfordMean _ m) = Just m +---------------------------------------------------------------- --- | Intermediate quantities to calculate the standard deviation.+-- | Incremental algorithms for calculation the standard deviation. 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) --- | 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>---+-- | 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) = Variance (n1+n2) (ta+tb) sumsq+ mappend (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)+ sumsq | n1 == 0 = sb+ | n2 == 0 = sa+ | otherwise = sa + sb + nom / ((na + nb) * na * nb) {-# INLINE mempty #-} {-# INLINE 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 #-}+ 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 _ _) = Nothing+ calcMean (Variance n s _) = Just (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 = Nothing+ | otherwise = Just $! s / fromIntegral (n - 1)+ calcVarianceML (Variance n _ s)+ | n < 1 = Nothing+ | otherwise = Just $! s / fromIntegral n +---------------------------------------------------------------- --- | 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)+-- | Calculate minimum of sample+newtype Min a = Min { calcMin :: Maybe a }+ deriving (Show,Eq,Ord,Typeable,Data,Generic) +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++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)++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++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 (Show,Typeable,Data,Generic)++instance Eq MinD where+ MinD a == MinD b+ | isNaN a && isNaN b = True+ | otherwise = a == b+ -- 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)+instance Monoid MinD where+ mempty = MinD (0/0)+ mappend (MinD x) (MinD y)+ | isNaN x = MinD y+ | isNaN y = MinD x+ | otherwise = MinD (min x y) {-# INLINE mempty #-}- {-# INLINE mappend #-} + {-# INLINE 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 (Show,Typeable,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)+instance Eq MaxD where+ MaxD a == MaxD b+ | isNaN a && isNaN b = True+ | otherwise = a == b++instance Monoid MaxD where+ mempty = MaxD (0/0)+ mappend (MaxD x) (MaxD y)+ | isNaN x = MaxD y+ | isNaN y = MaxD x+ | otherwise = MaxD (max x y) {-# INLINE mempty #-}- {-# INLINE mappend #-} + {-# INLINE mappend #-} -instance StatMonoid Max Double where- pappend !x m = mappend (Max x) m- {-# INLINE pappend #-}+instance a ~ Double => StatMonoid MaxD a where+ singletonMonoid = MaxD +---------------------------------------------------------------- +-- | Accumulator for binomial trials.+data BinomAcc = BinomAcc { binomAccSuccess :: !Int+ , binomAccTotal :: !Int+ }+ deriving (Show,Eq,Ord,Typeable,Data,Generic) +-- | Type restricted 'id'+asBinomAcc :: BinomAcc -> BinomAcc+asBinomAcc = id++instance Monoid BinomAcc where+ mempty = BinomAcc 0 0+ mappend (BinomAcc n1 m1) (BinomAcc n2 m2) = BinomAcc (n1+n2) (m1+m2)++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 ----------------------------------------------------------------- --- $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) - ---- | Statistics which could count number of elements in the sample+-- | Accumulator could be used to evaluate number of elements in+-- sample. class CalcCount m where -- | Number of elements in sample calcCount :: m -> Int --- | Statistics which could estimate mean of sample+-- | Monoids which could be used to calculate sample mean:+--+-- \[ \bar{x} = \frac{1}{N}\sum_{i=1}^N{x_i} \] class CalcMean m where- -- | Calculate esimate of mean of a sample- calcMean :: m -> Double- --- | Statistics which could estimate variance of sample+ -- | 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 biased estimate of variance- calcVariance :: m -> Double- -- | Calculate unbiased estimate of the variance, where the- -- denominator is $n-1$.- calcVarianceUnbiased :: m -> Double+ -- | 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 (biased estimator, $s$, where--- the denominator is $n-1$).-calcStddev :: CalcVariance m => m -> Double-calcStddev = sqrt . calcVariance-{-# INLINE calcStddev #-}+-- | 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 standard deviation of the sample--- (unbiased estimator, $\sigma$, where the denominator is $n$).-calcStddevUnbiased :: CalcVariance m => m -> Double-calcStddevUnbiased = sqrt . calcVarianceUnbiased-{-# INLINE calcStddevUnbiased #-}+-- | 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 ---------------------------------------------------------------- -- 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 "MeanKBN"+ [t| MeanKBN -> (Int,Double,Double) |]+ [| \(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) |]+ [| \(a,b,c) -> Variance a b c |]++derivingUnbox "MinD"+ [t| MinD -> Double |]+ [| calcMinD |]+ [| MinD |]++derivingUnbox "MaxD"+ [t| MaxD -> Double |]+ [| calcMaxD |]+ [| MaxD |]++-- $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.
monoid-statistics.cabal view
@@ -1,13 +1,12 @@-- Name: monoid-statistics-Version: 0.3.1-Cabal-Version: >= 1.6+Version: 1.0.0+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,39 @@ 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:+ 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.8 && <5+ , vector >=0.11 && <1+ , vector-th-unbox >=0.2.1.6+ , math-functions >=0.2.1.0 Exposed-modules: Data.Monoid.Statistics+ Data.Monoid.Statistics.Class Data.Monoid.Statistics.Numeric++test-suite 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.8 && <5+ , math-functions >=0.2.1+ , tasty >=0.11+ , tasty-quickcheck >=0.9+ , QuickCheck
+ tests/Main.hs view
@@ -0,0 +1,207 @@+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE ScopedTypeVariables #-}+--+{-# OPTIONS_GHC -fno-warn-orphans #-}+import Data.Monoid+import Data.Typeable+import Numeric.Sum+import Test.Tasty+import Test.Tasty.QuickCheck++import Data.Monoid.Statistics+++data T a = T++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++p_associativity+ :: forall m. (Monoid m, Arbitrary m, Show m, Eq m)+ => T m -> TestTree+p_associativity _+ = 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, Arbitrary m, Show m, Eq m)+ => T m -> TestTree+p_commutativity _+ = testProperty "commutativity" $ \(a :: m) b ->+ (a <> b) == (b <> a)++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 _ _+ = 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) ->+ 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++++----------------------------------------------------------------++testType :: forall m. Typeable m => T m -> [T m -> TestTree] -> TestTree+testType t props = testGroup (show (typeRep (Proxy :: Proxy m)))+ (fmap ($ t) props)+++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)+ ]+ ]++----------------------------------------------------------------++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 = Max <$> arbitrary++instance (Arbitrary a) => Arbitrary (Min a) where+ arbitrary = Min <$> arbitrary++instance Arbitrary MinD where+ arbitrary = frequency [ (1, pure mempty)+ , (4, MinD <$> arbitrary)+ ]++instance Arbitrary MaxD where+ arbitrary = frequency [ (1, pure mempty)+ , (4, MaxD <$> arbitrary)+ ]++instance Arbitrary BinomAcc where+ arbitrary = do+ NonNegative nSucc <- arbitrary+ NonNegative nFail <- arbitrary+ return $ BinomAcc nSucc (nFail + nSucc)++instance Arbitrary WelfordMean where+ arbitrary = arbitrary >>= \case+ NonNegative 0 -> return mempty+ NonNegative n -> do m <- arbitrary+ return (WelfordMean n m)++instance Arbitrary Variance where+ arbitrary = arbitrary >>= \case+ 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 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)