diff --git a/ChangeLog.md b/ChangeLog.md new file mode 100644 --- /dev/null +++ b/ChangeLog.md @@ -0,0 +1,10 @@ +0.5 +== + +- Switched to harpie, away from numhask-array + + +0.4.4 +=== + +- Added some common patterns diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) 2013, Tony Day + +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. diff --git a/mealy.cabal b/mealy.cabal --- a/mealy.cabal +++ b/mealy.cabal @@ -1,74 +1,71 @@ -cabal-version: 2.4 -name: mealy -version: 0.0.3 -synopsis: See readme.md - -description: See readme.md for description. -category: project +cabal-version: 3.0 +name: mealy +version: 0.5.1.1 +license: BSD-3-Clause +license-file: LICENSE +copyright: Tony Day (c) 2013 +category: algorithm author: Tony Day maintainer: tonyday567@gmail.com -copyright: Tony Day (c) AfterTimes -license: BSD-3-Clause homepage: https://github.com/tonyday567/mealy#readme bug-reports: https://github.com/tonyday567/mealy/issues +synopsis: Mealy machines for processing time-series and ordered data. +description: + @mealy@ reimagines statistics as a [mealy machine](https://en.wikipedia.org/wiki/Mealy_machine) processing data with some form of order such as time-series data. The 'Mealy', with the help of a decay function specifying the relative weights of recent values versus older value, can be treated as a compression or summary of the data stream into 'current state.' + Mealies are highly polymorphic, situated at a busy crossroad of theory and practice, and lend themselves to ergonmic, compact and realistic representations of a wide range of online phenomena. + + == Usage + + >>> import Mealy + + >>> fold ((,) <$> ma 0.9 <*> std 0.9) [1..100] + (91.00265621044142,9.472822805289121) + build-type: Simple +tested-with: + ghc ==9.10.3 + ghc ==9.12.2 + ghc ==9.14.1 + +extra-doc-files: + ChangeLog.md + readme.md + source-repository head type: git location: https://github.com/tonyday567/mealy -library - hs-source-dirs: - src - build-depends: - adjunctions >= 4.4, - base >=4.7 && <5, - containers >= 0.6, - folds, - generic-lens >= 2.0, - lens, - matrix >= 0.3.6 && < 0.3.7, - mwc-probability >= 2.3.1 && < 2.4, - numhask >= 0.7.1 && < 0.8, - numhask-array >= 0.8 && < 0.9, - primitive >= 0.7, - profunctors >= 5.5, - tdigest, - text, - vector, - vector-algorithms - exposed-modules: - Data.Mealy - Data.Mealy.Quantiles - Data.Mealy.Simulate - other-modules: - default-language: Haskell2010 - default-extensions: +common ghc-options-stanza ghc-options: -Wall -Wcompat + -Widentities -Wincomplete-record-updates -Wincomplete-uni-patterns + -Wpartial-fields -Wredundant-constraints - -fwrite-ide-info - -hiedir=.hie -test-suite test - type: exitcode-stdio-1.0 - main-is: test.hs - hs-source-dirs: - test +library + import: ghc-options-stanza + default-language: GHC2024 + hs-source-dirs: src build-depends: - base >=4.7 && <5, - doctest, - numhask >= 0.7 && < 0.8, - mealy - default-language: Haskell2010 - default-extensions: - ghc-options: - -Wall - -Wcompat - -Wincomplete-record-updates - -Wincomplete-uni-patterns - -Wredundant-constraints - -fwrite-ide-info - -hiedir=.hie + adjunctions >=4.0 && <4.5, + base >=4.14 && <5, + containers >=0.6 && <0.9, + harpie >=0.1 && <0.3, + harpie-numhask >=0.1 && <0.3, + mwc-probability >=2.3.1 && <2.4, + numhask >=0.11 && <0.14, + primitive >=0.7.2 && <0.10, + profunctors >=5.6.2 && <5.7, + tdigest >=0.2.1 && <0.4, + text >=1.2 && <2.2, + vector >=0.12.3 && <0.14, + vector-algorithms >=0.8.0 && <0.10, + + exposed-modules: + Data.Mealy + Data.Mealy.Quantiles + Data.Mealy.Simulate + diff --git a/readme.md b/readme.md new file mode 100644 --- /dev/null +++ b/readme.md @@ -0,0 +1,34 @@ +# mealy + +[![Hackage](https://img.shields.io/hackage/v/mealy.svg)](https://hackage.haskell.org/package/mealy) [![Build Status](https://github.com/tonyday567/mealy/actions/workflows/haskell-ci.yml/badge.svg)](https://github.com/tonyday567/mealy/actions/workflows/haskell-ci.yml) + +A 'Mealy' is a triple of functions + +- (a -> b) **inject**: Convert (initial) input into the (initial) state type. +- (b -> a -> b) **step**: Update state given prior state and (new) input. +- (c -> b) **extract**: Convert state to the output type. + +A sum, for example, looks like `M id (+) id` where the first id is the initial injection and the second id is the covariant extraction. + +This library provides support for computing statistics (such as an average or a standard deviation) as current state within a mealy context. + +## Usage + +Usage is to supply a decay function representing the relative weights of recent values versus older ones, in the manner of exponentially-weighted averages. The library attempts to be polymorphic in the statistic which can be combined in applicative style. + +```haskell +import Prelude +import Data.Mealy +``` + +```haskell +fold ((,) <$> ma 0.9 <*> std 0.9) [1..100::Double] +``` + +``` +(91.00265621044142,9.472822805289121) +``` + +## Reference + +[Finite State Transducers](https://stackoverflow.com/questions/27997155/finite-state-transducers-in-haskell) diff --git a/src/Data/Mealy.hs b/src/Data/Mealy.hs --- a/src/Data/Mealy.hs +++ b/src/Data/Mealy.hs @@ -1,32 +1,14 @@ -{-# 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 StrictData #-} -{-# LANGUAGE TupleSections #-} +{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE NoImplicitPrelude #-} -{-# OPTIONS_GHC -Wall #-} -{-# OPTIONS_GHC -Wno-incomplete-patterns #-} -{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-} -{-# OPTIONS_GHC -Wno-name-shadowing #-} -{-# OPTIONS_GHC -Wno-orphans #-} -{-# OPTIONS_GHC -Wno-type-defaults #-} --- | Online statistics for ordered data (such as time-series data), modelled as [mealy machines](https://en.wikipedia.org/wiki/Mealy_machine) +-- | A 'Mealy' is a polymorphic statistic reimagined as a [mealy machine](https://en.wikipedia.org/wiki/Mealy_machine) processing and summarising data with some form of order (such as time-series data), where recent data is less relevant than old data. module Data.Mealy ( -- * Types Mealy (..), + dipure, pattern M, scan, fold, @@ -37,7 +19,7 @@ online, -- * Statistics - -- $setup + -- $example-set ma, absma, sqma, @@ -53,84 +35,148 @@ reg, asum, aconst, + last, + maybeLast, delay1, delay, - depState, - Model1 (..), - zeroModel1, - depModel1, + window, + diff, + gdiff, + same, + countM, + sumM, + listify, -- * median Medianer (..), onlineL1, - onlineL1', maL1, - absmaL1, ) where -import Control.Lens hiding (Empty, Unwrapped, Wrapped, index, (:>), (|>)) -import Data.Fold hiding (M) +import Control.Category +import Control.Exception +import Data.Bifunctor import Data.Functor.Rep -import Data.Generics.Labels () - --- import qualified Numeric.LinearAlgebra as LA - -import qualified Data.Matrix as M -import qualified Data.Sequence as Seq -import qualified NumHask.Array.Fixed as F -import NumHask.Array.Shape (HasShape) -import NumHask.Prelude hiding (L1, State, StateT, asum, fold, get, replace, runState, runStateT, state) +import Data.List (scanl') +import Data.Map qualified as Map +import Data.Profunctor +import Data.Sequence (Seq) +import Data.Sequence qualified as Seq +import Data.Text (Text) +import GHC.TypeLits +import Harpie.Fixed qualified as F +import Harpie.NumHask qualified as N +import Harpie.Shape qualified as S +import NumHask.Prelude hiding (asum, diff, fold, id, last, (.)) -- $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.Mealy.Simulate +-- >>> import Harpie.Fixed qualified as F -- >>> 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 +-- $example-set +-- The doctest examples are composed from some random series generated with Data.Mealy.Simulate. -- --- * (a -> b) __inject__ Convert an input into the state type. --- * (b -> a -> b) __step__ Update state given prior state and (new) input. --- * (c -> b) __extract__ Convert state to the output type. +-- - xs0, xs1 & xs2 are samples from N(0,1) -- --- By adopting this order, a Mealy sum looks like: +-- - xsp is a pair of N(0,1)s with a correlation of 0.8 -- +-- >>> :set -XDataKinds +-- >>> import Data.Mealy.Simulate +-- >>> g <- create +-- >>> xs0 <- rvs g 10000 +-- >>> xs1 <- rvs g 10000 +-- >>> xs2 <- rvs g 10000 +-- >>> xsp <- rvsp g 10000 0.8 + +newtype MealyError = MealyError {mealyErrorMessage :: Text} + deriving (Show) + +instance Exception MealyError + +-- | A 'Mealy' a b is a triple of functions +-- +-- * (a -> s) __inject__ Convert an input into an internal state type. +-- * (s -> a -> s) __step__ Update state given prior state and (new) input. +-- * (s -> b) __extract__ Convert state to the output type. +-- +-- By adopting this order, a sum, for example, looks like: +-- -- > M id (+) id -- --- where the first id is the initial injection to a contravariant position, and the second id is the covriant extraction. +-- where the first id is the initial injection to a contravariant position, and the second id is the covariant extraction. -- --- __inject__ kicks off state on the initial element of the Foldable, but is otherwise be independent of __step__. +-- __inject__ kicks off state on the initial element of the Foldable, but is otherwise independent of __step__. -- --- > 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 +-- > scan (M i s e) (x : xs) = e <$> scanl' s (i x) xs +data Mealy a b = forall s. Mealy (a -> s) (s -> a -> s) (s -> b) --- | Pattern for a 'Mealy'. +-- | Strict Pair +data Pair' a b = Pair' !a !b deriving (Eq, Ord, Show, Read) + +instance (Semigroup a, Semigroup b) => Semigroup (Pair' a b) where + Pair' a b <> Pair' c d = Pair' (a <> c) (b <> d) + {-# INLINE (<>) #-} + +instance (Monoid a, Monoid b) => Monoid (Pair' a b) where + mempty = Pair' mempty mempty + +instance Functor (Mealy a) where + fmap f (Mealy z h k) = Mealy z h (f . k) + +instance Applicative (Mealy a) where + pure x = Mealy (const ()) (\() _ -> ()) (\() -> x) + Mealy zf hf kf <*> Mealy za ha ka = + Mealy + (\a -> Pair' (zf a) (za a)) + (\(Pair' x y) a -> Pair' (hf x a) (ha y a)) + (\(Pair' x y) -> kf x (ka y)) + +instance Category Mealy where + id = Mealy id (\_ a -> a) id + Mealy z h k . Mealy z' h' k' = Mealy z'' h'' (\(Pair' b _) -> k b) + where + z'' a = Pair' (z (k' b)) b where b = z' a + h'' (Pair' c d) a = Pair' (h c (k' d')) d' where d' = h' d a + +instance Profunctor Mealy where + dimap f g (Mealy z h k) = Mealy (z . f) (\a -> h a . f) (g . k) + lmap f (Mealy z h k) = Mealy (z . f) (\a -> h a . f) k + rmap g (Mealy z h k) = Mealy z h (g . k) + +instance Strong Mealy where + first' (M i s e) = M (first i) (\(cl, _) (al, ar) -> (s cl al, ar)) (first e) + +-- The right type for Choice would be something like: -- +-- left' :: p a b -> p (Either a c) (These b c) + +-- | Convenience pattern for a 'Mealy'. +-- -- @M extract step inject@ pattern M :: (a -> c) -> (c -> a -> c) -> (c -> b) -> Mealy a b -pattern M i s e = Mealy (L1 e s i) +pattern M i s e = Mealy i s e {-# COMPLETE M #-} +-- | Create a Mealy from a (pure) function +dipure :: (a -> a -> a) -> Mealy a a +dipure f = M id f id + -- | Fold a list through a 'Mealy'. -- -- > cosieve == fold fold :: Mealy a b -> [a] -> b -fold _ [] = panic "on the streets of Birmingham." +fold _ [] = throw (MealyError "empty list") fold (M i s e) (x : xs) = e $ foldl' s (i x) xs -- | Run a list through a 'Mealy' and return a list of values for every step @@ -175,9 +221,25 @@ 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 f g@ is a 'Mealy' where f is a transformation of the data and +-- g is a decay function (usually convergent to zero) applied at each step. -- -- > online id id == av +-- +-- @online@ is best understood by examining usage +-- to produce a moving average and standard deviation: +-- +-- An exponentially-weighted moving average with a decay rate of 0.9 +-- +-- > ma r == online id (*r) +-- +-- An exponentially-weighted moving average of the square. +-- +-- > sqma r = online (\x -> x * x) (* r) +-- +-- Applicative-style exponentially-weighted standard deviation computation: +-- +-- > std r = (\s ss -> sqrt (ss - s ** 2)) <$> ma r <*> sqma r online :: (Divisive b, Additive b) => (a -> b) -> (b -> b) -> Mealy a b online f g = M intract step av where @@ -188,18 +250,14 @@ -- | 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]) +-- >>> fold (ma 0) ([1..100]) -- 100.0 -- --- >>> fold (ma 1) (fromList [1..100]) +-- >>> fold (ma 1) ([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) +-- 9.713356299018187e-2 ma :: (Divisive a, Additive a) => a -> Mealy a a ma r = online id (* r) {-# INLINEABLE ma #-} @@ -207,8 +265,8 @@ -- | absolute average -- -- >>> fold (absma 1) xs0 --- 0.7894201075535578 -absma :: (Divisive a, Signed a) => a -> Mealy a a +-- 0.8075705557429647 +absma :: (Divisive a, Absolute a) => a -> Mealy a a absma r = online abs (* r) {-# INLINEABLE absma #-} @@ -237,17 +295,15 @@ -- 99.28328803163829 -- -- >>> fold (std 1) xs0 --- 0.9923523681261158 --- --- ![std chart](other/ex-std.svg) -std :: (Divisive a, ExpField a) => a -> Mealy a a +-- 1.0126438036262801 +std :: (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 +-- 0.7818936662586868 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 @@ -256,10 +312,11 @@ -- | correlation of a tuple, specialised to Guassian -- -- >>> fold (corrGauss 1) xsp --- 0.8020637696465039 +-- 0.7978347126677433 corrGauss :: (ExpField a) => a -> Mealy (a, a) a corrGauss r = - (\cov' stdx stdy -> cov' / (stdx * stdy)) <$> cov (ma r) + (\cov' stdx stdy -> cov' / (stdx * stdy)) + <$> cov (ma r) <*> lmap fst (std r) <*> lmap snd (std r) {-# INLINEABLE corrGauss #-} @@ -267,12 +324,13 @@ -- | a generalised version of correlation of a tuple -- -- >>> fold (corr (ma 1) (std 1)) xsp --- 0.8020637696465039 +-- 0.7978347126677433 -- -- > 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 + (\cov' stdx stdy -> cov' / (stdx * stdy)) + <$> cov central <*> lmap fst deviation <*> lmap snd deviation {-# INLINEABLE corr #-} @@ -290,10 +348,11 @@ -- \] -- -- >>> fold (beta1 (ma 1)) $ zipWith (\x y -> (y, x + y)) xs0 xs1 --- 0.9953875263096014 +-- 0.999747321294513 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 + (\xy x' y' x2 -> (xy - x' * y') / (x2 - x' * x')) + <$> lmap (uncurry (*)) m <*> lmap fst m <*> lmap snd m <*> lmap (\(x, _) -> x * x) m @@ -310,7 +369,7 @@ -- \] -- -- >>> fold (alpha1 (ma 1)) $ zipWith (\x y -> ((3+y), x + 0.5 * (3 + y))) xs0 xs1 --- 1.1880996822796197e-2 +-- 1.3680496627365146e-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 #-} @@ -318,7 +377,7 @@ -- | 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) +-- (1.3680496627365146e-2,0.4997473212944953) reg1 :: (ExpField a) => Mealy a a -> Mealy (a, a) (a, a) reg1 m = (,) <$> alpha1 m <*> beta1 m @@ -344,41 +403,24 @@ -- \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, KnownNat n, Fractional a, Eq a) => a -> Mealy (F.Array '[n] a, a) (F.Array '[n] a) +-- >>> let ys = zipWith3 (\x y z -> 0.1 * x + 0.5 * y + 1 * z) xs0 xs1 xs2 +-- >>> let zs = zip (zipWith (\x y -> F.array @'[2] [x,y]) xs1 xs2) ys +-- >>> fold (beta 0.99) zs +-- [0.6228820021456606,0.8461936860075405] +beta :: (ExpField a, KnownNat n, Eq a) => a -> Mealy (F.Array '[n] a, a) (F.Array '[n] a) beta r = M inject step extract where -- extract :: Averager (RegressionState n a) a -> (F.Array '[n] a) extract (A (RegressionState xx x xy y) c) = - (\a b -> inverse a `F.mult` b) - ((one / c) .* (xx - F.expand (*) x x)) - ((xy - (y .* x)) *. (one / c)) + (\a b -> recip a `N.mult` b) + ((one / c) *| (xx - F.expand (*) x x)) + ((xy - (y *| x)) |* (one / c)) step x (xs, y) = rsOnline r x (inject (xs, y)) -- inject :: (F.Array '[n] a, a) -> Averager (RegressionState n a) a inject (xs, y) = - A (RegressionState (F.expand (*) xs xs) xs (y .* xs) y) one + A (RegressionState (F.expand (*) xs xs) xs (y *| xs) y) one {-# INLINEABLE beta #-} -toMatrix :: (KnownNat n, KnownNat m) => F.Array [m, n] a -> M.Matrix a -toMatrix a = M.matrix m n (index a . (\(i, j) -> [i, j])) - where - (m : n : _) = F.shape a - -fromMatrix :: (KnownNat n, KnownNat m) => M.Matrix a -> F.Array [m, n] a -fromMatrix = fromList . M.toList - -data MatrixException = MatrixException - deriving (Show) - -instance Exception MatrixException - --- | The inverse of a square matrix. -inverse :: (KnownNat n, Fractional a, Eq a) => F.Array [n, n] a -> F.Array [n, n] a -inverse = either (const $ throw MatrixException) fromMatrix . M.inverse . toMatrix - 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') @@ -386,11 +428,11 @@ d s s' = r * s + s' -- | alpha in a multiple regression -alpha :: (ExpField a, KnownNat n, Fractional a, Eq a) => a -> Mealy (F.Array '[n] a, a) a +alpha :: (ExpField a, KnownNat n, Eq a) => 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 :: (S.KnownNats 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 @@ -399,11 +441,11 @@ -- | 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 :: (ExpField a, KnownNat n, Fractional a, Eq a) => a -> Mealy (F.Array '[n] a, a) (F.Array '[n] a, a) +-- >>> let ys = zipWith3 (\x y z -> 0.1 * x + 0.5 * y + 1 * z) xs0 xs1 xs2 +-- >>> let zs = zip (zipWith (\x y -> F.array @'[2] [x,y]) xs1 xs2) ys +-- >>> fold (reg 0.99) zs +-- ([0.6228820021456606,0.8461936860075405],2.536775201287266e-2) +reg :: (ExpField a, KnownNat n, Eq a) => a -> Mealy (F.Array '[n] a, a) (F.Array '[n] a, a) reg r = (,) <$> beta r <*> alpha r {-# INLINEABLE reg #-} @@ -415,6 +457,14 @@ aconst :: b -> Mealy a b aconst b = M (const ()) (\_ _ -> ()) (const b) +-- | most recent value +last :: Mealy a a +last = M id (\_ a -> a) id + +-- | most recent value if it exists, previous value otherwise. +maybeLast :: a -> Mealy (Maybe a) a +maybeLast def = M (fromMaybe def) fromMaybe id + -- | delay input values by 1 delay1 :: a -> Mealy a a delay1 x0 = M (x0,) (\(_, x) a -> (x, a)) fst @@ -441,102 +491,44 @@ where inject a = Seq.fromList x0 Seq.|> a extract :: Seq a -> a - extract Seq.Empty = panic "ACAB" + extract Seq.Empty = throw (MealyError "empty seq") extract (x Seq.:<| _) = x step :: Seq a -> a -> Seq a - step Seq.Empty _ = panic "ACAB" + step Seq.Empty _ = throw (MealyError "empty seq") 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 sInject sStep sExtract) = M inject step extract - where - inject a = (a, sInject a) - step (_, x) a = let a' = f a (sExtract x) in (a', sStep x a') - extract (a, _) = a +-- | a moving window of a's, most recent at the front of the sequence +window :: Int -> Mealy a (Seq.Seq a) +window n = M Seq.singleton (\xs x -> Seq.take n (x Seq.<| xs)) id +{-# INLINEABLE window #-} --- | 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) +-- | binomial operator applied to last and this value +diff :: (a -> a -> b) -> Mealy a b +diff f = f <$> id <*> delay1 undefined -zeroModel1 :: Model1 -zeroModel1 = Model1 0 0 0 0 0 0 +-- | generalised diff function. +gdiff :: (a -> b) -> (a -> a -> b) -> Mealy a b +gdiff d0 d = M (\a -> (d0 a, a)) (\(_, a') a -> (d a a', a)) fst --- | 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 +-- | Unchanged since last time. +same :: (Eq b) => (a -> b) -> Mealy a Bool +same b = M (\a -> (True, b a)) (\(s, x) a -> (s && b a == x, x)) fst + +-- | Count observed values +countM :: (Ord a) => Mealy a (Map.Map a Int) +countM = M (`Map.singleton` 1) (\m k -> Map.insertWith (+) k 1 m) id + +-- | Sum values of a key-value pair. +sumM :: (Ord a, Additive b) => Mealy (a, b) (Map.Map a b) +sumM = M (uncurry Map.singleton) (\m (k, v) -> Map.insertWith (+) k v m) id + +-- | Convert a Mealy to a Mealy operating on lists. +listify :: Mealy a b -> Mealy [a] [b] +listify (M sExtract sStep sInject) = M extract step inject 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) + extract = fmap sExtract + step = zipWith sStep + inject = fmap sInject -- | A rough Median. -- The average absolute value of the stat is used to callibrate estimate drift towards the median @@ -547,9 +539,9 @@ } -- | onlineL1' takes a function and turns it into a `Mealy` 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) -> Mealy a (b, b) -onlineL1' i d f g = M inject step extract +onlineL1 :: + (Ord b, Field b, Absolute b) => b -> b -> (a -> b) -> (b -> b) -> Mealy a b +onlineL1 i d f g = snd <$> M inject step extract where inject a = let s = abs (f a) in Medianer s one (i * s) step (Medianer s c m) a = @@ -567,29 +559,9 @@ | f a > m = one | f a < m = negate one | otherwise = zero -{-# 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) -> Mealy 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 -> Mealy a a +maL1 :: (Ord a, Field a, Absolute a) => a -> a -> a -> Mealy 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 -> Mealy a a -absmaL1 i d r = fst <$> onlineL1' i d id (* r) -{-# INLINEABLE absmaL1 #-} diff --git a/src/Data/Mealy/Quantiles.hs b/src/Data/Mealy/Quantiles.hs --- a/src/Data/Mealy/Quantiles.hs +++ b/src/Data/Mealy/Quantiles.hs @@ -1,21 +1,28 @@ {-# LANGUAGE DataKinds #-} {-# LANGUAGE RebindableSyntax #-} -{-# LANGUAGE StrictData #-} +-- | Mealy quantile statistics. module Data.Mealy.Quantiles ( median, quantiles, digitize, + signalize, + OnlineTDigest (..), + emptyOnlineTDigest, + onlineInsert, + onlineCompress, + onlineForceCompress, ) where +import Control.Monad.ST import Data.Mealy import Data.Ord import Data.TDigest hiding (median) import Data.TDigest.Internal 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 Data.Vector.Algorithms.Heap qualified as VHeap +import Data.Vector.Unboxed qualified as VU import NumHask.Prelude hiding (fold) data OnlineTDigest = OnlineTDigest @@ -28,7 +35,7 @@ emptyOnlineTDigest :: Double -> OnlineTDigest emptyOnlineTDigest = OnlineTDigest (emptyTDigest :: TDigest n) 0 --- | decaying quantiles based on the tdigest library +-- | Mealy quantiles based on the tdigest library quantiles :: Double -> [Double] -> Mealy Double [Double] quantiles r qs = M inject step extract where @@ -38,6 +45,9 @@ where (OnlineTDigest t _ _) = onlineForceCompress x +-- | Mealy median using 'tdigest' +-- +-- The tdigest algorithm works best at extremes and can be unreliable in the centre. median :: Double -> Mealy Double Double median r = M inject step extract where @@ -50,7 +60,7 @@ onlineInsert' :: Double -> OnlineTDigest -> OnlineTDigest onlineInsert' x (OnlineTDigest td' n r) = OnlineTDigest - (insertCentroid (x, r ^^ (- (fromIntegral $ n + 1) :: Integer)) td') + (insertCentroid (x, r ^^ (-(fromIntegral $ n + 1) :: Integer)) td') (n + 1) r @@ -62,7 +72,7 @@ onlineCompress otd@(OnlineTDigest t _ _) | Data.TDigest.Tree.Internal.size t > relMaxSize * compression && Data.TDigest.Tree.Internal.size t > absMaxSize = - onlineForceCompress otd + onlineForceCompress otd | otherwise = otd where compression = 25 @@ -83,6 +93,12 @@ VHeap.sortBy (comparing snd) v VU.unsafeFreeze v +-- | A mealy that computes the running quantile bucket. For example, +-- in a scan, @digitize 0.9 [0.5]@ returns: +-- +-- * 0 if the current value is less than the current mealy median. +-- +-- * 1 if the current value is greater than the current mealy median. digitize :: Double -> [Double] -> Mealy Double Int digitize r qs = M inject step extract where @@ -100,3 +116,8 @@ else 1 ) <$> xs + +-- | transform an input to a [0,1] signal, via digitalization. +signalize :: Double -> [Double] -> Mealy Double Double +signalize r qs' = + (\x -> fromIntegral x / fromIntegral (length qs' + 1)) <$> digitize r qs' diff --git a/src/Data/Mealy/Simulate.hs b/src/Data/Mealy/Simulate.hs --- a/src/Data/Mealy/Simulate.hs +++ b/src/Data/Mealy/Simulate.hs @@ -1,11 +1,6 @@ -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE ExtendedDefaultRules #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE RankNTypes #-} {-# LANGUAGE NoImplicitPrelude #-} -{-# OPTIONS_GHC -Wall #-} -{-# OPTIONS_GHC -fno-warn-type-defaults #-} +-- | simulation to support testing of Mealy's using mwc-probability module Data.Mealy.Simulate ( rvs, rvsp, @@ -21,33 +16,32 @@ -- >>> :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 -- +-- >>> import Data.Mealy +-- >>> import Data.Mealy.Simulate +-- >>> gen <- create +-- >>> rvs gen 3 +-- [1.8005943761746166e-2,0.36444481359059255,-1.2939898115295387] +-- -- >>> rs <- rvs gen 10000 -- >>> fold (ma 1) rs --- -1.735737734197327e-3 +-- 1.29805301109162e-2 -- -- >>> fold (std 1) rs --- 0.9923615647768976 +-- 1.0126527176272948 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 -- +-- >>> rvsp gen 3 0.8 +-- [(1.8005943761746166e-2,7.074509906249835e-2),(0.36444481359059255,-0.7073208451897444),(-1.2939898115295387,-0.643930709405127)] +-- -- >>> rsp <- rvsp gen 10000 0.8 -- >>> fold (corr (ma 1) (std 1)) rsp --- 0.7933213647252008 +-- 0.8050112742986588 rvsp :: Gen (PrimState IO) -> Int -> Double -> IO [(Double, Double)] rvsp gen n c = do s0 <- rvs gen n diff --git a/test/test.hs b/test/test.hs deleted file mode 100644 --- a/test/test.hs +++ /dev/null @@ -1,16 +0,0 @@ -{-# OPTIONS_GHC -Wall #-} -{-# OPTIONS_GHC -fno-warn-unused-imports #-} - -module Main where - -import NumHask.Prelude -import Test.DocTest -import Data.Mealy - -main :: IO () -main = - doctest - [ "src/Data/Mealy.hs", - "src/Data/Mealy/Quantiles.hs", - "src/Data/Mealy/Simulate.hs" - ]