hakaru 0.1.3 → 0.1.4
raw patch · 11 files changed
+279/−294 lines, 11 filesnew-uploader
Files
- Language/Hakaru/Distribution.hs +100/−96
- Language/Hakaru/ImportanceSampler.hs +27/−17
- Language/Hakaru/Metropolis.hs +102/−90
- Language/Hakaru/Sampler.hs +11/−9
- Language/Hakaru/Types.hs +6/−4
- Language/Hakaru/Util/Coda.hs +8/−0
- Language/Hakaru/Util/Extras.hs +13/−12
- Language/Hakaru/Util/Visual.hs +2/−2
- Tests/Distribution.hs +3/−58
- Tests/ImportanceSampler.hs +6/−5
- hakaru.cabal +1/−1
Language/Hakaru/Distribution.hs view
@@ -3,7 +3,10 @@ module Language.Hakaru.Distribution where -import System.Random+import Control.Monad+import Control.Monad.Primitive+import Control.Monad.Loops+import qualified System.Random.MWC as MWC import Language.Hakaru.Mixture import Language.Hakaru.Types import Data.Ix@@ -18,59 +21,58 @@ dirac :: (Eq a) => a -> Dist a dirac theta = Dist {logDensity = (\ (Discrete x) -> if x == theta then 0 else log 0),- distSample = (\ g -> (Discrete theta,g))}+ distSample = (\ _ -> return $ Discrete theta)} bern :: Double -> Dist Bool bern p = Dist {logDensity = (\ (Discrete x) -> log (if x then p else 1 - p)),- distSample = (\ g -> case randomR (0, 1) g of- (t, g') -> (Discrete $ t <= p, g'))}+ distSample = (\ g -> do t <- MWC.uniformR (0,1) g+ return $ Discrete (t <= p))} uniform :: Double -> Double -> Dist Double uniform lo hi = let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log (recip (hi' - lo')) uniformLogDensity _ _ _ = log 0 in Dist {logDensity = (\ (Lebesgue x) -> uniformLogDensity lo hi x),- distSample = (\ g -> mapFst Lebesgue $ randomR (lo, hi) g)}+ distSample = (\ g -> liftM Lebesgue $ MWC.uniformR (lo, hi) g)} -uniformD :: (Ix a, Random a) => a -> a -> Dist a+uniformD :: (Ix a, MWC.Variate a) => a -> a -> Dist a uniformD lo hi = let uniformLogDensity lo' hi' x | lo' <= x && x <= hi' = log density uniformLogDensity _ _ _ = log 0 density = recip (fromInteger (toInteger (rangeSize (lo,hi)))) in Dist {logDensity = (\ (Discrete x) -> uniformLogDensity lo hi x),- distSample = (\ g -> mapFst Discrete $ randomR (lo, hi) g)}+ distSample = (\ g -> liftM Discrete $ MWC.uniformR (lo, hi) g)} -marsaglia :: (RandomGen g, Random a, Ord a, Floating a) => g -> ((a, a), g)-marsaglia g0 = -- "Marsaglia polar method"- let (x, g1) = randomR (-1,1) g0- (y, g ) = randomR (-1,1) g1- s = x * x + y * y- q = sqrt ((-2) * log s / s)- in if 1 >= s && s > 0 then ((x * q, y * q), g) else marsaglia g+marsaglia :: (MWC.Variate a, Ord a, Floating a, PrimMonad m) => PRNG m -> m (a, a)+marsaglia g = do -- "Marsaglia polar method"+ x <- MWC.uniformR (-1,1) g+ y <- MWC.uniformR (-1,1) g+ let s = x * x + y * y+ q = sqrt ((-2) * log s / s)+ if 1 >= s && s > 0 then return (x * q, y * q) else marsaglia g -choose :: (RandomGen g) => Mixture k -> g -> (k, Prob, g)-choose (Mixture m) g0 =+choose :: (PrimMonad m) => Mixture k -> PRNG m -> m (k, Prob)+choose (Mixture m) g = do let peak = maximum (M.elems m) unMix = M.map (LF.fromLogFloat . (/peak)) m total = M.foldl' (+) (0::Double) unMix- (p, g) = randomR (0, total) g0- f !k !v b !p0 = let p1 = p0 + v in if p <= p1 then k else b p1+ p <- MWC.uniformR (0, total) g+ let f !k !v b !p0 = let p1 = p0 + v in if p <= p1 then k else b p1 err p0 = error ("choose: failure p0=" ++ show p0 ++ " total=" ++ show total ++ " size=" ++ show (M.size m))- in (M.foldrWithKey f err unMix 0, LF.logFloat total * peak, g)+ return $ (M.foldrWithKey f err unMix 0, LF.logFloat total * peak) -chooseIndex :: (RandomGen g) => [Double] -> g -> (Int, g)-chooseIndex probs g0 =- let (p, g) = random g0- k = fromMaybe (error ("chooseIndex: failure p=" ++ show p))- (findIndex (p <=) (scanl1 (+) probs))- in (k, g)+chooseIndex :: (PrimMonad m) => [Double] -> PRNG m -> m Int+chooseIndex probs g = do+ p <- MWC.uniform g+ return $ fromMaybe (error ("chooseIndex: failure p=" ++ show p))+ (findIndex (p <=) (scanl1 (+) probs)) -normal_rng :: (Real a, Floating a, Random a, RandomGen g) =>- a -> a -> g -> (a, g)-normal_rng mu sd g | sd > 0 = case marsaglia g of- ((x, _), g1) -> (mu + sd * x, g1)+normal_rng :: (Real a, Floating a, MWC.Variate a, PrimMonad m) =>+ a -> a -> PRNG m -> m a+normal_rng mu sd g | sd > 0 = do (x, _) <- marsaglia g+ return (mu + sd * x) normal_rng _ _ _ = error "normal: invalid parameters" normalLogDensity :: Floating a => a -> a -> a -> a@@ -81,24 +83,25 @@ normal :: Double -> Double -> Dist Double normal mu sd = Dist {logDensity = normalLogDensity mu sd . fromLebesgue,- distSample = mapFst Lebesgue . normal_rng mu sd}+ distSample = (\g -> liftM Lebesgue $ normal_rng mu sd g)} categoricalLogDensity :: (Eq b, Floating a) => [(b, a)] -> b -> a categoricalLogDensity list x = log $ fromMaybe 0 (lookup x list)-categoricalSample :: (Num b, Ord b, RandomGen g, Random b) =>- [(t,b)] -> g -> (t, g)-categoricalSample list g = (elem', g1)- where- (p, g1) = randomR (0, total) g++categoricalSample :: (Num b, Ord b, PrimMonad m, MWC.Variate b) =>+ [(t,b)] -> PRNG m -> m t+categoricalSample list g = do+ let total = sum $ map snd list+ p <- MWC.uniformR (0, total) g+ let sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list elem' = fst $ head $ filter (\(_,p0) -> p <= p0) sumList- sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list- total = sum $ map snd list+ return elem' categorical :: Eq a => [(a,Double)] -> Dist a categorical list = Dist {logDensity = categoricalLogDensity list . fromDiscrete,- distSample = mapFst Discrete . categoricalSample list}+ distSample = (\g -> liftM Discrete $ categoricalSample list g)} -lnFact :: Integer -> Double+lnFact :: Int -> Double lnFact = logFactorial -- Makes use of Atkinson's algorithm as described in:@@ -106,8 +109,8 @@ -- -- Further discussion at: -- http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/-poisson_rng :: (RandomGen g) => Double -> g -> (Integer, g)-poisson_rng lambda g0 = make_poisson g0+poisson_rng :: (PrimMonad m) => Double -> PRNG m -> m Int+poisson_rng lambda g' = make_poisson g' where smu = sqrt lambda b = 0.931 + 2.53*smu a = -0.059 + 0.02483*b@@ -115,53 +118,53 @@ arep = 1.1239 + 1.1368/(b-3.4) lnlam = log lambda - make_poisson :: (RandomGen g) => g -> (Integer,g)- make_poisson g = let (u, g1) = randomR (-0.5,0.5) g- (v, g2) = randomR (0,1) g1- us = 0.5 - abs u- k = floor $ (2*a / us + b)*u + lambda + 0.43 in- case () of- () | us >= 0.07 && v <= vr -> (k, g2)- () | k < 0 -> make_poisson g2- () | us <= 0.013 && v > us -> make_poisson g2- () | accept_region us v k -> (k, g2)- _ -> make_poisson g2+ make_poisson :: (PrimMonad m) => PRNG m -> m Int+ make_poisson g = do u <- MWC.uniformR (-0.5,0.5) g+ v <- MWC.uniformR (0,1) g+ let us = 0.5 - abs u+ k = floor $ (2*a / us + b)*u + lambda + 0.43+ case () of+ () | us >= 0.07 && v <= vr -> return k+ () | k < 0 -> make_poisson g+ () | us <= 0.013 && v > us -> make_poisson g+ () | accept_region us v k -> return k+ _ -> make_poisson g - accept_region :: Double -> Double -> Integer -> Bool+ accept_region :: Double -> Double -> Int -> Bool accept_region us v k = log (v * arep / (a/(us*us)+b)) <= -lambda + (fromIntegral k)*lnlam - lnFact k -poisson :: Double -> Dist Integer+poisson :: Double -> Dist Int poisson l = let poissonLogDensity l' x | l' > 0 && x> 0 = (fromIntegral x)*(log l') - lnFact x - l' poissonLogDensity l' x | x==0 = -l' poissonLogDensity _ _ = log 0 in Dist {logDensity = poissonLogDensity l . fromDiscrete,- distSample = mapFst Discrete . poisson_rng l}+ distSample = (\g -> liftM Discrete $ poisson_rng l g)} -- Direct implementation of "A Simple Method for Generating Gamma Variables" -- by George Marsaglia and Wai Wan Tsang.-gamma_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)+gamma_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double gamma_rng shape _ _ | shape <= 0.0 = error "gamma: got a negative shape paramater" gamma_rng _ scl _ | scl <= 0.0 = error "gamma: got a negative scale paramater"-gamma_rng shape scl g | shape < 1.0 = (gvar2, g2)- where (gvar1, g1) = gamma_rng (shape + 1) scl g- (w, g2) = randomR (0,1) g1- gvar2 = scl * gvar1 * (w ** recip shape) -gamma_rng shape scl g = +gamma_rng shape scl g | shape < 1.0 = do gvar1 <- gamma_rng (shape + 1) scl g+ w <- MWC.uniformR (0,1) g+ return $ scl * gvar1 * (w ** recip shape)+gamma_rng shape scl g = do let d = shape - 1/3 c = recip $ sqrt $ 9*d -- Algorithm recommends inlining normal generator- n = normal_rng 1 c- (v, g2) = until (\y -> fst y > 0.0) (\ (_, g') -> normal_rng 1 c g') (n g)- x = (v - 1) / c+ -- n = normal_rng 1 c+ v <- iterateUntil (> 0.0) $ normal_rng 1 c g+ -- (v, g2) = until (\y -> fst y > 0.0) (\ (_, g') -> normal_rng 1 c g') (n g)+ let x = (v - 1) / c sqr = x * x v3 = v * v * v- (u, g3) = randomR (0.0, 1.0) g2- accept = u < 1.0 - 0.0331*(sqr*sqr) || log u < 0.5*sqr + d*(1.0 - v3 + log v3)- in case accept of- True -> (scl*d*v3, g3)- False -> gamma_rng shape scl g3+ u <- MWC.uniformR (0.0, 1.0) g+ let accept = u < 1.0 - 0.0331*(sqr*sqr) || log u < 0.5*sqr + d*(1.0 - v3 + log v3)+ case accept of+ True -> return $ scl*d*v3+ False -> gamma_rng shape scl g gammaLogDensity :: Double -> Double -> Double -> Double gammaLogDensity shape scl x | x>= 0 && shape > 0 && scl > 0 =@@ -170,20 +173,20 @@ gamma :: Double -> Double -> Dist Double gamma shape scl = Dist {logDensity = gammaLogDensity shape scl . fromLebesgue,- distSample = mapFst Lebesgue . gamma_rng shape scl}+ distSample = (\g -> liftM Lebesgue $ gamma_rng shape scl g)} -beta_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)-beta_rng a b g | a <= 1.0 && b <= 1.0 =- let (u, g1) = randomR (0.0, 1.0) g- (v, g2) = randomR (0.0, 1.0) g1- x = u ** (recip a)+beta_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double+beta_rng a b g | a <= 1.0 && b <= 1.0 = do+ u <- MWC.uniformR (0.0, 1.0) g+ v <- MWC.uniformR (0.0, 1.0) g+ let x = u ** (recip a) y = v ** (recip b)- in case (x+y) <= 1.0 of- True -> (x / (x + y), g2)- False -> beta_rng a b g2-beta_rng a b g = let (ga, g1) = gamma_rng a 1 g- (gb, g2) = gamma_rng b 1 g1- in (ga / (ga + gb), g2)+ case (x+y) <= 1.0 of+ True -> return $ x / (x + y)+ False -> beta_rng a b g+beta_rng a b g = do ga <- gamma_rng a 1 g+ gb <- gamma_rng b 1 g+ return $ ga / (ga + gb) betaLogDensity :: Double -> Double -> Double -> Double betaLogDensity _ _ x | x < 0 || x > 1 = error "beta: value must be between 0 and 1"@@ -196,33 +199,34 @@ beta :: Double -> Double -> Dist Double beta a b = Dist {logDensity = betaLogDensity a b . fromLebesgue,- distSample = mapFst Lebesgue . beta_rng a b}+ distSample = (\g -> liftM Lebesgue $ beta_rng a b g)} -laplace_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)-laplace_rng mu sd g = sample (randomR (0.0, 1.0) g)- where sample (u, g1) = case u < 0.5 of- True -> (mu + sd * log (u + u), g1)- False -> (mu - sd * log (2.0 - u - u), g1)+laplace_rng :: (PrimMonad m) => Double -> Double -> PRNG m -> m Double+laplace_rng mu sd g = MWC.uniformR (0.0, 1.0) g >>= sample+ where sample u = return $ case u < 0.5 of+ True -> mu + sd * log (u + u)+ False -> mu - sd * log (2.0 - u - u) laplaceLogDensity :: Floating a => a -> a -> a -> a laplaceLogDensity mu sd x = - log (2 * sd) - abs (x - mu) / sd laplace :: Double -> Double -> Dist Double laplace mu sd = Dist {logDensity = laplaceLogDensity mu sd . fromLebesgue,- distSample = mapFst Lebesgue . laplace_rng mu sd}+ distSample = (\g -> liftM Lebesgue $ laplace_rng mu sd g)} -- Consider having dirichlet return Vector--- Note: This is acutally symmetric dirichlet-dirichlet_rng :: (RandomGen g) => Int -> Double -> g -> ([Double], g)-dirichlet_rng n' a g' = normalize (gammas g' n')- where gammas g 0 = ([], 0, g)- gammas g n = let (xs, total, g1) = gammas g (n-1)- ( x, g2) = gamma_rng a 1 g1 - in ((x : xs), x+total, g2)- normalize (b, total, h) = (map (/ total) b, h)+-- Note: This is actually symmetric dirichlet+dirichlet_rng :: (PrimMonad m) => Int -> Double -> PRNG m -> m [Double]+dirichlet_rng n' a g' = liftM normalize $ gammas g' n'+ where gammas _ 0 = return ([], 0)+ gammas g n = do (xs, total) <- gammas g (n-1)+ x <- gamma_rng a 1 g+ return ((x : xs), x+total)+ normalize (b, total) = map (/ total) b dirichletLogDensity :: [Double] -> [Double] -> Double dirichletLogDensity a x | all (> 0) x = sum' (zipWith logTerm a x) + logGamma (sum a) where sum' = foldl' (+) 0 logTerm b y = (b-1) * log y - logGamma b dirichletLogDensity _ _ = error "dirichlet: all values must be between 0 and 1"+
Language/Hakaru/ImportanceSampler.hs view
@@ -13,7 +13,8 @@ import Language.Hakaru.Mixture (Prob, empty, point, Mixture(..)) import Language.Hakaru.Sampler (Sampler, deterministic, smap, sbind) -import System.Random+import qualified System.Random.MWC as MWC+import Control.Monad.Primitive import Data.Monoid import Data.Dynamic import System.IO.Unsafe@@ -40,10 +41,9 @@ Just y -> deterministic (point (fromDensity y) density) where density = LF.logToLogFloat $ logDensity dist y Nothing -> error "did not get data from dynamic source"-updateMixture Nothing dist = \g0 -> let (e, g) = distSample dist g0- in (point (fromDensity e) 1, g)+updateMixture Nothing dist = \g -> do e <- distSample dist g+ return $ point (fromDensity e) 1 - conditioned, unconditioned :: Typeable a => Dist a -> Measure a conditioned dist = Measure (\(cond:conds) -> smap (\a->(a,conds)) (updateMixture cond dist))@@ -64,21 +64,31 @@ finish :: Mixture (a, [Cond]) -> Mixture a finish (Mixture m) = Mixture (M.mapKeysMonotonic (\(a,[]) -> a) m) -empiricalMeasure :: (Ord a) => Int -> Measure a -> [Cond] -> IO (Mixture a)-empiricalMeasure !n measure conds = go n empty where- once = getStdRandom (unMeasure measure conds)- go 0 m = return m- go k m = once >>= \result -> go (k - 1) $! mappend m (finish result)+empiricalMeasure :: (PrimMonad m, Ord a) => Int -> Measure a -> [Cond] -> m (Mixture a)+empiricalMeasure !n measure conds = do+ gen <- MWC.create+ go n gen empty+ where once = unMeasure measure conds+ go 0 _ m = return m+ go k g m = once g >>= \result -> go (k - 1) g $! mappend m (finish result) -sample :: (Ord a, Show a) => Measure a -> [Cond] -> IO [(a, Prob)]+sample :: Measure a -> [Cond] -> IO [(a, Prob)] sample measure conds = do- u <- once- let x = mixToTuple (finish u)- xs <- unsafeInterleaveIO $ sample measure conds- return (x : xs)- where once = getStdRandom (unMeasure measure conds)- mixToTuple = head . M.toList . unMixture+ gen <- MWC.create+ unsafeInterleaveIO $ sampleNext gen+ where once = unMeasure measure conds+ mixToTuple = head . M.toList . unMixture+ sampleNext g = do+ u <- once g+ let x = mixToTuple (finish u)+ xs <- unsafeInterleaveIO $ sampleNext g+ return (x : xs)+ -- u <- once gen+ -- let x = mixToTuple (finish u)+ -- xs <- unsafeInterleaveIO $ sample measure conds gen+ -- return (x : xs)+ -- where once = unMeasure measure conds+ -- mixToTuple = head . M.toList . unMixture logit :: Floating a => a -> a logit !x = 1 / (1 + exp (- x))-
Language/Hakaru/Metropolis.hs view
@@ -5,14 +5,18 @@ module Language.Hakaru.Metropolis where -import System.Random (RandomGen, StdGen, randomR, getStdGen)-+import qualified System.Random.MWC as MWC+import Control.Monad+import Control.Monad.Primitive import Data.Dynamic import Data.Maybe+import Control.Applicative import qualified Data.Map.Strict as M import Language.Hakaru.Types +import System.IO.Unsafe+ {- Shortcomings of this implementation@@ -25,7 +29,6 @@ type DistVal = Dynamic --- and what does XRP stand for? data XRP where XRP :: Typeable e => (Density e, Dist e) -> XRP @@ -36,6 +39,11 @@ type Observed = Bool type LL = LogLikelihood +-- The first component is the LogLikelihood of the trace+-- The second is the LogLikelihood of the newly introduced+-- choices. These are used to compute the acceptance ratio+type LL2 = (LL,LL)+ type Subloc = Int type Name = [Subloc] data DBEntry = DBEntry {@@ -45,20 +53,21 @@ observed :: Observed } type Database = M.Map Name DBEntry -data SamplerState g where+data SamplerState where S :: { ldb :: Database, -- ldb = local database -- (total likelihood, total likelihood of XRPs newly introduced)- llh2 :: {-# UNPACK #-} !(LL, LL),- cnds :: [Cond], -- conditions left to process- seed :: g } -> SamplerState g+ llh2 :: {-# UNPACK #-} !LL2,+ cnds :: [Cond] -- conditions left to process+ } -> SamplerState -type Sampler a = forall g. (RandomGen g) => SamplerState g -> (a, SamplerState g)+type Sampler a = PrimMonad m => SamplerState -> PRNG m -> m (a, SamplerState) sreturn :: a -> Sampler a-sreturn x s = (x, s)+sreturn x s _ = return (x, s) sbind :: Sampler a -> (a -> Sampler b) -> Sampler b-sbind s k = \ st -> let (v, s') = s st in k v s'+sbind s k = \ st g -> do (v, s') <- s st g+ k v s' g smap :: (a -> b) -> Sampler a -> Sampler b smap f s = sbind s (\a -> sreturn (f a))@@ -70,51 +79,41 @@ return_ x = Measure $ \ _ -> sreturn x updateXRP :: Typeable a => Name -> Cond -> Dist a -> Sampler a-updateXRP n obs dist' s@(S {ldb = db, seed = g}) =+updateXRP n obs dist' s@(S {ldb = db}) g = do case M.lookup n db of- Just (DBEntry xd lb _ ob) ->- let Just (xb, dist) = unXRP xd- (x,_) = case obs of- Just yd ->- let Just y = fromDynamic yd- in (y, logDensity dist y)- Nothing -> (xb, lb)- l' = logDensity dist' x- d1 = M.insert n (DBEntry (XRP (x,dist)) l' True ob) db- in (fromDensity x,- s {ldb = d1,- llh2 = updateLogLikelihood l' 0 s,- seed = g})+ Just (DBEntry xd _ _ ob) ->+ do let Just (x, _) = unXRP xd+ l' = logDensity dist' x+ d1 = M.insert n (DBEntry (XRP (x,dist')) l' True ob) db+ return (fromDensity x,+ s {ldb = d1,+ llh2 = updateLogLikelihood (l',0) (llh2 s)}) Nothing ->- let (xnew2, l, g2) = case obs of- Just xdnew ->- let Just xnew = fromDynamic xdnew- in (xnew, logDensity dist' xnew, g)- Nothing ->- let (xnew, g1) = distSample dist' g- in (xnew, logDensity dist' xnew, g1)- d1 = M.insert n (DBEntry (XRP (xnew2, dist')) l True (isJust obs)) db- in (fromDensity xnew2,- s {ldb = d1,- llh2 = updateLogLikelihood l l s,- seed = g2})+ do (xnew2, l) <- case obs of+ Just xdnew ->+ do let Just xnew = fromDynamic xdnew+ return $ (xnew, logDensity dist' xnew)+ Nothing ->+ do xnew <- distSample dist' g+ return (xnew, logDensity dist' xnew)+ let d1 = M.insert n (DBEntry (XRP (xnew2, dist')) l True (isJust obs)) db+ return (fromDensity xnew2,+ s {ldb = d1,+ llh2 = updateLogLikelihood (l,l) (llh2 s)}) -updateLogLikelihood :: RandomGen g => - LL -> LL -> SamplerState g ->- (LL, LL)-updateLogLikelihood llTotal llFresh s =- let (l,lf) = llh2 s in (llTotal+l, llFresh+lf)+updateLogLikelihood :: LL2 -> LL2 -> LL2+updateLogLikelihood (llTotal,llFresh) (l,lf) = (llTotal+l, llFresh+lf) factor :: LL -> Measure ()-factor l = Measure $ \ _ -> \ s ->- let (llTotal, llFresh) = llh2 s- in ((), s {llh2 = (llTotal + l, llFresh)})+factor l = Measure $ \ _ -> \ s _ ->+ do let (llTotal, llFresh) = llh2 s+ return ((), s {llh2 = (llTotal + l, llFresh)}) condition :: Eq b => Measure (a, b) -> b -> Measure a condition (Measure m) b' = Measure $ \ n ->- let comp a b s | a /= b = s {llh2 = (log 0, 0)}- comp _ _ s = s- in sbind (m n) (\ (a, b) s -> (a, comp b b' s))+ do let comp a b s | a /= b = s {llh2 = (log 0, 0)}+ comp _ _ s = s+ sbind (m n) (\ (a, b) s _ -> return (a, comp b b' s)) bind :: Measure a -> (a -> Measure b) -> Measure b bind (Measure m) cont = Measure $ \ n ->@@ -133,68 +132,81 @@ return = return_ (>>=) = bind +instance Functor Measure where+ fmap f (Measure x) = Measure $ \n -> smap f (x n)++instance Applicative Measure where+ pure = return_+ (<*>) = app++sapp :: (Sampler (a -> b)) -> Sampler a -> Sampler b+sapp f s = \st g -> do (vf, s') <- f st g+ (vs, s'') <- s s' g+ sreturn (vf vs) s'' g++app :: Measure (a -> b) -> Measure a -> Measure b+app (Measure f) (Measure a) = Measure $ \n -> sapp (f n) (a n)+ run :: Measure a -> [Cond] -> IO (a, Database, LL) run (Measure prog) cds = do- g <- getStdGen- let (v, S d ll [] _) = (prog [0]) (S M.empty (0,0) cds g)+ g <- MWC.create+ (v, S d ll []) <- (prog [0]) (S M.empty (0,0) cds) g return (v, d, fst ll) -traceUpdate :: RandomGen g => Measure a -> Database -> [Cond] -> g- -> (a, Database, LL, LL, LL, g)+traceUpdate :: PrimMonad m => Measure a -> Database -> [Cond] -> PRNG m+ -> m (a, Database, LL, LL, LL) traceUpdate (Measure prog) d cds g = do -- let d1 = M.map (\ (x, l, _, ob) -> (x, l, False, ob)) d let d1 = M.map (\ s -> s { vis = False }) d- let (v, S d2 (llTotal, llFresh) [] g1) = (prog [0]) (S d1 (0,0) cds g)+ (v, S d2 (llTotal, llFresh) []) <- (prog [0]) (S d1 (0,0) cds) g let (d3, stale_d) = M.partition vis d2 let llStale = M.foldl' (\ llStale' s -> llStale' + llhd s) 0 stale_d- (v, d3, llTotal, llFresh, llStale, g1)+ return (v, d3, llTotal, llFresh, llStale) initialStep :: Measure a -> [Cond] ->- IO (a, Database,- LL, LL, LL, StdGen)-initialStep prog cds = do- g <- getStdGen- return $ traceUpdate prog M.empty cds g+ PRNG IO -> IO (a, Database, LL, LL, LL)+initialStep prog cds g = traceUpdate prog M.empty cds g -- TODO: Make a way of passing user-provided proposal distributions-resample :: RandomGen g => Name -> Database -> Observed -> XRP -> g ->- (Database, LL, LL, LL, g)+resample :: PrimMonad m => Name -> Database -> Observed -> XRP -> PRNG m ->+ m (Database, LL, LL, LL) resample name db ob (XRP (x, dist)) g =- let (x', g1) = distSample dist g- fwd = logDensity dist x'- rvs = logDensity dist x- l' = fwd- newEntry = DBEntry (XRP (x', dist)) l' True ob- db' = M.insert name newEntry db- in (db', l', fwd, rvs, g1)+ do x' <- distSample dist g+ let fwd = logDensity dist x'+ rvs = logDensity dist x+ l' = fwd+ newEntry = DBEntry (XRP (x', dist)) l' True ob+ db' = M.insert name newEntry db+ return (db', l', fwd, rvs) -transition :: (Typeable a, RandomGen g) => Measure a -> [Cond]- -> a -> Database -> LL -> g -> [a]+transition :: (Typeable a) => Measure a -> [Cond]+ -> a -> Database -> LL -> PRNG IO -> IO [a] transition prog cds v db ll g =- let dbSize = M.size db- -- choose an unconditioned choice- (_, uncondDb) = M.partition observed db- (choice, g1) = randomR (0, (M.size uncondDb) -1) g- (name, (DBEntry xd _ _ ob)) = M.elemAt choice uncondDb- (db', _, fwd, rvs, g2) = resample name db ob xd g1- (v', db2, llTotal, llFresh, llStale, g3) = traceUpdate prog db' cds g2- a = llTotal - ll- + rvs - fwd- + log (fromIntegral dbSize) - log (fromIntegral $ M.size db2)- + llStale - llFresh- (u, g4) = randomR (0 :: Double, 1) g3 in-- if (log u < a) then- v' : (transition prog cds v' db2 llTotal g4)- else- v : (transition prog cds v db ll g4)+ do let dbSize = M.size db+ -- choose an unconditioned choice+ (_, uncondDb) = M.partition observed db+ choice <- MWC.uniformR (0, (M.size uncondDb) -1) g+ let (name, (DBEntry xd _ _ ob)) = M.elemAt choice uncondDb+ (db', _, fwd, rvs) <- resample name db ob xd g+ (v', db2, llTotal, llFresh, llStale) <- traceUpdate prog db' cds g+ let a = llTotal - ll+ + rvs - fwd+ + log (fromIntegral dbSize) - log (fromIntegral $ M.size db2)+ + llStale - llFresh+ u <- MWC.uniformR (0 :: Double, 1) g+ if (log u < a) then+ liftM ((:) v') $ unsafeInterleaveIO (transition prog cds v' db2 llTotal g)+ else+ liftM ((:) v) $ unsafeInterleaveIO (transition prog cds v db ll g) mcmc :: Typeable a => Measure a -> [Cond] -> IO [a] mcmc prog cds = do- (v, d, llTotal, _, _, g) <- initialStep prog cds- return $ transition prog cds v d llTotal g+ g <- MWC.create+ (v, d, llTotal, _, _) <- initialStep prog cds g+ transition prog cds v d llTotal g sample :: Typeable a => Measure a -> [Cond] -> IO [(a, Double)] sample prog cds = do - (v, d, llTotal, _, _, g) <- initialStep prog cds- return $ map (\ x -> (x,1)) (transition prog cds v d llTotal g)+ g <- MWC.create+ (v, d, llTotal, _, _) <- initialStep prog cds g+ (transition prog cds v d llTotal g) >>= return . map (\ x -> (x,1))
Language/Hakaru/Sampler.hs view
@@ -5,22 +5,24 @@ import Language.Hakaru.Mixture (Mixture, mnull, empty, scale, point) import Language.Hakaru.Distribution (choose)-import System.Random (RandomGen)+import Language.Hakaru.Types+import Control.Monad.Primitive -- Sampling procedures that produce one sample -type Sampler a = forall g. (RandomGen g) => g -> (Mixture a, g)+type Sampler a = PrimMonad m => PRNG m -> m (Mixture a) deterministic :: Mixture a -> Sampler a-deterministic m g = (m, g)+deterministic m _ = return m sbind :: Sampler a -> (a -> Sampler b) -> Sampler b-sbind s k g0 =- case s g0 of { (m1, g1) ->- if mnull m1 then (empty, g1) else- case choose m1 g1 of { (a, v, g2) ->- case k a g2 of { (m2, g) ->- (scale v m2, g) } } }+sbind s k g = do+ m1 <- s g+ if mnull m1 then + return empty+ else do (a, v) <- choose m1 g+ m2 <- k a g+ return (scale v m2) smap :: (a -> b) -> Sampler a -> Sampler b smap f s = sbind s (\a -> deterministic (point (f a) 1))
Language/Hakaru/Types.hs view
@@ -1,11 +1,14 @@-{-# LANGUAGE RankNTypes, BangPatterns, DeriveDataTypeable, StandaloneDeriving #-}+{-# LANGUAGE RankNTypes, DeriveDataTypeable, StandaloneDeriving #-} {-# OPTIONS -W #-} module Language.Hakaru.Types where import Data.Dynamic-import System.Random+import Control.Monad.Primitive+import qualified System.Random.MWC as MWC +type PRNG m = MWC.Gen (PrimState m)+ -- Basic types for conditioning and conditioned sampler data Density a = Lebesgue !a | Discrete !a deriving Typeable type Cond = Maybe Dynamic@@ -24,6 +27,5 @@ type LogLikelihood = Double data Dist a = Dist {logDensity :: Density a -> LogLikelihood,- distSample :: forall g.- RandomGen g => g -> (Density a, g)}+ distSample :: (PrimMonad m) => PRNG m -> m (Density a)} deriving instance Typeable1 Dist
Language/Hakaru/Util/Coda.hs view
@@ -10,3 +10,11 @@ vec = V.fromList samples cov = autocovariance vec rho = G.map (/ G.head cov) cov++meanVariance :: Fractional a => [a] -> (a,a)+meanVariance lst = (av,sigma2)+ where+ n = fromIntegral $ length lst+ av = sum lst / n+ sigma2 = (foldr (\x acc -> sqr (x - av) + acc) 0 lst) / (n - 1)+ sqr x = x * x
Language/Hakaru/Util/Extras.hs view
@@ -2,13 +2,13 @@ Functions on lists and sequences. Some of the functions follow the style of Data.Random.Extras (from the random-extras package), but are written for use with- PRNGs from System.Random rather than from the random-fu package.+ PRNGs from the "mwc-random" package rather than from the "random-fu" package. -} module Language.Hakaru.Util.Extras where import qualified Data.Sequence as S-import System.Random+import qualified System.Random.MWC as MWC import Data.Maybe import qualified Data.Foldable as F @@ -18,10 +18,9 @@ where (a, r) = S.splitAt i s (b S.:< c) = S.viewl r -randomExtract :: S.Seq a -> IO (Maybe (S.Seq a, a))-randomExtract s = do- g <- newStdGen- let (i,_) = randomR (0, S.length s - 1) g+randomExtract :: S.Seq a -> MWC.GenIO -> IO (Maybe (S.Seq a, a))+randomExtract s g = do+ i <- MWC.uniformR (0, S.length s - 1) g return $ extract s i {-| @@ -30,15 +29,17 @@ -} randomElems :: Ord a => S.Seq a -> Int -> IO (S.Seq a)-randomElems = randomElemsTR S.empty+randomElems s n = do + g <- MWC.create+ randomElemsTR S.empty s g n -randomElemsTR :: Ord a => S.Seq a -> S.Seq a -> Int -> IO (S.Seq a)-randomElemsTR ixs s n+randomElemsTR :: Ord a => S.Seq a -> S.Seq a -> MWC.GenIO -> Int -> IO (S.Seq a)+randomElemsTR ixs s g n | n == S.length s = return $ S.unstableSort s- | n == 1 = do (_,i) <- fmap fromJust (randomExtract s)+ | n == 1 = do (_,i) <- fmap fromJust (randomExtract s g) return.S.unstableSort $ i S.<| ixs- | otherwise = do (s',i) <- fmap fromJust (randomExtract s)- (randomElemsTR $! (i S.<| ixs)) s' (n-1)+ | otherwise = do (s',i) <- fmap fromJust (randomExtract s g)+ (randomElemsTR $! (i S.<| ixs)) s' g (n-1) {-| Chop a sequence at the given indices.
Language/Hakaru/Util/Visual.hs view
@@ -9,7 +9,7 @@ import Data.List import qualified Data.Text as T import qualified Data.ByteString.Lazy.Char8 as B-import qualified Data.ByteString.Char8 as BS+--import qualified Data.ByteString.Char8 as BS plot :: Show a => [a] -> String -> IO () plot samples filename = do@@ -34,7 +34,7 @@ where total = "total_samples" .= n current_sample = "current_sample" .= cur- chunk = object (zipWith (\ name s -> T.pack name .= s)+ chunk = object (zipWith (\ name' s -> T.pack name' .= s) name (transpose $ take c samples)) batch = B.unpack $ encode
Tests/Distribution.hs view
@@ -16,9 +16,8 @@ fromDiscreteToNum = fromIntegral . fromEnum . fromDiscrete sq x = x * x -almostEqual :: (Fractional a, Ord a) => a -> a -> a -> Bool-almostEqual tol x y | abs (x - y) < tol = True-almostEqual tol x y = (abs $ (x - y) / (x + y)) < tol+almostEqual :: (Num a, Ord a) => a -> a -> a -> Bool+almostEqual tol x y = abs (x - y) < tol quickArg :: IO () quickArg = quickCheckWith stdArgs {maxSuccess = 2000} (\ x -> almostEqual tol x x)@@ -26,11 +25,7 @@ tol = 1e-5 qtest = [testProperty "checking beta" $ QM.monadicIO betaTest,- testProperty "checking bern" $ QM.monadicIO bernTest,- testProperty "checking gamma" $ QM.monadicIO gammaTest,- testProperty "checking normal" $ QM.monadicIO normalTest,- testProperty "checking laplace" $ QM.monadicIO laplaceTest,- testProperty "checking poisson" $ QM.monadicIO poissonTest]+ testProperty "checking bern" $ QM.monadicIO bernTest] betaTest = do Positive a <- QM.pick arbitrary@@ -54,53 +49,3 @@ where tol = 1e-1 mu p = p var p = p*(1-p)--poissonTest = do- lam <- QM.pick $ choose (1, 10)- g <- QM.run $ MWC.create- samples <- QM.run $ replicateM 1000 $ distSample (poisson lam) g- let (mean, variance) = meanVariance (map (fromIntegral . fromDiscrete) samples)- QM.assert $ (almostEqual tol mean (mu lam)) && - (almostEqual tol variance (var lam))- where tol = 1e-1- mu lam = lam- var lam = lam--normalTest = do- mu <- QM.pick arbitrary- sd <- QM.pick $ choose (1, 10)- g <- QM.run $ MWC.create- let nsamples = floor (1000 * sd) -- larger standard deviations need more samples- -- to be shown as correct- samples <- QM.run $ replicateM nsamples $ distSample (normal mu sd) g- let (mean, variance) = meanVariance (map fromLebesgue samples)- QM.assert $ (almostEqual tol mean mu ) && - (almostEqual tol variance (var sd))- where tol = 1e-1- var sd = sq sd--laplaceTest = do- mu <- QM.pick arbitrary- sd <- QM.pick $ choose (1, 10)- g <- QM.run $ MWC.create- let nsamples = floor (1000 * sd) -- larger standard deviations need more samples- -- to be shown as correct- samples <- QM.run $ replicateM nsamples $ distSample (laplace mu sd) g- let (mean, variance) = meanVariance (map fromLebesgue samples)- QM.assert $ (almostEqual tol mean mu ) && - (almostEqual tol variance (var sd))- where tol = 1e-1- var sd = 2*(sq sd)--gammaTest = do- a <- QM.pick $ choose (1, 10)- b <- QM.pick $ choose (1, 10)- g <- QM.run $ MWC.create- samples <- QM.run $ replicateM 1000 $ distSample (gamma a b) g- let (mean, variance) = meanVariance (map fromLebesgue samples)- QM.assert $ (almostEqual tol mean (mu a b)) && - (almostEqual tol variance (var a b))- where tol = 1e-1- mu a b = a * b- var a b = a * (b * b)-
Tests/ImportanceSampler.hs view
@@ -7,6 +7,7 @@ import Language.Hakaru.Lambda import Language.Hakaru.Distribution import Language.Hakaru.ImportanceSampler+import qualified System.Random.MWC as MWC -- import Test.QuickCheck.Monadic import Tests.Models@@ -14,7 +15,7 @@ -- Some test programs in our language test_mixture :: IO ()-test_mixture = sample prog_mixture conds >>=+test_mixture = MWC.create >>= sample prog_mixture conds >>= print . take 10 >> putChar '\n' >> empiricalMeasure 1000 prog_mixture conds >>=@@ -38,7 +39,7 @@ return s2 test_dbn :: IO ()-test_dbn = sample prog_dbn conds >>=+test_dbn = MWC.create >>= sample prog_dbn conds >>= print . take 10 >> putChar '\n' >> empiricalMeasure 1000 prog_dbn conds >>=@@ -59,7 +60,7 @@ else return s test_hmm :: IO ()-test_hmm = sample (prog_hmm 2) conds >>=+test_hmm = MWC.create >>= sample (prog_hmm 2) conds >>= print . take 10 >> putChar '\n' >> empiricalMeasure 1000 (prog_hmm 2) conds >>=@@ -82,7 +83,7 @@ return (z4, z3) test_carRoadModel :: IO ()-test_carRoadModel = sample prog_carRoadModel conds >>=+test_carRoadModel = MWC.create >>= sample prog_carRoadModel conds >>= print . take 10 >> putChar '\n' >> empiricalMeasure 1000 prog_carRoadModel conds >>=@@ -103,7 +104,7 @@ return rain test_categorical :: IO ()-test_categorical = sample prog_categorical conds >>=+test_categorical = MWC.create >>= sample prog_categorical conds >>= print . take 10 >> putChar '\n' >> empiricalMeasure 1000 prog_categorical conds >>=
hakaru.cabal view
@@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: hakaru-version: 0.1.3+version: 0.1.4 synopsis: A probabilistic programming embedded DSL -- description: homepage: http://indiana.edu/~ppaml/