diff --git a/Data/Monoid/Statistics.hs b/Data/Monoid/Statistics.hs --- a/Data/Monoid/Statistics.hs +++ b/Data/Monoid/Statistics.hs @@ -1,7 +1,3 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE DeriveDataTypeable #-} -- | -- Module : Data.Monoid.Statistics -- Copyright : Copyright (c) 2010, Alexey Khudyakov @@ -9,136 +5,10 @@ -- Maintainer : Alexey Khudyakov -- 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 diff --git a/Data/Monoid/Statistics/Class.hs b/Data/Monoid/Statistics/Class.hs new file mode 100644 --- /dev/null +++ b/Data/Monoid/Statistics/Class.hs @@ -0,0 +1,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 +-- License : BSD3 +-- Maintainer : Alexey Khudyakov +-- 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 |] diff --git a/Data/Monoid/Statistics/Numeric.hs b/Data/Monoid/Statistics/Numeric.hs --- a/Data/Monoid/Statistics/Numeric.hs +++ b/Data/Monoid/Statistics/Numeric.hs @@ -1,256 +1,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. --- --- --- +-- | 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. +-- +-- * [Chan1979] Chan, Tony F.; Golub, Gene H.; LeVeque, Randall +-- J. (1979), Updating Formulae and a Pairwise Algorithm for +-- Computing Sample Variances., Technical Report STAN-CS-79-773, +-- Department of Computer Science, Stanford University. Page 4. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,7 @@ +# monoid-statistics parallelizable constant space estimators + +[![Build Status](https://travis-ci.org/Shimuuar/monoid-statistics.png?branch=master)](https://travis-ci.org/Shimuuar/monoid-statistics) + +Monoids for calculation of statistics of sample. This approach allows to +calculate many statistics in one pass over data and possibility to parallelize +calculations. However not all statistics could be calculated this way. diff --git a/monoid-statistics.cabal b/monoid-statistics.cabal --- a/monoid-statistics.cabal +++ b/monoid-statistics.cabal @@ -1,13 +1,12 @@ - - Name: monoid-statistics -Version: 0.3.1 -Cabal-Version: >= 1.6 +Version: 1.0.0 +Cabal-Version: >= 1.10 License: BSD3 License-File: LICENSE Author: Alexey Khudyakov Maintainer: Alexey Khudyakov -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 diff --git a/tests/Main.hs b/tests/Main.hs new file mode 100644 --- /dev/null +++ b/tests/Main.hs @@ -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)