diff --git a/mealy.cabal b/mealy.cabal
--- a/mealy.cabal
+++ b/mealy.cabal
@@ -1,74 +1,58 @@
 cabal-version: 2.4
 name:          mealy
-version:       0.0.3
-synopsis: See readme.md
+version:       0.1.0
+license:       BSD-3-Clause
+copyright:     Tony Day (c) AfterTimes
+maintainer:    tonyday567@gmail.com
+author:        Tony Day
+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@ provides support for computing statistics (such as an average or a standard deviation)
+     as current state. 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.
+    .
+    == Usage
+    .
+    >>> import Mealy
+    .
+    >>> fold ((,) <$> ma 0.9 <*> std 0.9) [1..100]
+    (91.00265621044142,9.472822805289121)
 
-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/mealy#readme
-bug-reports: https://github.com/tonyday567/mealy/issues
-build-type: Simple
+category:      folding
+build-type:    Simple
+tested-with:   GHC ==8.8.4 || ==8.10.4 || ==9.0.1 || ==9.2.0.20210821
+
 source-repository head
-  type: git
+  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:
+
+  hs-source-dirs:   src
   default-language: Haskell2010
-  default-extensions:
   ghc-options:
-    -Wall
-    -Wcompat
-    -Wincomplete-record-updates
-    -Wincomplete-uni-patterns
-    -Wredundant-constraints
-    -fwrite-ide-info
-    -hiedir=.hie
+    -Wall -Wcompat -Wincomplete-record-updates
+    -Wincomplete-uni-patterns -Wredundant-constraints -fwrite-ide-info
+    -hiedir=.hie -Wunused-packages
 
-test-suite test
-  type: exitcode-stdio-1.0
-  main-is: test.hs
-  hs-source-dirs:
-    test
   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.4
+    , base               >=4.13 && <5
+    , containers         ^>=0.6.2
+    , folds              ^>=0.7.6
+    , generic-lens       ^>=2.2.0
+    , lens               ^>=5.0.1
+    , matrix             ^>=0.3.6
+    , mwc-probability    ^>=2.3.1
+    , numhask            ^>=0.8.1
+    , numhask-array      ^>=0.9.1
+    , primitive          ^>=0.7.2
+    , profunctors        ^>=5.6.2
+    , tdigest            ^>=0.2.1
+    , text               ^>=1.2.4
+    , vector             ^>=0.12.3
+    , vector-algorithms  ^>=0.8.0
diff --git a/src/Data/Mealy.hs b/src/Data/Mealy.hs
--- a/src/Data/Mealy.hs
+++ b/src/Data/Mealy.hs
@@ -37,7 +37,7 @@
     online,
 
     -- * Statistics
-    -- $setup
+    -- $example-set
     ma,
     absma,
     sqma,
@@ -69,26 +69,25 @@
   )
 where
 
+import Control.Category
+import Control.Exception
 import Control.Lens hiding (Empty, Unwrapped, Wrapped, index, (:>), (|>))
 import Data.Fold hiding (M)
 import Data.Functor.Rep
 import Data.Generics.Labels ()
-
--- import qualified Numeric.LinearAlgebra as LA
-
+import Data.List (scanl')
 import qualified Data.Matrix as M
+import Data.Sequence (Seq)
 import qualified Data.Sequence as Seq
+import Data.Text (Text)
+import Data.Typeable (Typeable)
+import GHC.TypeLits
 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 NumHask.Prelude hiding (L1, asum, fold, id, (.))
 
 -- $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
@@ -99,6 +98,26 @@
 -- >>> xs2 <- rvs g 10000
 -- >>> xsp <- rvsp g 10000 0.8
 
+-- $example-set
+-- The doctest examples are composed from some random series generated with Data.Mealy.Simulate.
+--
+-- - 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 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, Typeable)
+
+instance Exception MealyError
+
 -- | A 'Mealy' is a triple of functions
 --
 -- * (a -> b) __inject__ Convert an input into the state type.
@@ -130,7 +149,7 @@
 --
 -- > 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 +194,26 @@
 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 +224,15 @@
 
 -- | 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.
+-- 9.713356299018187e-2
 --
--- ![ma chart](other/ex-ma.svg)
 ma :: (Divisive a, Additive a) => a -> Mealy a a
 ma r = online id (* r)
 {-# INLINEABLE ma #-}
@@ -207,7 +240,7 @@
 -- | absolute average
 --
 -- >>> fold (absma 1) xs0
--- 0.7894201075535578
+-- 0.8075705557429647
 absma :: (Divisive a, Signed a) => a -> Mealy a a
 absma r = online abs (* r)
 {-# INLINEABLE absma #-}
@@ -237,9 +270,8 @@
 -- 99.28328803163829
 --
 -- >>> fold (std 1) xs0
--- 0.9923523681261158
+-- 1.0126438036262801
 --
--- ![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 #-}
@@ -247,7 +279,7 @@
 -- | 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,7 +288,7 @@
 -- | 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)
@@ -267,7 +299,7 @@
 -- | 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
@@ -290,7 +322,7 @@
 -- \]
 --
 -- >>> 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
@@ -310,7 +342,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 +350,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
 
@@ -441,10 +473,10 @@
   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.
@@ -477,17 +509,8 @@
 -- >>> -- 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)
+-- 0.10000000000000009
 --
--- ![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
@@ -513,14 +536,16 @@
   }
   deriving (Eq, Show, Generic)
 
+-- | zeroised Model1
 zeroModel1 :: Model1
 zeroModel1 = Model1 0 0 0 0 0 0
 
 -- | Apply a model1 relationship using a single decay factor.
 --
 -- >>> :set -XOverloadedLabels
+-- >>> import Control.Lens
 -- >>> fold (depModel1 0.01 (zeroModel1 & #betaMa2X .~ 0.1)) xs0
--- -0.47228537123218206
+-- -0.4591515493154126
 depModel1 :: Double -> Model1 -> Mealy Double Double
 depModel1 r m1 =
   depState fX st
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
@@ -2,6 +2,7 @@
 {-# LANGUAGE RebindableSyntax #-}
 {-# LANGUAGE StrictData #-}
 
+-- | Mealy quantile statistics.
 module Data.Mealy.Quantiles
   ( median,
     quantiles,
@@ -9,6 +10,7 @@
   )
 where
 
+import Control.Monad.ST
 import Data.Mealy
 import Data.Ord
 import Data.TDigest hiding (median)
@@ -28,7 +30,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 +40,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 +55,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
 
@@ -83,6 +88,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,0.5,1]@ 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
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
@@ -6,6 +6,7 @@
 {-# 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 +22,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"
-  ]