streamly-statistics (empty) → 0.1.0
raw patch · 9 files changed
+2098/−0 lines, 9 filesdep +QuickCheckdep +basedep +containerssetup-changed
Dependencies added: QuickCheck, base, containers, deepseq, deque, fusion-plugin, hspec, hspec-core, mwc-random, random, statistics, streamly-core, streamly-statistics, tasty, tasty-bench, vector
Files
- CHANGELOG.md +5/−0
- LICENSE +177/−0
- NOTICE +5/−0
- README.md +18/−0
- Setup.hs +2/−0
- benchmark/Main.hs +274/−0
- src/Streamly/Statistics.hs +1155/−0
- streamly-statistics.cabal +157/−0
- test/Main.hs +305/−0
+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Changelog++## 0.1.0 (Apr 2023)++* Initial version
+ LICENSE view
@@ -0,0 +1,177 @@++ Apache License+ Version 2.0, January 2004+ http://www.apache.org/licenses/++ TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION++ 1. Definitions.++ "License" shall mean the terms and conditions for use, reproduction,+ and distribution as defined by Sections 1 through 9 of this document.++ "Licensor" shall mean the copyright owner or entity authorized by+ the copyright owner that is granting the License.++ "Legal Entity" shall mean the union of the acting entity and all+ other entities that control, are controlled by, or are under common+ control with that entity. For the purposes of this definition,+ "control" means (i) the power, direct or indirect, to cause the+ direction or management of such entity, whether by contract or+ otherwise, or (ii) ownership of fifty percent (50%) or more of the+ outstanding shares, or (iii) beneficial ownership of such entity.++ "You" (or "Your") shall mean an individual or Legal Entity+ exercising permissions granted by this License.++ "Source" form shall mean the preferred form for making modifications,+ including but not limited to software source code, documentation+ source, and configuration files.++ "Object" form shall mean any form resulting from mechanical+ transformation or translation of a Source form, including but+ not limited to compiled object code, generated documentation,+ and conversions to other media types.++ "Work" shall mean the work of authorship, whether in Source or+ Object form, made available under the License, as indicated by a+ copyright notice that is included in or attached to the work+ (an example is provided in the Appendix below).++ "Derivative Works" shall mean any work, whether in Source or Object+ form, that is based on (or derived from) the Work and for which the+ editorial revisions, annotations, elaborations, or other modifications+ represent, as a whole, an original work of authorship. For the purposes+ of this License, Derivative Works shall not include works that remain+ separable from, or merely link (or bind by name) to the interfaces of,+ the Work and Derivative Works thereof.++ "Contribution" shall mean any work of authorship, including+ the original version of the Work and any modifications or additions+ to that Work or Derivative Works thereof, that is intentionally+ submitted to Licensor for inclusion in the Work by the copyright owner+ or by an individual or Legal Entity authorized to submit on behalf of+ the copyright owner. For the purposes of this definition, "submitted"+ means any form of electronic, verbal, or written communication sent+ to the Licensor or its representatives, including but not limited to+ communication on electronic mailing lists, source code control systems,+ and issue tracking systems that are managed by, or on behalf of, the+ Licensor for the purpose of discussing and improving the Work, but+ excluding communication that is conspicuously marked or otherwise+ designated in writing by the copyright owner as "Not a Contribution."++ "Contributor" shall mean Licensor and any individual or Legal Entity+ on behalf of whom a Contribution has been received by Licensor and+ subsequently incorporated within the Work.++ 2. Grant of Copyright License. Subject to the terms and conditions of+ this License, each Contributor hereby grants to You a perpetual,+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable+ copyright license to reproduce, prepare Derivative Works of,+ publicly display, publicly perform, sublicense, and distribute the+ Work and such Derivative Works in Source or Object form.++ 3. Grant of Patent License. Subject to the terms and conditions of+ this License, each Contributor hereby grants to You a perpetual,+ worldwide, non-exclusive, no-charge, royalty-free, irrevocable+ (except as stated in this section) patent license to make, have made,+ use, offer to sell, sell, import, and otherwise transfer the Work,+ where such license applies only to those patent claims licensable+ by such Contributor that are necessarily infringed by their+ Contribution(s) alone or by combination of their Contribution(s)+ with the Work to which such Contribution(s) was submitted. If You+ institute patent litigation against any entity (including a+ cross-claim or counterclaim in a lawsuit) alleging that the Work+ or a Contribution incorporated within the Work constitutes direct+ or contributory patent infringement, then any patent licenses+ granted to You under this License for that Work shall terminate+ as of the date such litigation is filed.++ 4. Redistribution. You may reproduce and distribute copies of the+ Work or Derivative Works thereof in any medium, with or without+ modifications, and in Source or Object form, provided that You+ meet the following conditions:++ (a) You must give any other recipients of the Work or+ Derivative Works a copy of this License; and++ (b) You must cause any modified files to carry prominent notices+ stating that You changed the files; and++ (c) You must retain, in the Source form of any Derivative Works+ that You distribute, all copyright, patent, trademark, and+ attribution notices from the Source form of the Work,+ excluding those notices that do not pertain to any part of+ the Derivative Works; and++ (d) If the Work includes a "NOTICE" text file as part of its+ distribution, then any Derivative Works that You distribute must+ include a readable copy of the attribution notices contained+ within such NOTICE file, excluding those notices that do not+ pertain to any part of the Derivative Works, in at least one+ of the following places: within a NOTICE text file distributed+ as part of the Derivative Works; within the Source form or+ documentation, if provided along with the Derivative Works; or,+ within a display generated by the Derivative Works, if and+ wherever such third-party notices normally appear. The contents+ of the NOTICE file are for informational purposes only and+ do not modify the License. You may add Your own attribution+ notices within Derivative Works that You distribute, alongside+ or as an addendum to the NOTICE text from the Work, provided+ that such additional attribution notices cannot be construed+ as modifying the License.++ You may add Your own copyright statement to Your modifications and+ may provide additional or different license terms and conditions+ for use, reproduction, or distribution of Your modifications, or+ for any such Derivative Works as a whole, provided Your use,+ reproduction, and distribution of the Work otherwise complies with+ the conditions stated in this License.++ 5. Submission of Contributions. Unless You explicitly state otherwise,+ any Contribution intentionally submitted for inclusion in the Work+ by You to the Licensor shall be under the terms and conditions of+ this License, without any additional terms or conditions.+ Notwithstanding the above, nothing herein shall supersede or modify+ the terms of any separate license agreement you may have executed+ with Licensor regarding such Contributions.++ 6. Trademarks. This License does not grant permission to use the trade+ names, trademarks, service marks, or product names of the Licensor,+ except as required for reasonable and customary use in describing the+ origin of the Work and reproducing the content of the NOTICE file.++ 7. Disclaimer of Warranty. Unless required by applicable law or+ agreed to in writing, Licensor provides the Work (and each+ Contributor provides its Contributions) on an "AS IS" BASIS,+ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or+ implied, including, without limitation, any warranties or conditions+ of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A+ PARTICULAR PURPOSE. You are solely responsible for determining the+ appropriateness of using or redistributing the Work and assume any+ risks associated with Your exercise of permissions under this License.++ 8. Limitation of Liability. In no event and under no legal theory,+ whether in tort (including negligence), contract, or otherwise,+ unless required by applicable law (such as deliberate and grossly+ negligent acts) or agreed to in writing, shall any Contributor be+ liable to You for damages, including any direct, indirect, special,+ incidental, or consequential damages of any character arising as a+ result of this License or out of the use or inability to use the+ Work (including but not limited to damages for loss of goodwill,+ work stoppage, computer failure or malfunction, or any and all+ other commercial damages or losses), even if such Contributor+ has been advised of the possibility of such damages.++ 9. Accepting Warranty or Additional Liability. While redistributing+ the Work or Derivative Works thereof, You may choose to offer,+ and charge a fee for, acceptance of support, warranty, indemnity,+ or other liability obligations and/or rights consistent with this+ License. However, in accepting such obligations, You may act only+ on Your own behalf and on Your sole responsibility, not on behalf+ of any other Contributor, and only if You agree to indemnify,+ defend, and hold each Contributor harmless for any liability+ incurred by, or claims asserted against, such Contributor by reason+ of your accepting any such warranty or additional liability.++ END OF TERMS AND CONDITIONS
+ NOTICE view
@@ -0,0 +1,5 @@+streamly-statistics+Copyright 2021 Composewell Technologies++This product includes software developed at+Composewell Technologies (http://www.composewell.com).
+ README.md view
@@ -0,0 +1,18 @@+# streamly-statistics++Statistical measures for finite or infinite data streams.++All operations use numerically stable floating point arithmetic. Measurements+can be performed over the entire input stream or on a sliding window of fixed+or variable size. Where possible, measures are computed online without+buffering the input stream.++Includes:++* Summary: length, sum, powerSum+* Location: minimum, maximum, rawMoments, means, exponential smoothing+* Spread: range, variance, deviations+* Shape: skewness, kurtosis+* Sample statistics, resampling+* Probablity distribution: frequency, mode, histograms+* Transforms: Fast fourier transform
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ benchmark/Main.hs view
@@ -0,0 +1,274 @@+{-# LANGUAGE TupleSections #-}++import Control.DeepSeq (NFData)+import Streamly.Data.Fold (Fold)+import Streamly.Data.Stream (Stream)+import System.Random (randomRIO)++import qualified Streamly.Data.Fold as Fold+import qualified Streamly.Data.Stream as Stream+import qualified Streamly.Data.Array as Array+import qualified Streamly.Internal.Data.Ring.Unboxed as Ring+import qualified Streamly.Statistics as Statistics++import Gauge++{-# INLINE source #-}+source :: (Monad m, Num a, Stream.Enumerable a) => Int -> a -> Stream m a+source len from =+ Stream.enumerateFromThenTo from (from + 1) (from + fromIntegral len)++{-# INLINE sourceDescending #-}+sourceDescending :: (Monad m, Num a, Stream.Enumerable a) =>+ Int -> a -> Stream m a+sourceDescending len from =+ Stream.enumerateFromThenTo+ (from + fromIntegral len)+ (from + fromIntegral (len - 1))+ from++{-# INLINE sourceDescendingInt #-}+sourceDescendingInt :: Monad m => Int -> Int -> Stream m Int+sourceDescendingInt = sourceDescending++{-# INLINE benchWith #-}+benchWith :: (Num a, NFData a) =>+ (Int -> a -> Stream IO a) -> Int -> String -> Fold IO a a -> Benchmark+benchWith src len name f =+ bench name+ $ nfIO+ $ randomRIO (1, 1 :: Int) >>= Stream.fold f . src len . fromIntegral++{-# INLINE benchWithFold #-}+benchWithFold :: Int -> String -> Fold IO Double Double -> Benchmark+benchWithFold len name f = benchWith source len name f++{-# INLINE benchWithFoldInt #-}+benchWithFoldInt :: Int -> String -> Fold IO Int Int -> Benchmark+benchWithFoldInt len name f = benchWith source len name f++{-# INLINE benchWithPostscan #-}+benchWithPostscan :: Int -> String -> Fold IO Double Double -> Benchmark+benchWithPostscan len name f =+ bench name $ nfIO $ randomRIO (1, 1) >>=+ Stream.fold Fold.drain . Stream.postscan f . source len++{-# INLINE benchWithResample #-}+benchWithResample :: Int -> String -> Benchmark+benchWithResample len name = bench name $ nfIO $ do+ i <- randomRIO (1, 1)+ arr <- Stream.fold Array.write (source len i :: Stream IO Double)+ Stream.fold Fold.drain $ Stream.unfold Statistics.resample arr++{-# INLINE benchWithFoldResamples #-}+benchWithFoldResamples :: Int -> String -> Fold IO Double Double -> Benchmark+benchWithFoldResamples len name f = bench name $ nfIO $ do+ i <- randomRIO (1, 1)+ arr <- Stream.fold Array.write (source len i :: Stream IO Double)+ Stream.fold Fold.drain $ Statistics.foldResamples len arr f++{-# INLINE numElements #-}+numElements :: Int+numElements = 100000++main :: IO ()+main =+ defaultMain+ [ bgroup+ "fold"+ [ benchWithFold numElements "minimum (window size 100)"+ (Ring.slidingWindow 100 Statistics.minimum)+ , benchWithFold numElements "minimum (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.minimum)+ , benchWith sourceDescendingInt numElements+ "minimum descending (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.minimum)++ , benchWithFold numElements "maximum (window size 100)"+ (Ring.slidingWindow 100 Statistics.maximum)+ , benchWithFold numElements "maximum (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.maximum)+ , benchWith sourceDescendingInt numElements+ "maximum descending (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.maximum)++ , benchWithFold numElements "range (window size 100)"+ (Ring.slidingWindow 100 Statistics.range)+ , benchWithFold numElements "range (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.range)++ , benchWithFoldInt numElements "sumInt (window size 100)"+ (Ring.slidingWindow 100 Statistics.sumInt)+ , benchWithFoldInt numElements "sum for Int (window size 100)"+ (Ring.slidingWindow 100 Statistics.sum)++ , benchWithFold numElements "sum (window size 100)"+ (Ring.slidingWindow 100 Statistics.sum)+ , benchWithFold numElements "sum (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.sum)+ , benchWithFold numElements "sum (entire stream)"+ (Statistics.cumulative Statistics.sum)+ , benchWithFold numElements "sum (Data.Fold)"+ (Fold.sum)++ , benchWithFold numElements "mean (window size 100)"+ (Ring.slidingWindow 100 Statistics.mean)+ , benchWithFold numElements "mean (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.mean)+ , benchWithFold numElements "mean (entire stream)"+ (Statistics.cumulative Statistics.mean)+ , benchWithFold numElements "mean (Data.Fold)"+ (Fold.mean)++ , benchWithFold+ numElements+ "welfordMean (window size 100)"+ (Ring.slidingWindow 100 Statistics.welfordMean)+ , benchWithFold+ numElements+ "welfordMean (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.welfordMean)+ , benchWithFold+ numElements+ "welfordMean (entire stream)"+ (Statistics.cumulative Statistics.welfordMean)++ , benchWithFold numElements "geometricMean (window size 100)"+ (Ring.slidingWindow 100 Statistics.geometricMean)+ , benchWithFold numElements "geometricMean (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.geometricMean)+ , benchWithFold numElements "geometricMean (entire stream)"+ (Statistics.cumulative Statistics.geometricMean)++ , benchWithFold numElements "harmonicMean (window size 100)"+ (Ring.slidingWindow 100 Statistics.harmonicMean)+ , benchWithFold numElements "harmonicMean (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.harmonicMean)+ , benchWithFold numElements "harmonicMean (entire stream)"+ (Statistics.cumulative Statistics.harmonicMean)++ , benchWithFold numElements "quadraticMean (window size 100)"+ (Ring.slidingWindow 100 Statistics.quadraticMean)+ , benchWithFold numElements "quadraticMean (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.quadraticMean)+ , benchWithFold numElements "quadraticMean (entire stream)"+ (Statistics.cumulative Statistics.quadraticMean)++ , benchWithFold numElements "powerSum 2 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerSum 2))+ , benchWithFold numElements "powerSum 2 (entire stream)"+ (Statistics.cumulative (Statistics.powerSum 2))++ , benchWithFold numElements "rawMoment 2 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerSum 2))+ , benchWithFold numElements "rawMoment 2 (entire stream)"+ (Statistics.cumulative (Statistics.rawMoment 2))++ , benchWithFold numElements "powerMean 1 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerMean 1))+ , benchWithFold numElements "powerMean 2 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerMean 2))+ , benchWithFold numElements "powerMean 10 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerMean 10))++ , benchWithFold numElements "powerMeanFrac (-1) (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerMeanFrac (-1)))+ , benchWithFold numElements "powerMeanFrac 1 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerMeanFrac 1))+ , benchWithFold numElements "powerMeanFrac 2 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerMeanFrac 2))+ , benchWithFold numElements "powerMeanFrac 10 (window size 100)"+ (Ring.slidingWindow 100 (Statistics.powerMeanFrac 10))++ , benchWithFold numElements "ewma (entire stream)"+ (Statistics.ewma 0.5)+ , benchWithFold numElements "ewmaAfterMean (entire stream)"+ (Statistics.ewmaAfterMean 10 0.5)+ , benchWithFold numElements "ewmaRampUpSmoothing (entire stream)"+ (Statistics.ewmaRampUpSmoothing 0.5 0.5)++ , benchWithFold numElements "variance (window size 100)"+ (Ring.slidingWindow 100 (Statistics.variance))+ , benchWithFold numElements "variance (entire stream)"+ (Statistics.cumulative (Statistics.variance))+ -- , benchWithFold numElements "variance (Data.Fold)"+ -- (Fold.variance)++ , benchWithFold numElements "sampleVariance (window size 100)"+ (Ring.slidingWindow 100 (Statistics.sampleVariance))+ , benchWithFold numElements "sampleVariance (entire stream)"+ (Statistics.cumulative (Statistics.sampleVariance))++ , benchWithFold numElements "stdDev (window size 100)"+ (Ring.slidingWindow 100 (Statistics.stdDev))+ , benchWithFold numElements "stdDev (entire stream)"+ (Statistics.cumulative (Statistics.stdDev))+ -- , benchWithFold numElements "stdDev (Data.Fold)"+ -- (Fold.stdDev)++ , benchWithFold numElements "sampleStdDev (window size 100)"+ (Ring.slidingWindow 100 (Statistics.sampleStdDev))+ , benchWithFold numElements "sampleStdDev (entire stream)"+ (Statistics.cumulative (Statistics.sampleStdDev))++ , benchWithFold numElements "stdErrMean (window size 100)"+ (Ring.slidingWindow 100 (Statistics.stdErrMean))+ , benchWithFold numElements "stdErrMean (entire stream)"+ (Statistics.cumulative (Statistics.stdErrMean))++-- These benchmarks take a lot of time/memory with fusion-plugin possibly+-- because of the use of Tee.+#ifndef FUSION_PLUGIN+ , benchWithFold numElements "skewness (window size 100)"+ (Ring.slidingWindow 100 (Statistics.skewness))+ , benchWithFold numElements "skewness (entire stream)"+ (Statistics.cumulative (Statistics.skewness))++ , benchWithFold numElements "kurtosis (window size 100)"+ (Ring.slidingWindow 100 (Statistics.kurtosis))+ , benchWithFold numElements "kurtosis (entire stream)"+ (Statistics.cumulative (Statistics.kurtosis))+#endif+ , benchWithFold numElements "md (window size 100)"+ (Ring.slidingWindowWith 100 Statistics.md)+ ]+ , bgroup+ "scan"+ [ benchWithPostscan numElements "minimum (window size 100)"+ (Ring.slidingWindow 100 Statistics.minimum)+ , benchWithPostscan numElements "minimum (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.minimum)+ , benchWithPostscan numElements "maximum (window size 100)"+ (Ring.slidingWindow 100 Statistics.maximum)+ , benchWithPostscan numElements "maximum (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.maximum)+ , benchWithPostscan numElements "range (window size 100)"+ (Ring.slidingWindow 100 Statistics.range)+ , benchWithPostscan numElements "range (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.range)+ , benchWithPostscan numElements "sum (window size 100)"+ (Ring.slidingWindow 100 Statistics.sum)+ , benchWithPostscan numElements "sum (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.sum)+ , benchWithPostscan numElements "mean (window size 100)"+ (Ring.slidingWindow 100 Statistics.mean)+ , benchWithPostscan numElements "mean (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.mean)+ , benchWithPostscan+ numElements+ "welfordMean (window size 100)"+ (Ring.slidingWindow 100 Statistics.welfordMean)+ , benchWithPostscan+ numElements+ "welfordMean (window size 1000)"+ (Ring.slidingWindow 1000 Statistics.welfordMean)+ , benchWithPostscan+ numElements+ "md (window size 100)"+ (Ring.slidingWindowWith 100 Statistics.md)+ -- XXX These benchmarks measure the cost of creating the array as well,+ -- we can do that outside the benchmark.+ , benchWithResample numElements "Resample"+ , benchWithFoldResamples 316 "FoldResamples 316" Fold.mean+ ]+ ]
+ src/Streamly/Statistics.hs view
@@ -0,0 +1,1155 @@+-- |+-- Module : Streamly.Statistics+-- Copyright : (c) 2020 Composewell Technologies+-- License : Apache-2.0+-- Maintainer : streamly@composewell.com+-- Stability : experimental+-- Portability : GHC+--+-- Statistical measures over a stream of data. All operations use numerically+-- stable floating point arithmetic.+--+-- Measurements can be performed over the entire input stream or on a sliding+-- window of fixed or variable size. Where possible, measures are computed+-- online without buffering the input stream.+--+-- Currently there is no overflow detection.+--+-- References:+--+-- * https://en.wikipedia.org/wiki/Statistics+-- * https://mathworld.wolfram.com/topics/ProbabilityandStatistics.html++-- Resources:+--+-- This may be another useful resource for incremental (non-windowed)+-- computation:+--+-- https://www.researchgate.net/publication/287152055_Incremental_Statistical_Measures+--+-- Sample Statistics+--+-- Terms+--+-- Population: the complete data set from which statistical samples are taken.+--+-- Sample: a subset of the population.+--+-- https://en.wikipedia.org/wiki/Sample_(statistics)+--+-- Estimator:+--+-- Statistical measures can be computed either from the actual population+-- or from samples. Statistical measures computed from the samples provide an+-- estimate of the actual measures of the entire population. Measures computed+-- from samples may not truly reflect the actual measures and may have to be+-- adjusted for biases or errors.+--+-- An "estimator" is a method or function to compute a statistical measure from+-- sampled data. For example, the sample variance is an esitmator of the+-- population variance.+--+-- https://en.wikipedia.org/wiki/Estimator+--+-- Bias:+--+-- The result computed by an estimator may not be centered at the true value as+-- determined by computing the measure for the actual population. Such an+-- estimator is called a biased estimator. For example, notice how+-- 'sampleVariance' is adjusted for bias.+--+-- https://en.wikipedia.org/wiki/Bias_of_an_estimator+--+-- Consistency:+--+-- https://en.wikipedia.org/wiki/Consistent_estimator++{-# LANGUAGE ScopedTypeVariables #-}+module Streamly.Statistics+ (+ -- * Incremental Folds+ -- | Folds of type @Fold m (a, Maybe a) b@ are incremental sliding window+ -- folds. An input of type @(a, Nothing)@ indicates that the input element+ -- @a@ is being inserted in the window without ejecting an old value+ -- increasing the window size by 1. An input of type @(a, Just a)@+ -- indicates that the first element is being inserted in the window and the+ -- second element is being removed from the window, the window size remains+ -- the same. The window size can only increase and never decrease.+ --+ -- You can compute the statistics over the entire stream using sliding+ -- window folds by keeping the second element of the input tuple as+ -- @Nothing@.+ --+ Window.lmap+ , Window.cumulative++ -- * Summary Statistics+ -- | See https://en.wikipedia.org/wiki/Summary_statistics .++ -- ** Sums+ , Window.length+ , Window.sum+ , Window.sumInt+ , Window.powerSum++ -- ** Location+ -- | See https://en.wikipedia.org/wiki/Location_parameter .+ --+ -- See https://en.wikipedia.org/wiki/Central_tendency .+ , minimum+ , maximum+ , rawMoment+ , rawMomentFrac++ -- Pythagorean means (https://en.wikipedia.org/wiki/Pythagorean_means)+ , mean+ , welfordMean+ , geometricMean+ , harmonicMean++ , quadraticMean++ -- Generalized mean+ , powerMean+ , powerMeanFrac++ -- ** Weighted Means+ -- | Exponential Smoothing.+ , ewma+ , ewmaAfterMean+ , ewmaRampUpSmoothing++ -- ** Spread+ -- | Second order central moment is a statistical measure of dispersion.+ -- The \(k\)th moment about the mean (or \(k\)th central moment) is defined+ -- as:+ --+ -- \(\mu_k = \frac{1}{n}\sum_{i=1}^n {(x_{i}-\mu)}^k\)+ --+ -- See https://mathworld.wolfram.com/CentralMoment.html .+ --+ -- See https://en.wikipedia.org/wiki/Statistical_dispersion .+ , range+ , md+ , variance+ , stdDev++ -- ** Shape+ -- | Third and fourth order central moments are a measure of shape.+ --+ -- See https://en.wikipedia.org/wiki/Shape_parameter .+ --+ -- See https://en.wikipedia.org/wiki/Standardized_moment .+ , skewness+ , kurtosis++ -- XXX Move to Statistics.Sample or Statistics.Estimation module?+ -- ** Estimation+ , sampleVariance+ , sampleStdDev+ , stdErrMean++ -- ** Resampling+ , resample+ , foldResamples+ , jackKnifeMean+ , jackKnifeVariance+ , jackKnifeStdDev++ -- ** Probability Distribution+ , frequency+ , frequency'+ , mode++ -- Histograms+ , HistBin (..)+ , binOffsetSize+ , binFromSizeN+ , binFromToN+ , binBoundaries+ , histogram++ -- * Transforms+ , fft+ )+where++import Control.Exception (assert)+import Control.Monad (when)+import Control.Monad.IO.Class (MonadIO(..))+import Data.Bits (Bits(complement, shiftL, shiftR, (.&.), (.|.)))+import Data.Complex (Complex ((:+)))+import Data.Function ((&))+import Data.Functor.Identity (runIdentity, Identity)+import Data.Map.Strict (Map)+import Data.Maybe (fromMaybe)+import Streamly.Data.Array (Array, length, Unbox)+import Streamly.Data.Fold (Tee(..))+import Streamly.Data.Stream (Stream)+import Streamly.Internal.Data.Array.Type (unsafeIndexIO)+import Streamly.Internal.Data.Fold.Type (Fold(..), Step(..))+import Streamly.Internal.Data.Stream.StreamD.Step (Step(..))+import Streamly.Internal.Data.Tuple.Strict (Tuple'(..), Tuple3'(..))+import Streamly.Internal.Data.Unfold.Type (Unfold(..))+import System.Random.MWC (createSystemRandom, uniformRM)++import qualified Data.Map.Strict as Map+import qualified Deque.Strict as Deque+import qualified Streamly.Data.Fold as Fold+import qualified Streamly.Data.Array as Array hiding (read)+import qualified Streamly.Internal.Data.Array as Array (read)+import qualified Streamly.Data.MutArray as MA+import qualified Streamly.Internal.Data.Array.Mut as MA+ (getIndexUnsafe, putIndexUnsafe, unsafeSwapIndices)+import qualified Streamly.Internal.Data.Fold.Window as Window+import qualified Streamly.Data.Stream as Stream++import Prelude hiding (length, sum, minimum, maximum)++-- TODO: Overflow checks. Would be good if we can directly replace the+-- operations with overflow checked operations.+--+-- See https://hackage.haskell.org/package/safe-numeric+-- See https://hackage.haskell.org/package/safeint+--+-- TODO We have many of these functions in Streamly.Data.Fold as well. Need to+-- think about deduplication.++-------------------------------------------------------------------------------+-- Transforms+-------------------------------------------------------------------------------++-- XXX These utility functions can be moved to streamly-numeric++-- | Test if the given integer value is a power of 2.+{-# INLINE isPower2 #-}+isPower2 :: Int -> Bool+isPower2 n = n .&. (n - 1) == 0++-- | Create a power of 2+--+-- Argument must be less than 64 assuming 64-bit Int size.+--+{-# INLINE _power2 #-}+_power2 :: Int -> Int+_power2 n = shiftL 1 n++-- | Create a bit mask with lower n bits 0 and the rest as 1.+--+-- Argument must be less than 64 assuming 64-bit Int size.+--+{-# INLINE maskLowerN #-}+maskLowerN :: Int -> Int+maskLowerN n = complement (shiftL 1 n - 1)++-- | Compute the base 2 logarithm of the given value.+--+-- Assumes the Int size to be 64-bit.+--+{-# INLINE logBase2 #-}+logBase2 :: Int -> Int+logBase2 v0+ | v0 <= 0 = error $ "logBase2: input must be greater than 0 " ++ show v0+ | otherwise = go 32 0 v0++ where++ go !bits !result !v+ | bits == 0 = result+ | v .&. maskLowerN bits /= 0 =+ go (bits `shiftR` 1) (result .|. bits) (v `shiftR` bits)+ | otherwise = go (bits `shiftR` 1) result v++-- Algo translated from the statistics library.+--+-- XXX We can use a wrapper API that takes an array of Double input instead of+-- array of Complex.+--+-- | Compute fast fourier transform of an array of 'Complex' values.+--+-- Array length must be power of 2.+--+{-# INLINE fft #-}+fft :: MonadIO m => MA.MutArray (Complex Double) -> m ()+fft marr+ | isPower2 len = bitReverse 0 0+ | otherwise = error "fft: Array length must be power of 2"++ where++ len = MA.length marr++ halve x = x `shiftR` 1++ twice x = x `shiftL` 1++ inner i j k+ | k <= j = inner i (j - k) (halve k)+ | otherwise = bitReverse (i + 1) (j + k)++ bitReverse i j+ | i == len - 1 = stage 0 1+ | otherwise = do+ when (i < j) $ MA.unsafeSwapIndices i j marr+ inner i j (halve len)++ log2len = logBase2 len++ stage l !l1+ | l == log2len = return ()+ | otherwise = do+ let !l2 = twice l1+ !e = -6.283185307179586/fromIntegral l2+ flight j !a | j == l1 = stage (l + 1) l2+ | otherwise = do+ let butterfly i | i >= len = flight (j + 1) (a + e)+ | otherwise = do+ let i1 = i + l1+ xi1 :+ yi1 <- MA.getIndexUnsafe i1 marr+ let !c = cos a+ !s = sin a+ d = (c * xi1 - s * yi1) :+ (s * xi1 + c * yi1)+ ci <- MA.getIndexUnsafe i marr+ MA.putIndexUnsafe i1 marr (ci - d)+ MA.putIndexUnsafe i marr (ci + d)+ butterfly (i + l2)+ butterfly j+ flight 0 0++-------------------------------------------------------------------------------+-- Location+-------------------------------------------------------------------------------++-- Theoretically, we can approximate minimum in a rolling window by using a+-- 'powerMean' with sufficiently large negative power.+--+-- XXX If we need to know the minimum in the window only once in a while then+-- we can use linear search when it is extracted and not pay the cost all the+-- time.+--+-- | The minimum element in a rolling window.+--+-- For smaller window sizes (< 30) Streamly.Data.Fold.Window.minimum performs+-- better. If you want to compute the minimum of the entire stream Fold.min+-- from streamly package would be much faster.+--+-- /Time/: \(\mathcal{O}(n*w)\) where \(w\) is the window size.+--+{-# INLINE minimum #-}+minimum :: (Monad m, Ord a) => Fold m (a, Maybe a) a+minimum = Fold step initial extract++ where++ initial =+ return+ $ Partial+ $ Tuple3' (0 :: Int) (0 :: Int) (mempty :: Deque.Deque (Int, a))++ step (Tuple3' i w q) (a, ma) =+ case ma of+ Nothing ->+ return+ $ Partial+ $ Tuple3'+ (i + 1)+ (w + 1)+ (headCheck i q (w + 1) & dqloop (i, a))+ Just _ ->+ return+ $ Partial+ $ Tuple3' (i + 1) w (headCheck i q w & dqloop (i,a))++ {-# INLINE headCheck #-}+ headCheck i q w =+ case Deque.uncons q of+ Nothing -> q+ Just (ia', q') ->+ if fst ia' <= i - w+ then q'+ else q++ dqloop ia q =+ case Deque.unsnoc q of+ Nothing -> Deque.snoc ia q+ -- XXX This can be improved for the case of `=`+ Just (ia', q') ->+ if snd ia <= snd ia'+ then dqloop ia q'+ else Deque.snoc ia q++ extract (Tuple3' _ _ q) =+ return+ $ snd+ $ fromMaybe (0, error "minimum: Empty stream")+ $ Deque.head q++-- Theoretically, we can approximate maximum in a rolling window by using a+-- 'powerMean' with sufficiently large positive power.+--+-- | The maximum element in a rolling window.+--+-- For smaller window sizes (< 30) Streamly.Data.Fold.Window.maximum performs+-- better. If you want to compute the maximum of the entire stream+-- Streamly.Data.Fold.maximum from streamly package would be much faster.+--+-- /Time/: \(\mathcal{O}(n*w)\) where \(w\) is the window size.+--+{-# INLINE maximum #-}+maximum :: (Monad m, Ord a) => Fold m (a, Maybe a) a+maximum = Fold step initial extract++ where++ initial =+ return+ $ Partial+ $ Tuple3' (0 :: Int) (0 :: Int) (mempty :: Deque.Deque (Int, a))++ step (Tuple3' i w q) (a, ma) =+ case ma of+ Nothing ->+ return+ $ Partial+ $ Tuple3'+ (i + 1)+ (w + 1)+ (headCheck i q (w + 1) & dqloop (i, a))+ Just _ ->+ return+ $ Partial+ $ Tuple3' (i + 1) w (headCheck i q w & dqloop (i,a))++ {-# INLINE headCheck #-}+ headCheck i q w =+ case Deque.uncons q of+ Nothing -> q+ Just (ia', q') ->+ if fst ia' <= i - w+ then q'+ else q++ dqloop ia q =+ case Deque.unsnoc q of+ Nothing -> Deque.snoc ia q+ -- XXX This can be improved for the case of `=`+ Just (ia', q') ->+ if snd ia >= snd ia'+ then dqloop ia q'+ else Deque.snoc ia q++ extract (Tuple3' _ _ q) =+ return+ $ snd+ $ fromMaybe (0, error "maximum: Empty stream")+ $ Deque.head q++-------------------------------------------------------------------------------+-- Mean+-------------------------------------------------------------------------------++-- | Arithmetic mean of elements in a sliding window:+--+-- \(\mu = \frac{\sum_{i=1}^n x_{i}}{n}\)+--+-- This is also known as the Simple Moving Average (SMA) when used in the+-- sliding window and Cumulative Moving Avergae (CMA) when used on the entire+-- stream.+--+-- Mean is the same as the first raw moment.+--+-- \(\mu = \mu'_1\)+--+-- >>> mean = rawMoment 1+-- >>> mean = powerMean 1+-- >>> mean = Fold.teeWith (/) sum length+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE mean #-}+mean :: forall m a. (Monad m, Fractional a) => Fold m (a, Maybe a) a+mean = Window.mean++-- | Recompute mean from old mean when an item is removed from the sample.+{-# INLINE _meanSubtract #-}+_meanSubtract :: Fractional a => Int -> a -> a -> a+_meanSubtract n oldMean oldItem =+ let delta = (oldItem - oldMean) / fromIntegral (n - 1)+ in oldMean - delta++-- | Recompute mean from old mean when an item is added to the sample.+{-# INLINE meanAdd #-}+meanAdd :: Fractional a => Int -> a -> a -> a+meanAdd n oldMean newItem =+ let delta = (newItem - oldMean) / fromIntegral (n + 1)+ in oldMean + delta++-- We do not carry rounding errors, therefore, this would be less numerically+-- stable than the kbn mean.+--+-- | Recompute mean from old mean when an item in the sample is replaced.+{-# INLINE meanReplace #-}+meanReplace :: Fractional a => Int -> a -> a -> a -> a+meanReplace n oldMean oldItem newItem =+ let n1 = fromIntegral n+ -- Compute two deltas instead of a single (newItem - oldItem) because+ -- the latter would be too small causing rounding errors.+ delta1 = (newItem - oldMean) / n1+ delta2 = (oldItem - oldMean) / n1+ in (oldMean + delta1) - delta2++-- | Same as 'mean' but uses Welford's algorithm to compute the mean+-- incrementally.+--+-- It maintains a running mean instead of a running sum and adjusts the mean+-- based on a new value. This is slower than 'mean' because of using the+-- division operation on each step and it is numerically unstable (as of now).+-- The advantage over 'mean' could be no overflow if the numbers are large,+-- because we do not maintain a sum, but that is a highly unlikely corner case.+--+-- /Internal/+{-# INLINE welfordMean #-}+welfordMean :: forall m a. (Monad m, Fractional a) => Fold m (a, Maybe a) a+welfordMean = Fold step initial extract++ where++ initial =+ return+ $ Partial+ $ Tuple'+ (0 :: a) -- running mean+ (0 :: Int) -- count of items in the window++ step (Tuple' oldMean w) (new, mOld) =+ return+ $ Partial+ $ case mOld of+ Nothing -> Tuple' (meanAdd w oldMean new) (w + 1)+ Just old -> Tuple' (meanReplace w oldMean old new) w++ extract (Tuple' x _) = return x++-------------------------------------------------------------------------------+-- Moments+-------------------------------------------------------------------------------++-- XXX We may have chances of overflow if the powers are high or the numbers+-- are large. A limited mitigation could be to use welford style avg+-- computation. Do we need an overflow detection?+--+-- | Raw moment is the moment about 0. The \(k\)th raw moment is defined as:+--+-- \(\mu'_k = \frac{\sum_{i=1}^n x_{i}^k}{n}\)+--+-- >>> rawMoment k = Fold.teeWith (/) (powerSum p) length+--+-- See https://en.wikipedia.org/wiki/Moment_(mathematics) .+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE rawMoment #-}+rawMoment :: (Monad m, Fractional a) => Int -> Fold m (a, Maybe a) a+rawMoment k = Fold.teeWith (/) (Window.powerSum k) Window.length++-- | Like 'rawMoment' but powers can be negative or fractional. This is+-- slower than 'rawMoment' for positive intergal powers.+--+-- >>> rawMomentFrac p = Fold.teeWith (/) (powerSumFrac p) length+--+{-# INLINE rawMomentFrac #-}+rawMomentFrac :: (Monad m, Floating a) => a -> Fold m (a, Maybe a) a+rawMomentFrac k = Fold.teeWith (/) (Window.powerSumFrac k) Window.length++-- XXX Overflow can happen when large powers or large numbers are used. We can+-- keep a running mean instead of running sum but that won't mitigate the+-- overflow possibility by much. The overflow can still happen when computing+-- the mean incrementally.++-- | The \(k\)th power mean of numbers \(x_1, x_2, \ldots, x_n\) is:+--+-- \(M_k = \left( \frac{1}{n} \sum_{i=1}^n x_i^k \right)^{\frac{1}{k}}\)+--+-- \(powerMean(k) = (rawMoment(k))^\frac{1}{k}\)+--+-- >>> powerMean k = (** (1 / fromIntegral k)) <$> rawMoment k+--+-- All other means can be expressed in terms of power mean. It is also known as+-- the generalized mean.+--+-- See https://en.wikipedia.org/wiki/Generalized_mean+--+{-# INLINE powerMean #-}+powerMean :: (Monad m, Floating a) => Int -> Fold m (a, Maybe a) a+powerMean k = (** (1 / fromIntegral k)) <$> rawMoment k++-- | Like 'powerMean' but powers can be negative or fractional. This is+-- slower than 'powerMean' for positive intergal powers.+--+-- >>> powerMeanFrac k = (** (1 / k)) <$> rawMomentFrac k+--+{-# INLINE powerMeanFrac #-}+powerMeanFrac :: (Monad m, Floating a) => a -> Fold m (a, Maybe a) a+powerMeanFrac k = (** (1 / k)) <$> rawMomentFrac k++-- | The harmonic mean of the positive numbers \(x_1, x_2, \ldots, x_n\) is+-- defined as:+--+-- \(HM = \frac{n}{\frac1{x_1} + \frac1{x_2} + \cdots + \frac1{x_n}}\)+--+-- \(HM = \left(\frac{\sum\limits_{i=1}^n x_i^{-1}}{n}\right)^{-1}\)+--+-- >>> harmonicMean = Fold.teeWith (/) length (lmap recip sum)+-- >>> harmonicMean = powerMeanFrac (-1)+--+-- See https://en.wikipedia.org/wiki/Harmonic_mean .+--+{-# INLINE harmonicMean #-}+harmonicMean :: (Monad m, Fractional a) => Fold m (a, Maybe a) a+harmonicMean = Fold.teeWith (/) Window.length (Window.lmap recip Window.sum)++-- | Geometric mean, defined as:+--+-- \(GM = \sqrt[n]{x_1 x_2 \cdots x_n}\)+--+-- \(GM = \left(\prod_{i=1}^n x_i\right)^\frac{1}{n}\)+--+-- or, equivalently, as the arithmetic mean in log space:+--+-- \(GM = e ^{{\frac{\sum_{i=1}^{n}\ln a_i}{n}}}\)+--+-- >>> geometricMean = exp <$> lmap log mean+--+-- See https://en.wikipedia.org/wiki/Geometric_mean .+{-# INLINE geometricMean #-}+geometricMean :: (Monad m, Floating a) => Fold m (a, Maybe a) a+geometricMean = exp <$> Window.lmap log mean++-- | The quadratic mean or root mean square (rms) of the numbers+-- \(x_1, x_2, \ldots, x_n\) is defined as:+--+-- \(RMS = \sqrt{ \frac{1}{n} \left( x_1^2 + x_2^2 + \cdots + x_n^2 \right) }.\)+--+-- >>> quadraticMean = powerMean 2+--+-- See https://en.wikipedia.org/wiki/Root_mean_square .+--+{-# INLINE quadraticMean #-}+quadraticMean :: (Monad m, Floating a) => Fold m (a, Maybe a) a+quadraticMean = powerMean 2++-------------------------------------------------------------------------------+-- Weighted Means+-------------------------------------------------------------------------------++-- XXX Is this numerically stable? We can use the kbn summation here.+-- | ewmaStep smoothing-factor old-value new-value+{-# INLINE ewmaStep #-}+ewmaStep :: Double -> Double -> Double -> Double+ewmaStep k x0 x1 = (1 - k) * x0 + k * x1++-- XXX Compute this in a sliding window?+--+-- | @ewma smoothingFactor@.+--+-- @ewma@ of an empty stream is 0.+--+-- Exponential weighted moving average, \(s_n\), of \(n\) values,+-- \(x_1,\ldots,x_n\), is defined recursively as:+--+-- \(\begin{align} s_0& = x_0\\ s_n & = \alpha x_{n} + (1-\alpha)s_{n-1},\quad n>0 \end{align}\)+--+-- If we expand the recursive term it becomes an exponential series:+--+-- \(s_n = \alpha \left[x_n + (1-\alpha)x_{n-1} + (1-\alpha)^2 x_{n-2} + \cdots + (1-\alpha)^{n-1} x_1 \right] + (1-\alpha)^n x_0\)+--+-- where \(\alpha\), the smoothing factor, is in the range \(0 <\alpha < 1\).+-- More the value of \(\alpha\), the more weight is given to newer values. As+-- a special case if it is 0 then the weighted sum would always be the same as+-- the oldest value, if it is 1 then the sum would always be the same as the+-- newest value.+--+-- See https://en.wikipedia.org/wiki/Moving_average+--+-- See https://en.wikipedia.org/wiki/Exponential_smoothing+--+{-# INLINE ewma #-}+ewma :: Monad m => Double -> Fold m Double Double+ewma k = extract <$> Fold.foldl' step (Tuple' 0 1)++ where++ step (Tuple' x0 k1) x = Tuple' (ewmaStep k1 x0 x) k++ extract (Tuple' x _) = x++-- XXX It can perhaps perform better if implemented as a custom fold?+--+-- | @ewma n k@ is like 'ewma' but uses the mean of the first @n@ values and+-- then uses that as the initial value for the @ewma@ of the rest of the+-- values.+--+-- This can be used to reduce the effect of volatility of the initial value+-- when k is too small.+--+{-# INLINE ewmaAfterMean #-}+ewmaAfterMean :: Monad m => Int -> Double -> Fold m Double Double+ewmaAfterMean n k =+ Fold.concatMap (\i -> (Fold.foldl' (ewmaStep k) i)) (Fold.take n Fold.mean)++-- | @ewma n k@ is like 'ewma' but uses 1 as the initial smoothing factor and+-- then exponentially smooths it to @k@ using @n@ as the smoothing factor.+--+-- This is significantly faster than 'ewmaAfterMean'.+--+{-# INLINE ewmaRampUpSmoothing #-}+ewmaRampUpSmoothing :: Monad m => Double -> Double -> Fold m Double Double+ewmaRampUpSmoothing n k1 = extract <$> Fold.foldl' step initial++ where++ initial = Tuple' 0 1++ step (Tuple' x0 k0) x1 =+ let x = ewmaStep k0 x0 x1+ k = ewmaStep n k0 k1+ in Tuple' x k++ extract (Tuple' x _) = x++-------------------------------------------------------------------------------+-- Spread/Dispersion+-------------------------------------------------------------------------------++-- | The difference between the maximum and minimum elements of a rolling window.+--+-- >>> range = Fold.teeWith (-) maximum minimum+--+-- If you want to compute the range of the entire stream @Fold.teeWith (-)+-- Fold.maximum Fold.minimum@ from the streamly package would be much faster.+--+-- /Space/: \(\mathcal{O}(n)\) where @n@ is the window size.+--+-- /Time/: \(\mathcal{O}(n*w)\) where \(w\) is the window size.+--+{-# INLINE range #-}+range :: (Monad m, Num a, Ord a) => Fold m (a, Maybe a) a+range = Fold.teeWith (-) maximum minimum++-- | @md n@ computes the mean absolute deviation (or mean deviation) in a+-- sliding window of last @n@ elements in the stream.+--+-- The mean absolute deviation of the numbers \(x_1, x_2, \ldots, x_n\) is:+--+-- \(MD = \frac{1}{n}\sum_{i=1}^n |x_i-\mu|\)+--+-- Note: It is expensive to compute MD in a sliding window. We need to+-- maintain a ring buffer of last n elements and maintain a running mean, when+-- the result is extracted we need to compute the difference of all elements+-- from the mean and get the average. Using standard deviation may be+-- computationally cheaper.+--+-- See https://en.wikipedia.org/wiki/Average_absolute_deviation .+--+-- /Pre-release/+{-# INLINE md #-}+md :: MonadIO m => Fold m ((Double, Maybe Double), m (MA.MutArray Double)) Double+md =+ Fold.rmapM computeMD+ $ Fold.tee (Fold.lmap fst mean) (Fold.lmap snd Fold.latest)++ where++ computeMD (mn, rng) =+ case rng of+ Just action -> do+ arr <- action+ Stream.fold Fold.mean+ $ fmap (\a -> abs (mn - a))+ $ Stream.unfold MA.reader arr+ Nothing -> return 0.0++-- | The variance \(\sigma^2\) of a population of \(n\) equally likely values+-- is defined as the average of the squares of deviations from the mean+-- \(\mu\). In other words, second moment about the mean:+--+-- \(\sigma^2 = \frac{1}{n}\sum_{i=1}^n {(x_{i}-\mu)}^2\)+--+-- \(\sigma^2 = rawMoment(2) - \mu^2\)+--+-- \(\mu_2 = -(\mu'_1)^2 + \mu'_2\)+--+-- Note that the variance would be biased if applied to estimate the population+-- variance from a sample of the population. See 'sampleVariance'.+--+-- See https://en.wikipedia.org/wiki/Variance.+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE variance #-}+variance :: (Monad m, Fractional a) => Fold m (a, Maybe a) a+variance = Fold.teeWith (\p2 m -> p2 - m ^ (2 :: Int)) (rawMoment 2) mean++-- | Standard deviation \(\sigma\) is the square root of 'variance'.+--+-- This is the population standard deviation or uncorrected sample standard+-- deviation.+--+-- >>> stdDev = sqrt <$> variance+--+-- See https://en.wikipedia.org/wiki/Standard_deviation .+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE stdDev #-}+stdDev :: (Monad m, Floating a) => Fold m (a, Maybe a) a+stdDev = sqrt <$> variance++-- | Skewness \(\gamma\) is the standardized third central moment defined as:+--+-- \(\tilde{\mu}_3 = \frac{\mu_3}{\sigma^3}\)+--+-- The third central moment can be computed in terms of raw moments:+--+-- \(\mu_3 = 2(\mu'_1)^3 - 3\mu'_1\mu'_2 + \mu'_3\)+--+-- Substituting \(\mu'_1 = \mu\), and \(\mu'_2 = \mu^2 + \sigma^2\):+--+-- \(\mu_3 = -\mu^3 - 3\mu\sigma^2 + \mu'_3\)+--+-- Skewness is a measure of symmetry of the probability distribution. It is 0+-- for a symmetric distribution, negative for a distribution that is skewed+-- towards left, positive for a distribution skewed towards right.+--+-- For a normal like distribution the median can be found around+-- \(\mu - \frac{\gamma\sigma}{6}\) and the mode can be found around+-- \(\mu - \frac{\gamma \sigma}{2}\).+--+-- See https://en.wikipedia.org/wiki/Skewness .+--+{-# INLINE skewness #-}+skewness :: (Monad m, Floating a) => Fold m (a, Maybe a) a+skewness =+ unTee+ $ (\rm3 sd mu ->+ rm3 / sd ^ (3 :: Int) - 3 * (mu / sd) - (mu / sd) ^ (3 :: Int)+ )+ <$> Tee (rawMoment 3)+ <*> Tee stdDev+ <*> Tee mean++-- XXX We can compute the 2nd, 3rd, 4th raw moments by repeatedly multiplying+-- instead of computing the powers every time.+--+-- | Kurtosis \(\kappa\) is the standardized fourth central moment, defined as:+--+-- \(\tilde{\mu}_4 = \frac{\mu_4}{\sigma^4}\)+--+-- The fourth central moment can be computed in terms of raw moments:+--+-- \(\mu_4 = -3(\mu'_1)^4 + 6(\mu'_1)^2\mu'_2 - 4\mu'_1\mu'_3\ + \mu'_4\)+--+-- Substituting \(\mu'_1 = \mu\), and \(\mu'_2 = \mu^2 + \sigma^2\):+--+-- \(\mu_4 = 3\mu^4 + 6\mu^2\sigma^2 - 4\mu\mu'_3 + \mu'_4\)+--+-- It is always non-negative. It is 0 for a point distribution, low for light+-- tailed (platykurtic) distributions and high for heavy tailed (leptokurtic)+-- distributions.+--+-- \(\kappa >= \gamma^2 + 1\)+--+-- For a normal distribution \(\kappa = 3\sigma^4\).+--+-- See https://en.wikipedia.org/wiki/Kurtosis .+--+{-# INLINE kurtosis #-}+kurtosis :: (Monad m, Floating a) => Fold m (a, Maybe a) a+kurtosis =+ unTee+ $ (\rm4 rm3 sd mu ->+ ( 3 * mu ^ (4 :: Int)+ + 6 * mu ^ (2 :: Int) * sd ^ (2 :: Int)+ - 4 * mu * rm3+ + rm4) / (sd ^ (4 :: Int))+ )+ <$> Tee (rawMoment 4)+ <*> Tee (rawMoment 3)+ <*> Tee stdDev+ <*> Tee mean++-------------------------------------------------------------------------------+-- Estimation+-------------------------------------------------------------------------------++-- | Unbiased sample variance i.e. the variance of a sample corrected to+-- better estimate the variance of the population, defined as:+--+-- \(s^2 = \frac{1}{n - 1}\sum_{i=1}^n {(x_{i}-\mu)}^2\)+--+-- \(s^2 = \frac{n}{n - 1} \times \sigma^2\).+--+-- See https://en.wikipedia.org/wiki/Bessel%27s_correction.+--+{-# INLINE sampleVariance #-}+sampleVariance :: (Monad m, Fractional a) => Fold m (a, Maybe a) a+sampleVariance = Fold.teeWith (\n s2 -> n * s2 / (n - 1)) Window.length variance++-- | Sample standard deviation:+--+-- \(s = \sqrt{sampleVariance}\)+--+-- >>> sampleStdDev = sqrt <$> sampleVariance+--+-- See https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation+-- .+--+{-# INLINE sampleStdDev #-}+sampleStdDev :: (Monad m, Floating a) => Fold m (a, Maybe a) a+sampleStdDev = sqrt <$> sampleVariance++-- | Standard error of the sample mean (SEM), defined as:+--+-- \( SEM = \frac{sampleStdDev}{\sqrt{n}} \)+--+-- See https://en.wikipedia.org/wiki/Standard_error .+--+-- /Space/: \(\mathcal{O}(1)\)+--+-- /Time/: \(\mathcal{O}(n)\)+{-# INLINE stdErrMean #-}+stdErrMean :: (Monad m, Floating a) => Fold m (a, Maybe a) a+stdErrMean = Fold.teeWith (\sd n -> sd / sqrt n) sampleStdDev Window.length++-------------------------------------------------------------------------------+-- Resampling+-------------------------------------------------------------------------------++{-# INLINE foldArray #-}+foldArray :: Unbox a => Fold Identity a b -> Array a -> b+foldArray f = runIdentity . Stream.fold f . Array.read++-- XXX Is this numerically stable? Should we keep the rounding error in the sum+-- and take it into account when subtracting?+--+-- | Given an array of @n@ items, compute mean of @(n - 1)@ items at a time,+-- producing a stream of all possible mean values omitting a different item+-- every time.+--+{-# INLINE jackKnifeMean #-}+jackKnifeMean :: (Monad m, Fractional a, Unbox a) => Array a -> Stream m a+jackKnifeMean arr = do+ let len = fromIntegral (length arr - 1)+ s = foldArray Fold.sum arr+ in fmap (\b -> (s - b) / len) $ Array.read arr++-- | Given an array of @n@ items, compute variance of @(n - 1)@ items at a time,+-- producing a stream of all possible variance values omitting a different item+-- every time.+--+{-# INLINE jackKnifeVariance #-}+jackKnifeVariance :: (Monad m, Fractional a, Unbox a) =>+ Array a -> Stream m a+jackKnifeVariance arr = do+ let len = fromIntegral $ length arr - 1+ foldSums (s, s2) x = (s + x, s2 + x ^ (2 :: Int))+ (sum, sum2) = foldArray (Fold.foldl' foldSums (0.0, 0.0)) arr+ var x = (sum2 - x ^ (2 :: Int)) / len - ((sum - x) / len) ^ (2::Int)+ in fmap var $ Array.read arr++-- | Standard deviation computed from 'jackKnifeVariance'.+--+{-# INLINE jackKnifeStdDev #-}+jackKnifeStdDev :: (Monad m, Unbox a, Floating a) =>+ Array a -> Stream m a+jackKnifeStdDev = fmap sqrt . jackKnifeVariance++-- XXX This can be made more modular if the replicateM unfold can take count+-- from the seed.+--+-- | Randomly select elements from an array, with replacement, producing+-- a stream of the same size as the original array.+{-# INLINE resample #-}+resample :: (MonadIO m, Unbox a) => Unfold m (Array a) a+resample = Unfold step inject++ where++ inject arr = liftIO $ do+ g <- createSystemRandom+ return $ (g, arr, length arr, 0)++ chooseOne g arr len = do+ i <- uniformRM (0, len - 1) g+ unsafeIndexIO i arr++ step (g, arr, len, idx) = liftIO $ do+ if idx >= len+ then return Stop+ else do+ e <- chooseOne g arr len+ return $ Yield e (g, arr, len, idx + 1)++-- XXX Use concurrent combinators++-- | Resample an array multiple times and run the supplied fold on each+-- resampled stream, producing a stream of fold results. The fold is usually an+-- estimator fold.+{-# INLINE foldResamples #-}+foldResamples :: (MonadIO m, Unbox a) =>+ Int -- ^ Number of resamples to compute.+ -> Array a -- ^ Original sample.+ -> Fold m a b -- ^ Estimator fold+ -> Stream m b+foldResamples n arr fld =+ Stream.sequence+ $ Stream.replicate n (Stream.fold fld $ Stream.unfold resample arr)++-------------------------------------------------------------------------------+-- Probability Distribution+-------------------------------------------------------------------------------++-- XXX We can use a Windowed classifyWith operation, that will allow us to+-- express windowed frequency, mode, histograms etc idiomatically.++-- | Count the frequency of elements in a sliding window.+--+-- >>> input = Stream.fromList [1,1,3,4,4::Int]+-- >>> f = Ring.slidingWindow 4 Statistics.frequency+-- >>> Stream.fold f input+-- fromList [(1,1),(3,1),(4,2)]+--+{-# INLINE frequency #-}+frequency :: (Monad m, Ord a) => Fold m (a, Maybe a) (Map a Int)+frequency = Fold.foldl' step Map.empty++ where++ decrement v =+ if v == 1+ then Nothing+ else Just (v - 1)++ step refCountMap (new, mOld) =+ let m1 = Map.insertWith (+) new 1 refCountMap+ in case mOld of+ Just k -> Map.update decrement k m1+ Nothing -> m1++-- XXX Check if the performance of window frequency is the same as this in the+-- full case, if so remove this.+-- XXX This is available in the streamly package as well.++-- | Determine the frequency of each element in the stream.+--+{-# INLINE frequency' #-}+frequency' :: (Monad m, Ord a) => Fold m a (Map a Int)+frequency' = Fold.toMap id Fold.length++-- | Find out the most frequently ocurring element in the stream and its+-- frequency.+--+{-# INLINE mode #-}+mode :: (Monad m, Ord a) => Fold m a (Maybe (a, Int))+mode = Fold.rmapM findMax frequency'++ where++ fmax k v Nothing = Just (k, v)+ fmax k v old@(Just (_, v1))+ | v > v1 = Just (k, v)+ | otherwise = old++ findMax = return . Map.foldrWithKey fmax Nothing++-------------------------------------------------------------------------------+-- Histograms+-------------------------------------------------------------------------------++-- | @binOffsetSize offset binSize input@. Given an integral input value,+-- return its bin index provided that each bin contains @binSize@ items and the+-- bins are aligned such that the 0 index bin starts at @offset@ from 0. If+-- offset = 0 then the bin with index 0 would have values from 0 to binSize -+-- 1.+--+-- This API does not put a bound on the number of bins, therefore, the number+-- of bins could be potentially large depending on the range of values.+--+{-# INLINE binOffsetSize #-}+binOffsetSize :: Integral a => a -> a -> a -> a+binOffsetSize offset binSize x = (x - offset) `div` binSize++data HistBin a = BelowRange | InRange a | AboveRange deriving (Eq, Show)++instance (Ord a) => Ord (HistBin a) where+ compare BelowRange BelowRange = EQ+ compare BelowRange (InRange _) = LT+ compare BelowRange AboveRange = LT++ compare (InRange _) BelowRange = GT+ compare (InRange x) (InRange y)= x `compare` y+ compare (InRange _) AboveRange = LT++ compare AboveRange BelowRange = GT+ compare AboveRange (InRange _) = GT+ compare AboveRange AboveRange = EQ++-- | @binFromSizeN low binSize nbins input@. Classify @input@ into bins+-- specified by a @low@ limit, @binSize@ and @nbins@. Inputs below the lower+-- limit are classified into 'BelowRange' and inputs above the highest bin are+-- classified into 'AboveRange'. 'InRange' inputs are classified into bins+-- starting from bin index 0.+--+{-# INLINE binFromSizeN #-}+binFromSizeN :: Integral a => a -> a -> a -> a -> HistBin a+binFromSizeN low binSize nbins x =+ let high = low + binSize * nbins+ in if x < low+ then BelowRange+ else if x >= high+ then AboveRange+ else InRange ((x - low) `div` binSize)++-- | @binFromToN low high nbins input@. Like @binFromSizeN@ except that a range+-- of lower and higher limit is specified. @binSize@ is computed using the+-- range and @nbins@. @nbins@ is rounded to the range @0 < nbins < (high - low+-- + 1)@. @high >= low@ must hold.+--+{-# INLINE binFromToN #-}+binFromToN :: Integral a => a -> a -> a -> a -> HistBin a+binFromToN low high n x =+ let count = high - low + 1+ n1 = max n 1+ n2 = min n1 count+ binSize = count `div` n2+ nbins =+ if binSize * n2 < count+ then n2 + 1+ else n2+ in assert (high >= low) (binFromSizeN low binSize nbins x)++-- Use binary search to find the bin+--+-- | Classify an input value to bins using the bin boundaries specified in an+-- array.+--+-- /Unimplemented/+--+{-# INLINE binBoundaries #-}+binBoundaries :: -- Integral a =>+ Array.Array a -> a -> HistBin a+binBoundaries = undefined++-- | Given a bin classifier function and a stream of values, generate a+-- histogram map from indices of bins to the number of items in the bin.+--+-- >>> Stream.fold (histogram (binOffsetSize 0 3)) $ Stream.fromList [1..15]+-- fromList [(0,2),(1,3),(2,3),(3,3),(4,3),(5,1)]+--+{-# INLINE histogram #-}+histogram :: (Monad m, Ord k) => (a -> k) -> Fold m a (Map k Int)+histogram bin = Fold.toMap bin Fold.length
+ streamly-statistics.cabal view
@@ -0,0 +1,157 @@+cabal-version: 2.4+name: streamly-statistics+version: 0.1.0+synopsis:+ Statistical measures for finite or infinite data streams.+description:+ Statistical measures for finite or infinite data streams.+ .+ All operations use numerically stable floating point arithmetic.+ Measurements can be performed over the entire input stream or on a sliding+ window of fixed or variable size. Where possible, measures are computed+ online without buffering the input stream.+ .+ Includes\:+ .+ * Summary: length, sum, powerSum+ * Location: minimum, maximum, rawMoments, means, exponential smoothing+ * Spread: range, variance, deviations+ * Shape: skewness, kurtosis+ * Sample statistics, resampling+ * Probablity distribution: frequency, mode, histograms+ * Transforms: Fast fourier transform+homepage: https://streamly.composewell.com+bug-reports: https://github.com/composewell/streamly-statistics/issues+license: Apache-2.0+license-file: LICENSE+tested-with:+ GHC==8.10.7+ , GHC==9.0.2+ , GHC==9.2.2+ , GHC==9.4.4+author: Composewell Technologies+maintainer: streamly@composewell.com+copyright: 2019 Composewell Technologies+category: Streamly, Statistics++extra-source-files:+ CHANGELOG.md+ , NOTICE+ , README.md++source-repository head+ type: git+ location: https://github.com/composewell/streamly-statistics++flag fusion-plugin+ description: Use fusion plugin for benchmarks+ manual: True+ default: True++common default-extensions+ default-extensions:+ BangPatterns+ CApiFFI+ CPP+ ConstraintKinds+ DeriveDataTypeable+ DeriveGeneric+ DeriveTraversable+ ExistentialQuantification+ FlexibleContexts+ FlexibleInstances+ GeneralizedNewtypeDeriving+ InstanceSigs+ KindSignatures+ LambdaCase+ MagicHash+ MultiParamTypeClasses+ PatternSynonyms+ RankNTypes+ RecordWildCards+ ScopedTypeVariables+ TupleSections+ TypeApplications+ TypeFamilies+ ViewPatterns++ -- MonoLocalBinds, enabled by TypeFamilies, causes performance+ -- regressions. Disable it. This must come after TypeFamilies,+ -- otherwise TypeFamilies will enable it again.+ NoMonoLocalBinds++ -- UndecidableInstances -- Does not show any perf impact+ -- UnboxedTuples -- interferes with (#.)++common compile-options+ default-language: Haskell2010+ ghc-options: -Wall+ -Wcompat+ -Wunrecognised-warning-flags+ -Widentities+ -Wincomplete-record-updates+ -Wincomplete-uni-patterns+ -Wredundant-constraints+ -Wnoncanonical-monad-instances+ -Rghc-timing++common optimization-options+ ghc-options: -O2+ -fdicts-strict+ -fspec-constr-recursive=16+ -fmax-worker-args=16+ -fsimpl-tick-factor=200++common ghc-options+ import: default-extensions, compile-options, optimization-options++library+ import: ghc-options+ exposed-modules: Streamly.Statistics+ build-depends: base >= 4.9 && < 5+ , streamly-core == 0.1.0+ , containers >= 0.5 && < 0.7+ , random >= 1.2 && < 1.3+ , mwc-random >= 0.15 && < 0.16+ , deque >= 0.4.4 && < 0.4.5+ hs-source-dirs: src++test-suite test+ import: ghc-options+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: Main.hs+ build-depends: streamly-statistics+ , streamly-core == 0.1.0+ , base >= 4.9 && < 5+ , QuickCheck >= 2.10 && < 2.15+ , hspec >= 2.0 && < 3+ , hspec-core >= 2.0 && < 3+ , random >= 1.0.0 && < 2+ , containers >= 0.5 && < 0.7+ -- XXX Should remove these dependencies+ , vector >= 0.11 && < 0.14+ , statistics >= 0.15 && < 0.17++benchmark benchmark+ import: ghc-options+ ghc-options: +RTS -M3G -RTS+ type: exitcode-stdio-1.0+ hs-source-dirs: benchmark+ main-is: Main.hs+ build-depends: streamly-statistics+ , streamly-core == 0.1.0+ , base >= 4.9 && < 5+ , random >= 1.0.0 && < 2+ , deepseq >= 1.4.1 && < 1.5+ , tasty-bench >= 0.3 && < 0.4+ , tasty >= 1.4.1 && < 1.5+ mixins: tasty-bench+ (Test.Tasty.Bench as Gauge+ ,Test.Tasty.Bench as Gauge.Main+ )+ if flag(fusion-plugin) && !impl(ghcjs) && !impl(ghc < 8.6)+ cpp-options: -DFUSION_PLUGIN+ ghc-options: -fplugin Fusion.Plugin+ build-depends:+ fusion-plugin >= 0.2 && < 0.3
+ test/Main.hs view
@@ -0,0 +1,305 @@+{-# LANGUAGE TupleSections #-}++import Control.Monad.IO.Class (liftIO)+import Data.Complex (Complex ((:+)))+import Data.Functor.Classes (liftEq2)+import Streamly.Data.Array (Unbox)+import Streamly.Data.Stream (Stream)+import Test.Hspec.Core.Spec (SpecM)+import Test.Hspec.QuickCheck (prop)+import Test.QuickCheck+ (elements, chooseInt, choose, forAll, Property, vectorOf)+import Test.QuickCheck.Monadic (monadicIO, assert)++import qualified Data.Map.Strict as Map+import qualified Data.Set as Set+import qualified Data.Vector as V+import qualified Statistics.Sample.Powers as STAT+import qualified Statistics.Transform as STAT+import qualified Streamly.Data.Array as Array+import qualified Streamly.Data.Fold as Fold+import qualified Streamly.Data.MutArray as MA+import qualified Streamly.Internal.Data.Ring.Unboxed as Ring+import qualified Streamly.Data.Stream as Stream+import qualified Streamly.Data.Stream as S++import Prelude hiding (sum, maximum, minimum)++import Streamly.Statistics+import Test.Hspec++tolerance :: Double+tolerance = 0.00001++validate :: Double -> Bool+validate delta = delta < tolerance++jackKnifeInput :: [Double]+jackKnifeInput = [1.0::Double, 2.0, 3.0, 4.0]++jackMeanRes :: [Double]+jackMeanRes = [3.0, 2.6666666666666665, 2.3333333333333335, 2.0]++jackVarianceRes :: [Double]+jackVarianceRes =+ [ 0.6666666666666661+ , 1.5555555555555554+ , 1.5555555555555545+ , 0.666666666666667+ ]++jackStdDevRes :: [Double]+jackStdDevRes =+ [ 0.8164965809277257+ , 1.247219128924647+ , 1.2472191289246466+ , 0.8164965809277263+ ]++testDistributions+ :: (STAT.Powers -> Double)+ -> Fold.Fold IO (Double, Maybe Double) Double+ -> Property+testDistributions func fld =+ forAll (chooseInt (1, 1000)) $ \list_length ->+ forAll (vectorOf list_length (choose (-50.0 :: Double, 100.0)))+ $ \ls ->+ monadicIO $ do+ let var2 = func . STAT.powers 2 $ V.fromList ls+ strm = S.fromList ls+ var1 <-+ liftIO $ S.fold (Ring.slidingWindow list_length fld) strm+ assert (validate $ abs (var1 - var2))++testVariance :: Property+testVariance = testDistributions STAT.variance variance++testStdDev :: Property+testStdDev = testDistributions STAT.stdDev stdDev++testFuncMD ::+ Fold.Fold IO ((Double, Maybe Double), IO (MA.MutArray Double)) Double -> Spec+testFuncMD f = do+ let c = S.fromList [10.0, 11.0, 12.0, 14.0]+ a1 <- runIO $ S.fold (Ring.slidingWindowWith 2 f) c+ a2 <- runIO $ S.fold (Ring.slidingWindowWith 3 f) c+ a3 <- runIO $ S.fold (Ring.slidingWindowWith 4 f) c+ it ("MD should be 1.0 , 1.1111111111111114 , 1.25 but actual is "+ ++ show a1 ++ " " ++ show a2 ++ " " ++ show a3)+ ( validate (abs (a1 - 1.0))+ && validate (abs (a2 - 1.1111111))+ && validate (abs (a3 - 1.25))+ )++testFuncKurt :: Spec+testFuncKurt = do+ let c = S.fromList+ [21.3 :: Double, 38.4, 12.7, 41.6]+ krt <- runIO $ S.fold (Ring.slidingWindow 4 kurtosis) c+ it ( "kurtosis should be 1.2762447351370185 Actual is " +++ show krt+ )++ (validate $ abs (krt - 1.2762447))++testJackKnife :: (Show a, Eq a, Unbox a) =>+ (Array.Array a -> Stream (SpecM ()) a)+ -> [a]+ -> [a]+ -> Spec+testJackKnife f ls expRes = do+ let arr = Array.fromList ls+ res <- Stream.fold Fold.toList $ f arr+ it ("testJackKnife result should be ="+ ++ show expRes+ ++ " Actual is = " ++show res+ )+ (res == expRes)++testFuncHistogram :: Spec+testFuncHistogram = do+ let strm = S.fromList [1..15]+ res <- runIO $+ S.fold (histogram (binOffsetSize (0::Int) (3::Int))) strm+ let expected = Map.fromList+ [ (0::Int, 2::Int)+ , (1, 3)+ , (2, 3)+ , (3, 3)+ , (4, 3)+ , (5, 1)+ ]++ it ("Map should be = "+ ++ show expected+ ++ " Actual is = "+ ++ show res) (expected == res)++testFuncbinFromSizeN :: Int -> Int -> Int -> Int -> HistBin Int -> SpecWith (Arg Bool)+testFuncbinFromSizeN low binSize nbins x exp0 = do+ let res = binFromSizeN low binSize nbins x+ it ("Bin should be = "+ ++ show exp0+ ++ " Actual is = "+ ++ show res) (res == exp0)++testFuncbinFromToN :: Int -> Int -> Int -> Int -> HistBin Int -> SpecWith ()+testFuncbinFromToN low high n x exp0 = do+ let res = binFromToN low high n x+ it ("Bin should be = "+ ++ show exp0+ ++ " Actual is = "+ ++ show res) (res == exp0)++testFrequency :: [Int] -> Map.Map Int Int -> Spec+testFrequency inputList result = do+ freq <- S.fold frequency' $ S.fromList inputList+ it ("Frequency " ++ show freq) $ liftEq2 (==) (==) freq result++testMode :: [Int] -> Maybe (Int, Int) -> Spec+testMode inputList res = do+ mode0 <- S.fold mode $ S.fromList inputList+ it ("Mode " ++ show mode0) $ mode0 == res++testFFT :: Property+testFFT = do+ let lengths = [2, 4, 8, 16]+ forAll (elements lengths) $ \list_length ->+ forAll (vectorOf list_length (choose (-50.0 :: Double, 100.0)))+ $ \ls ->+ monadicIO $ do+ let tc = map (\x -> x :+ 0) ls+ let vr = V.toList (STAT.fft (V.fromList tc)+ :: V.Vector STAT.CD)+ marr <- MA.fromList tc+ fft marr+ res <- MA.toList marr+ assert (vr == res)++sampleList :: [Double]+sampleList = [1.0, 2.0, 3.0, 4.0, 5.0]++testResample :: [Double] -> Spec+testResample sample = do+ let sampleArr = Array.fromList sample+ sampleSet = Set.fromList sample+ resampleList <- runIO $ S.fold Fold.toList $ S.unfold resample sampleArr+ let resampleSet = Set.fromList resampleList+ sub = Set.isSubsetOf resampleSet sampleSet+ -- XXX We should not use dynamic output in test description+ it ("resample " ++ show resampleList)+ (Prelude.length resampleList == Array.length sampleArr && sub)++testFoldResamples :: Int -> [Double] -> Spec+testFoldResamples n sample = do+ let arr = Array.fromList sample+ a <- runIO $ S.fold Fold.toList $ foldResamples n arr Fold.mean+ -- XXX We should not use dynamic output in test description+ it ("foldResamples " ++ show a) (Prelude.length a == n)++main :: IO ()+main = hspec $ do+ describe "Numerical stability while streaming" $ do+ let numElem = 80000+ winSize = 800+ testCaseChunk = [9007199254740992, 1, 1.0 :: Double,+ 9007199254740992, 1, 1, 1, 9007199254740992]+ testCase = take numElem $ cycle testCaseChunk+ deviationLimit = 1+ testFunc f = do+ let c = S.fromList testCase+ a <- runIO $ S.fold (Ring.slidingWindow winSize f) c+ b <- runIO $ S.fold f $ S.drop (numElem - winSize)+ $ fmap (, Nothing) c+ let c1 = a - b+ it ("should not deviate more than " ++ show deviationLimit)+ $ c1 >= -1 * deviationLimit && c1 <= deviationLimit++ describe "Sum" $ testFunc sum+ describe "mean" $ testFunc mean+ describe "welfordMean" $ testFunc welfordMean++ describe "Correctness" $ do+ let winSize = 3+ testCase1 = [31, 41, 59, 26, 53, 58, 97] :: [Double]+ testCase2 = replicate 5 1.0 ++ [7.0]++ testFunc tc f sI sW = do+ let c = S.fromList tc+ a <- runIO $ S.fold Fold.toList $ S.postscan f $ fmap (, Nothing) c+ b <- runIO $ S.fold Fold.toList $ S.postscan+ (Ring.slidingWindow winSize f) c+ it "Infinite" $ a == sI+ it ("Finite " ++ show winSize) $ b == sW++ -- Resampling+ describe "JackKnife Mean" $+ testJackKnife jackKnifeMean jackKnifeInput jackMeanRes+ describe "JackKnife Variance" $ do+ testJackKnife jackKnifeVariance jackKnifeInput jackVarianceRes+ describe "JackKnife StdDev" $+ testJackKnife jackKnifeStdDev jackKnifeInput jackStdDevRes++ describe "resample" $ do+ testResample sampleList+ describe "foldResamples 4" $ do+ testFoldResamples 4 sampleList+ describe "foldResamples 6" $ do+ testFoldResamples 6 sampleList++ -- Spread/Mean+ describe "MD" $ testFuncMD md+ describe "Kurt" testFuncKurt+ prop "fft" testFFT+ describe "minimum" $ do+ let scanInf = [31, 31, 31, 26, 26, 26, 26] :: [Double]+ scanWin = [31, 31, 31, 26, 26, 26, 53] :: [Double]+ testFunc testCase1 minimum scanInf scanWin+ describe "maximum" $ do+ let scanInf = [31, 41, 59, 59, 59, 59, 97] :: [Double]+ scanWin = [31, 41, 59, 59, 59, 58, 97] :: [Double]+ testFunc testCase1 maximum scanInf scanWin+ describe "range" $ do+ let scanInf = [0, 10, 28, 33, 33, 33, 71] :: [Double]+ scanWin = [0, 10, 28, 33, 33, 32, 44] :: [Double]+ testFunc testCase1 range scanInf scanWin+ describe "sum" $ do+ let scanInf = [1, 2, 3, 4, 5, 12] :: [Double]+ scanWin = [1, 2, 3, 3, 3, 9] :: [Double]+ testFunc testCase2 sum scanInf scanWin+ describe "mean" $ do+ let scanInf = [1, 1, 1, 1, 1, 2] :: [Double]+ scanWin = [1, 1, 1, 1, 1, 3] :: [Double]+ testFunc testCase2 mean scanInf scanWin+ describe "welfordMean" $ do+ let scanInf = [1, 1, 1, 1, 1, 2] :: [Double]+ scanWin = [1, 1, 1, 1, 1, 3] :: [Double]+ testFunc testCase2 welfordMean scanInf scanWin++ -- Probability Distribution+ describe "frequency"+ $ testFrequency+ [1::Int, 1, 2, 3, 3, 3]+ (Map.fromList [(1, 2), (2, 1), (3, 3)])+ describe "Mode" $ testMode [1::Int, 1, 2, 3, 3, 3] (Just (3, 3))+ describe "Mode Empty " $ testMode ([]::[Int]) Nothing++ describe "histogram" testFuncHistogram+ describe "binFromSizeN AboveRange" $+ testFuncbinFromSizeN (0::Int) 2 10 55 AboveRange+ describe "binFromSizeN BelowRange" $+ testFuncbinFromSizeN (0::Int) 2 10 (-1) BelowRange+ describe "binFromSizeN InRange" $+ testFuncbinFromSizeN (0::Int) 2 10 19 (InRange 9)+ describe "binFromSizeN AboveRange" $+ testFuncbinFromSizeN (0::Int) 2 10 20 AboveRange+ describe "binFromToN AboveRange" $+ testFuncbinFromToN (0::Int) 49 10 55 AboveRange+ describe "binFromToN BelowRange" $+ testFuncbinFromToN (0::Int) 49 10 (-1) BelowRange+ describe "binFromToN InRange" $+ testFuncbinFromToN (0::Int) 49 10 19 (InRange 3)+ describe "binFromToN AboveRange" $+ testFuncbinFromToN (0::Int) 50 10 20 (InRange 4)+ prop "variance" testVariance+ prop "stdDev" testStdDev