commit793d0bba901f116754b8602bf3e042662c60d3bfparentd8a2aaaa24130c8c7bf2a7a257dd4e4818d9f4c8Author:Jared Tobin <jared@jtobin.ca>Date:Mon, 21 Oct 2013 14:06:22 +1300 Add fromMassFunctionT, categorical example.Diffstat:

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

M | src/Measurable/Generic.hs | | | 27 | +++++++++++++++++++++++++++ |

2 files changed, 76 insertions(+), 4 deletions(-)diff --git a/src/Examples.hs b/src/Examples.hs@@ -7,6 +7,7 @@ import Data.HashMap.Strict (HashMap) import qualified Data.HashMap.Strict as HashMap import Data.Vector (singleton) import Measurable.Generic +import Numeric.SpecFunctions import Statistics.Distribution import System.Random.MWC import System.Random.MWC.Distributions @@ -32,8 +33,9 @@ genNormalSamples n m t g = replicateM n $ normal m (1 / t) g -- | Normal-gamma model. Note the resulting type is a probability measure on -- tuples. -- --- X | t ~ N(mu, 1/(lambda * t)) --- t ~ gamma(a, b) +-- X | t ~ N(mu, 1/(lambda * t)) +-- t ~ gamma(a, b) +-- (X, t) ~ NormalGamma(mu, lambda, a, b) normalGammaMeasure :: (Fractional r, PrimMonad m) => Int @@ -105,10 +107,44 @@ altNormalNormalGammaMeasure n a b mu lambda g = do normalSamples <- lift $ genNormalSamples n m t g fromObservationsT normalSamples +-- | A binomial density (with respect to counting measure). +binom :: Double -> Int -> Int -> Double +binom p n k + | n <= 0 = 0 + | k < 0 = 0 + | n < k = 0 + | otherwise = n `choose` k * p ^ k * (1 - p) ^ (n - k) + +-- | Generate a binomial measure from its mass function. +binomMeasure + :: (Applicative m, Monad m) + => Int + -> Double + -> MeasureT Double m Int +binomMeasure n p = fromMassFunctionT (return . binom p n) [0..n] + +-- | Note that we can handle all sorts of things that are 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). +data Group = A | B | C deriving (Enum, Eq, Ord, Show) + +categoricalOnGroupDensity :: Fractional a => Group -> a +categoricalOnGroupDensity g + | g == A = 0.3 + | g == B = 0.6 + | g == C = 0.1 + +-- | Here's a measure defined on the Group data type. +categoricalOnGroupMeasure + :: (Applicative m, Monad m, Fractional a) + => MeasureT a m Group +categoricalOnGroupMeasure = + fromMassFunctionT (return . categoricalOnGroupDensity) [A, B, C] + main :: IO () main = do - let nng = normalNormalGammaMeasure 100 2 6 1 0.5 - anng = altNormalNormalGammaMeasure 100 2 6 1 0.5 + let nng = normalNormalGammaMeasure 30 2 6 1 0.5 + anng = altNormalNormalGammaMeasure 30 2 6 1 0.5 m0 <- withSystemRandom . asGenIO $ \g -> expectationT id $ nng g @@ -121,9 +157,18 @@ main = do p1 <- withSystemRandom . asGenIO $ \g -> expectationT id $ 2 `to` 3 <$> anng g + binomialMean <- expectationT fromIntegral (binomMeasure 10 0.5) + + groupProbBC <- expectationT id + (containing [B, C] <$> categoricalOnGroupMeasure) + print m0 print m1 print p0 print p1 + print binomialMean + + print groupProbBC +diff --git a/src/Measurable/Generic.hs b/src/Measurable/Generic.hs@@ -5,7 +5,12 @@ module Measurable.Generic where import Control.Applicative import Control.Monad import Control.Monad.Trans.Cont +import qualified Data.Foldable as Foldable import Data.List +import Data.Monoid +import qualified Data.Set as Set +import Data.Traversable hiding (mapM) +import Numeric.Integration.TanhSinh type MeasureT r m a = ContT r m a @@ -17,6 +22,21 @@ measureT = runContT fromObservationsT :: (Monad m, Fractional r) => [a] -> ContT r m a fromObservationsT 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. +fromMassFunctionT + :: (Num r, Applicative f) + => (a -> f r) + -> [a] + -> ContT r f a +fromMassFunctionT p support = ContT $ \f -> + fmap 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`. @@ -46,6 +66,12 @@ to a b x | x >= a && x <= b = 1 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 => [a] -> a average xs = fst $ foldl' @@ -62,3 +88,4 @@ weightedAverageM :: (Fractional c, Monad m) => (a -> m c) -> [a] -> m c weightedAverageM f = liftM average . mapM f {-# INLINE weightedAverageM #-} +