welford-online-mean-variance (empty) → 0.1.0.0

raw patch · 9 files changed

+334/−0 lines, 9 filesdep +QuickCheckdep +basedep +cerealsetup-changed

Dependencies added: QuickCheck, base, cereal, deepseq, tasty, tasty-discover, tasty-quickcheck, vector, welford-online-mean-variance

Files

• CHANGELOG.md +11/−0
• LICENSE +30/−0
• README.md +1/−0
• Setup.hs +2/−0
• src/Lib.hs +6/−0
• src/Statistics/Sample/WelfordOnlineMeanVariance.hs +167/−0
• test/Spec.hs +1/−0
• test/Statistics/Sample/WelfordOnlineTest.hs +53/−0
• welford-online-mean-variance.cabal +63/−0

@@ -0,0 +1,11 @@+# Changelog for `welford`++All notable changes to this project will be documented in this file.++The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),+and this project adheres to the+[Haskell Package Versioning Policy](https://pvp.haskell.org/).++## Unreleased++## 0.1.0.0 - YYYY-MM-DD

@@ -0,0 +1,30 @@+Copyright Author name here (c) 2022++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Author name here nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

@@ -0,0 +1,1 @@+# Welford: Online mean and variance computation

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,6 @@+module Lib+    ( someFunc+    ) where++someFunc :: IO ()+someFunc = putStrLn "someFunc"

@@ -0,0 +1,167 @@+{-# LANGUAGE DeriveAnyClass    #-}+{-# LANGUAGE DeriveGeneric     #-}+{-# LANGUAGE FlexibleInstances #-}+module Statistics.Sample.WelfordOnlineMeanVariance+  ( WelfordExistingAggregate(..)+  , addValue+  , finalize+  , nextValue+  , newWelfordAggregate+  , Mean+  , Variance+  , SampleVariance+  ) where++import           Control.DeepSeq+import           Data.Serialize+import qualified Data.Vector          as VB+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Unboxed  as VU+import           GHC.Generics++type Mean a = a+type Variance a = a+type SampleVariance a = a+++-- | For the storage of required information.+data WelfordExistingAggregate a+  = WelfordExistingAggregateEmpty+  | WelfordExistingAggregate+      { welfordCount :: !Int+      , welfordMean  :: !a+      , welfordM2    :: !a+      }+  deriving (Eq, Show, Read, Generic, NFData, Serialize)+++-- | Create a new empty Aggreate for the calculation.+newWelfordAggregate :: WelfordExistingAggregate a+newWelfordAggregate = WelfordExistingAggregateEmpty+++-- | Class for all data strucutres that can be used to computer the Welford approximation. For instance, this can be used to compute the Welford algorithm on a `Vector`s of `Fractional`, while only+-- requiring to handle one `WelfordExistingAggregate`.+class WelfordOnline a where+  plus :: a -> a -> a+  minus :: a -> a -> a+  multiply :: a -> a -> a+  divideInt :: a -> Int -> a+++instance (WelfordOnline a) => WelfordOnline (VB.Vector a) where+  plus          = VB.zipWith plus+  {-# INLINE plus #-}+  minus         = VB.zipWith minus+  {-# INLINE minus #-}+  multiply      = VB.zipWith multiply+  {-# INLINE multiply #-}+  divideInt x i = VB.map (`divideInt` i) x+  {-# INLINE divideInt #-}++instance (WelfordOnline a, VS.Storable a) => WelfordOnline (VS.Vector a) where+  plus          = VS.zipWith plus+  {-# INLINE plus #-}+  minus         = VS.zipWith minus+  {-# INLINE minus #-}+  multiply      = VS.zipWith multiply+  {-# INLINE multiply #-}+  divideInt x i = VS.map (`divideInt` i) x+  {-# INLINE divideInt #-}+++instance (WelfordOnline a, VU.Unbox a) => WelfordOnline (VU.Vector a) where+  plus          = VU.zipWith plus+  {-# INLINE plus #-}+  minus         = VU.zipWith minus+  {-# INLINE minus #-}+  multiply      = VU.zipWith multiply+  {-# INLINE multiply #-}+  divideInt x i = VU.map (`divideInt` i) x+  {-# INLINE divideInt #-}+++instance WelfordOnline Double where+  plus = (+)+  {-# INLINE plus #-}+  minus = (-)+  {-# INLINE minus #-}+  multiply = (*)+  {-# INLINE multiply #-}+  divideInt x i = x / fromIntegral i+  {-# INLINE divideInt #-}++instance WelfordOnline Float where+  plus = (+)+  {-# INLINE plus #-}+  minus = (-)+  {-# INLINE minus #-}+  multiply = (*)+  {-# INLINE multiply #-}+  divideInt x i = x / fromIntegral i+  {-# INLINE divideInt #-}++instance WelfordOnline Rational where+  plus = (+)+  {-# INLINE plus #-}+  minus = (-)+  {-# INLINE minus #-}+  multiply = (*)+  {-# INLINE multiply #-}+  divideInt x i = x / fromIntegral i+  {-# INLINE divideInt #-}+++-- | Add one value to the current aggregate.+addValue :: (WelfordOnline a) => WelfordExistingAggregate a -> a -> WelfordExistingAggregate a+addValue (WelfordExistingAggregate count mean m2) val =+  let count' = count + 1+      delta = val `minus` mean+      mean' = mean `plus` (delta `divideInt` count')+      delta2 = val `minus` mean'+      m2' = m2 `plus` (delta `multiply` delta2)+   in WelfordExistingAggregate count' mean' m2'+addValue WelfordExistingAggregateEmpty val =+  let count' = 1+      delta' = val+      mean' = delta' `divideInt` count'+      delta2' = val `minus` mean'+      m2' = delta' `multiply` delta2'+   in WelfordExistingAggregate count' mean' m2'+++-- | Calculate mean, variance and sample variance from aggregate.+finalize :: (WelfordOnline a) => WelfordExistingAggregate a -> (Mean a, Variance a, SampleVariance a)+finalize (WelfordExistingAggregate count mean m2)+  | count < 2 = (mean, m2, m2)+  | otherwise = (mean, m2 `divideInt` count, m2 `divideInt` (count - 1))+finalize WelfordExistingAggregateEmpty = error "finalize: Emtpy Welford Online Aggreate. Add data first!"++-- | Add a new sample to the aggregate and compute mean and variances.+nextValue :: (WelfordOnline a) => WelfordExistingAggregate a -> a -> (WelfordExistingAggregate a, (Mean a, Variance a, SampleVariance a))+nextValue agg val =+  let agg' = addValue agg val+   in (agg', finalize agg')+++-- # For a new value newValue, compute the new count, new mean, the new M2.+-- # mean accumulates the mean of the entire dataset+-- # M2 aggregates the squared distance from the mean+-- # count aggregates the number of samples seen so far+-- def update(existingAggregate, newValue):+--     (count, mean, M2) = existingAggregate+--     count' += 1+--     delta = newValue - mean+--     mean' += delta / count'+--     delta2 = newValue - mean'+--     M2' += delta * delta2+--     return (count', mean', M2')++-- # Retrieve the mean, variance and sample variance from an aggregate+-- def finalize(existingAggregate):+--     (count, mean, M2) = existingAggregate+--     if count < 2:+--         return float("nan")+--     else:+--         (mean, variance, sampleVariance) = (mean, M2 / count, M2 / (count - 1))+--         return (mean, variance, sampleVariance)

@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF tasty-discover #-}

@@ -0,0 +1,53 @@+module Statistics.Sample.WelfordOnlineTest where++import           Control.Monad+import           Data.List                                   (foldl')+import qualified Data.Vector                                 as VB+import           Debug.Trace+import           Test.Tasty.QuickCheck++import           Statistics.Sample.WelfordOnlineMeanVariance++main :: IO ()+main = putStrLn "Test suite not yet implemented"+++-- | Recursively generate candles+generateNValues :: Int -> Gen [Double]+generateNValues n+  | n <= 0 = return []+  | otherwise =  do+      x <- arbitrary+      (x:) <$> generateNValues (n - 1)++epsEqWith :: (Ord a, Fractional a, Show a) => a -> a -> a -> Property+epsEqWith eps a b+  | abs (a - b) <= eps = a === a+  | otherwise = a === b+++prop_MeanAndVariance :: Gen Property+prop_MeanAndVariance = do+  vals <- sized $ \n -> generateNValues (10 + 5 * n)+  let n = fromIntegral (length vals)+      mean = sum vals / n+      var = sum (map (\x -> (x - mean) ^ 2) vals) / (n - 1)+      (wMean, _, wVarSample) = snd $ foldl' (nextValue . fst) (WelfordExistingAggregateEmpty, (0, 0, 0)) vals+      eps = min 0.01 $ max 0.001 (0.001 * mean)+  return $ epsEqWith eps wMean mean .&&. epsEqWith eps wVarSample var+++prop_MeanAndVarianceVector :: Gen Property+prop_MeanAndVarianceVector = do+  len <- chooseInt (1, 100)+  vecs <- sized $ \n -> replicateM len (VB.fromList <$> generateNValues (10 + 5 * n))+  let mkRes vals =+        let n = fromIntegral (VB.length vals)+            mean = VB.sum vals / n+            var = VB.sum (VB.map (\x -> (x - mean) ^ 2) vals) / (n - 1)+            (wMean, _, wVarSample) = snd $ foldl' (nextValue . fst) (WelfordExistingAggregateEmpty, (0, 0, 0)) (VB.toList vals)+            eps = min 0.01 $ max 0.001 (0.001 * mean)+        in (eps, wMean, mean, wVarSample, var)+      ress = map mkRes vecs+  return $+    foldl1 (.&&.) $ map (\(eps, wMean, mean, wVarSample, var) -> epsEqWith eps wMean mean .&&. epsEqWith eps wVarSample var) ress

@@ -0,0 +1,63 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.34.4.+--+-- see: https://github.com/sol/hpack++name:           welford-online-mean-variance+version:        0.1.0.0+synopsis:       Online computation of mean and variance using the Welford algorithm.+description:    Please see the README on GitHub at <https://github.com/githubuser/welford-online-mean-variance#readme>+category:       Statistics+homepage:       https://github.com/schnecki/welford-online-mean-variance#readme+bug-reports:    https://github.com/schnecki/welford-online-mean-variance/issues+author:         Manuel Schneckenreither+maintainer:     manuel.schnecki@gmail.com+copyright:      2022 Manuel Schneckenreither+license:        BSD3+license-file:   LICENSE+build-type:     Simple+extra-source-files:+    README.md+    CHANGELOG.md++source-repository head+  type: git+  location: https://github.com/schnecki/welford-online-mean-variance++library+  exposed-modules:+      Lib+      Statistics.Sample.WelfordOnlineMeanVariance+  other-modules:+      Paths_welford_online_mean_variance+  hs-source-dirs:+      src+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints+  build-depends:+      base >=4.7 && <5+    , cereal+    , deepseq+    , vector+  default-language: Haskell2010++test-suite welford-online-mean-variance-test+  type: exitcode-stdio-1.0+  main-is: Spec.hs+  other-modules:+      Statistics.Sample.WelfordOnlineTest+      Paths_welford_online_mean_variance+  hs-source-dirs:+      test+  ghc-options: -Wall -Wcompat -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns -Wmissing-export-lists -Wmissing-home-modules -Wpartial-fields -Wredundant-constraints -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      QuickCheck+    , base >=4.7 && <5+    , cereal+    , deepseq+    , tasty+    , tasty-discover+    , tasty-quickcheck+    , vector+    , welford-online-mean-variance+  default-language: Haskell2010