commitba51e4e625f4c04d5b8cb65b7c7b17502409e274parent66b12a36f8dde6457ef3ee487d465dac11aed6cdAuthor:Jared Tobin <jared@jtobin.ca>Date:Thu, 2 Apr 2015 18:12:28 +1000 Add docs, Util module.Diffstat:

M | measurable.cabal | | | 4 | +++- |

M | src/Measurable/Core.hs | | | 137 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------- |

A | src/Measurable/Util.hs | | | 36 | ++++++++++++++++++++++++++++++++++++ |

3 files changed, 141 insertions(+), 36 deletions(-)diff --git a/measurable.cabal b/measurable.cabal@@ -38,9 +38,11 @@ source-repository head location: git://github.com/jtobin/measurable.git library - exposed-modules: Measurable.Core hs-source-dirs: src default-language: Haskell2010 + exposed-modules: + Measurable.Core + , Measurable.Util other-extensions: BangPatternsdiff --git a/src/Measurable/Core.hs b/src/Measurable/Core.hs@@ -6,18 +6,19 @@ module Measurable.Core where -import Control.Arrow import Control.Applicative -import Control.Foldl (Fold) -import qualified Control.Foldl as Foldl import Control.Monad import Control.Monad.Trans import Data.Foldable (Foldable) import qualified Data.Foldable as Foldable import Data.Functor.Identity import Data.Traversable +import Measurable.Util import Numeric.Integration.TanhSinh +-- | A hand-rolled continuation type. Exactly like the standard one you'd find +-- in @Control.Monad.Trans.Cont@, but without the supporting functions like +-- @callCC@, etc. included in that module. newtype ContT r m a = ContT { runContT :: (a -> m r) -> m r } type Cont r = ContT r Identity @@ -28,6 +29,16 @@ runCont m k = runIdentity $ runContT m (Identity . k) cont :: ((a -> r) -> r) -> Cont r a cont f = ContT $ \c -> Identity $ f (runIdentity . c) +-- | A measure can be represented by nothing more than a continuation with a +-- restricted output type corresponding to the reals. +-- +-- A @Functor@ instance implements pushforward or image measures - merely +-- @fmap@ a measurable function over a measure to create one. +-- +-- An @Applicative@ instance adds measure convolution, subtraction, and +-- multiplication by enabling a @Num@ instance via 'liftA2' and an implicit +-- marginalizing effect. A @Monad@ instance lumps the ability to create +-- measures from graphs of measures on top of that. type Measure a = Cont Double a type MeasureT m a = ContT Double m a @@ -49,6 +60,15 @@ instance Monad (ContT r m) where instance MonadTrans (ContT r) where lift m = ContT (m >>=) +-- | The 'integrate' function is just 'runCont' with its arguments reversed +-- in order to resemble the conventional mathematical notation, in which one +-- integrates a measurable function against a measure. +-- +-- >>> let mu = fromSamples [-1, 0, 1] +-- >>> expectation mu +-- 0.0 +-- >>> expectation mu +-- 1.0 integrate :: (a -> Double) -> Measure a -> Double integrate = flip runCont @@ -63,6 +83,16 @@ instance (Applicative m, Num a) => Num (ContT Double m a) where signum = fmap signum fromInteger = pure . fromInteger +-- | Create a 'Measure' from a probability mass function and its support, +-- provided as a foldable container. +-- +-- The requirement to supply the entire support is restrictive but necessary; +-- for approximations, consider using 'fromSamples' or +-- 'fromSamplingFunction'. +-- +-- >>> let mu = fromMassFunction (binomialPmf 10 0.2) [0..10] +-- >>> integrate fromIntegral mu +-- 2.0 fromMassFunction :: Foldable f => (a -> Double) -> f a -> Measure a fromMassFunction f support = cont $ \g -> weightedAverage (g /* f) support @@ -73,10 +103,28 @@ fromMassFunctionT :: (Applicative m, Traversable t) fromMassFunctionT f support = ContT $ \g -> fmap Foldable.sum . traverse (g //* (pure . f)) $ support +-- | Create a 'Measure' from a probability density function. +-- +-- Note that queries on measures constructed with @fromDensityFunction@ are +-- subject to numerical error due to the underlying dependency on quadrature! +-- +-- >>> let f x = 1 / (sqrt (2 * pi)) * exp (- (x ^ 2) / 2) +-- >>> let mu = fromDensityFunction f +-- >>> expectation mu +-- 0.0 +-- >>> variance mu +-- 1.0000000000000002 fromDensityFunction :: (Double -> Double) -> Measure Double fromDensityFunction d = cont $ \f -> quadratureTanhSinh $ f /* d where quadratureTanhSinh = result . last . everywhere trap +-- | Create a measure from a collection of observations. +-- +-- Useful for creating general purpose empirical measures. +-- +-- >>> let mu = fromSamples [(1, 2), (3, 4)] +-- >>> integrate (uncurry (+)) mu +-- 5.0 fromSamples :: Foldable f => f a -> Measure a fromSamples = cont . flip weightedAverage @@ -86,6 +134,8 @@ fromSamplesT -> MeasureT m a fromSamplesT = ContT . flip weightedAverageM +-- | Create a measure from a sampling function. Runs the sampling function +-- the provided number of times and runs 'fromSamples' on the result. fromSamplingFunction :: (Monad m, Applicative m) => (t -> m b) @@ -94,43 +144,44 @@ fromSamplingFunction -> MeasureT m b fromSamplingFunction f n g = (lift $ replicateM n (f g)) >>= fromSamplesT +-- | A simple alias for @fmap@. push :: (a -> b) -> Measure a -> Measure b push = fmap pushT :: Monad m => (a -> b) -> MeasureT m a -> MeasureT m b pushT = fmap --- | Expectation is integration against the identity function. +-- | The expectation of a measure is typically understood to be its expected +-- value, which is found by integrating it against the identity function. expectation :: Measure Double -> Double expectation = integrate id expectationT :: Applicative m => MeasureT m Double -> m Double expectationT = integrateT id +-- | The variance of a measure, as per the usual formula +-- @var X = E^2 X - EX^2@. variance :: Measure Double -> Double variance mu = integrate (^ 2) mu - expectation mu ^ 2 varianceT :: Applicative m => MeasureT m Double -> m Double varianceT mu = liftA2 (-) (integrateT (^ 2) mu) ((^ 2) <$> expectationT mu) -meanVariance :: Measure Double -> (Double, Double) -meanVariance = expectation &&& variance - -meanVarianceT :: Applicative m => MeasureT m Double -> m (Double, Double) -meanVarianceT mu = liftA2 (,) (expectationT mu) (varianceT mu) - +-- | The @nth@ raw moment of a 'Measure'. rawMoment :: Int -> Measure Double -> Double rawMoment n = integrate (^ n) rawMomentT :: (Applicative m, Monad m) => Int -> MeasureT m Double -> m Double rawMomentT n = integrateT (^ n) +-- | All raw moments of a 'Measure'. rawMoments :: Measure Double -> [Double] rawMoments mu = (`rawMoment` mu) <$> [1..] rawMomentsT :: (Applicative m, Monad m) => MeasureT m Double -> Int -> m [Double] rawMomentsT mu n = traverse (`rawMomentT` mu) $ take n [1..] +-- | The @nth@ central moment of a 'Measure'. centralMoment :: Int -> Measure Double -> Double centralMoment n mu = integrate (\x -> (x - rm) ^ n) $ mu where rm = rawMoment 1 mu @@ -140,62 +191,78 @@ centralMomentT n mu = integrateT (^ n) $ do rm <- lift $ rawMomentT 1 mu (subtract rm) <$> mu +-- | All central moments of a 'Measure'. centralMoments :: Measure Double -> [Double] centralMoments mu = (`centralMoment` mu) <$> [1..] centralMomentsT :: (Applicative m, Monad m) => MeasureT m Double -> Int -> m [Double] centralMomentsT mu n = traverse (`centralMomentT` mu) $ take n [1..] +-- | The moment generating function corresponding to a 'Measure'. +-- +-- >>> let mu = fromSamples [1..10] +-- >>> let mgfMu = momentGeneratingFunction mu +-- >>> fmap mgfMu [0, 0.5, 1] +-- [1.0,37.4649671547254,3484.377384533132] momentGeneratingFunction :: Measure Double -> Double -> Double momentGeneratingFunction mu t = integrate (exp . (* t)) mu +-- | The cumulant generating function corresponding to a 'Measure'. +-- +-- >>> let mu = fromSamples [1..10] +-- >>> let cgfMu = cumulantGeneratingFunction mu +-- >>> fmap cgfMu [0, 0.5, 1] +-- [0.0,3.6234062871236543,8.156044651432666] cumulantGeneratingFunction :: Measure Double -> Double -> Double cumulantGeneratingFunction mu = log . momentGeneratingFunction mu +-- | Calculates the volume of a 'Measure' over its entire space. Trivially 1 +-- for any probability measure. +-- +-- >>> let mu = fromSamples [1..10] +-- >>> volume mu +-- 1.0 volume :: Measure a -> Double volume = integrate $ const 1 volumeT :: Applicative m => MeasureT m a -> m Double volumeT = integrateT $ const 1 +-- | The cumulative distribution function corresponding to a 'Measure' +-- +-- >>> let mu = fromSamples [1..10] +-- >>> let cdfMu = cdf mu +-- >>> fmap cdfMu [0..10] +-- [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0] cdf :: Measure Double -> Double -> Double cdf mu x = expectation $ negativeInfinity `to` x <$> mu cdfT :: Applicative m => MeasureT m Double -> Double -> m Double cdfT mu x = expectationT $ negativeInfinity `to` x <$> mu +-- | A helpful utility for calculating the volume of a region in a measure +-- space. +-- +-- >>> let mu = fromSamples [1..10] +-- >>> integrate (2 `to` 8) mu +-- 0.7 to :: (Num a, Ord a) => a -> a -> a -> a to a b x | x >= a && x <= b = 1 | otherwise = 0 +-- | An analogue of 'to' for measures defined over non-ordered domains. +-- +-- >>> data Group = A | B | C deriving Eq +-- >>> let mu = fromSamples [A, A, A, B, A, B, C] +-- >>> integrate (containing [B]) mu +-- 0.2857142857142857 +-- >>> integrate (containing [A,C]) mu +-- 0.7142857142857143 +-- >>> integrate (containing [A,B,C]) mu +-- 1.0 containing :: (Num a, Eq b) => [b] -> b -> a containing xs x | x `elem` xs = 1 | otherwise = 0 -negativeInfinity :: Fractional a => a -negativeInfinity = negate $ 1 / 0 - -weightedAverage :: (Foldable f, Fractional r) => (a -> r) -> f a -> r -weightedAverage f = Foldl.fold (weightedAverageFold f) - -weightedAverageM - :: (Traversable t, Applicative f, Fractional r) - => (a -> f r) - -> t a - -> f r -weightedAverageM f = fmap (Foldl.fold averageFold) . traverse f - -weightedAverageFold :: Fractional r => (a -> r) -> Fold a r -weightedAverageFold f = Foldl.premap f averageFold - -averageFold :: Fractional a => Fold a a -averageFold = (/) <$> Foldl.sum <*> Foldl.genericLength - -(/*) :: (Num c, Applicative f) => f c -> f c -> f c -(/*) = liftA2 (*) - -(//*) :: (Num c, Applicative f, Applicative g) => f (g c) -> f (g c) -> f (g c) -(//*) = liftA2 (/*) -diff --git a/src/Measurable/Util.hs b/src/Measurable/Util.hs@@ -0,0 +1,36 @@ + +module Measurable.Util where + +import Control.Applicative +import Control.Foldl +import qualified Control.Foldl as Foldl +import Data.Foldable (Foldable) +import Data.Traversable + +negativeInfinity :: Fractional a => a +negativeInfinity = negate $ 1 / 0 + +weightedAverage :: (Foldable f, Fractional r) => (a -> r) -> f a -> r +weightedAverage f = Foldl.fold (weightedAverageFold f) + +weightedAverageM + :: (Traversable t, Applicative f, Fractional r) + => (a -> f r) + -> t a + -> f r +weightedAverageM f = fmap (Foldl.fold averageFold) . traverse f + +weightedAverageFold :: Fractional r => (a -> r) -> Fold a r +weightedAverageFold f = Foldl.premap f averageFold + +averageFold :: Fractional a => Fold a a +averageFold = (/) <$> Foldl.sum <*> Foldl.genericLength + +-- | Lifted multiplication. +(/*) :: (Num c, Applicative f) => f c -> f c -> f c +(/*) = liftA2 (*) + +-- | Doubly-lifted multiplication. +(//*) :: (Num c, Applicative f, Applicative g) => f (g c) -> f (g c) -> f (g c) +(//*) = liftA2 (/*) +