commiteb72ec258bc003a528a62f780c10612a0580f71bparentca72b1b4923b49b606c3746707b6dc021ef95b49Author:Jared Tobin <jared@jtobin.ca>Date:Thu, 24 Oct 2013 20:57:45 +1300 Create Measure/MeasureT types and lump them in Core.Diffstat:

M | src/Examples.hs | | | 125 | +++++++++++++++++++++++++++---------------------------------------------------- |

D | src/Examples/ChineseRestaurantProcess.hs | | | 82 | ------------------------------------------------------------------------------- |

D | src/Measurable.hs | | | 164 | ------------------------------------------------------------------------------- |

A | src/Measurable/Core.hs | | | 153 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |

D | src/Measurable/Generic.hs | | | 98 | ------------------------------------------------------------------------------- |

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

6 files changed, 204 insertions(+), 439 deletions(-)diff --git a/src/Examples.hs b/src/Examples.hs@@ -13,12 +13,35 @@ import Data.IntMap.Strict (IntMap) import qualified Data.IntMap.Strict as IntMap import Data.Vector (singleton) import qualified Data.Traversable as Traversable -import Measurable.Generic +import Measurable.Core import Numeric.SpecFunctions import Statistics.Distribution +import Statistics.Distribution.Normal +import Statistics.Distribution.Beta +import Statistics.Distribution.ChiSquared import System.Random.MWC import System.Random.MWC.Distributions +-- | Some workhorse densities (with respect to Lebesgue measure). +genNormal m v = density $ normalDistr m v +genBeta a b = density $ betaDistr a b +genChiSq d = density $ chiSquared d + +-- | Measures created from densities. Notice that the binomial measure has to +-- be treated differently than the measures absolutely continuous WRT Lebesgue +-- measure. +normalMeasure m v = fromDensityLebesgue $ genNormal m v +betaMeasure a b = fromDensityLebesgue $ genBeta a b +chiSqMeasure d = fromDensityLebesgue $ genChiSq d + +-- | A standard beta-binomial conjugate model. Notice how naturally it's +-- expressed using do-notation! +betaBinomialConjugate :: Double -> Double -> Int -> Measure Double Int +betaBinomialConjugate a b n = do + p <- betaMeasure a b + binomMeasure n p + +-- | Observe a gamma distribution. genGammaSamples :: (Applicative m, PrimMonad m) => Int @@ -28,6 +51,7 @@ genGammaSamples -> m [Double] genGammaSamples n a b g = replicateM n $ gamma a b g +-- | Observe a normal distributions. genNormalSamples :: (Applicative m, PrimMonad m) => Int @@ -54,10 +78,10 @@ normalGammaMeasure -> MeasureT r m (Double, Double) normalGammaMeasure n a b mu lambda g = do gammaSamples <- lift $ genGammaSamples n a b g - precision <- fromObservations gammaSamples + precision <- fromObservationsT gammaSamples normalSamples <- lift $ genNormalSamples n mu (lambda * precision) g - location <- fromObservations normalSamples + location <- fromObservationsT normalSamples return (location, precision) @@ -75,13 +99,14 @@ altNormalGammaMeasure -> MeasureT r m (HashMap String Double) altNormalGammaMeasure n a b mu lambda g = do gammaSamples <- lift $ genGammaSamples n a b g - precision <- fromObservations gammaSamples + precision <- fromObservationsT gammaSamples normalSamples <- lift $ genNormalSamples n mu (lambda * precision) g - location <- fromObservations normalSamples + location <- fromObservationsT normalSamples return $ HashMap.fromList [("location", location), ("precision", precision)] +-- | A normal-normal gamma conjugate model normalNormalGammaMeasure :: (Fractional r, Applicative m, PrimMonad m) => Int @@ -94,8 +119,9 @@ normalNormalGammaMeasure normalNormalGammaMeasure n a b mu lambda g = do (m, t) <- normalGammaMeasure n a b mu lambda g normalSamples <- lift $ genNormalSamples n m t g - fromObservations normalSamples + fromObservationsT normalSamples +-- | Alternate normal-normal gamma conjugate model. altNormalNormalGammaMeasure :: (Fractional r, Applicative m, PrimMonad m) => Int @@ -112,9 +138,9 @@ altNormalNormalGammaMeasure n a b mu lambda g = do t = fromMaybe (error "no precision!") $ HashMap.lookup "precision" parameterHash normalSamples <- lift $ genNormalSamples n m t g - fromObservations normalSamples + fromObservationsT normalSamples --- | A binomial density (with respect to counting measure). +-- | The binomial density. binom :: Double -> Int -> Int -> Double binom p n k | n <= 0 = 0 @@ -122,19 +148,21 @@ binom p n k | n < k = 0 | otherwise = n `choose` k * p ^ k * (1 - p) ^ (n - k) --- | Generate a binomial measure from its mass function. +-- | Generate a measure from the binomial density. binomMeasure :: (Applicative m, Monad m) => Int -> Double -> MeasureT Double m Int -binomMeasure n p = fromMassFunction (return . binom p n) [0..n] +binomMeasure n p = fromDensityCountingT (binom p n) [0..n] -- | Note that we can handle all sorts of things that have densities w/respect --- to counting measure. They don't necessarily have to have integral --- domains (or even have Ordered domains, though that's the case here). +-- to counting measure. They don't necessarily have to have domains that +-- are instances of Num (or even have Ordered domains, though that's the case +-- here). data Group = A | B | C deriving (Enum, Eq, Ord, Show) +-- | Density of a categorical measure. categoricalOnGroupDensity :: Fractional a => Group -> a categoricalOnGroupDensity g | g == A = 0.3 @@ -146,7 +174,7 @@ categoricalOnGroupMeasure :: (Applicative m, Monad m, Fractional r) => MeasureT r m Group categoricalOnGroupMeasure = - fromMassFunction (return . categoricalOnGroupDensity) [A, B, C] + fromDensityCountingT categoricalOnGroupDensity [A, B, C] -- | A gaussian mixture model, with mixing probabilities based on observed -- groups. Again, note that Group is not an instance of Num! We can compose @@ -166,78 +194,11 @@ gaussianMixtureModel -> Gen (PrimState m) -> MeasureT r m Double gaussianMixtureModel n observed g = do - s <- fromObservations observed + s <- fromObservationsT observed samples <- case s of A -> lift $ genNormalSamples n (-2) 1 g B -> lift $ genNormalSamples n 0 1 g C -> lift $ genNormalSamples n 1 1 g - fromObservations samples - --- | A bizarre random measure. -weirdMeasure - :: Fractional r - => [Group] - -> [Bool] - -> MeasureT r IO Bool -weirdMeasure [] acc = fromObservations acc -weirdMeasure (m:ms) acc - | m == A = do - j <- lift $ withSystemRandom . asGenIO $ uniform - k <- lift $ withSystemRandom . asGenIO $ uniform - if j - then weirdMeasure ms (k : acc) - else weirdMeasure ms acc - | otherwise = weirdMeasure ms acc - -main :: IO () -main = do - let nng = normalNormalGammaMeasure 30 2 6 1 0.5 - anng = altNormalNormalGammaMeasure 30 2 6 1 0.5 - m0 <- withSystemRandom . asGenIO $ \g -> - expectation id $ nng g - - m1 <- withSystemRandom . asGenIO $ \g -> - expectation id $ anng g - - p0 <- withSystemRandom . asGenIO $ \g -> - expectation id $ 2 `to` 3 <$> nng g - - p1 <- withSystemRandom . asGenIO $ \g -> - expectation id $ 2 `to` 3 <$> anng g - - binomialMean <- expectation fromIntegral (binomMeasure 10 0.5) - - groupProbBC <- expectation id - (containing [B, C] <$> categoricalOnGroupMeasure) - - let groupsObserved = [A, A, A, B, A, B, B, A, C, B, B, A, A, B, C, A, A] - categoricalFromObservationsMeasure = fromObservations groupsObserved - - groupsObservedProbBC <- expectation id - (containing [B, C] <$> categoricalFromObservationsMeasure) - - mixtureModelMean <- withSystemRandom . asGenIO $ \g -> - expectation id (gaussianMixtureModel 30 groupsObserved g) - - print m0 - print m1 - - print p0 - print p1 - - print binomialMean - - print groupProbBC - print groupsObservedProbBC - - print mixtureModelMean - - weirdProbabilityOfTrue <- expectation id $ - containing [True] <$> weirdMeasure [B, B, A, C, A, C, C, A, B, A] [] - - print weirdProbabilityOfTrue - - - + fromObservationsT samplesdiff --git a/src/Examples/ChineseRestaurantProcess.hs b/src/Examples/ChineseRestaurantProcess.hs@@ -1,82 +0,0 @@ -{-# LANGUAGE TemplateHaskell #-} - -import Control.Applicative -import Control.Lens -import Control.Monad.Trans -import Data.IntMap.Strict (IntMap) -import qualified Data.IntMap.Strict as IntMap -import Measurable.Generic - -data Table = Table { - _number :: {-# UNPACK #-} !Int - , _people :: {-# UNPACK #-} !Int - } deriving (Eq, Show) - -instance Ord Table where - t1 < t2 = _people t1 < _people t2 - -makeLenses ''Table - --- | Mass function for a given table. It's dependent on the state of the --- restaurant via 'n' and 'newestTable'. -tableMass n a newestTable table - | table^.number == newestTable^.number = a / (fromIntegral n + a) - | otherwise = fromIntegral (table^.people) / (fromIntegral n + a) - --- | A dependent measure over tables. -tableMeasure n a newestTable = - fromMassFunction (return . tableMass n a newestTable) - --- | A probability measure over restaurants, represented by IntMaps. -restaurantMeasure a restaurant = do - let numberOfCustomers = sumOf (traverse.people) restaurant - numberOfTables = lengthOf traverse restaurant - nextTableNum = succ numberOfTables - possibleTable = Table nextTableNum 1 - possibleRestaurant = IntMap.insert nextTableNum possibleTable restaurant - - table <- tableMeasure numberOfCustomers a possibleTable possibleRestaurant - - let newTable | table^.number == possibleTable^.number = table - | otherwise = table&people %~ succ - - return $ IntMap.insert (newTable^.number) newTable restaurant - --- | The Chinese Restaurant process. --- --- This implementation is dismally inefficient as-is, but appears to be --- correct. I think I need to look at doing memoization under the hood. -chineseRestaurantProcess n a = go n IntMap.empty - where go 0 restaurant = return restaurant - go j restaurant = restaurantMeasure a restaurant >>= go (pred j) - -main :: IO () -main = do - -- Easily verified by hand. - meanTinyRestaurant <- expectation (fromIntegral . lengthOf traverse) - (chineseRestaurantProcess 2 1) - - -- Easily verified by hand. - meanBiggerRestaurant <- expectation (fromIntegral . lengthOf traverse) - (chineseRestaurantProcess 3 1) - - meanBigRestaurant <- expectation (fromIntegral . lengthOf traverse) - (chineseRestaurantProcess 9 1) - - meanBigRestaurantAntisocial <- expectation (fromIntegral . lengthOf traverse) - (chineseRestaurantProcess 9 3) - - -- We can answer other questions by changing our transformation function. - -- Trivial example: the expected number of customers for a CRP observed for n - -- epochs is always n. - - differentQuestionSameMeasure <- - expectation (fromIntegral . sumOf (traverse.people)) - (chineseRestaurantProcess 9 3) - - print meanTinyRestaurant - print meanBiggerRestaurant - print meanBigRestaurant - print meanBigRestaurantAntisocial - print differentQuestionSameMeasure -diff --git a/src/Measurable.hs b/src/Measurable.hs@@ -1,164 +0,0 @@ -{-# LANGUAGE BangPatterns #-} - -module Measurable where - -import Control.Applicative -import Control.Monad -import Data.List -import Data.Monoid -import Numeric.Integration.TanhSinh - --- | A measure is a set function from some sigma-algebra to the extended real --- line. In practical terms we define probability in terms of measures; for a --- probability measure /P/ and some measurable set /A/, the measure of the set --- is defined to be that set's probability. --- --- For any sigma field, there is a one-to-one correspondence between measures --- and increasing linear functionals on its associated space of positive --- measurable functions. That is, --- --- P(A) = f(I_A) --- --- For A a measurable set and I_A its indicator function. So we can generally --- abuse notation to write something like P(I_A), even though we don't --- directly apply the measure to a function. --- --- So we can actually deal with measures in terms of measurable functions, --- 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. --- --- 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., just function application. --- --- 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. A measure is actually represented as a --- *program* that, given a measurable function, integrates that function. --- It's thus completely equivalent to the Continuation monad, albeit with --- a restricted result type. --- --- The Functor instance provides the ability to create pushforward/image --- measures. That's handled by good ol' fmap. The strength of the Monad --- instance is that it allows us to do Bayesian inference. That is, --- --- priorMeasure >>= likelihoodMeasure == posteriorPredictiveMeasure --- --- or, using arguably more modern terminology: --- --- parameterModel >>= dataModel == model --- --- Ex, given 'betaMeasure a b' and 'binomMeasure n p' functions that create --- the obvious measures, we can express a (Bayesian) beta-binomial model like --- so: --- --- betaBinomialConjugate :: Double -> Double -> Int -> Measure Double --- betaBinomialConjugate a b n = do --- p <- betaMeasure a b --- binomMeasure n p --- - -newtype Measure a = Measure { measure :: (a -> Double) -> Double } - -instance Num a => Num (Measure a) where - (+) = liftA2 (+) - (-) = liftA2 (-) - (*) = liftA2 (*) - abs = id - signum mu = error "fromInteger: not supported for Measures" - fromInteger = error "fromInteger: not supported for Measures" - -instance Fractional a => Monoid (Measure a) where - mempty = identityMeasure - mappend = (+) - -instance Functor Measure where - fmap f mu = Measure $ \g -> measure mu $ g . f -- pushforward/image measure - -instance Applicative Measure where - pure = return - (<*>) = ap - -instance Monad Measure where - return x = Measure (\f -> f x) - mu >>= f = Measure $ \d -> - measure mu $ \g -> - measure (f g) d - --- | The volume is obtained by integrating against a constant. This is '1' for --- any probability measure. -volume :: Measure a -> Double -volume mu = measure mu (const 1) - --- | The expectation is obtained by integrating against the identity function. -expectation :: Measure Double -> Double -expectation mu = measure mu id - --- | The variance is obtained by integrating against the usual function. -variance :: Measure Double -> Double -variance mu = measure mu (^ 2) - expectation mu ^ 2 - --- | Create a measure from a collection of observations from some distribution. -fromObservations :: Fractional a => [a] -> Measure a -fromObservations xs = Measure (`weightedAverage` xs) - --- | Create a measure from a density function. -fromDensity :: (Double -> Double) -> Measure Double -fromDensity d = Measure $ \f -> quadratureTanhSinh $ liftA2 (*) f d - where quadratureTanhSinh = result . last . everywhere trap - --- | Create a measure from a mass function. -fromMassFunction :: (a -> Double) -> [a] -> Measure a -fromMassFunction p support = Measure $ \f -> - sum . map (liftA2 (*) f p) $ support - --- | The (sum) identity measure. -identityMeasure :: Fractional a => Measure a -identityMeasure = fromObservations [] - --- | 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 average #-} - --- | Weighted average. -weightedAverage :: Fractional c => (a -> c) -> [a] -> c -weightedAverage f = average . map f -{-# INLINE weightedAverage #-} - --- | Integrate from a to b. -to :: (Num a, Ord a) => a -> a -> a -> a -to a b x | x >= a && x <= b = 1 - | otherwise = 0 - --- | Cumulative distribution function for a measure. --- --- Really cool; works perfectly in both discrete and continuous cases. --- --- > let f = cdf (fromObservations [1..10]) --- > cdf 0 --- 0.0 --- > cdf 1 --- 0.1 --- > cdf 10 --- 1.0 --- > cdf 11 --- 1.0 --- --- > let g = cdf (fromDensity standardNormal) --- > cdf 0 --- 0.504 --- > cdf 1 --- 0.838 --- > cdf (1 / 0) --- 0.999 -cdf :: Measure Double -> Double -> Double -cdf mu b = expectation $ negate (1 / 0) `to` b <$> mu -diff --git a/src/Measurable/Core.hs b/src/Measurable/Core.hs@@ -0,0 +1,153 @@ +{-# LANGUAGE BangPatterns #-} + +module Measurable.Core where + +import Control.Applicative +import Control.Monad +import Control.Monad.Trans.Cont +import Data.Foldable (Foldable) +import qualified Data.Foldable as Foldable +import Data.Monoid +import qualified Data.Set as Set +import Data.Traversable hiding (mapM) +import Numeric.Integration.TanhSinh + +-- | A measure is represented as a continuation. +type Measure r a = Cont r a +type MeasureT r m a = ContT r m a + +-- | A more measure-y alias for runCont. +measure :: Measure r a -> (a -> r) -> r +measure = runCont + +measureT :: MeasureT r m a -> (a -> m r) -> m r +measureT = runContT + +-- | Things like convolution are trivially expressed by lifted arithmetic +-- operators. +instance (Monad m, Num a) => Num (ContT r m a) where + (+) = liftA2 (+) + (-) = liftA2 (-) + (*) = liftA2 (*) + abs = id + signum = error "signum: not supported for Measures" + fromInteger = error "fromInteger: not supported for measures" + +-- | Create a measure from a density w/respect to counting measure. +fromDensityCounting + :: (Num r, Functor f, Foldable f) + => (a -> r) + -> f a + -> Measure r a +fromDensityCounting f support = cont $ \g -> + Foldable.sum . fmap (liftA2 (*) g f) $ support + +fromDensityCountingT + :: (Num r, Applicative m, Traversable t, Monad m) + => (a -> r) + -> t a + -> MeasureT r m a +fromDensityCountingT p support = ContT $ \f -> + fmap Foldable.sum . traverse (liftA2 (liftA2 (*)) f (return . p)) $ support + +-- | Create a measure from a density w/respect to Lebesgue measure. +-- +-- NOTE The quality of this implementation depends entirely on the underlying +-- quadrature routine. As we're presently using the +-- Numeric.Integration.TanhSinh module, the types are also constricted +-- to Doubles. +-- +-- This is included moreso for interest's sake, and doesn't produce +-- particularly accurate results. +fromDensityLebesgue :: (Double -> Double) -> Measure Double Double +fromDensityLebesgue d = cont $ \f -> quadratureTanhSinh $ liftA2 (*) f d + where quadratureTanhSinh = result . last . everywhere trap + +-- | Create a measure from observations sampled from some distribution. +fromObservations + :: (Functor f, Foldable f, Fractional r) + => f a + -> Measure r a +fromObservations = cont . flip weightedAverage + +fromObservationsT + :: (Applicative m, Monad m, Fractional r, Traversable f) + => f a + -> MeasureT r m a +fromObservationsT = ContT . flip weightedAverageM + +-- | Integration of a function against a measure is just running a +-- continuation. +integrate :: (a -> r) -> Measure r a -> r +integrate = flip measure + +integrateT :: Monad m => (a -> r) -> MeasureT r m a -> m r +integrateT f = (`measureT` (return . f)) + +-- | Expectation is integration against the identity function. +expectation :: Measure r r -> r +expectation = integrate id + +expectationT :: Monad m => MeasureT r m r -> m r +expectationT = integrateT id + +-- | The variance is obtained by integrating against the usual function. +variance :: Num r => Measure r r -> r +variance mu = measure mu (^ 2) - expectation mu ^ 2 + +varianceT :: (Monad m, Num r) => MeasureT r m r -> m r +varianceT mu = liftM2 (-) (measureT mu (return . (^ 2))) + (liftM (^ 2) (expectationT mu)) + +-- | The measure of the underlying space. This is trivially 1 for any +-- probability measure. +volume :: Num r => Measure r r -> r +volume = integrate (const 1) + +volumeT :: (Num r, Monad m) => MeasureT r m r -> m r +volumeT = integrateT (const 1) + +-- | Cumulative distribution function. Only makes sense for Fractional/Ord +-- inputs. +cdf :: (Fractional r, Ord r) => Measure r r -> r -> r +cdf mu x = expectation $ negativeInfinity `to` x <$> mu + +cdfT :: (Fractional r, Ord r, Monad m) => MeasureT r m r -> r -> m r +cdfT mu x = expectationT $ negativeInfinity `to` x <$> mu + +-- | Indicator function for the interval a <= x <= b. Useful for integrating +-- from a to b. +to :: (Num a, Ord a) => a -> a -> a -> a +to a b x | x >= a && x <= b = 1 + | otherwise = 0 + +-- | Integrate over an ordered, discrete set. +containing :: (Num a, Ord b) => [b] -> b -> a +containing xs x | x `Set.member` set = 1 + | otherwise = 0 + where set = Set.fromList xs + +-- | End of the line. +negativeInfinity :: Fractional a => a +negativeInfinity = negate (1 / 0) + +-- | Simple average. +average :: (Fractional a, Foldable f) => f a -> a +average xs = fst $ Foldable.foldl' + (\(!m, !n) x -> (m + (x - m) / fromIntegral (n + 1), n + 1)) (0, 0) xs +{-# INLINE average #-} + +-- | Weighted average. +weightedAverage :: (Functor f, Foldable f, Fractional c) => (a -> c) -> f a -> c +weightedAverage f = average . fmap f +{-# INLINE weightedAverage #-} + +-- | Monadic weighted average. +weightedAverageM + :: (Fractional c, Traversable f, Monad m, Applicative m) + => (a -> m c) + -> f a + -> m c +weightedAverageM f = liftM average . traverse f +{-# INLINE weightedAverageM #-} +diff --git a/src/Measurable/Generic.hs b/src/Measurable/Generic.hs@@ -1,98 +0,0 @@ -{-# LANGUAGE BangPatterns #-} - -module Measurable.Generic where - -import Control.Applicative -import Control.Monad -import Control.Monad.Trans.Cont -import Data.Foldable (Foldable) -import qualified Data.Foldable as Foldable -import Data.List -import Data.Monoid -import qualified Data.Set as Set -import Data.Traversable hiding (mapM) - -type MeasureT r m a = ContT r m a - --- | A more measure-y alias for runContT. -measureT :: MeasureT r m a -> (a -> m r) -> m r -measureT = runContT - --- | Create a measure from observations (samples) from some distribution. -fromObservations - :: (Applicative m, Monad m, Fractional r, Traversable f) - => f a - -> MeasureT r m a -fromObservations xs = ContT (`weightedAverageM` xs) - --- A mass function is close to universal when dealing with discrete objects, but --- the problem is that we need to create it over the entire support. In terms --- of getting answers out, that breaks down immediately for something as trivial --- as the natural numbers. --- --- Maybe we can use something like an 'observed support'. You can probably get --- inspiration from how the Dirichlet process is handled in practice. -fromMassFunction - :: (Num r, Applicative f, Traversable t) - => (a -> f r) - -> t a - -> MeasureT r f a -fromMassFunction p support = ContT $ \f -> - fmap Foldable.sum . traverse (liftA2 (liftA2 (*)) f p) $ support - --- | Expectation is obtained by integrating against the identity function. We --- provide an additional function for mapping the input type to the output --- type -- if these are equivalent, this would just `id`. --- --- NOTE should we have this transformation handled elsewhere? I.e. make fmap --- responsible for transforming the type? -expectation :: Monad m => (a -> r) -> MeasureT r m a -> m r -expectation f = (`measureT` (return . f)) - --- | The volume is obtained by integrating against a constant. This is '1' for --- any probability measure. -volume :: (Num r, Monad m) => MeasureT r m r -> m r -volume mu = measureT mu (return . const 1) - --- | Cumulative distribution function. Only makes sense for Fractional/Ord --- inputs. Lots of potentially interesting cases where this isn't necessarily --- true. -cdf :: (Fractional r, Ord r, Monad m) => MeasureT r m r -> r -> m r -cdf mu x = expectation id $ (negativeInfinity `to` x) <$> mu - --- | Integrate from a to b. -to :: (Num a, Ord a) => a -> a -> a -> a -to a b x | x >= a && x <= b = 1 - | otherwise = 0 - --- | End of the line. -negativeInfinity :: Fractional a => a -negativeInfinity = negate (1 / 0) - --- | Integrate over an ordered, discrete set. -containing :: (Num a, Ord b) => [b] -> b -> a -containing xs x | x `Set.member` set = 1 - | otherwise = 0 - where set = Set.fromList xs - --- | Simple average. -average :: (Fractional a, Foldable f) => f a -> a -average xs = fst $ Foldable.foldl' - (\(!m, !n) x -> (m + (x - m) / fromIntegral (n + 1), n + 1)) (0, 0) xs -{-# INLINE average #-} - --- | Weighted average. -weightedAverage :: (Functor f, Foldable f, Fractional c) => (a -> c) -> f a -> c -weightedAverage f = average . fmap f -{-# INLINE weightedAverage #-} - --- | Monadic weighted average. -weightedAverageM - :: (Fractional c, Traversable f, Monad m, Applicative m) - => (a -> m c) - -> f a - -> m c -weightedAverageM f = liftM average . traverse f -{-# INLINE weightedAverageM #-} - -diff --git a/tests/Test.hs b/tests/Test.hs@@ -3,7 +3,7 @@ import Control.Applicative import Control.Monad import Data.Vector (singleton) -import Measurable +import Measurable.Core import Numeric.SpecFunctions import Statistics.Distribution hiding (mean, variance) import Statistics.Distribution.Normal @@ -27,11 +27,10 @@ binom p n k -- | Measures created from densities. Notice that the binomial measure has to -- be treated differently than the measures absolutely continuous WRT Lebesgue -- measure. -normalMeasure m v = fromDensity $ genNormal m v -betaMeasure a b = fromDensity $ genBeta a b -chiSqMeasure d = fromDensity $ genChiSq d -binomMeasure n p = fromMassFunction (binom p n . truncate) - (fromIntegral <$> [0..n] :: [Double]) +normalMeasure m v = fromDensityLebesgue $ genNormal m v +betaMeasure a b = fromDensityLebesgue $ genBeta a b +chiSqMeasure d = fromDensityLebesgue $ genChiSq d +binomMeasure n p = fromDensityCounting (binom p n) [0..n] -- | Sampling functions. generateExpSamples n l g = replicateM n (exponential l g) @@ -39,7 +38,7 @@ generateNormalSamples n m v g = replicateM n (normal m v g) -- | A standard beta-binomial conjugate model. Notice how naturally it's -- expressed using do-notation! -betaBinomialConjugate :: Double -> Double -> Int -> Measure Double +betaBinomialConjugate :: Double -> Double -> Int -> Measure Double Int betaBinomialConjugate a b n = do p <- betaMeasure a b binomMeasure n p @@ -114,13 +113,9 @@ main = do putStrLn "X | p ~ binomial(10, p)" putStrLn "p ~ beta(1, 4)" - putStrLn "" - putStrLn "(integrate out p..)" - putStrLn "" - - let phi = betaBinomialConjugate 1 4 10 + let phi = fromIntegral <$> betaBinomialConjugate 1 4 10 - putStrLn $ "E(X): " + putStrLn $ "E(X) (unconditional): " ++ show (expectation phi) putStrLn $ "P(X == 5): "