commit5481d0f68053ba386c59244228d69503d85fc3b9parent61b5bd98535154a88821ad37e70412d6747674f6Author:Jared Tobin <jared@jtobin.ca>Date:Fri, 18 Oct 2013 11:18:47 +1300 Generalize, add instances, remove moment stuff.Diffstat:

M | src/Measurable.hs | | | 141 | +++++++++++++++++++++++++++++++++---------------------------------------------- |

M | tests/Test.hs | | | 18 | +++++++++++++----- |

2 files changed, 71 insertions(+), 88 deletions(-)diff --git a/src/Measurable.hs b/src/Measurable.hs@@ -1,5 +1,10 @@ +{-# LANGUAGE BangPatterns #-} + module Measurable where +import Control.Monad +import Data.List +import Data.Monoid import Control.Applicative import Numeric.Integration.TanhSinh @@ -22,106 +27,76 @@ import Numeric.Integration.TanhSinh -- rather than via the sigma algebra or measurable sets. Instead of taking -- a set and returning a probability, we can take a function and return a -- probability. - -type Measure a = (a -> Double) -> Double - --- | Once we have a measure, we use it by integrating against it. Take a +-- +-- Once we have a measure, we use it by integrating against it. Take a -- real-valued random variable (i.e., measurable function) /X/. The mean of X -- is E X = int_R X d(Px), for Px the distribution (image measure) of X. -- -- We can generalize this to work with arbitrary measurable functions and -- measures. Expectation can be defined by taking a measurable function and --- applying a measure to it - i.e., it's just function application. --- --- expectation :: Measure Double -> (Double -> Double) -> Double --- expectation --- :: ((Double -> Double) -> Double) -- (a -> b) --- -> (Double -> Double) -> Double -- a -> b --- --- ($) :: (a -> b) -> a -> b +-- applying a measure to it - i.e., just function application. -- --- expectation == ($) +-- So really, a Measure in this sense is an expression of a particular +-- computational process - expectation. We leave a callback to be +-- plugged in - a measurable function - and from there, we can finish the +-- computation and return a value. -- --- So there's no need to express expectation as its own operator. We can just --- define particular expectations that are interesting. Moments, for example. +-- This is equivalent to the type that forms the Continuation monad, albeit +-- with constant (Double) result type. The functor and monad instances follow +-- suit. The strength of the Monad instance is that it allows us to meld +-- together measures with different input types. --- NOTE want to add cumulants +newtype Measure a = Measure { measure :: (a -> Double) -> Double } --- | The nth raw moment. -rawMoment :: Integral a => Measure Double -> a -> Double -rawMoment mu n = mu (^^ n) +instance Functor Measure where + fmap = push --- | All positive raw moments. -rawMoments :: Measure Double -> [Double] -rawMoments mu = map (rawMoment mu) [1..] +instance Fractional a => Monoid (Measure a) where + mempty = identityMeasure + mappend = convolute --- | Alias for first raw moment. -mean :: Measure Double -> Double -mean mu = rawMoment mu 1 +instance Monad Measure where + return x = Measure (\f -> f x) + mu >>= f = Measure $ \d -> + measure mu $ \g -> + measure (f g) d --- | The nth central moment. --- --- NOTE slow-as in ghci. Need to memoize or something, or this might just --- disappear when compiling. -centralMoment :: Integral a => Measure Double -> a -> Double -centralMoment mu n = mu $ (^^ n) . \x -> x - rawMoment mu 1 +-- | The pushforward measure is obtained by 'pushing' a function onto an +-- existing measure. +push :: (a -> b) -> Measure a -> Measure b +push f mu = Measure pushforward + where pushforward g = measure mu $ g . f --- | All positive central moments. -centralMoments :: Measure Double -> [Double] -centralMoments mu = map (centralMoment mu) [1..] +-- | The mean is obtained by integrating against the identity function. +mean :: Measure Double -> Double +mean mu = measure mu id --- | Alias for second central moment. +-- | The variance is obtained by integrating against the usual function. variance :: Measure Double -> Double -variance mu = centralMoment mu 2 - --- | The nth normalized moment. -normalizedMoment :: Integral a => Measure Double -> a -> Double -normalizedMoment mu n = (/ (sd ^ n)) $ centralMoment mu n - where sd = sqrt $ centralMoment mu 2 +variance mu = measure mu (^ 2) - mean mu ^ 2 --- | All normalized moments. -normalizedMoments :: Measure Double -> [Double] -normalizedMoments mu = map (normalizedMoment mu) [1..] +-- | Create a measure from a collection of observations from some distribution. +fromObservations :: Fractional a => [a] -> Measure a +fromObservations xs = Measure $ \f -> average . map f $ xs --- | The moment generating function about a point. -momentGeneratingFunction :: Double -> Measure Double -> Double -momentGeneratingFunction t mu = mu $ exp . (* t) . id - --- | The cumulant generating function about a point. -cumulantGeneratingFunction :: Double -> Measure Double -> Double -cumulantGeneratingFunction t mu = log $ momentGeneratingFunction t mu - --- | We want two ways to create measures; empirically (i.e. from observations) --- or directly from some integrable function (i.e. a density). - --- | Construct an empirical measure from observations. -fromObservations :: [Double] -> Measure Double -fromObservations xs f = normalize . sum . map f $ xs - where normalize = (/ fromIntegral (length xs)) - --- | Construct a measure from a density function. -fromDensity :: (Double -> Double) -> (Double -> Double) -> Double -fromDensity d f = quadratureTanhSinh $ liftA2 (*) f d +-- | Create a measure from a density function. +fromDensity :: (Double -> Double) -> Measure Double +fromDensity d = Measure $ \f -> quadratureTanhSinh $ liftM2 (*) f d where quadratureTanhSinh = result . last . everywhere trap --- | For a random variable X, measurable function f, and measure P, we can --- construct the image (pushforward) measure P_(f X). -push :: (Double -> Double) -> Measure Double -> Measure Double -push f mu g = mu $ g . f - --- | Measure composition is convolution. This allows us to compose measures, --- independent of how they were constructed. --- --- This is (.) on the category of measures. -convolute :: Measure Double -> Measure Double -> Measure Double -convolute mu nu f = nu $ \y -> mu $ \x -> f $ x + y - --- | The identity measure. --- --- This is 'id' on the category of measures. -identity :: Measure Double -identity = fromObservations [0] - --- instance Functor Measure where --- fmap f mu = push f mu +-- | Measure addition is convolution. +convolute :: Num a => Measure a -> Measure a -> Measure a +convolute mu nu = Measure $ \f -> measure nu + $ \y -> measure mu + $ \x -> f (x + y) + +-- | The (sum) identity measure. +identityMeasure :: Fractional a => Measure a +identityMeasure = fromObservations [0] + +-- | Simple average. +average :: Fractional a => [a] -> a +average xs = fst $ foldl' + (\(!m, !n) x -> (m + (x - m) / fromIntegral (n + 1), n + 1)) (0, 0) xs +{-# INLINE mean #-}diff --git a/tests/Test.hs b/tests/Test.hs@@ -5,7 +5,7 @@ -- -- Push trigonometric functions on each and convolute the results. -- --- Push an exponential function on that, and calculate the mean of that +-- Push an exponential function on that, and calculate the mean of the resulting -- distribution. -- @@ -19,13 +19,21 @@ import System.Random.MWC.Distributions standardNormal = density $ normalDistr 0 1 main = do - g <- create - samples <- replicateM 30 $ exponential 1 g + expSamples <- withSystemRandom . asGenIO $ \g -> + replicateM 100 $ exponential 1 g + + normSamples <- withSystemRandom . asGenIO $ \g -> + replicateM 100 $ normal 0 1 g let mu = fromDensity standardNormal - nu = fromObservations samples + nu = fromObservations expSamples rho = convolute (push cos mu) (push sin nu) eta = push exp rho - print $ mean eta + putStrLn $ "mean of normal samples (should be around 0): " ++ + show (mean . fromObservations $ normSamples) + putStrLn $ "variance of normal samples (should be around 1): " ++ + show (variance . fromObservations $ normSamples) + putStrLn $ "let X ~ N(0, 1), Y ~ observed. mean of exp(cos X + sin Y): " ++ + show (mean eta)