online 0.5.0 → 0.6.0

raw patch · 12 files changed

+874/−588 lines, 12 filesdep +adjunctionsdep +containersdep +foldsdep −tastydep ~backpropsetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: adjunctions, containers, folds, generic-lens, hmatrix, lens, mwc-probability, mwc-random, numhask, numhask-array, online, primitive, profunctors, text

Dependencies removed: tasty

Dependency ranges changed: backprop

API changes (from Hackage documentation)

Files

• LICENSE +0/−30
• Setup.hs +0/−2
• online.cabal +65/−33
• src/Data/Mealy.hs +596/−0
• src/Data/Quantiles.hs +145/−0
• src/Data/Simulate.hs +59/−0
• src/Online.hs +0/−12
• src/Online/Averages.hs +0/−169
• src/Online/AveragesB.hs +0/−64
• src/Online/Medians.hs +0/−125
• src/Online/Quantiles.hs +0/−141
• test/test.hs +9/−12

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

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

@@ -1,48 +1,80 @@-cabal-version: 3.0-name:           online-version:        0.5.0-synopsis:       online statistics-description:    transformation of statistics to online algorithms-category:       statistics-homepage:       https://github.com/tonyday567/online#readme-bug-reports:    https://github.com/tonyday567/online/issues-author:         Tony Day-maintainer:     tonyday567@gmail.com-copyright:      Tony Day-license:        BSD-3-Clause-license-file:   LICENSE-build-type:     Simple-+cabal-version: 2.4+name:          online+version:       0.6.0+synopsis: See readme.md+description: See readme.md for description.+category: project+author: Tony Day+maintainer: tonyday567@gmail.com+copyright: Tony Day (c) AfterTimes+license: BSD-3-Clause+homepage: https://github.com/tonyday567/online#readme+bug-reports: https://github.com/tonyday567/online/issues+build-type: Simple source-repository head   type: git   location: https://github.com/tonyday567/online  library-  exposed-modules:-      Online-      Online.Averages-      Online.AveragesB-      Online.Medians-      Online.Quantiles   hs-source-dirs:-      src-  ghc-options: -Wall -Wcompat -Wincomplete-record-updates -Wincomplete-uni-patterns -Wredundant-constraints+    src   build-depends:-      base >=4.7 && <5-    , foldl-    , tdigest-    , vector-    , vector-algorithms-    , backprop+    adjunctions >= 4.4,+    backprop >= 0.2.6.4 && < 0.3,+    base >=4.7 && <5,+    containers >= 0.6,+    foldl,+    folds,+    generic-lens >= 2.0,+    hmatrix >= 0.20,+    lens,+    mwc-probability,+    mwc-random,+    numhask >= 0.6 && < 0.7,+    numhask-array >= 0.7 && < 0.8,+    primitive >= 0.7,+    profunctors >= 5.5,+    tdigest,+    text,+    vector,+    vector-algorithms+  exposed-modules:+    Data.Mealy+    Data.Quantiles+    Data.Simulate+  other-modules:   default-language: Haskell2010+  default-extensions:+    NoImplicitPrelude+    NegativeLiterals+    OverloadedStrings+    UnicodeSyntax+  ghc-options:+    -Wall+    -Wcompat+    -Wincomplete-record-updates+    -Wincomplete-uni-patterns+    -Wredundant-constraints  test-suite test   type: exitcode-stdio-1.0   main-is: test.hs   hs-source-dirs:-      test+    test   build-depends:-      base >=4.7 && <5-    , doctest-    , tasty+    base >=4.7 && <5,+    doctest,+    numhask >= 0.6 && < 0.7,+    online   default-language: Haskell2010+  default-extensions:+    NoImplicitPrelude+    NegativeLiterals+    OverloadedStrings+    UnicodeSyntax+  ghc-options:+    -Wall+    -Wcompat+    -Wincomplete-record-updates+    -Wincomplete-uni-patterns+    -Wredundant-constraints

@@ -0,0 +1,596 @@+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DerivingVia #-}+{-# LANGUAGE DuplicateRecordFields #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE OverloadedLabels #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# OPTIONS_GHC -Wall #-}+{-# OPTIONS_GHC -Wno-name-shadowing #-}+{-# OPTIONS_GHC -Wno-orphans #-}+{-# OPTIONS_GHC -Wno-type-defaults #-}+{-# OPTIONS_GHC -Wno-incomplete-patterns #-}++-- | Online statistics for ordered data (such as time-series data), modelled as [mealy machines](https://en.wikipedia.org/wiki/Mealy_machine)+module Data.Mealy+  ( -- * Types+    Mealy (..),+    pattern M,+    scan,+    fold,+    Averager (..),+    pattern A,+    av,+    av_,+    online,++    -- * Statistics+    -- $setup+    ma,+    absma,+    sqma,+    std,+    cov,+    corrGauss,+    corr,+    beta1,+    alpha1,+    reg1,+    beta,+    alpha,+    reg,+    asum,+    aconst,+    delay1,+    delay,+    depState,+    Model1 (..),+    zeroModel1,+    depModel1,++    -- * conversion+    fromFoldl,+    foldB,+    maB,++    -- * median+    Medianer (..),+    onlineL1,+    onlineL1',+    maL1,+    absmaL1,+  )+where++import qualified Control.Foldl as L+import Control.Lens hiding ((:>), Empty, Unwrapped, Wrapped, index, (|>))+import Data.Fold hiding (M)+import Data.Functor.Rep+import Data.Generics.Labels ()+import qualified Data.Sequence as Seq+import qualified NumHask.Array.Fixed as F+import qualified NumHask.Array.HMatrix as HM+import NumHask.Array.Shape (HasShape)+import NumHask.Prelude hiding (L1, State, StateT, asum, fold, get, replace, runState, runStateT, state)+import qualified Numeric.Backprop as B+import Numeric.Backprop (BVar, Reifies, W)+import qualified Numeric.LinearAlgebra as LA+import qualified Prelude.Backprop as PB+import qualified Prelude as P++-- $setup+-- Generate some random variates for the examples.+--+-- xs0, xs1 & xs2 are samples from N(0,1)+--+-- xsp is a pair of N(0,1)s with a correlation of 0.8+--+-- >>> :set -XDataKinds+-- >>> import Control.Category ((>>>))+-- >>> import Data.List+-- >>> import Data.Simulate+-- >>> g <- create+-- >>> xs0 <- rvs g 10000+-- >>> xs1 <- rvs g 10000+-- >>> xs2 <- rvs g 10000+-- >>> xsp <- rvsp g 10000 0.8++-- | A 'Mealy' is a triple of functions+--+-- * (c -> b) __extract__ Convert state to the output type.+-- * (c -> a -> c) __step__ Update state given prior state and (new) input.+-- * (a -> c) __inject__ Convert an input into the state type.+--+-- The type is a newtype wrapper around 'L1' in 'Data.Fold'.+--+-- __inject__ is necessary to kick off state on a 'fold' or 'scan', rather than a state existing prior to the fold or scan (this is a Moore machine or 'M1' in 'Data.Fold').+--+-- > scan (M e s i) (x : xs) = e <$> scanl' s (i x) xs+newtype Mealy a b = Mealy {l1 :: L1 a b}+  deriving (Profunctor, Category) via L1+  deriving (Functor, Applicative) via L1 a++-- | Pattern for a 'Mealy'.+--+-- @M extract step inject@+pattern M :: (c -> b) -> (c -> a -> c) -> (a -> c) -> Mealy a b+pattern M e s i = Mealy (L1 e s i)++{-# COMPLETE M #-}++-- | Fold a list through a 'Mealy'.+--+-- > cosieve == fold+fold :: Mealy a b -> [a] -> b+fold _ [] = panic "on the streets of Birmingham."+fold (M e s i) (x : xs) = e $ foldl' s (i x) xs++-- | Run a list through a 'Mealy' and return a list of values for every step+--+-- > length (scan _ xs) == length xs+scan :: Mealy a b -> [a] -> [b]+scan _ [] = []+scan (M e s i) (x : xs) = fromList (e <$> scanl' s (i x) xs)++-- | Most common statistics are averages, which are some sort of aggregation of values (sum) and some sort of sample size (count).+newtype Averager a b+  = Averager+      { sumCount :: (a, b)+      }+  deriving (Eq, Show)++-- | Pattern for an 'Averager'.+--+-- @A sum count@+pattern A :: a -> b -> Averager a b+pattern A s c = Averager (s, c)++{-# COMPLETE A #-}++instance (Additive a, Additive b) => Semigroup (Averager a b) where+  (<>) (A s c) (A s' c') = A (s + s') (c + c')++-- |+-- > av mempty == nan+instance (Additive a, Additive b) => Monoid (Averager a b) where+  mempty = A zero zero+  mappend = (<>)++-- | extract the average from an 'Averager'+--+-- av gives NaN on zero divide+av :: (Divisive a) => Averager a a -> a+av (A s c) = s / c++-- | substitute a default value on zero-divide+--+-- > av_ (Averager (0,0)) x == x+av_ :: (Eq a, Additive a, Divisive a) => Averager a a -> a -> a+av_ (A s c) def = bool def (s / c) (c == zero)++-- | @online f g@ is a 'Mealy' where f is a transformation of the data and g is a decay function (convergent tozero) applied at each step.+--+-- > online id id == av+online :: (Divisive b, Additive b) => (a -> b) -> (b -> b) -> Mealy a b+online f g = M av step intract+  where+    intract a = A (f a) one+    step (A s c) a =+      let (A s' c') = intract a+       in A (g s + s') (g c + c')++-- | A moving average using a decay rate of r. r=1 represents the simple average, and r=0 represents the latest value.+--+-- >>> fold (ma 0) (fromList [1..100])+-- 100.0+--+-- >>> fold (ma 1) (fromList [1..100])+-- 50.5+--+-- >>> fold (ma 0.99) xs0+-- -4.292501077490672e-2+--+-- A change in the underlying mean at n=10000 in the chart below highlights the trade-off between stability of the statistic and response to non-stationarity.+--+-- ![ma chart](other/ex-ma.svg)+ma :: (Divisive a, Additive a) => a -> Mealy a a+ma r = online id (* r)+{-# INLINEABLE ma #-}++-- | absolute average+--+-- >>> fold (absma 1) xs0+-- 0.7894201075535578+absma :: (Divisive a, Additive a, Signed a) => a -> Mealy a a+absma r = online abs (* r)+{-# INLINEABLE absma #-}++-- | average square+--+-- > fold (ma r) . fmap (**2) == fold (sqma r)+sqma :: (Divisive a, Additive a) => a -> Mealy a a+sqma r = online (\x -> x * x) (* r)+{-# INLINEABLE sqma #-}++-- | standard deviation+--+-- The construction of standard deviation, using the Applicative instance of a 'Mealy':+--+-- > (\s ss -> sqrt (ss - s ** (one+one))) <$> ma r <*> sqma r+--+-- The average deviation of the numbers 1..1000 is about 1 / sqrt 12 * 1000+-- <https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)#Standard_uniform>+--+-- >>> fold (std 1) [0..1000]+-- 288.9636655359978+--+-- The average deviation with a decay of 0.99+--+-- >>> fold (std 0.99) [0..1000]+-- 99.28328803163829+--+-- >>> fold (std 1) xs0+-- 0.9923523681261158+--+-- ![std chart](other/ex-std.svg)+std :: (Divisive a, ExpField a) => a -> Mealy a a+std r = (\s ss -> sqrt (ss - s ** (one + one))) <$> ma r <*> sqma r+{-# INLINEABLE std #-}++-- | The covariance of a tuple given an underlying central tendency fold.+--+-- >>> fold (cov (ma 1)) xsp+-- 0.8011368250045314+cov :: (Field a) => Mealy a a -> Mealy (a, a) a+cov m =+  (\xy x' y' -> xy - x' * y') <$> lmap (uncurry (*)) m <*> lmap fst m <*> lmap snd m+{-# INLINEABLE cov #-}++-- | correlation of a tuple, specialised to Guassian+--+-- >>> fold (corrGauss 1) xsp+-- 0.8020637696465039+corrGauss :: (ExpField a) => a -> Mealy (a, a) a+corrGauss r =+  (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> cov (ma r)+    <*> lmap fst (std r)+    <*> lmap snd (std r)+{-# INLINEABLE corrGauss #-}++-- | a generalised version of correlation of a tuple+--+-- >>> fold (corr (ma 1) (std 1)) xsp+-- 0.8020637696465039+--+-- > corr (ma r) (std r) == corrGauss r+corr :: (ExpField a) => Mealy a a -> Mealy a a -> Mealy (a, a) a+corr central deviation =+  (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> cov central+    <*> lmap fst deviation+    <*> lmap snd deviation+{-# INLINEABLE corr #-}++-- | The beta in a simple linear regression of an (independent variable, single dependent variable) tuple given an underlying central tendency fold.+--+-- This is a generalisation of the classical regression formula, where averages are replaced by 'Mealy' statistics.+--+-- \[+-- \begin{align}+-- \beta & = \frac{n\sum xy - \sum x \sum y}{n\sum x^2 - (\sum x)^2} \\+--     & = \frac{n^2 \overline{xy} - n^2 \bar{x} \bar{y}}{n^2 \overline{x^2} - n^2 \bar{x}^2} \\+--     & = \frac{\overline{xy} - \bar{x} \bar{y}}{\overline{x^2} - \bar{x}^2} \\+-- \end{align}+-- \]+--+-- >>> fold (beta1 (ma 1)) $ zipWith (\x y -> (y, x + y)) xs0 xs1+-- 0.9953875263096014+beta1 :: (ExpField a) => Mealy a a -> Mealy (a, a) a+beta1 m =+  (\xy x' y' x2 -> (xy - x' * y') / (x2 - x' * x')) <$> lmap (uncurry (*)) m+    <*> lmap fst m+    <*> lmap snd m+    <*> lmap (\(x, _) -> x * x) m+{-# INLINEABLE beta1 #-}++-- | The alpha in a simple linear regression of an (independent variable, single dependent variable) tuple given an underlying central tendency fold.+--+-- \[+-- \begin{align}+-- \alpha & = \frac{\sum y \sum x^2 - \sum x \sum xy}{n\sum x^2 - (\sum x)^2} \\+--     & = \frac{n^2 \bar{y} \overline{x^2} - n^2 \bar{x} \overline{xy}}{n^2 \overline{x^2} - n^2 \bar{x}^2} \\+--     & = \frac{\bar{y} \overline{x^2} - \bar{x} \overline{xy}}{\overline{x^2} - \bar{x}^2} \\+-- \end{align}+-- \]+--+-- >>> fold (alpha1 (ma 1)) $ zipWith (\x y -> ((3+y), x + 0.5 * (3 + y))) xs0 xs1+-- 1.1880996822796197e-2+alpha1 :: (ExpField a) => Mealy a a -> Mealy (a, a) a+alpha1 m = (\x b y -> y - b * x) <$> lmap fst m <*> beta1 m <*> lmap snd m+{-# INLINEABLE alpha1 #-}++-- | The (alpha, beta) tuple in a simple linear regression of an (independent variable, single dependent variable) tuple given an underlying central tendency fold.+--+-- >>> fold (reg1 (ma 1)) $ zipWith (\x y -> ((3+y), x + 0.5 * (3 + y))) xs0 xs1+-- (1.1880996822796197e-2,0.49538752630956845)+reg1 :: (ExpField a) => Mealy a a -> Mealy (a, a) (a, a)+reg1 m = (,) <$> alpha1 m <*> beta1 m++data RegressionState (n :: Nat) a+  = RegressionState+      { _xx :: F.Array '[n, n] a,+        _x :: F.Array '[n] a,+        _xy :: F.Array '[n] a,+        _y :: a+      }+  deriving (Functor)++-- | multiple regression+--+-- \[+-- \begin{align}+-- {\hat  {{\mathbf  {B}}}}=({\mathbf  {X}}^{{{\rm {T}}}}{\mathbf  {X}})^{{ -1}}{\mathbf  {X}}^{{{\rm {T}}}}{\mathbf  {Y}}+-- \end{align}+-- \]+--+-- \[+-- \begin{align}+-- {\mathbf  {X}}={\begin{bmatrix}{\mathbf  {x}}_{1}^{{{\rm {T}}}}\\{\mathbf  {x}}_{2}^{{{\rm {T}}}}\\\vdots \\{\mathbf  {x}}_{n}^{{{\rm {T}}}}\end{bmatrix}}={\begin{bmatrix}x_{{1,1}}&\cdots &x_{{1,k}}\\x_{{2,1}}&\cdots &x_{{2,k}}\\\vdots &\ddots &\vdots \\x_{{n,1}}&\cdots &x_{{n,k}}\end{bmatrix}}+-- \end{align}+-- \]+--+-- > let ys = zipWith3 (\x y z -> 0.1 * x + 0.5 * y + 1 * z) xs0 xs1 xs2+-- > let zs = zip (zipWith (\x y -> fromList [x,y] :: F.Array '[2] Double) xs1 xs2) ys+-- > fold (beta 0.99) zs+-- [0.4982692361226971, 1.038192474255091]+beta :: (Field a, LA.Field a, KnownNat n) => a -> Mealy (F.Array '[n] a, a) (F.Array '[n] a)+beta r = M extract step inject+  where+    extract (A (RegressionState xx x xy y) c) =+      liftHM2+        (\a b -> LA.pinv a LA.<> LA.tr b)+        ((one / c) *. (xx - F.expand (*) x x))+        ((xy - (y *. x)) .* (one / c))+    step x (xs, y) = rsOnline r x (inject (xs, y))+    inject (xs, y) =+      A (RegressionState (F.expand (*) xs xs) xs (y *. xs) y) one+{-# INLINEABLE beta #-}++liftHM2 :: (LA.Field a, HasShape s, HasShape s', HasShape s'') => (LA.Matrix a -> LA.Matrix a -> LA.Matrix a) -> F.Array s a -> F.Array s' a -> F.Array s'' a+liftHM2 f a b = HM.toFixed $ HM.Array $ f (HM.unArray . HM.fromFixed $ a) (HM.unArray . HM.fromFixed $ b)++rsOnline :: (Field a, KnownNat n) => a -> Averager (RegressionState n a) a -> Averager (RegressionState n a) a -> Averager (RegressionState n a) a+rsOnline r (A (RegressionState xx x xy y) c) (A (RegressionState xx' x' xy' y') c') =+  A (RegressionState (liftR2 d xx xx') (liftR2 d x x') (liftR2 d xy xy') (d y y')) (d c c')+  where+    d s s' = r * s + s'++-- | alpha in a multiple regression+alpha :: (LA.Field a, ExpField a, KnownNat n) => a -> Mealy (F.Array '[n] a, a) a+alpha r = (\xs b y -> y - sum (liftR2 (*) b xs)) <$> lmap fst (arrayify $ ma r) <*> beta r <*> lmap snd (ma r)+{-# INLINEABLE alpha #-}++arrayify :: (HasShape s) => Mealy a b -> Mealy (F.Array s a) (F.Array s b)+arrayify (M sExtract sStep sInject) = M extract step inject+  where+    extract = fmap sExtract+    step = liftR2 sStep+    inject = fmap sInject++-- | multiple regression+--+-- > let ys = zipWith3 (\x y z -> 0.1 * x + 0.5 * y + 1 * z) xs0 xs1 xs2+-- > let zs = zip (zipWith (\x y -> fromList [x,y] :: F.Array '[2] Double) xs1 xs2) ys+-- > fold (reg 0.99) zs+-- ([0.4982692361226971, 1.038192474255091],2.087160803386695e-3)+reg :: (LA.Field a, ExpField a, KnownNat n) => a -> Mealy (F.Array '[n] a, a) (F.Array '[n] a, a)+reg r = (,) <$> beta r <*> alpha r+{-# INLINEABLE reg #-}++-- | accumulated sum+asum :: (Additive a) => Mealy a a+asum = M id (+) id++-- | constant Mealy+aconst :: b -> Mealy a b+aconst a = M (const a) (\_ _ -> ()) (const ())++-- | delay input values by 1+delay1 :: a -> Mealy a a+delay1 x0 = M fst (\(_, x) a -> (x, a)) (x0,)++-- | delays values by n steps+--+-- delay [0] == delay1 0+--+-- delay [] == id+--+-- delay [1,2] = delay1 2 . delay1 1+--+-- >>> scan (delay [-2,-1]) [0..3]+-- [-2,-1,0,1]+--+-- Autocorrelation example:+--+-- > scan (((,) <$> id <*> delay [0]) >>> beta (ma 0.99)) xs0+delay ::+  -- | initial statistical values, delay equals length+  [a] ->+  Mealy a a+delay x0 = M extract step inject+  where+    inject a = Seq.fromList x0 Seq.|> a+    extract :: Seq a -> a+    extract Seq.Empty = panic "ACAB"+    extract (x Seq.:<| _) = x+    step :: Seq a -> a -> Seq a+    step Seq.Empty _ = panic "ACAB"+    step (_ Seq.:<| xs) a = xs Seq.|> a++-- | Add a state dependency to a series.+--+-- Typical regression analytics tend to assume that moments of a distributional assumption are unconditional with respect to prior instantiations of the stochastics being studied.+--+-- For time series analytics, a major preoccupation is estimation of the current moments given what has happened in the past.+--+-- IID:+--+-- \[+-- \begin{align}+-- x_{t+1} & = alpha_t^x + s_{t+1}\\+-- s_{t+1} & = alpha_t^s * N(0,1)+-- \end{align}+-- \]+--+-- Example: including a linear dependency on moving average history:+--+-- \[+-- \begin{align}+-- x_{t+1} & = (alpha_t^x + beta_t^{x->x} * ma_t^x) + s_{t+1}\\+-- s_{t+1} & = alpha_t^s * N(0,1)+-- \end{align}+-- \]+--+-- >>> let xs' = scan (depState (\a m -> a + 0.1 * m) (ma 0.99)) xs0+-- >>> let ma' = scan ((ma (1 - 0.01)) >>> delay [0]) xs'+-- >>> let xsb = fold (beta1 (ma (1 - 0.001))) $ drop 1 $ zip ma' xs'+-- >>> -- beta measurement if beta of ma was, in reality, zero.+-- >>> let xsb0 = fold (beta1 (ma (1 - 0.001))) $ drop 1 $ zip ma' xs0+-- >>> xsb - xsb0+-- 9.999999999999976e-2+--+-- This simple model of relationship between a series and it's historical average shows how fragile the evidence can be.+--+-- ![madep](other/ex-madep.svg)+--+-- In unravelling the drivers of this result, the standard deviation of a moving average scan seems well behaved for r > 0.01, but increases substantively for values less than this.  This result seems to occur for wide beta values. For high r, the standard deviation of the moving average seems to be proprtional to r**0.5, and equal to around (0.5*r)**0.5.+--+-- > fold (std 1) (scan (ma (1 - 0.01)) xs0)+--+-- ![stdma](other/ex-stdma.svg)+depState :: (a -> b -> a) -> Mealy a b -> Mealy a a+depState f (M sExtract sStep sInject) = M fst step inject+  where+    inject a = (a, sInject a)+    step (_, m) a = let a' = f a (sExtract m) in (a', sStep m a')++-- | a linear model of state dependencies for the first two moments+--+-- \[+-- \begin{align}+-- x_{t+1} & = (alpha_t^x + beta_t^{x->x} * ma_t^x + beta_t^{s->x} * std_t^x) + s_{t+1}\\+-- s_{t+1} & = (alpha_t^s + beta_t^{x->s} * ma_t^x + beta_t^{s->s} * std_t^x) * N(0,1)+-- \end{align}+-- \]+data Model1+  = Model1+      { alphaX :: Double,+        alphaS :: Double,+        betaMa2X :: Double,+        betaMa2S :: Double,+        betaStd2X :: Double,+        betaStd2S :: Double+      }+  deriving (Eq, Show, Generic)++zeroModel1 :: Model1+zeroModel1 = Model1 0 0 0 0 0 0++-- | Apply a model1 relationship using a single decay factor.+--+-- >>> :set -XOverloadedLabels+-- >>> fold (depModel1 0.01 (zeroModel1 & #betaMa2X .~ 0.1)) xs0+-- -0.47228537123218206+depModel1 :: Double -> Model1 -> Mealy Double Double+depModel1 r m1 =+  depState fX st+  where+    st = (,) <$> ma (1 - r) <*> std (1 - r)+    fX a (m, s) =+      a+        * ( (1 + m1 ^. #alphaS)+              + (m1 ^. #betaMa2S) * m+              + (m1 ^. #betaStd2S) * (s - 1)+          )+        + m1 ^. #alphaX+        + (m1 ^. #betaMa2X) * m+        + (m1 ^. #betaStd2X) * (s - 1)++fromFoldl :: L.Fold a b -> Mealy a b+fromFoldl (L.Fold step begin e) = M e step (step begin)++foldB :: (Reifies s W) => (BVar s Double -> BVar s Double) -> BVar s Double -> BVar s [Double] -> BVar s Double+foldB f r xs = divide (PB.foldl' (step' f r) (B.T2 0 0) xs)+  where+    step' f' r' (B.T2 s c) a = uncurry B.T2 ((r P.*) $ s P.+ f' a, (r' P.*) $ c P.+ 1)+    divide (B.T2 s c) = s P./ c++maB :: Reifies s W => BVar s Double -> BVar s [Double] -> BVar s Double+maB = foldB id++-- | A rough Median.+-- The average absolute value of the stat is used to callibrate estimate drift towards the median+data Medianer a b+  = Medianer+      { medAbsSum :: a,+        medCount :: b,+        medianEst :: a+      }++-- | onlineL1' takes a function and turns it into a `Control.Foldl.Fold` where the step is an incremental update of an (isomorphic) median statistic.+onlineL1' ::+  (Ord b, Field b, Signed b) => b -> b -> (a -> b) -> (b -> b) -> L.Fold a (b, b)+onlineL1' i d f g = L.Fold step begin extract+  where+    begin = Medianer zero zero zero+    step (Medianer s c m) a =+      Medianer+        (g $ s + abs (f a))+        (g $ c + one)+        ((one - d) * (m + s' * i * s / c') + d * f a)+      where+        c' =+          if c == zero+            then one+            else c+        s'+          | f a > m = one+          | f a < m = negate one+          | otherwise = zero+    extract (Medianer s c m) = (s / c, m)+{-# INLINEABLE onlineL1' #-}++-- | onlineL1 takes a function and turns it into a `Control.Foldl.Fold` where the step is an incremental update of an (isomorphic) median statistic.+onlineL1 :: (Ord b, Field b, Signed b) => b -> b -> (a -> b) -> (b -> b) -> L.Fold a b+onlineL1 i d f g = snd <$> onlineL1' i d f g+{-# INLINEABLE onlineL1 #-}++-- $setup+--+-- >>> import qualified Control.Foldl as L+-- >>> let n = 100+-- >>> let inc = 0.1+-- >>> let d = 0+-- >>> let r = 0.9++-- | moving median+-- > L.fold (maL1 inc d r) [1..n]+-- 93.92822312742108+maL1 :: (Ord a, Field a, Signed a) => a -> a -> a -> L.Fold a a+maL1 i d r = onlineL1 i d id (* r)+{-# INLINEABLE maL1 #-}++-- | moving absolute deviation+absmaL1 :: (Ord a, Field a, Signed a) => a -> a -> a -> L.Fold a a+absmaL1 i d r = fst <$> onlineL1' i d id (* r)+{-# INLINEABLE absmaL1 #-}

@@ -0,0 +1,145 @@+{-# LANGUAGE DataKinds #-}++module Data.Quantiles+  ( tDigest,+    tDigestQuantiles,+    tDigestHist,+    OnlineTDigest (..),+    onlineQuantiles,+    Data.Quantiles.median,+    onlineDigitize,+    onlineDigestHist,+  )+where++import qualified Control.Foldl as L+import Control.Monad.ST (runST)+import Data.Foldable+import Data.List.NonEmpty (NonEmpty)+import Data.Maybe+import Data.Ord+import Data.TDigest+import Data.TDigest.Internal+import Data.TDigest.Postprocess (HistBin, histogram)+import Data.TDigest.Tree.Internal (TDigest (..), absMaxSize, emptyTDigest, insertCentroid, relMaxSize, size, toMVector)+import qualified Data.Vector.Algorithms.Heap as VHeap+import qualified Data.Vector.Unboxed as VU+import Prelude++-- | a raw non-online tdigest fold+tDigest :: L.Fold Double (TDigest 25)+tDigest = L.Fold step begin done+  where+    step x a = insert a x+    begin = tdigest ([] :: [Double]) :: TDigest 25+    done = id++-- | non-online version+tDigestQuantiles :: [Double] -> L.Fold Double [Double]+tDigestQuantiles qs = L.Fold step begin done+  where+    step x a = insert a x+    begin = tdigest ([] :: [Double]) :: TDigest 25+    done x = fromMaybe (0 / 0) . (`quantile` compress x) <$> qs++-- | non-online version+tDigestHist :: L.Fold Double (Maybe (NonEmpty HistBin))+tDigestHist = L.Fold step begin done+  where+    step x a = insert a x+    begin = tdigest ([] :: [Double]) :: TDigest 25+    done = histogram . compress++data OnlineTDigest+  = OnlineTDigest+      { td :: TDigest 25,+        tdN :: Int,+        tdRate :: Double+      }+  deriving (Show)++emptyOnlineTDigest :: Double -> OnlineTDigest+emptyOnlineTDigest = OnlineTDigest (emptyTDigest :: TDigest n) 0++-- | decaying quantiles based on the tdigest library+onlineQuantiles :: Double -> [Double] -> L.Fold Double [Double]+onlineQuantiles r qs = L.Fold step begin done+  where+    step x a = onlineInsert a x+    begin = emptyOnlineTDigest r+    done x = fromMaybe (0 / 0) . (`quantile` t) <$> qs+      where+        (OnlineTDigest t _ _) = onlineForceCompress x++median :: Double -> L.Fold Double Double+median r = L.Fold step begin done+  where+    step x a = onlineInsert a x+    begin = emptyOnlineTDigest r+    done x = fromMaybe (0 / 0) (quantile 0.5 t)+      where+        (OnlineTDigest t _ _) = onlineForceCompress x++onlineInsert' :: Double -> OnlineTDigest -> OnlineTDigest+onlineInsert' x (OnlineTDigest td' n r) =+  OnlineTDigest+    (insertCentroid (x, r ^^ (- (fromIntegral $ n + 1) :: Integer)) td')+    (n + 1)+    r++onlineInsert :: Double -> OnlineTDigest -> OnlineTDigest+onlineInsert x otd = onlineCompress (onlineInsert' x otd)++onlineCompress :: OnlineTDigest -> OnlineTDigest+onlineCompress otd@(OnlineTDigest Nil _ _) = otd+onlineCompress otd@(OnlineTDigest t _ _)+  | Data.TDigest.Tree.Internal.size t > relMaxSize * compression+      && Data.TDigest.Tree.Internal.size t > absMaxSize =+    onlineForceCompress otd+  | otherwise = otd+  where+    compression = 25++onlineForceCompress :: OnlineTDigest -> OnlineTDigest+onlineForceCompress otd@(OnlineTDigest Nil _ _) = otd+onlineForceCompress (OnlineTDigest t n r) = OnlineTDigest t' 0 r+  where+    t' =+      foldl' (flip insertCentroid) emptyTDigest $+        (\(m, w) -> (m, w * (r ^^ n))) . fst <$> VU.toList centroids+    -- Centroids are shuffled based on space+    centroids :: VU.Vector (Centroid, Double)+    centroids =+      runST $ do+        v <- toMVector t+        -- sort by cumulative weight+        VHeap.sortBy (comparing snd) v+        VU.unsafeFreeze v++onlineDigitize :: Double -> [Double] -> L.Fold Double Int+onlineDigitize r qs = L.Fold step begin done+  where+    step (x, _) a = (onlineInsert a x, a)+    begin = (emptyOnlineTDigest r, 0 / 0)+    done (x, l) = bucket' qs' l+      where+        qs' = fromMaybe (0 / 0) . (`quantile` t) <$> qs+        (OnlineTDigest t _ _) = onlineForceCompress x+        bucket' xs l' =+          L.fold L.sum $+            ( \x' ->+                if x' > l'+                  then 0+                  else 1+            )+              <$> xs++-- | decaying histogram based on the tdigest library+onlineDigestHist :: Double -> L.Fold Double (Maybe (NonEmpty HistBin))+onlineDigestHist r = L.Fold step begin done+  where+    step x a = onlineInsert a x+    begin = emptyOnlineTDigest r+    done x = histogram . compress $ t+      where+        (OnlineTDigest t _ _) = onlineForceCompress x

@@ -0,0 +1,59 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE ExtendedDefaultRules #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE NoImplicitPrelude #-}+{-# OPTIONS_GHC -Wall #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Data.Simulate+  ( rvs,+    rvsp,+    create+  )+where++import Control.Monad.Primitive (PrimState)+import NumHask.Prelude hiding (fold)+import System.Random.MWC+import System.Random.MWC.Probability hiding (beta)++-- $setup+-- >>> :set -XFlexibleContexts+-- >>> import Data.Mealy+-- >>> gen <- create+-- >>> let n = 3+-- >>> let eq' a b = all nearZero $ zipWith (-) a b+-- >>> let eq'p a b = all (\x -> x) $ zipWith (\(x0,x1) (y0,y1) -> nearZero (x0-y0) && nearZero (x1-y1)) a b++-- | rvs creates a list of standard normal random variates.+-- >>> t <- rvs gen n+-- >>> t `eq'` [-0.8077385934202513,-1.3423948150518445,-0.4900206084002882]+-- True+--+-- >>> rs <- rvs gen 10000+-- >>> fold (ma 1) rs+-- -1.735737734197327e-3+--+-- >>> fold (std 1) rs+-- 0.9923615647768976+rvs :: Gen (PrimState IO) -> Int -> IO [Double]+rvs gen n = samples n standardNormal gen++-- | rvsPair generates a list of correlated random variate tuples+-- |+-- >>> t <- rvsp gen 3 0.8+-- >>> t `eq'p` [(-0.8077385934202513,-1.4591410449385904),(-1.3423948150518445,-0.6046212701237168),(-0.4900206084002882,0.923007518547542)]+-- True+--+-- >>> rsp <- rvsp gen 10000 0.8+-- >>> fold (corr (ma 1) (std 1)) rsp+-- 0.7933213647252008+rvsp :: Gen (PrimState IO) -> Int -> Double -> IO [(Double, Double)]+rvsp gen n c = do+  s0 <- rvs gen n+  s1 <- rvs gen n+  let s1' = zipWith (\x y -> c * x + sqrt (1 - c * c) * y) s0 s1+  pure $ zip s0 s1'

@@ -1,12 +0,0 @@--- | online library-module Online-  ( module Online.Averages-  , module Online.AveragesB-  , module Online.Medians-  , module Online.Quantiles-  ) where--import Online.Quantiles-import Online.Averages-import Online.AveragesB-import Online.Medians

@@ -1,169 +0,0 @@-{-# LANGUAGE CPP #-}---- | online statistics based on a moving average-module Online.Averages-  ( Averager-  , online-    -- * online statistics-  , ma-  , absma-  , sqma-  , std-  , cov-  , corr-  , corrGauss-  , beta-  , alpha-  , autocorr-  , mconst-  ) where--import qualified Control.Foldl as L-import Control.Foldl (Fold(..))-import Prelude---- | Most common statistics are averages.-newtype Averager a b = Averager-  { _averager :: (a, b)-  }--instance (Semigroup a, Semigroup b) => Semigroup (Averager a b) where-  (Averager (s, c)) <> (Averager (s', c')) =-    Averager (s <> s', c <> c')--instance (Semigroup a, Semigroup b, Monoid a, Monoid b) => Monoid (Averager a b) where-  mempty = Averager (mempty, mempty)-  mappend = (<>)---- | online takes a function and turns it into a `Control.Foldl.Fold` where the step is an incremental update of the (isomorphic) statistic.-online :: (Fractional b) => (a -> b) -> (b -> b) -> Fold a b-online f g = Fold step begin extract-  where-    begin = Averager (0, 0)-    step (Averager (s, c)) a = Averager (g $ s + f a, g $ c + 1)-    extract (Averager (s, c)) = s / c-{-# INLINABLE online #-}---- $setup------ >>> import qualified Control.Foldl as L--- >>> let n = 100--- >>> let r = 0.9---- | moving average with a decay rate--- --- so 'ma 1' is the simple average (no decay in the statistic), and 'ma 0.00001' is the last value (insta-decay)------ >>> L.fold (ma 1) [0..100]--- 50.0------ >>> L.fold (ma 1e-12) [0..100]--- 99.999999999999------ >>> L.fold (ma 0.9) [0..100]--- 91.00241448887785----ma :: (Fractional a) => a -> Fold a a-ma r = online id (* r)-{-# INLINABLE ma #-}---- | absolute average-absma :: (Fractional a) => a -> Fold a a-absma r = online abs (* r)-{-# INLINABLE absma #-}---- | average square-sqma :: (Fractional a) => a -> Fold a a-sqma r = online (\x -> x * x) (* r)-{-# INLINABLE sqma #-}---- | standard deviation------ The formulae for standard deviation, expressed in online terminology, highlights how this statistic is composed of averages:------ > (\s ss -> sqrt (ss - s ** (one+one))) <$> ma r <*> sqma r------ The average deviation of the numbers 1..1000 is about 1 / sqrt 12 * 1000 (see <<https://en.wikipedia.org/wiki/Uniform_distribution_(continuous)#Standard_uniform wiki>>)------ >>> L.fold (std 1) [0..1000]--- 288.9636655359978------ The average deviation with a decay of 0.99------ >>> L.fold (std 0.99) [0..1000]--- 99.28328803164005-std :: (Fractional a, Floating a) => a -> Fold a a-std r = (\s ss -> sqrt (ss - s ** 2)) <$> ma r <*> sqma r-{-# INLINABLE std #-}---- | the covariance of a tuple--- given an underlying central tendency fold-cov :: (Num a) => Fold a a -> Fold (a, a) a-cov m =-  (\xy x' y' -> xy - x' * y') <$> L.premap (uncurry (*)) m <*> L.premap fst m <*>-  L.premap snd m-{-# INLINABLE cov #-}---- | correlation of a tuple, specialised to Guassian-corrGauss :: (Floating a) => a -> Fold (a, a) a-corrGauss r =-  (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> cov (ma r) <*>-  L.premap fst (std r) <*>-  L.premap snd (std r)-{-# INLINABLE corrGauss #-}---- | a generalised version of correlation of a tuple-corr :: (Floating a) => Fold a a -> Fold a a -> Fold (a, a) a-corr central deviation =-  (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> cov central <*>-  L.premap fst deviation <*>-  L.premap snd deviation-{-# INLINABLE corr #-}---- | the beta in a simple linear regression of a tuple--- given an underlying central tendency fold-beta :: (Fractional a) => Fold a a -> Fold (a, a) a-beta m =-  (\xy x' y' x2 -> (xy - x' * y') / (x2 - x' * x')) <$> L.premap (uncurry (*)) m <*>-  L.premap fst m <*>-  L.premap snd m <*>-  L.premap (\(x, _) -> x * x) m-{-# INLINABLE beta #-}---- | the alpha of a tuple-alpha :: (Fractional a) => Fold a a -> Fold (a, a) a-alpha m = (\y b x -> y - b * x) <$> L.premap fst m <*> beta m <*> L.premap snd m-{-# INLINABLE alpha #-}--{-| autocorrelation is a slippery concept.  This method starts with the concept that there is an underlying random error process (e), and autocorrelation is a process on top of that ie for a one-step correlation relationship.--value@t = e@t + k * e@t-1--where k is the autocorrelation.--There are thus two online rates needed: one for the average being considered to be the dependent variable, and one for the online of the correlation calculation between the most recent value and the moving average. For example,--> L.fold (autocorr zero one)--would estimate the one-step autocorrelation relationship of the previous value and the current value over the entire sample set.---}-autocorr :: (RealFloat a) => Fold a a -> Fold (a, a) a -> Fold a a-autocorr central corrf =-  case central of-    (Fold mStep mBegin mDone) ->-      case corrf of-        (Fold dStep dBegin dDone) ->-          let begin = (mBegin, dBegin)-              step (mAcc, dAcc) a =-                ( mStep mAcc a-                , if isNaN (mDone mAcc)-                    then dAcc-                    else dStep dAcc (mDone mAcc, a))-              done = dDone . snd-          in Fold step begin done-{-# INLINABLE autocorr #-}---- | a constant fold-mconst :: a -> L.Fold a a-mconst a = L.Fold (\() _ -> ()) () (const a)

@@ -1,64 +0,0 @@-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE CPP #-}-{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}---- | online statistics based on a moving average-module Online.AveragesB-  ( foldB,-    foldB',-    maB,-    absmaB,-    sqmaB,-    stdB,-    std'-  ) where--import Prelude-import Numeric.Backprop as B-import qualified Prelude.Backprop as PB-import Control.Foldl (Fold(..))--foldB' :: (Backprop a, Backprop b, Reifies s W, Fractional b) => (BVar s a -> BVar s b) -> BVar s b -> BVar s [a] -> BVar s b-foldB' f r xs = divide (PB.foldl' (step' f r) (T2 0 0) xs) where-  step' f' r' (T2 s c) a = uncurry T2 ((r'*) $ s + f' a, (r'*) $ c + 1)  -  divide (T2 s c) = s / c---online :: (Reifies s W, Fractional b) => (BVar s a -> BVar s b) -> (BVar s b -> BVar s b) -> Fold (BVar s a) (BVar s b)-online f g = Fold step begin extract-  where-    begin = (0, 0)-    step (s, c) a = (g $ s + f a, g $ c + 1)-    extract (s, c) = s / c--ma' :: (Reifies s W, Fractional b) => BVar s b -> Fold (BVar s b) (BVar s b)-ma' r = online id (*r)--sqma' :: (Reifies s W, Fractional b) => BVar s b -> Fold (BVar s b) (BVar s b)-sqma' r = online (\x -> x * x) (*r)--std' :: (Reifies s W, Floating b) => BVar s b -> Fold (BVar s b) (BVar s b)-std' r = (\s ss -> sqrt (ss - s ** 2)) <$> ma' r <*> sqma' r---- coerceVar--foldB :: (Reifies s W) => (BVar s Double -> BVar s Double) -> BVar s Double -> BVar s [Double] -> BVar s Double-foldB f r xs = divide (PB.foldl' (step' f r) (T2 0 0) xs) where-  step' f' r' (T2 s c) a = uncurry T2 ((r*) $ s + f' a, (r'*) $ c + 1)  -  divide (T2 s c) = s / c--stdB :: Reifies s W => BVar s Double -> BVar s [Double] -> BVar s Double-stdB r xs = (\ss s -> sqrt (ss - s ** 2)) (foldB id r xs) (foldB (\x -> x * x) r xs)---- stdB' :: Reifies s W => BVar s Double -> BVar s [Double] -> BVar s Double--- stdB' r xs = (\(T2 ss s) -> sqrt (ss - s ** 2)) (foldB' (\x -> (T2 x (x*x))) (T2 r r) xs)--maB :: Reifies s W => BVar s Double -> BVar s [Double] -> BVar s Double-maB r = foldB id r--absmaB :: Reifies s W => BVar s Double -> BVar s [Double] -> BVar s Double-absmaB r = foldB abs r--sqmaB :: Reifies s W => BVar s Double -> BVar s [Double] -> BVar s Double-sqmaB r = foldB (\x -> x * x) r-

@@ -1,125 +0,0 @@-module Online.Medians-  ( -- * convert a statistic to an online median stat equivalent to L1-    Medianer(..)-  , onlineL1-  , onlineL1'--    -- * online statistics-  , maL1-  , absmaL1-  , covL1-  , corrL1-  , betaL1-  , alphaL1-  , autocorrL1-  ) where--import qualified Control.Foldl as L-import Control.Foldl (Fold(..))-import Prelude---- | A rough Median.--- The average absolute value of the stat is used to callibrate estimate drift towards the median-data Medianer a b = Medianer-  { medAbsSum :: a-  , medCount :: b-  , medianEst :: a-  }---- | onlineL1' takes a function and turns it into a `Control.Foldl.Fold` where the step is an incremental update of an (isomorphic) median statistic.-onlineL1' ::-     (Ord b, Fractional b) => b -> b -> (a -> b) -> (b -> b) -> Fold a (b, b)-onlineL1' i d f g = Fold step begin extract-  where-    begin = Medianer 0 0 0-    step (Medianer s c m) a =-      Medianer-        (g $ s + abs (f a))-        (g $ c + 1)-        ((1 - d) * (m + s' * i * s / c') + d * f a)-      where-        c' =-          if c == 0-            then 1-            else c-        s'-          | f a > m = 1-          | f a < m = -1-          | otherwise = 0-    extract (Medianer s c m) = (s / c, m)-{-# INLINABLE onlineL1' #-}---- | onlineL1 takes a function and turns it into a `Control.Foldl.Fold` where the step is an incremental update of an (isomorphic) median statistic.-onlineL1 :: (Ord b, Fractional b) => b -> b -> (a -> b) -> (b -> b) -> Fold a b-onlineL1 i d f g = snd <$> onlineL1' i d f g-{-# INLINABLE onlineL1 #-}---- $setup------ >>> import qualified Control.Foldl as L--- >>> let n = 100--- >>> let inc = 0.1--- >>> let d = 0--- >>> let r = 0.9---- | moving median--- >>> L.fold (maL1 inc d r) [1..n]--- 93.92822312742108----maL1 :: (Ord a, Fractional a) => a -> a -> a -> Fold a a-maL1 i d r = onlineL1 i d id (* r)-{-# INLINABLE maL1 #-}---- | moving absolute deviation-absmaL1 :: (Ord a, Fractional a) => a -> a -> a -> Fold a a-absmaL1 i d r = fst <$> onlineL1' i d id (* r)-{-# INLINABLE absmaL1 #-}---- | covariance of a tuple-covL1 :: (Ord a, Fractional a) => a -> a -> a -> Fold (a, a) a-covL1 i d r =-  (\xy xbar ybar -> xy - xbar * ybar) <$> onlineL1 i d (uncurry (*)) (* r) <*>-  onlineL1 i d fst (* r) <*>-  onlineL1 i d snd (* r)-{-# INLINABLE covL1 #-}---- | correlation of a tuple-corrL1 :: (Ord a, Floating a) => a -> a -> a -> Fold (a, a) a-corrL1 i d r =-  (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> covL1 i d r <*>-  L.premap fst (absmaL1 i d r) <*>-  L.premap snd (absmaL1 i d r)-{-# INLINABLE corrL1 #-}---- | the beta in a simple linear regression of a tuple-betaL1 :: (Ord a, Floating a) => a -> a -> a -> Fold (a, a) a-betaL1 i d r =-  (\xy x' y' x2 -> (xy - x' * y') / (x2 - x' * x')) <$>-  L.premap (uncurry (*)) (maL1 i d r) <*>-  L.premap fst (maL1 i d r) <*>-  L.premap snd (maL1 i d r) <*>-  L.premap (\(x, _) -> x * x) (maL1 i d r)-{-# INLINABLE betaL1 #-}---- | the alpha in a simple linear regression of `snd` on `fst`-alphaL1 :: (Ord a, Floating a) => a -> a -> a -> Fold (a, a) a-alphaL1 i d r =-  (\y b x -> y - b * x) <$> L.premap fst (maL1 i d r) <*> betaL1 i d r <*>-  L.premap snd (maL1 i d r)-{-# INLINABLE alphaL1 #-}--autocorrL1 :: (Floating a, RealFloat a) => a -> a -> a -> a -> Fold a a-autocorrL1 i d maR corrR =-  case maL1 i d maR of-    (Fold maStep maBegin maDone) ->-      case corrL1 i d corrR of-        (Fold corrStep corrBegin corrDone) ->-          let begin = (maBegin, corrBegin)-              step (maAcc, corrAcc) a =-                ( maStep maAcc a-                , if isNaN (maDone maAcc)-                    then corrAcc-                    else corrStep corrAcc (maDone maAcc, a))-              done = corrDone . snd-          in Fold step begin done-{-# INLINABLE autocorrL1 #-}

@@ -1,141 +0,0 @@-{-# LANGUAGE DataKinds #-}--module Online.Quantiles-  ( tDigest-  , tDigestQuantiles-  , tDigestHist-  , OnlineTDigest(..)-  , onlineQuantiles-  , Online.Quantiles.median-  , onlineDigitize-  , onlineDigestHist-  )-where--import qualified Control.Foldl as L-import Data.List.NonEmpty (NonEmpty)-import Data.TDigest-import Data.TDigest.Internal-import Data.TDigest.Tree.Internal (TDigest(..), size, emptyTDigest, insertCentroid, relMaxSize, absMaxSize, toMVector)-import Data.TDigest.Postprocess (HistBin, histogram)-import qualified Data.Vector.Algorithms.Heap as VHeap-import qualified Data.Vector.Unboxed as VU-import Prelude-import Control.Monad.ST (runST)-import Data.Maybe-import Data.Ord-import Data.Foldable---- | a raw non-online tdigest fold-tDigest :: L.Fold Double (TDigest 25)-tDigest = L.Fold step begin done-  where-    step x a = insert a x-    begin = tdigest ([] :: [Double]) :: TDigest 25-    done = id---- | non-online version-tDigestQuantiles :: [Double] -> L.Fold Double [Double]-tDigestQuantiles qs = L.Fold step begin done-  where-    step x a = insert a x-    begin = tdigest ([] :: [Double]) :: TDigest 25-    done x = fromMaybe (0/0) . (`quantile` compress x) <$> qs---- | non-online version-tDigestHist :: L.Fold Double (Maybe (NonEmpty HistBin))-tDigestHist = L.Fold step begin done-  where-    step x a = insert a x-    begin = tdigest ([] :: [Double]) :: TDigest 25-    done = histogram . compress--data OnlineTDigest = OnlineTDigest-  { td :: TDigest 25-  , tdN :: Int-  , tdRate :: Double-  } deriving (Show)--emptyOnlineTDigest :: Double -> OnlineTDigest-emptyOnlineTDigest = OnlineTDigest (emptyTDigest :: TDigest n) 0---- | decaying quantiles based on the tdigest library-onlineQuantiles :: Double -> [Double] -> L.Fold Double [Double]-onlineQuantiles r qs = L.Fold step begin done-  where-    step x a = onlineInsert a x-    begin = emptyOnlineTDigest r-    done x = fromMaybe (0/0) . (`quantile` t) <$> qs-      where-        (OnlineTDigest t _ _) = onlineForceCompress x--median :: Double -> L.Fold Double Double-median r = L.Fold step begin done-  where-    step x a = onlineInsert a x-    begin = emptyOnlineTDigest r-    done x = fromMaybe (0/0) (quantile 0.5 t)-      where-        (OnlineTDigest t _ _) = onlineForceCompress x--onlineInsert' :: Double -> OnlineTDigest -> OnlineTDigest-onlineInsert' x (OnlineTDigest td' n r) =-  OnlineTDigest-    (insertCentroid (x, r ^^ (-(fromIntegral $ n + 1) :: Integer)) td')-    (n + 1)-    r--onlineInsert :: Double -> OnlineTDigest -> OnlineTDigest-onlineInsert x otd = onlineCompress (onlineInsert' x otd)--onlineCompress :: OnlineTDigest -> OnlineTDigest-onlineCompress otd@(OnlineTDigest Nil _ _) = otd-onlineCompress otd@(OnlineTDigest t _ _)-  | Data.TDigest.Tree.Internal.size t > relMaxSize * compression &&-      Data.TDigest.Tree.Internal.size t > absMaxSize = onlineForceCompress otd-  | otherwise = otd-  where-    compression = 25--onlineForceCompress :: OnlineTDigest -> OnlineTDigest-onlineForceCompress otd@(OnlineTDigest Nil _ _) = otd-onlineForceCompress (OnlineTDigest t n r) = OnlineTDigest t' 0 r-  where-    t' =-      foldl' (flip insertCentroid) emptyTDigest $-      (\(m, w) -> (m, w * (r ^^ n))) . fst <$> VU.toList centroids-    -- Centroids are shuffled based on space-    centroids :: VU.Vector (Centroid, Double)-    centroids =-      runST $ do-        v <- toMVector t-        -- sort by cumulative weight-        VHeap.sortBy (comparing snd) v-        VU.unsafeFreeze v--onlineDigitize :: Double -> [Double] -> L.Fold Double Int-onlineDigitize r qs = L.Fold step begin done-  where-    step (x, _) a = (onlineInsert a x, a)-    begin = (emptyOnlineTDigest r, 0/0)-    done (x, l) = bucket' qs' l-      where-        qs' = fromMaybe (0/0) . (`quantile` t) <$> qs-        (OnlineTDigest t _ _) = onlineForceCompress x-        bucket' xs l' =-          L.fold L.sum $-          (\x' ->-             if x' > l'-               then 0-               else 1) <$>-          xs---- | decaying histogram based on the tdigest library-onlineDigestHist :: Double -> L.Fold Double (Maybe (NonEmpty HistBin))-onlineDigestHist r = L.Fold step begin done-  where-    step x a = onlineInsert a x-    begin = emptyOnlineTDigest r-    done x = histogram . compress $ t-      where-        (OnlineTDigest t _ _) = onlineForceCompress x

@@ -1,19 +1,16 @@ {-# OPTIONS_GHC -Wall #-}+{-# OPTIONS_GHC -fno-warn-unused-imports #-}  module Main where -import Prelude-import Test.Tasty (TestTree, testGroup, defaultMain)+import NumHask.Prelude import Test.DocTest+import Data.Mealy  main :: IO ()-main = do-    doctest ["src/Online/Averages.hs", "src/Online/Medians.hs"]-    defaultMain tests--tests :: TestTree-tests =-    testGroup ""-    [-    ]-+main =+  doctest+  [ "src/Data/Mealy.hs",+    "src/Data/Quantiles.hs",+    "src/Data/Simulate.hs"+  ]