-- |
-- Module      : Streamly.Statistics.Scanl
-- Copyright   : (c) 2024 Composewell Technologies
-- License     : Apache-2.0
-- Maintainer  : streamly@composewell.com
-- Stability   : experimental
-- Portability : GHC
--
-- See "Streamly.Statistics" for general information. This module provides
-- scans instead of folds.

{-# LANGUAGE ScopedTypeVariables #-}
module Streamly.Statistics.Scanl
    (
    -- * Incremental Scans
    -- | Scans of type @Scanl m (a, Maybe a) b@ are incremental sliding window
    -- scans. 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@.
    --
    -- Also see "Streamly.Data.Scanl" for some basic window scans.

    -- * Summary Statistics
    -- | See https://en.wikipedia.org/wiki/Summary_statistics .

    -- ** Location
    -- | See https://en.wikipedia.org/wiki/Location_parameter .
    --
    -- See https://en.wikipedia.org/wiki/Central_tendency .
      incrMinimum
    , incrMaximum
    , incrRawMoment
    , incrRawMomentFrac

    -- Pythagorean means (https://en.wikipedia.org/wiki/Pythagorean_means)
    , incrWelfordMean
    , incrGeometricMean
    , incrHarmonicMean

    , incrQuadraticMean

    -- Generalized mean
    , incrPowerMean
    , incrPowerMeanFrac

    -- ** Weighted Means
    -- | Exponential Smoothing.
    , ewma
    , ewmaRampUpSmoothing
    , incrEwma

    -- ** 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 .
    , incrRange
    , incrMd
    , incrVariance
    , incrStdDev

    -- ** 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 .
    , incrSkewness
    , incrKurtosis

    -- XXX Move to Statistics.Sample or Statistics.Estimation module?
    -- ** Estimation
    , incrSampleVariance
    , incrSampleStdDev
    , incrStdErrMean

    -- ** Probability Distribution
    , incrFrequency
    )
where

import Control.Exception (assert)
import Control.Monad (when)
import Control.Monad.IO.Class (MonadIO(..))
import Data.Function ((&))
import Data.Map.Strict (Map)
import Data.Maybe (fromMaybe)
import Streamly.Internal.Data.Fold (Step(..))
import Streamly.Internal.Data.Scanl (Scanl(..), Incr(..))
import Streamly.Internal.Data.Tuple.Strict (Tuple'(..), Tuple3'(..))

import qualified Data.Map.Strict as Map
import qualified Deque.Strict as Deque
import qualified Streamly.Data.Fold as Fold
import qualified Streamly.Internal.Data.RingArray as Ring
import qualified Streamly.Internal.Data.Scanl as Scanl
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.

-------------------------------------------------------------------------------
-- 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.Internal.Data.Fold.windowMinimum 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 incrMinimum #-}
incrMinimum :: (Monad m, Ord a) => Scanl m (Incr a) a
incrMinimum :: forall (m :: * -> *) a. (Monad m, Ord a) => Scanl m (Incr a) a
incrMinimum = (Tuple3' Int Int (Deque (Int, a))
 -> Incr a -> m (Step (Tuple3' Int Int (Deque (Int, a))) a))
-> m (Step (Tuple3' Int Int (Deque (Int, a))) a)
-> (Tuple3' Int Int (Deque (Int, a)) -> m a)
-> (Tuple3' Int Int (Deque (Int, a)) -> m a)
-> Scanl m (Incr a) a
forall (m :: * -> *) a b s.
(s -> a -> m (Step s b))
-> m (Step s b) -> (s -> m b) -> (s -> m b) -> Scanl m a b
Scanl Tuple3' Int Int (Deque (Int, a))
-> Incr a -> m (Step (Tuple3' Int Int (Deque (Int, a))) a)
forall {m :: * -> *} {a} {b} {b}.
(Monad m, Num a, Ord a, Ord b) =>
Tuple3' a a (Deque (a, b))
-> Incr b -> m (Step (Tuple3' a a (Deque (a, b))) b)
step m (Step (Tuple3' Int Int (Deque (Int, a))) a)
forall {a} {b}. m (Step (Tuple3' Int Int (Deque (Int, a))) b)
initial Tuple3' Int Int (Deque (Int, a)) -> m a
forall {m :: * -> *} {a} {a} {b} {a}.
(Monad m, Num a) =>
Tuple3' a b (Deque (a, a)) -> m a
extract Tuple3' Int Int (Deque (Int, a)) -> m a
forall {m :: * -> *} {a} {a} {b} {a}.
(Monad m, Num a) =>
Tuple3' a b (Deque (a, a)) -> m a
extract

    where

    initial :: m (Step (Tuple3' Int Int (Deque (Int, a))) b)
initial =
        Step (Tuple3' Int Int (Deque (Int, a))) b
-> m (Step (Tuple3' Int Int (Deque (Int, a))) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (Step (Tuple3' Int Int (Deque (Int, a))) b
 -> m (Step (Tuple3' Int Int (Deque (Int, a))) b))
-> Step (Tuple3' Int Int (Deque (Int, a))) b
-> m (Step (Tuple3' Int Int (Deque (Int, a))) b)
forall a b. (a -> b) -> a -> b
$ Tuple3' Int Int (Deque (Int, a))
-> Step (Tuple3' Int Int (Deque (Int, a))) b
forall s b. s -> Step s b
Partial
            (Tuple3' Int Int (Deque (Int, a))
 -> Step (Tuple3' Int Int (Deque (Int, a))) b)
-> Tuple3' Int Int (Deque (Int, a))
-> Step (Tuple3' Int Int (Deque (Int, a))) b
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Deque (Int, a) -> Tuple3' Int Int (Deque (Int, a))
forall a b c. a -> b -> c -> Tuple3' a b c
Tuple3' (Int
0 :: Int) (Int
0 :: Int) (Deque (Int, a)
forall {a}. Deque (Int, a)
forall a. Monoid a => a
mempty :: Deque.Deque (Int, a))

    step :: Tuple3' a a (Deque (a, b))
-> Incr b -> m (Step (Tuple3' a a (Deque (a, b))) b)
step (Tuple3' a
i a
w Deque (a, b)
q) (Insert b
a) =
                Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
                    (Step (Tuple3' a a (Deque (a, b))) b
 -> m (Step (Tuple3' a a (Deque (a, b))) b))
-> Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a b. (a -> b) -> a -> b
$ Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b
forall s b. s -> Step s b
Partial
                    (Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b)
-> Tuple3' a a (Deque (a, b))
-> Step (Tuple3' a a (Deque (a, b))) b
forall a b. (a -> b) -> a -> b
$ a -> a -> Deque (a, b) -> Tuple3' a a (Deque (a, b))
forall a b c. a -> b -> c -> Tuple3' a b c
Tuple3'
                        (a
i a -> a -> a
forall a. Num a => a -> a -> a
+ a
1)
                        (a
w a -> a -> a
forall a. Num a => a -> a -> a
+ a
1)
                        (a -> Deque (a, b) -> a -> Deque (a, b)
forall {a} {b}.
(Ord a, Num a) =>
a -> Deque (a, b) -> a -> Deque (a, b)
headCheck a
i Deque (a, b)
q (a
w a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) Deque (a, b) -> (Deque (a, b) -> Deque (a, b)) -> Deque (a, b)
forall a b. a -> (a -> b) -> b
& (a, b) -> Deque (a, b) -> Deque (a, b)
forall {a} {a}. Ord a => (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a
i, b
a))

    step (Tuple3' a
i a
w Deque (a, b)
q) (Replace b
_ b
new) =
        Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (Step (Tuple3' a a (Deque (a, b))) b
 -> m (Step (Tuple3' a a (Deque (a, b))) b))
-> Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a b. (a -> b) -> a -> b
$ Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b
forall s b. s -> Step s b
Partial
            (Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b)
-> Tuple3' a a (Deque (a, b))
-> Step (Tuple3' a a (Deque (a, b))) b
forall a b. (a -> b) -> a -> b
$ a -> a -> Deque (a, b) -> Tuple3' a a (Deque (a, b))
forall a b c. a -> b -> c -> Tuple3' a b c
Tuple3' (a
i a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) a
w (a -> Deque (a, b) -> a -> Deque (a, b)
forall {a} {b}.
(Ord a, Num a) =>
a -> Deque (a, b) -> a -> Deque (a, b)
headCheck a
i Deque (a, b)
q a
w Deque (a, b) -> (Deque (a, b) -> Deque (a, b)) -> Deque (a, b)
forall a b. a -> (a -> b) -> b
& (a, b) -> Deque (a, b) -> Deque (a, b)
forall {a} {a}. Ord a => (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a
i, b
new))

    {-# INLINE headCheck #-}
    headCheck :: a -> Deque (a, b) -> a -> Deque (a, b)
headCheck a
i Deque (a, b)
q a
w =
        case Deque (a, b) -> Maybe ((a, b), Deque (a, b))
forall a. Deque a -> Maybe (a, Deque a)
Deque.uncons Deque (a, b)
q of
            Maybe ((a, b), Deque (a, b))
Nothing -> Deque (a, b)
q
            Just ((a, b)
ia', Deque (a, b)
q') ->
                if (a, b) -> a
forall a b. (a, b) -> a
fst (a, b)
ia' a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
i a -> a -> a
forall a. Num a => a -> a -> a
- a
w
                then Deque (a, b)
q'
                else Deque (a, b)
q

    dqloop :: (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a, a)
ia Deque (a, a)
q =
        case Deque (a, a) -> Maybe ((a, a), Deque (a, a))
forall a. Deque a -> Maybe (a, Deque a)
Deque.unsnoc Deque (a, a)
q of
            Maybe ((a, a), Deque (a, a))
Nothing -> (a, a) -> Deque (a, a) -> Deque (a, a)
forall a. a -> Deque a -> Deque a
Deque.snoc (a, a)
ia Deque (a, a)
q
            -- XXX This can be improved for the case of `=`
            Just ((a, a)
ia', Deque (a, a)
q') ->
                if (a, a) -> a
forall a b. (a, b) -> b
snd (a, a)
ia a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= (a, a) -> a
forall a b. (a, b) -> b
snd (a, a)
ia'
                then (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a, a)
ia Deque (a, a)
q'
                else (a, a) -> Deque (a, a) -> Deque (a, a)
forall a. a -> Deque a -> Deque a
Deque.snoc (a, a)
ia Deque (a, a)
q

    extract :: Tuple3' a b (Deque (a, a)) -> m a
extract (Tuple3' a
_ b
_ Deque (a, a)
q) =
        a -> m a
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (a -> m a) -> a -> m a
forall a b. (a -> b) -> a -> b
$ (a, a) -> a
forall a b. (a, b) -> b
snd
            ((a, a) -> a) -> (a, a) -> a
forall a b. (a -> b) -> a -> b
$ (a, a) -> Maybe (a, a) -> (a, a)
forall a. a -> Maybe a -> a
fromMaybe (a
0, [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"minimum: Empty stream")
            (Maybe (a, a) -> (a, a)) -> Maybe (a, a) -> (a, a)
forall a b. (a -> b) -> a -> b
$ Deque (a, a) -> Maybe (a, a)
forall a. Deque a -> Maybe a
Deque.head Deque (a, a)
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.Internal.Data.Fold.windowMaximum
-- 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 incrMaximum #-}
incrMaximum :: (Monad m, Ord a) => Scanl m (Incr a) a
incrMaximum :: forall (m :: * -> *) a. (Monad m, Ord a) => Scanl m (Incr a) a
incrMaximum = (Tuple3' Int Int (Deque (Int, a))
 -> Incr a -> m (Step (Tuple3' Int Int (Deque (Int, a))) a))
-> m (Step (Tuple3' Int Int (Deque (Int, a))) a)
-> (Tuple3' Int Int (Deque (Int, a)) -> m a)
-> (Tuple3' Int Int (Deque (Int, a)) -> m a)
-> Scanl m (Incr a) a
forall (m :: * -> *) a b s.
(s -> a -> m (Step s b))
-> m (Step s b) -> (s -> m b) -> (s -> m b) -> Scanl m a b
Scanl Tuple3' Int Int (Deque (Int, a))
-> Incr a -> m (Step (Tuple3' Int Int (Deque (Int, a))) a)
forall {m :: * -> *} {a} {b} {b}.
(Monad m, Num a, Ord a, Ord b) =>
Tuple3' a a (Deque (a, b))
-> Incr b -> m (Step (Tuple3' a a (Deque (a, b))) b)
step m (Step (Tuple3' Int Int (Deque (Int, a))) a)
forall {a} {b}. m (Step (Tuple3' Int Int (Deque (Int, a))) b)
initial Tuple3' Int Int (Deque (Int, a)) -> m a
forall {m :: * -> *} {a} {a} {b} {a}.
(Monad m, Num a) =>
Tuple3' a b (Deque (a, a)) -> m a
extract Tuple3' Int Int (Deque (Int, a)) -> m a
forall {m :: * -> *} {a} {a} {b} {a}.
(Monad m, Num a) =>
Tuple3' a b (Deque (a, a)) -> m a
extract

    where

    initial :: m (Step (Tuple3' Int Int (Deque (Int, a))) b)
initial =
        Step (Tuple3' Int Int (Deque (Int, a))) b
-> m (Step (Tuple3' Int Int (Deque (Int, a))) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (Step (Tuple3' Int Int (Deque (Int, a))) b
 -> m (Step (Tuple3' Int Int (Deque (Int, a))) b))
-> Step (Tuple3' Int Int (Deque (Int, a))) b
-> m (Step (Tuple3' Int Int (Deque (Int, a))) b)
forall a b. (a -> b) -> a -> b
$ Tuple3' Int Int (Deque (Int, a))
-> Step (Tuple3' Int Int (Deque (Int, a))) b
forall s b. s -> Step s b
Partial
            (Tuple3' Int Int (Deque (Int, a))
 -> Step (Tuple3' Int Int (Deque (Int, a))) b)
-> Tuple3' Int Int (Deque (Int, a))
-> Step (Tuple3' Int Int (Deque (Int, a))) b
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Deque (Int, a) -> Tuple3' Int Int (Deque (Int, a))
forall a b c. a -> b -> c -> Tuple3' a b c
Tuple3' (Int
0 :: Int) (Int
0 :: Int) (Deque (Int, a)
forall {a}. Deque (Int, a)
forall a. Monoid a => a
mempty :: Deque.Deque (Int, a))

    step :: Tuple3' a a (Deque (a, b))
-> Incr b -> m (Step (Tuple3' a a (Deque (a, b))) b)
step (Tuple3' a
i a
w Deque (a, b)
q) (Insert b
a) =
        Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (Step (Tuple3' a a (Deque (a, b))) b
 -> m (Step (Tuple3' a a (Deque (a, b))) b))
-> Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a b. (a -> b) -> a -> b
$ Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b
forall s b. s -> Step s b
Partial
            (Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b)
-> Tuple3' a a (Deque (a, b))
-> Step (Tuple3' a a (Deque (a, b))) b
forall a b. (a -> b) -> a -> b
$ a -> a -> Deque (a, b) -> Tuple3' a a (Deque (a, b))
forall a b c. a -> b -> c -> Tuple3' a b c
Tuple3'
                (a
i a -> a -> a
forall a. Num a => a -> a -> a
+ a
1)
                (a
w a -> a -> a
forall a. Num a => a -> a -> a
+ a
1)
                (a -> Deque (a, b) -> a -> Deque (a, b)
forall {a} {b}.
(Ord a, Num a) =>
a -> Deque (a, b) -> a -> Deque (a, b)
headCheck a
i Deque (a, b)
q (a
w a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) Deque (a, b) -> (Deque (a, b) -> Deque (a, b)) -> Deque (a, b)
forall a b. a -> (a -> b) -> b
& (a, b) -> Deque (a, b) -> Deque (a, b)
forall {a} {a}. Ord a => (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a
i, b
a))

    step (Tuple3' a
i a
w Deque (a, b)
q) (Replace b
_ b
new) =
        Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (Step (Tuple3' a a (Deque (a, b))) b
 -> m (Step (Tuple3' a a (Deque (a, b))) b))
-> Step (Tuple3' a a (Deque (a, b))) b
-> m (Step (Tuple3' a a (Deque (a, b))) b)
forall a b. (a -> b) -> a -> b
$ Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b
forall s b. s -> Step s b
Partial
            (Tuple3' a a (Deque (a, b)) -> Step (Tuple3' a a (Deque (a, b))) b)
-> Tuple3' a a (Deque (a, b))
-> Step (Tuple3' a a (Deque (a, b))) b
forall a b. (a -> b) -> a -> b
$ a -> a -> Deque (a, b) -> Tuple3' a a (Deque (a, b))
forall a b c. a -> b -> c -> Tuple3' a b c
Tuple3' (a
i a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) a
w (a -> Deque (a, b) -> a -> Deque (a, b)
forall {a} {b}.
(Ord a, Num a) =>
a -> Deque (a, b) -> a -> Deque (a, b)
headCheck a
i Deque (a, b)
q a
w Deque (a, b) -> (Deque (a, b) -> Deque (a, b)) -> Deque (a, b)
forall a b. a -> (a -> b) -> b
& (a, b) -> Deque (a, b) -> Deque (a, b)
forall {a} {a}. Ord a => (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a
i, b
new))

    {-# INLINE headCheck #-}
    headCheck :: a -> Deque (a, b) -> a -> Deque (a, b)
headCheck a
i Deque (a, b)
q a
w =
        case Deque (a, b) -> Maybe ((a, b), Deque (a, b))
forall a. Deque a -> Maybe (a, Deque a)
Deque.uncons Deque (a, b)
q of
            Maybe ((a, b), Deque (a, b))
Nothing -> Deque (a, b)
q
            Just ((a, b)
ia', Deque (a, b)
q') ->
                if (a, b) -> a
forall a b. (a, b) -> a
fst (a, b)
ia' a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
i a -> a -> a
forall a. Num a => a -> a -> a
- a
w
                then Deque (a, b)
q'
                else Deque (a, b)
q

    dqloop :: (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a, a)
ia Deque (a, a)
q =
        case Deque (a, a) -> Maybe ((a, a), Deque (a, a))
forall a. Deque a -> Maybe (a, Deque a)
Deque.unsnoc Deque (a, a)
q of
            Maybe ((a, a), Deque (a, a))
Nothing -> (a, a) -> Deque (a, a) -> Deque (a, a)
forall a. a -> Deque a -> Deque a
Deque.snoc (a, a)
ia Deque (a, a)
q
            -- XXX This can be improved for the case of `=`
            Just ((a, a)
ia', Deque (a, a)
q') ->
                if (a, a) -> a
forall a b. (a, b) -> b
snd (a, a)
ia a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= (a, a) -> a
forall a b. (a, b) -> b
snd (a, a)
ia'
                then (a, a) -> Deque (a, a) -> Deque (a, a)
dqloop (a, a)
ia Deque (a, a)
q'
                else (a, a) -> Deque (a, a) -> Deque (a, a)
forall a. a -> Deque a -> Deque a
Deque.snoc (a, a)
ia Deque (a, a)
q

    extract :: Tuple3' a b (Deque (a, a)) -> m a
extract (Tuple3' a
_ b
_ Deque (a, a)
q) =
        a -> m a
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (a -> m a) -> a -> m a
forall a b. (a -> b) -> a -> b
$ (a, a) -> a
forall a b. (a, b) -> b
snd
            ((a, a) -> a) -> (a, a) -> a
forall a b. (a -> b) -> a -> b
$ (a, a) -> Maybe (a, a) -> (a, a)
forall a. a -> Maybe a -> a
fromMaybe (a
0, [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"maximum: Empty stream")
            (Maybe (a, a) -> (a, a)) -> Maybe (a, a) -> (a, a)
forall a b. (a -> b) -> a -> b
$ Deque (a, a) -> Maybe (a, a)
forall a. Deque a -> Maybe a
Deque.head Deque (a, a)
q

-------------------------------------------------------------------------------
-- Mean
-------------------------------------------------------------------------------

-- | Recompute mean from old mean when an item is added to the sample.
{-# INLINE meanAdd #-}
meanAdd :: Fractional a => Int -> a -> a -> a
meanAdd :: forall a. Fractional a => Int -> a -> a -> a
meanAdd Int
n a
oldMean a
newItem =
    let delta :: a
delta = (a
newItem a -> a -> a
forall a. Num a => a -> a -> a
- a
oldMean) a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
     in a
oldMean a -> a -> a
forall a. Num a => a -> a -> a
+ a
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 :: forall a. Fractional a => Int -> a -> a -> a -> a
meanReplace Int
n a
oldMean a
oldItem a
newItem =
    let n1 :: a
n1 = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
        -- Compute two deltas instead of a single (newItem - oldItem) because
        -- the latter would be too small causing rounding errors.
        delta1 :: a
delta1 = (a
newItem a -> a -> a
forall a. Num a => a -> a -> a
- a
oldMean) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
n1
        delta2 :: a
delta2 = (a
oldItem a -> a -> a
forall a. Num a => a -> a -> a
- a
oldMean) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
n1
     in (a
oldMean a -> a -> a
forall a. Num a => a -> a -> a
+ a
delta1) a -> a -> a
forall a. Num a => a -> a -> a
- a
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 incrWelfordMean #-}
incrWelfordMean :: forall m a. (Monad m, Fractional a) => Scanl m (Incr a) a
incrWelfordMean :: forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
incrWelfordMean = (Tuple' a Int -> Incr a -> m (Step (Tuple' a Int) a))
-> m (Step (Tuple' a Int) a)
-> (Tuple' a Int -> m a)
-> (Tuple' a Int -> m a)
-> Scanl m (Incr a) a
forall (m :: * -> *) a b s.
(s -> a -> m (Step s b))
-> m (Step s b) -> (s -> m b) -> (s -> m b) -> Scanl m a b
Scanl Tuple' a Int -> Incr a -> m (Step (Tuple' a Int) a)
forall {m :: * -> *} {a} {b}.
(Monad m, Fractional a) =>
Tuple' a Int -> Incr a -> m (Step (Tuple' a Int) b)
step m (Step (Tuple' a Int) a)
forall {b}. m (Step (Tuple' a Int) b)
initial Tuple' a Int -> m a
forall {m :: * -> *} {a} {b}. Monad m => Tuple' a b -> m a
extract Tuple' a Int -> m a
forall {m :: * -> *} {a} {b}. Monad m => Tuple' a b -> m a
extract

    where

    initial :: m (Step (Tuple' a Int) b)
initial =
        Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return
            (Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b))
-> Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b)
forall a b. (a -> b) -> a -> b
$ Tuple' a Int -> Step (Tuple' a Int) b
forall s b. s -> Step s b
Partial
            (Tuple' a Int -> Step (Tuple' a Int) b)
-> Tuple' a Int -> Step (Tuple' a Int) b
forall a b. (a -> b) -> a -> b
$ a -> Int -> Tuple' a Int
forall a b. a -> b -> Tuple' a b
Tuple'
                (a
0 :: a)   -- running mean
                (Int
0 :: Int) -- count of items in the window

    step :: Tuple' a Int -> Incr a -> m (Step (Tuple' a Int) b)
step (Tuple' a
oldMean Int
w) (Insert a
new) =
        Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return (Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b))
-> Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b)
forall a b. (a -> b) -> a -> b
$ Tuple' a Int -> Step (Tuple' a Int) b
forall s b. s -> Step s b
Partial (Tuple' a Int -> Step (Tuple' a Int) b)
-> Tuple' a Int -> Step (Tuple' a Int) b
forall a b. (a -> b) -> a -> b
$ a -> Int -> Tuple' a Int
forall a b. a -> b -> Tuple' a b
Tuple' (Int -> a -> a -> a
forall a. Fractional a => Int -> a -> a -> a
meanAdd Int
w a
oldMean a
new) (Int
w Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)

    step (Tuple' a
oldMean Int
w) (Replace a
old a
new) =
        Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return (Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b))
-> Step (Tuple' a Int) b -> m (Step (Tuple' a Int) b)
forall a b. (a -> b) -> a -> b
$ Tuple' a Int -> Step (Tuple' a Int) b
forall s b. s -> Step s b
Partial (Tuple' a Int -> Step (Tuple' a Int) b)
-> Tuple' a Int -> Step (Tuple' a Int) b
forall a b. (a -> b) -> a -> b
$ a -> Int -> Tuple' a Int
forall a b. a -> b -> Tuple' a b
Tuple' (Int -> a -> a -> a -> a
forall a. Fractional a => Int -> a -> a -> a -> a
meanReplace Int
w a
oldMean a
old a
new) Int
w

    extract :: Tuple' a b -> m a
extract (Tuple' a
x b
_) = a -> m a
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return a
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 (/) (Fold.windowPowerSum p) Fold.windowLength
--
-- See https://en.wikipedia.org/wiki/Moment_(mathematics) .
--
-- /Space/: \(\mathcal{O}(1)\)
--
-- /Time/: \(\mathcal{O}(n)\)
{-# INLINE incrRawMoment #-}
incrRawMoment :: (Monad m, Fractional a) => Int -> Scanl m (Incr a) a
incrRawMoment :: forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Int -> Scanl m (Incr a) a
incrRawMoment Int
k =
    (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (m :: * -> *) b c d a.
Monad m =>
(b -> c -> d) -> Scanl m a b -> Scanl m a c -> Scanl m a d
Scanl.teeWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) (Int -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Num a) =>
Int -> Scanl m (Incr a) a
Scanl.incrPowerSum Int
k) Scanl m (Incr a) a
forall (m :: * -> *) b a. (Monad m, Num b) => Scanl m (Incr a) b
Scanl.incrCount

-- | Like 'rawMoment' but powers can be negative or fractional. This is
-- slower than 'rawMoment' for positive intergal powers.
--
-- >>> rawMomentFrac p = Fold.teeWith (/) (Fold.windowPowerSumFrac p) Fold.windowLength
--
{-# INLINE incrRawMomentFrac #-}
incrRawMomentFrac :: (Monad m, Floating a) => a -> Scanl m (Incr a) a
incrRawMomentFrac :: forall (m :: * -> *) a.
(Monad m, Floating a) =>
a -> Scanl m (Incr a) a
incrRawMomentFrac a
k =
    (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (m :: * -> *) b c d a.
Monad m =>
(b -> c -> d) -> Scanl m a b -> Scanl m a c -> Scanl m a d
Scanl.teeWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) (a -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Floating a) =>
a -> Scanl m (Incr a) a
Scanl.incrPowerSumFrac a
k) Scanl m (Incr a) a
forall (m :: * -> *) b a. (Monad m, Num b) => Scanl m (Incr a) b
Scanl.incrCount

-- 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 incrPowerMean #-}
incrPowerMean :: (Monad m, Floating a) => Int -> Scanl m (Incr a) a
incrPowerMean :: forall (m :: * -> *) a.
(Monad m, Floating a) =>
Int -> Scanl m (Incr a) a
incrPowerMean Int
k = (a -> a -> a
forall a. Floating a => a -> a -> a
** (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k)) (a -> a) -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Int -> Scanl m (Incr a) a
incrRawMoment Int
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 incrPowerMeanFrac #-}
incrPowerMeanFrac :: (Monad m, Floating a) => a -> Scanl m (Incr a) a
incrPowerMeanFrac :: forall (m :: * -> *) a.
(Monad m, Floating a) =>
a -> Scanl m (Incr a) a
incrPowerMeanFrac a
k = (a -> a -> a
forall a. Floating a => a -> a -> a
** (a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
k)) (a -> a) -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> a -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Floating a) =>
a -> Scanl m (Incr a) a
incrRawMomentFrac a
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 incrHarmonicMean #-}
incrHarmonicMean :: (Monad m, Fractional a) => Scanl m (Incr a) a
incrHarmonicMean :: forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
incrHarmonicMean =
    (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (m :: * -> *) b c d a.
Monad m =>
(b -> c -> d) -> Scanl m a b -> Scanl m a c -> Scanl m a d
Scanl.teeWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/)
        Scanl m (Incr a) a
forall (m :: * -> *) b a. (Monad m, Num b) => Scanl m (Incr a) b
Scanl.incrCount ((Incr a -> Incr a) -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall a b (m :: * -> *) r. (a -> b) -> Scanl m b r -> Scanl m a r
Scanl.lmap ((a -> a) -> Incr a -> Incr a
forall a b. (a -> b) -> Incr a -> Incr b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
forall a. Fractional a => a -> a
recip) Scanl m (Incr a) a
forall (m :: * -> *) a. (Monad m, Num a) => Scanl m (Incr a) a
Scanl.incrSum)

-- | 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 incrGeometricMean #-}
incrGeometricMean :: (Monad m, Floating a) => Scanl m (Incr a) a
incrGeometricMean :: forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrGeometricMean = a -> a
forall a. Floating a => a -> a
exp (a -> a) -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Incr a -> Incr a) -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall a b (m :: * -> *) r. (a -> b) -> Scanl m b r -> Scanl m a r
Scanl.lmap ((a -> a) -> Incr a -> Incr a
forall a b. (a -> b) -> Incr a -> Incr b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap a -> a
forall a. Floating a => a -> a
log) Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
Scanl.incrMean

-- | 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 incrQuadraticMean #-}
incrQuadraticMean :: (Monad m, Floating a) => Scanl m (Incr a) a
incrQuadraticMean :: forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrQuadraticMean = Int -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Floating a) =>
Int -> Scanl m (Incr a) a
incrPowerMean Int
2

-------------------------------------------------------------------------------
-- Weighted Means
-------------------------------------------------------------------------------

-- XXX Is this numerically stable? We can use the kbn summation here.
-- XXX change the signature to use the newer value first.

-- | ewmaStep smoothing-factor old-value new-value
--
-- >>> ewmaStep k x0 x1 = k * x1 + (1 - k) * x0
--
{-# INLINE ewmaStep #-}
ewmaStep :: Double -> Double -> Double -> Double
ewmaStep :: Double -> Double -> Double -> Double
ewmaStep Double
k Double
x0 Double
x1 = (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x1

-- XXX Compute this in a sliding window?

-- | Exponentially weighted moving average is given by @ewma w@ where @w@ is
-- the smoothing factor varying between 0 and 1. To compute the moving average
-- the newest input is given a weight of @w@ and the running ewma is given a
-- weight of @1 - w@.
--
-- For an empty stream this API returns 0. For a non-empty stream the first
-- value in the stream is used as the initial ewma.
--
-- The higher the smoothing factor the more weight is given to the new value.
-- Consider some interesting special cases, when @w@ is 1 the ewma is always
-- the newest value, when @w@ is 0 then the ewma is always the oldest value. If
-- @w@ is @0.5@ then the new inputs get exponentially weighted by weights
-- @1\/2, 1\/4, 1\/8, ...@, see below for details.
--
-- Mathematically, we define ewma \(s_n\), of \(n\) values, \(x_1,\ldots,x_n\),
-- recursively as follows:
--
-- \(\begin{align} s_0& = x_0, \quad \text{initial value}\\ s_n & = \alpha x_{n} + (1-\alpha)s_{n-1},\quad n>0 \end{align}\)
--
-- If we expand the recursive term it reveals 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 -> Scanl m Double Double
ewma :: forall (m :: * -> *). Monad m => Double -> Scanl m Double Double
ewma Double
k = Tuple' Double Double -> Double
forall {a} {b}. Tuple' a b -> a
extract (Tuple' Double Double -> Double)
-> Scanl m Double (Tuple' Double Double) -> Scanl m Double Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Tuple' Double Double -> Double -> Tuple' Double Double)
-> Tuple' Double Double -> Scanl m Double (Tuple' Double Double)
forall (m :: * -> *) b a.
Monad m =>
(b -> a -> b) -> b -> Scanl m a b
Scanl.mkScanl Tuple' Double Double -> Double -> Tuple' Double Double
step (Double -> Double -> Tuple' Double Double
forall a b. a -> b -> Tuple' a b
Tuple' Double
0 Double
1)

    where

    step :: Tuple' Double Double -> Double -> Tuple' Double Double
step (Tuple' Double
x0 Double
k1) Double
x = Double -> Double -> Tuple' Double Double
forall a b. a -> b -> Tuple' a b
Tuple' (Double -> Double -> Double -> Double
ewmaStep Double
k1 Double
x0 Double
x) Double
k

    extract :: Tuple' a b -> a
extract (Tuple' a
x b
_) = a
x

-- | @ewma n k@ is like 'ewma' but the smoothing factor used is itself
-- exponentially smoothened starting from @1@ to the final value @k@ using @n@
-- as its smoothing factor. In other words, the smoothing factor is derived by
-- smoothing the series @1,k,k,k,...@ using @n@ as the smoothing factor. As
-- time goes on, the smoothing factor gets closer to k.
--
{-# INLINE ewmaRampUpSmoothing #-}
ewmaRampUpSmoothing :: Monad m => Double -> Double -> Scanl m Double Double
ewmaRampUpSmoothing :: forall (m :: * -> *).
Monad m =>
Double -> Double -> Scanl m Double Double
ewmaRampUpSmoothing Double
n Double
k1 = Tuple' Double Double -> Double
forall {a} {b}. Tuple' a b -> a
extract (Tuple' Double Double -> Double)
-> Scanl m Double (Tuple' Double Double) -> Scanl m Double Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Tuple' Double Double -> Double -> Tuple' Double Double)
-> Tuple' Double Double -> Scanl m Double (Tuple' Double Double)
forall (m :: * -> *) b a.
Monad m =>
(b -> a -> b) -> b -> Scanl m a b
Scanl.mkScanl Tuple' Double Double -> Double -> Tuple' Double Double
step Tuple' Double Double
initial

    where

    initial :: Tuple' Double Double
initial = Double -> Double -> Tuple' Double Double
forall a b. a -> b -> Tuple' a b
Tuple' Double
0 Double
1

    step :: Tuple' Double Double -> Double -> Tuple' Double Double
step (Tuple' Double
x0 Double
k0) Double
x1 =
        let x :: Double
x = Double -> Double -> Double -> Double
ewmaStep Double
k0 Double
x0 Double
x1
            k :: Double
k = Double -> Double -> Double -> Double
ewmaStep Double
n Double
k0 Double
k1
        in Double -> Double -> Tuple' Double Double
forall a b. a -> b -> Tuple' a b
Tuple' Double
x Double
k

    extract :: Tuple' a b -> a
extract (Tuple' a
x b
_) = a
x

data Ewma = EwmaInit | EwmaGo !Double !Double

-- | @incrEwma w@ computes the ewma of the elements incrementally in a window
-- using @w@ as the smoothing factor.
--
-- ewma can be calculated incrementally as follows. If \(x_i\) is the value
-- entering the window, n is the size of window, \(x_1\) is the second last value
-- in the old window and \(x_0\) is the value exiting the window:
--
-- \(s_{new} = \alpha x_i + (1-\alpha)s_{old} + (1-\alpha)^{n} (x_1 - x_0\)\)
--
-- /Unimplemented/
--
{-# INLINE incrEwma #-}
incrEwma :: MonadIO m =>
    Double -> Scanl m (Incr Double, Ring.RingArray Double) Double
incrEwma :: forall (m :: * -> *).
MonadIO m =>
Double -> Scanl m (Incr Double, RingArray Double) Double
incrEwma Double
w = (Ewma -> (Incr Double, RingArray Double) -> m (Step Ewma Double))
-> m (Step Ewma Double)
-> (Ewma -> m Double)
-> (Ewma -> m Double)
-> Scanl m (Incr Double, RingArray Double) Double
forall (m :: * -> *) a b s.
(s -> a -> m (Step s b))
-> m (Step s b) -> (s -> m b) -> (s -> m b) -> Scanl m a b
Scanl Ewma -> (Incr Double, RingArray Double) -> m (Step Ewma Double)
forall {m :: * -> *} {b}.
MonadIO m =>
Ewma -> (Incr Double, RingArray Double) -> m (Step Ewma b)
step m (Step Ewma Double)
forall {b}. m (Step Ewma b)
initial Ewma -> m Double
forall {m :: * -> *}. Monad m => Ewma -> m Double
extract Ewma -> m Double
forall {m :: * -> *}. Monad m => Ewma -> m Double
extract

    where

    initial :: m (Step Ewma b)
initial = do
        Bool -> m () -> m ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
w Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1)
            (m () -> m ()) -> m () -> m ()
forall a b. (a -> b) -> a -> b
$ [Char] -> m ()
forall a. HasCallStack => [Char] -> a
error [Char]
"incrEwma: weight must be >= 0 and <= 1"
        Step Ewma b -> m (Step Ewma b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return (Step Ewma b -> m (Step Ewma b)) -> Step Ewma b -> m (Step Ewma b)
forall a b. (a -> b) -> a -> b
$ Ewma -> Step Ewma b
forall s b. s -> Step s b
Partial Ewma
EwmaInit

    step :: Ewma -> (Incr Double, RingArray Double) -> m (Step Ewma b)
step Ewma
EwmaInit (Insert Double
x, RingArray Double
rng) = do
        let len :: Int
len = RingArray Double -> Int
forall a. Unbox a => RingArray a -> Int
Ring.length RingArray Double
rng
        Bool -> m () -> m ()
forall a. HasCallStack => Bool -> a -> a
assert (Int
len Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0) (() -> m ()
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return ())
        Step Ewma b -> m (Step Ewma b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return (Step Ewma b -> m (Step Ewma b)) -> Step Ewma b -> m (Step Ewma b)
forall a b. (a -> b) -> a -> b
$ Ewma -> Step Ewma b
forall s b. s -> Step s b
Partial (Double -> Double -> Ewma
EwmaGo Double
x Double
1)

    step Ewma
EwmaInit (Incr Double, RingArray Double)
_ = [Char] -> m (Step Ewma b)
forall a. HasCallStack => [Char] -> a
error [Char]
"incrEwma: the first window operation must be Insert"

    step (EwmaGo Double
s Double
k) (Insert Double
x, RingArray Double
_) = do
        let s1 :: Double
s1 = Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s
        Step Ewma b -> m (Step Ewma b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return (Step Ewma b -> m (Step Ewma b)) -> Step Ewma b -> m (Step Ewma b)
forall a b. (a -> b) -> a -> b
$ Ewma -> Step Ewma b
forall s b. s -> Step s b
Partial (Double -> Double -> Ewma
EwmaGo Double
s1 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w)))

    step (EwmaGo Double
s Double
k) (Replace Double
x0 Double
x, RingArray Double
rng) = do
        Double
x1 <- Int -> RingArray Double -> m Double
forall (m :: * -> *) a.
(MonadIO m, Unbox a) =>
Int -> RingArray a -> m a
Ring.unsafeGetIndex Int
0 RingArray Double
rng
        let s1 :: Double
s1 = Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
w) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
x1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x0)
        Step Ewma b -> m (Step Ewma b)
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return (Step Ewma b -> m (Step Ewma b)) -> Step Ewma b -> m (Step Ewma b)
forall a b. (a -> b) -> a -> b
$ Ewma -> Step Ewma b
forall s b. s -> Step s b
Partial (Double -> Double -> Ewma
EwmaGo Double
s1 Double
k)

    extract :: Ewma -> m Double
extract Ewma
EwmaInit = Double -> m Double
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
0
    extract (EwmaGo Double
x Double
_) = Double -> m Double
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return Double
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 incrRange #-}
incrRange :: (Monad m, Num a, Ord a) => Scanl m (Incr a) a
incrRange :: forall (m :: * -> *) a.
(Monad m, Num a, Ord a) =>
Scanl m (Incr a) a
incrRange = (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (m :: * -> *) b c d a.
Monad m =>
(b -> c -> d) -> Scanl m a b -> Scanl m a c -> Scanl m a d
Scanl.teeWith (-) Scanl m (Incr a) a
forall (m :: * -> *) a. (Monad m, Ord a) => Scanl m (Incr a) a
incrMaximum Scanl m (Incr a) a
forall (m :: * -> *) a. (Monad m, Ord a) => Scanl m (Incr a) a
incrMinimum

-- | @md n@ computes the mean absolute deviation (or mean deviation) in a
-- sliding window of last @n@ elements in the stream.
--
-- The input of the scan is (incr, ring), where incr is the incremental window
-- operation and ring is the contents of the entire window in a ring array.
--
-- 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 incrMd #-}
incrMd ::  MonadIO m =>
    Scanl m (Incr Double, Ring.RingArray Double) Double
incrMd :: forall (m :: * -> *).
MonadIO m =>
Scanl m (Incr Double, RingArray Double) Double
incrMd =
    ((Double, Maybe (RingArray Double)) -> m Double)
-> Scanl
     m
     (Incr Double, RingArray Double)
     (Double, Maybe (RingArray Double))
-> Scanl m (Incr Double, RingArray Double) Double
forall (m :: * -> *) b c a.
Monad m =>
(b -> m c) -> Scanl m a b -> Scanl m a c
Scanl.rmapM (Double, Maybe (RingArray Double)) -> m Double
forall {m :: * -> *} {b}.
(Fractional b, MonadIO m, Unbox b) =>
(b, Maybe (RingArray b)) -> m b
computeMD
        (Scanl
   m
   (Incr Double, RingArray Double)
   (Double, Maybe (RingArray Double))
 -> Scanl m (Incr Double, RingArray Double) Double)
-> Scanl
     m
     (Incr Double, RingArray Double)
     (Double, Maybe (RingArray Double))
-> Scanl m (Incr Double, RingArray Double) Double
forall a b. (a -> b) -> a -> b
$ Scanl m (Incr Double, RingArray Double) Double
-> Scanl
     m (Incr Double, RingArray Double) (Maybe (RingArray Double))
-> Scanl
     m
     (Incr Double, RingArray Double)
     (Double, Maybe (RingArray Double))
forall (m :: * -> *) a b c.
Monad m =>
Scanl m a b -> Scanl m a c -> Scanl m a (b, c)
Scanl.tee
            (((Incr Double, RingArray Double) -> Incr Double)
-> Scanl m (Incr Double) Double
-> Scanl m (Incr Double, RingArray Double) Double
forall a b (m :: * -> *) r. (a -> b) -> Scanl m b r -> Scanl m a r
Scanl.lmap (Incr Double, RingArray Double) -> Incr Double
forall a b. (a, b) -> a
fst Scanl m (Incr Double) Double
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
Scanl.incrMean) (((Incr Double, RingArray Double) -> RingArray Double)
-> Scanl m (RingArray Double) (Maybe (RingArray Double))
-> Scanl
     m (Incr Double, RingArray Double) (Maybe (RingArray Double))
forall a b (m :: * -> *) r. (a -> b) -> Scanl m b r -> Scanl m a r
Scanl.lmap (Incr Double, RingArray Double) -> RingArray Double
forall a b. (a, b) -> b
snd Scanl m (RingArray Double) (Maybe (RingArray Double))
forall (m :: * -> *) a. Monad m => Scanl m a (Maybe a)
Scanl.latest)

    where

    computeMD :: (b, Maybe (RingArray b)) -> m b
computeMD (b
mn, Maybe (RingArray b)
mRng) =
        case Maybe (RingArray b)
mRng of
            Just RingArray b
rng -> do
                Fold m b b -> Stream m b -> m b
forall (m :: * -> *) a b.
Monad m =>
Fold m a b -> Stream m a -> m b
Stream.fold Fold m b b
forall (m :: * -> *) a. (Monad m, Fractional a) => Fold m a a
Fold.mean
                    (Stream m b -> m b) -> Stream m b -> m b
forall a b. (a -> b) -> a -> b
$ (b -> b) -> Stream m b -> Stream m b
forall a b. (a -> b) -> Stream m a -> Stream m b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (\b
a -> b -> b
forall a. Num a => a -> a
abs (b
mn b -> b -> b
forall a. Num a => a -> a -> a
- b
a))
                    (Stream m b -> Stream m b) -> Stream m b -> Stream m b
forall a b. (a -> b) -> a -> b
$ RingArray b -> Stream m b
forall (m :: * -> *) a.
(MonadIO m, Unbox a) =>
RingArray a -> Stream m a
Ring.read RingArray b
rng
            Maybe (RingArray b)
Nothing -> b -> m b
forall a. a -> m a
forall (m :: * -> *) a. Monad m => a -> m a
return b
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 incrVariance #-}
incrVariance :: (Monad m, Fractional a) => Scanl m (Incr a) a
incrVariance :: forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
incrVariance =
    (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (m :: * -> *) b c d a.
Monad m =>
(b -> c -> d) -> Scanl m a b -> Scanl m a c -> Scanl m a d
Scanl.teeWith
        (\a
p2 a
m -> a
p2 a -> a -> a
forall a. Num a => a -> a -> a
- a
m a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)) (Int -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Int -> Scanl m (Incr a) a
incrRawMoment Int
2) Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
Scanl.incrMean

-- | 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 incrStdDev #-}
incrStdDev :: (Monad m, Floating a) => Scanl m (Incr a) a
incrStdDev :: forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrStdDev = a -> a
forall a. Floating a => a -> a
sqrt (a -> a) -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
incrVariance

-- XXX Need a tee3 operation for better performance.

-- | 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 incrSkewness #-}
incrSkewness :: (Monad m, Floating a) => Scanl m (Incr a) a
incrSkewness :: forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrSkewness =
          (\a
rm3 a
sd a
mu ->
            a
rm3 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
sd a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int) a -> a -> a
forall a. Num a => a -> a -> a
- a
3 a -> a -> a
forall a. Num a => a -> a -> a
* (a
mu a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
sd) a -> a -> a
forall a. Num a => a -> a -> a
- (a
mu a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
sd) a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
3 :: Int)
          )
        (a -> a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) (a -> a -> a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Int -> Scanl m (Incr a) a
incrRawMoment Int
3
        Scanl m (Incr a) (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) (a -> a)
forall a b.
Scanl m (Incr a) (a -> b)
-> Scanl m (Incr a) a -> Scanl m (Incr a) b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Scanl m (Incr a) a
forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrStdDev
        Scanl m (Incr a) (a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall a b.
Scanl m (Incr a) (a -> b)
-> Scanl m (Incr a) a -> Scanl m (Incr a) b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
Scanl.incrMean

-- XXX We can compute the 2nd, 3rd, 4th raw moments by repeatedly multiplying
-- instead of computing the powers every time.
-- XXX Need a tee4 operation for better performance.

-- | 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 incrKurtosis #-}
incrKurtosis :: (Monad m, Floating a) => Scanl m (Incr a) a
incrKurtosis :: forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrKurtosis =
          (\a
rm4 a
rm3 a
sd a
mu ->
             ( a
3 a -> a -> a
forall a. Num a => a -> a -> a
* a
mu a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
4 :: Int)
            a -> a -> a
forall a. Num a => a -> a -> a
+ a
6 a -> a -> a
forall a. Num a => a -> a -> a
* a
mu a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int) a -> a -> a
forall a. Num a => a -> a -> a
* a
sd a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
2 :: Int)
            a -> a -> a
forall a. Num a => a -> a -> a
- a
4 a -> a -> a
forall a. Num a => a -> a -> a
* a
mu a -> a -> a
forall a. Num a => a -> a -> a
* a
rm3
            a -> a -> a
forall a. Num a => a -> a -> a
+ a
rm4) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
sd a -> Int -> a
forall a b. (Num a, Integral b) => a -> b -> a
^ (Int
4 :: Int))
          )
        (a -> a -> a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) (a -> a -> a -> a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Int -> Scanl m (Incr a) a
incrRawMoment Int
4
        Scanl m (Incr a) (a -> a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) (a -> a -> a)
forall a b.
Scanl m (Incr a) (a -> b)
-> Scanl m (Incr a) a -> Scanl m (Incr a) b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int -> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Int -> Scanl m (Incr a) a
incrRawMoment Int
3
        Scanl m (Incr a) (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) (a -> a)
forall a b.
Scanl m (Incr a) (a -> b)
-> Scanl m (Incr a) a -> Scanl m (Incr a) b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Scanl m (Incr a) a
forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrStdDev
        Scanl m (Incr a) (a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall a b.
Scanl m (Incr a) (a -> b)
-> Scanl m (Incr a) a -> Scanl m (Incr a) b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
Scanl.incrMean

-------------------------------------------------------------------------------
-- 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 incrSampleVariance #-}
incrSampleVariance :: (Monad m, Fractional a) => Scanl m (Incr a) a
incrSampleVariance :: forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
incrSampleVariance =
    (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (m :: * -> *) b c d a.
Monad m =>
(b -> c -> d) -> Scanl m a b -> Scanl m a c -> Scanl m a d
Scanl.teeWith (\a
n a
s2 -> a
n a -> a -> a
forall a. Num a => a -> a -> a
* a
s2 a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
n a -> a -> a
forall a. Num a => a -> a -> a
- a
1)) Scanl m (Incr a) a
forall (m :: * -> *) b a. (Monad m, Num b) => Scanl m (Incr a) b
Scanl.incrCount Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
incrVariance

-- | Sample standard deviation:
--
-- \(s = \sqrt{sampleVariance}\)
--
-- >>> sampleStdDev = sqrt <$> sampleVariance
--
-- See https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation
-- .
--
{-# INLINE incrSampleStdDev #-}
incrSampleStdDev :: (Monad m, Floating a) => Scanl m (Incr a) a
incrSampleStdDev :: forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrSampleStdDev = a -> a
forall a. Floating a => a -> a
sqrt (a -> a) -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Scanl m (Incr a) a
forall (m :: * -> *) a.
(Monad m, Fractional a) =>
Scanl m (Incr a) a
incrSampleVariance

-- | 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 incrStdErrMean #-}
incrStdErrMean :: (Monad m, Floating a) => Scanl m (Incr a) a
incrStdErrMean :: forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrStdErrMean =
    (a -> a -> a)
-> Scanl m (Incr a) a -> Scanl m (Incr a) a -> Scanl m (Incr a) a
forall (m :: * -> *) b c d a.
Monad m =>
(b -> c -> d) -> Scanl m a b -> Scanl m a c -> Scanl m a d
Scanl.teeWith (\a
sd a
n -> a
sd a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt a
n) Scanl m (Incr a) a
forall (m :: * -> *) a. (Monad m, Floating a) => Scanl m (Incr a) a
incrSampleStdDev Scanl m (Incr a) a
forall (m :: * -> *) b a. (Monad m, Num b) => Scanl m (Incr a) b
Scanl.incrCount

-------------------------------------------------------------------------------
-- 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 incrFrequency #-}
incrFrequency :: (Monad m, Ord a) => Scanl m (Incr a) (Map a Int)
incrFrequency :: forall (m :: * -> *) a.
(Monad m, Ord a) =>
Scanl m (Incr a) (Map a Int)
incrFrequency = (Map a Int -> Incr a -> Map a Int)
-> Map a Int -> Scanl m (Incr a) (Map a Int)
forall (m :: * -> *) b a.
Monad m =>
(b -> a -> b) -> b -> Scanl m a b
Scanl.mkScanl Map a Int -> Incr a -> Map a Int
forall {k} {a}.
(Ord k, Num a, Eq a) =>
Map k a -> Incr k -> Map k a
step Map a Int
forall k a. Map k a
Map.empty

    where

    decrement :: a -> Maybe a
decrement a
v =
        if a
v a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1
        then Maybe a
forall a. Maybe a
Nothing
        else a -> Maybe a
forall a. a -> Maybe a
Just (a
v a -> a -> a
forall a. Num a => a -> a -> a
- a
1)

    step :: Map k a -> Incr k -> Map k a
step Map k a
refCountMap (Insert k
new) =
        (a -> a -> a) -> k -> a -> Map k a -> Map k a
forall k a. Ord k => (a -> a -> a) -> k -> a -> Map k a -> Map k a
Map.insertWith a -> a -> a
forall a. Num a => a -> a -> a
(+) k
new a
1 Map k a
refCountMap

    step Map k a
refCountMap (Replace k
old k
new) =
        let m1 :: Map k a
m1 = (a -> a -> a) -> k -> a -> Map k a -> Map k a
forall k a. Ord k => (a -> a -> a) -> k -> a -> Map k a -> Map k a
Map.insertWith a -> a -> a
forall a. Num a => a -> a -> a
(+) k
new a
1 Map k a
refCountMap
         in (a -> Maybe a) -> k -> Map k a -> Map k a
forall k a. Ord k => (a -> Maybe a) -> k -> Map k a -> Map k a
Map.update a -> Maybe a
forall {a}. (Eq a, Num a) => a -> Maybe a
decrement k
old Map k a
m1