diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,5 @@ +# Changelog + +## 0.1.0 (Apr 2023) + +* Initial version diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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 diff --git a/NOTICE b/NOTICE new file mode 100644 --- /dev/null +++ b/NOTICE @@ -0,0 +1,5 @@ +streamly-statistics +Copyright 2021 Composewell Technologies + +This product includes software developed at +Composewell Technologies (http://www.composewell.com). diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/benchmark/Main.hs b/benchmark/Main.hs new file mode 100644 --- /dev/null +++ b/benchmark/Main.hs @@ -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 + ] + ] diff --git a/src/Streamly/Statistics.hs b/src/Streamly/Statistics.hs new file mode 100644 --- /dev/null +++ b/src/Streamly/Statistics.hs @@ -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 diff --git a/streamly-statistics.cabal b/streamly-statistics.cabal new file mode 100644 --- /dev/null +++ b/streamly-statistics.cabal @@ -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 diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -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