hakaru 0.1.2 → 0.1.3
raw patch · 30 files changed
+1308/−1072 lines, 30 filesdep +Cabaldep +HUnitdep +QuickCheckdep ~basedep ~statistics
Dependencies added: Cabal, HUnit, QuickCheck, array, criterion, deepseq, directory, ghc-prim, hakaru, integration, monad-loops, mwc-random, parallel, primitive, test-framework, test-framework-hunit, test-framework-quickcheck2, transformers
Dependency ranges changed: base, statistics
Files
- Bench/Bench.hs +20/−0
- Examples/Tests.hs +0/−106
- LICENSE +30/−0
- Language/Hakaru/Arrow.hs +6/−0
- Language/Hakaru/Distribution.hs +228/−0
- Language/Hakaru/ImportanceSampler.hs +37/−127
- Language/Hakaru/Lambda.hs +22/−0
- Language/Hakaru/Metropolis.hs +114/−208
- Language/Hakaru/Mixture.hs +61/−0
- Language/Hakaru/Sampler.hs +26/−0
- Language/Hakaru/Symbolic.hs +161/−54
- Language/Hakaru/Syntax.hs +0/−185
- Language/Hakaru/Types.hs +29/−0
- Language/Hakaru/Util/Coda.hs +12/−0
- Language/Hakaru/Util/Csv.hs +40/−0
- Language/Hakaru/Util/Extras.hs +78/−0
- Language/Hakaru/Util/Visual.hs +43/−0
- Mixture.hs +0/−61
- RandomChoice.hs +0/−145
- Sampler.hs +0/−26
- Tests/Distribution.hs +106/−0
- Tests/ImportanceSampler.hs +118/−0
- Tests/Metropolis.hs +45/−0
- Tests/Models.hs +26/−0
- Tests/Tests.hs +24/−0
- Types.hs +0/−13
- Util/Coda.hs +0/−12
- Util/Extras.hs +0/−84
- Visual.hs +0/−43
- hakaru.cabal +82/−8
+ Bench/Bench.hs view
@@ -0,0 +1,20 @@+import Criterion.Main++import Control.Monad++import Language.Hakaru.Distribution+import qualified Language.Hakaru.Metropolis as MH+import qualified Language.Hakaru.ImportanceSampler as IS++giveLast :: IO [(a,b)] -> IO (a,b)+giveLast samples = do s <- samples+ return $ last (take 100 s)++main = defaultMain [+ bcompare [+ bench "is normal 10" $ whnfIO $ giveLast (IS.sample (replicateM 10 (IS.unconditioned (normal 0 10))) [])+ , bench "is normal 20" $ whnfIO $ giveLast (IS.sample (replicateM 20 (IS.unconditioned (normal 0 10))) [])+ , bench "mh normal 10" $ whnfIO $ giveLast (MH.sample (replicateM 10 (MH.unconditioned (normal 0 10))) [])+ , bench "mh normal 20" $ whnfIO $ giveLast (MH.sample (replicateM 20 (MH.unconditioned (normal 0 10))) [])+ ]+ ]
− Examples/Tests.hs
@@ -1,106 +0,0 @@-{-# LANGUAGE RankNTypes, NoMonomorphismRestriction, BangPatterns #-}--module Examples.Tests where--import Types-import Data.Dynamic-import Language.Hakaru.ImportanceSampler---- Some example/test programs in our language-test :: Measure Bool-test = do- c <- unconditioned (bern 0.5)- _ <- conditioned (ifThenElse c (normal (lit (1 :: Double)) (lit 1))- (uniformC (lit 0) (lit 3)))- return c--test_dup :: Measure (Bool, Bool)-test_dup = do- let c = unconditioned (bern 0.5)- x <- c- y <- c- return (x,y)--test_dbn :: Measure Bool-test_dbn = do- s0 <- unconditioned (bern 0.75)- s1 <- unconditioned (if s0 then bern 0.75 else bern 0.25)- _ <- conditioned (if s1 then bern 0.9 else bern 0.1)- s2 <- unconditioned (if s1 then bern 0.75 else bern 0.25)- _ <- conditioned (if s2 then bern 0.9 else bern 0.1)- return s2--test_hmm :: Integer -> Measure Bool-test_hmm n = do- s <- unconditioned (bern 0.75) - loop_hmm n s--loop_hmm :: Integer -> (Bool -> Measure Bool)-loop_hmm !numLoops s = do- _ <- conditioned (if s then bern 0.9 else bern 0.1)- u <- unconditioned (if s then bern 0.75 else bern 0.25)- if (numLoops > 1) then loop_hmm (numLoops - 1) u - else return s--test_carRoadModel :: Measure (Double, Double)-test_carRoadModel = do- speed <- unconditioned (uniformC (lit (5::Double)) (lit (15::Double)))- let z0 = lit 0 - _ <- conditioned (normal (z0 :: Double) (lit 1))- z1 <- unconditioned (normal (z0 + speed) (lit 1))- _ <- conditioned (normal z1 (lit 1))- z2 <- unconditioned (normal (z1 + speed) (lit 1)) - _ <- conditioned (normal z2 (lit 1))- z3 <- unconditioned (normal (z2 + speed) (lit 1)) - _ <- conditioned (normal z3 (lit 1))- z4 <- unconditioned (normal (z3 + speed) (lit 1)) - return (z4, z3)--test_categorical :: Measure Bool-test_categorical = do - rain <- unconditioned (categorical [(lit True, 0.2), (lit False, 0.8)]) - sprinkler <- unconditioned (if rain then bern 0.01 else bern 0.4)- _ <- conditioned (if rain then (if sprinkler then bern 0.99 else bern 0.8)- else (if sprinkler then bern 0.9 else bern 0.1))- return rain---- printing test results-main :: IO ()-main = sample_ 3 test conds >>- putChar '\n' >>- sample 1000 test conds >>=- print- where conds = [Lebesgue (toDyn (2 :: Double))]--main_dbn :: IO ()-main_dbn = sample_ 10 test_dbn conds >>- putChar '\n' >>- sample 1000 test_dbn conds >>=- print - where conds = [Discrete (toDyn (True :: Bool)),- Discrete (toDyn (True :: Bool))]--main_hmm :: IO ()-main_hmm = sample_ 10 (test_hmm 2) conds >>- putChar '\n' >>- sample 1000 (test_hmm 2) conds >>=- print - where conds = [Discrete (toDyn (True :: Bool)),- Discrete (toDyn (True :: Bool))]--main_carRoadModel :: IO ()-main_carRoadModel = sample_ 10 test_carRoadModel conds >>- putChar '\n' >>- sample 1000 test_carRoadModel conds >>=- print - where conds = [Lebesgue (toDyn (0 :: Double)),- Lebesgue (toDyn (11 :: Double)), - Lebesgue (toDyn (19 :: Double)),- Lebesgue (toDyn (33 :: Double))]--main_categorical :: IO ()-main_categorical = sample_ 10 test_categorical conds >>- putChar '\n' >>- sample 1000 test_categorical conds >>=- print - where conds = [Discrete (toDyn (True :: Bool))]
LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, The Hakaru Team++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of The Hakaru Team nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Language/Hakaru/Arrow.hs view
@@ -0,0 +1,6 @@+{-# LANGUAGE TypeOperators #-}+module Language.Hakaru.Arrow where++import Language.Hakaru.Types (Dist)++type a ~~> b = a -> Dist b
+ Language/Hakaru/Distribution.hs view
@@ -0,0 +1,228 @@+{-# LANGUAGE RankNTypes, BangPatterns, GADTs #-}+{-# OPTIONS -Wall #-}++module Language.Hakaru.Distribution where++import System.Random+import Language.Hakaru.Mixture+import Language.Hakaru.Types+import Data.Ix+import Data.Maybe (fromMaybe)+import Data.List (findIndex, foldl')+import Numeric.SpecFunctions+import qualified Data.Map.Strict as M+import qualified Data.Number.LogFloat as LF++mapFst :: (t -> s) -> (t, u) -> (s, u)+mapFst f (a,b) = (f a, b)++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))}++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'))}++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)}++uniformD :: (Ix a, Random 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)}++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++choose :: (RandomGen g) => Mixture k -> g -> (k, Prob, g)+choose (Mixture m) g0 =+ 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+ 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)++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)++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 _ _ _ = error "normal: invalid parameters"++normalLogDensity :: Floating a => a -> a -> a -> a+normalLogDensity mu sd x = (-tau * square (x - mu)+ + log (tau / pi / 2)) / 2+ where square y = y * y+ tau = 1 / square sd++normal :: Double -> Double -> Dist Double +normal mu sd = Dist {logDensity = normalLogDensity mu sd . fromLebesgue,+ distSample = mapFst Lebesgue . normal_rng mu sd}++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+ elem' = fst $ head $ filter (\(_,p0) -> p <= p0) sumList+ sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list+ total = sum $ map snd list++categorical :: Eq a => [(a,Double)] -> Dist a+categorical list = Dist {logDensity = categoricalLogDensity list . fromDiscrete,+ distSample = mapFst Discrete . categoricalSample list}++lnFact :: Integer -> Double+lnFact = logFactorial++-- Makes use of Atkinson's algorithm as described in:+-- Monte Carlo Statistical Methods pg. 55+--+-- 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+ where smu = sqrt lambda+ b = 0.931 + 2.53*smu+ a = -0.059 + 0.02483*b+ vr = 0.9277 - 3.6224/(b - 2)+ 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++ accept_region :: Double -> Double -> Integer -> Bool+ accept_region us v k = log (v * arep / (a/(us*us)+b)) <=+ -lambda + (fromIntegral k)*lnlam - lnFact k++poisson :: Double -> Dist Integer+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}++-- 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 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 = + 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+ 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++gammaLogDensity :: Double -> Double -> Double -> Double+gammaLogDensity shape scl x | x>= 0 && shape > 0 && scl > 0 =+ scl * log shape - scl * x + (shape - 1) * log x - logGamma shape+gammaLogDensity _ _ _ = log 0++gamma :: Double -> Double -> Dist Double+gamma shape scl = Dist {logDensity = gammaLogDensity shape scl . fromLebesgue,+ distSample = mapFst Lebesgue . gamma_rng shape scl}++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)+ 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)++betaLogDensity :: Double -> Double -> Double -> Double+betaLogDensity _ _ x | x < 0 || x > 1 = error "beta: value must be between 0 and 1"+betaLogDensity a b _ | a <= 0 || b <= 0 = error "beta: parameters must be positve" +betaLogDensity a b x = (logGamma (a + b)+ - logGamma a+ - logGamma b+ + (a - 1) * log x+ + (b - 1) * log (1 - x))++beta :: Double -> Double -> Dist Double+beta a b = Dist {logDensity = betaLogDensity a b . fromLebesgue,+ distSample = mapFst Lebesgue . beta_rng a b}++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)++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}++-- 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)++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
@@ -9,111 +9,18 @@ -- inputs. In exchange, we get to make Measure an instance of Monad, and we -- can express models whose number of observations is unknown at compile time. -import Types (Cond(..), CSampler(CSampler))-import RandomChoice (normal_rng, chooseIndex)-import Mixture (Prob, empty, point, Mixture(..))-import Sampler (Sampler, deterministic, smap, sbind)+import Language.Hakaru.Types+import Language.Hakaru.Mixture (Prob, empty, point, Mixture(..))+import Language.Hakaru.Sampler (Sampler, deterministic, smap, sbind) import System.Random import Data.Monoid-import Data.Ix import Data.Dynamic-import Data.List-import Control.Monad+import System.IO.Unsafe import qualified Data.Map.Strict as M import qualified Data.Number.LogFloat as LF -dirac :: (Eq a, Typeable a) => a -> CSampler a-dirac x = CSampler c where- c Unconditioned = deterministic (point x 1)- c (Discrete y) = case fromDynamic y of- Just y -> deterministic (if x == y then point x 1 else empty)- Nothing -> error "dirac: did not get data from dynamic source"- c _ = error "dirac: got a non-discrete sampler"--bern :: Double -> CSampler Bool-bern theta | 0 <= theta && theta <= 1 = CSampler c where- c Unconditioned = \g0 -> case randomR (0, 1) g0 of- (x, g) -> (point (x <= theta) 1, g)- c (Discrete y) = case fromDynamic y of- Just y -> deterministic (point y (LF.logFloat (if y then theta else 1 - theta)))- Nothing -> error "bern: did not get data from dynamic source"- c _ = error "bern: got a non-discrete sampler"-bern theta = error ("bernoulli: invalid parameter " ++ show theta)--uniformC :: (Fractional a, Real a, Random a, Typeable a) => a -> a -> CSampler a-uniformC lo hi | lo < hi = CSampler c where- c Unconditioned = \g0 -> case randomR (lo,hi) g0 of- (x, g) -> (point x 1, g)- c (Lebesgue y) = case fromDynamic y of- Just y -> deterministic (if lo < y && y < hi then point y density else empty)- Nothing -> error "uniformC: did not get data from dynamic source"- c _ = error "uniformC: got a discrete sampler"- density = fromRational (toRational (recip (hi - lo)))-uniformC _ _ = error "uniformC: invalid parameters"--uniformD :: (Ix a, Random a, Typeable a) => a -> a -> CSampler a-uniformD lo hi | lo <= hi = CSampler c where- c Unconditioned = \g0 -> case randomR (lo,hi) g0 of- (x, g) -> (point x 1, g)- c (Discrete y) = case fromDynamic y of- Just y -> deterministic (if lo <= y && y <= hi then point y density else empty)- Nothing -> error "uniformD: did not get data from dynamic source"- c _ = error "uniformD: got a non-discrete sampler"- density = recip (fromInteger (toInteger (rangeSize (lo,hi))))-uniformD _ _ = error "uniformD: invalid parameters"--poisson :: (Integral a, Typeable a) => Double -> CSampler a-poisson !l | 0 <= l = CSampler c where- c Unconditioned = \g0 ->- let probs = exp (-l) : zipWith (\k p -> p * l / k) [1..] probs- (k, g) = chooseIndex probs g0- in (point (fromInteger (toInteger k)) 1, g)- c (Discrete k) = case fromDynamic k of- Just k ->- deterministic- (if 0 <= k then point k (LF.logToLogFloat (-l)- * LF.logFloat l ^ k- / product (map fromIntegral [1..k]))- else empty)- Nothing -> error "poisson: did not get data from dynamic source"- c _ = error "poisson: got a non-discrete sampler"-poisson _ = error "poisson: invalid parameter"--normal :: (Real a, Floating a, Random a, Typeable a) => a -> a -> CSampler a-normal !mean !std | std > 0 = CSampler c where- c Unconditioned = \g0 -> let (x, g) = normal_rng mean std g0- in (point (mean + std * x) 1, g)- c (Lebesgue y) = case fromDynamic y of- Just y ->- let density = exp (square ((y - mean) / std) / (-2)) / std / sqrt (2 * pi)- square y = y * y- in deterministic (point y (fromRational (toRational density))) -- TODO: use log-density and LogFloat directly- Nothing -> error "normal: did not get data from dynamic source"- c _ = error "normal: got a discrete sampler"-normal _ _ = error "normal: invalid parameters"--categorical :: (Typeable a, Eq a) => [(a, Prob)] -> CSampler a-categorical list = CSampler c where- peak :: LF.LogFloat- peak = maximum (map snd list)- total :: Double- (total, list') = mapAccumL f 0 list- where f acc (a,b) = (acc', (a, (b', acc')))- where b' = b/peak- acc' :: Double- acc' = acc + LF.fromLogFloat b'- c Unconditioned =- \g0 -> let (p, g) = randomR (0, total) g0- (elem, _) : _ = filter (\(_,(_,p0)) -> p <= p0) list' in- (point elem 1, g)- c (Discrete y) = case fromDynamic y of- Just y -> deterministic (maybe empty (point y . (/ LF.logFloat total) . fst)- (lookup y list'))- Nothing -> error "categorical: did not get data from dynamic source"- c _ = error "categorical: got a non-discrete sampler"- -- Conditioned sampling newtype Measure a = Measure { unMeasure :: [Cond] -> Sampler (a, [Cond]) } @@ -121,53 +28,56 @@ bind measure continuation = Measure (\conds -> sbind (unMeasure measure conds)- (\(a,conds) -> unMeasure (continuation a) conds))+ (\(a,cds) -> unMeasure (continuation a) cds)) instance Monad Measure where return x = Measure (\conds -> deterministic (point (x,conds) 1)) (>>=) = bind -conditioned, unconditioned :: CSampler a -> Measure a-conditioned (CSampler f) = Measure (\(cond:conds) -> smap (\a->(a,conds)) (f cond ))-unconditioned (CSampler f) = Measure (\ conds -> smap (\a->(a,conds)) (f Unconditioned))+updateMixture :: Typeable a => Cond -> Dist a -> Sampler a+updateMixture (Just cond) dist =+ case fromDynamic cond of+ 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)+ +conditioned, unconditioned :: Typeable a => Dist a -> Measure a+conditioned dist = Measure (\(cond:conds) -> smap (\a->(a,conds))+ (updateMixture cond dist))+unconditioned dist = Measure (\ conds -> smap (\a->(a,conds))+ (updateMixture Nothing dist))+ factor :: Prob -> Measure () factor p = Measure (\conds -> deterministic (point ((), conds) p)) --- Our language also includes the usual goodies of a lambda calculus-var :: a -> a-var = id--lit :: a -> a-lit = id--lam :: (a -> b) -> (a -> b)-lam f = f--app :: (a -> b) -> a -> b-app f x = f x--fix :: ((a -> b) -> (a -> b)) -> (a -> b)-fix g = f where f = g f--ifThenElse :: Bool -> a -> a -> a-ifThenElse True t _ = t-ifThenElse False _ e = e+condition :: Eq b => Measure (a, b) -> b -> Measure a+condition m b' =+ Measure (\ conds -> + sbind (unMeasure m conds)+ (\ ((a,b), cds) ->+ deterministic (if b==b' then point (a,cds) 1 else empty))) -- Drivers for testing finish :: Mixture (a, [Cond]) -> Mixture a finish (Mixture m) = Mixture (M.mapKeysMonotonic (\(a,[]) -> a) m) -sample :: (Ord a) => Int -> Measure a -> [Cond] -> IO (Mixture a)-sample !n measure conds = go n empty where+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 n m = once >>= \result -> go (n - 1) $! mappend m (finish result)+ go k m = once >>= \result -> go (k - 1) $! mappend m (finish result) -sample_ :: (Ord a, Show a) => Int -> Measure a -> [Cond] -> IO ()-sample_ !n measure conds = replicateM_ n (once >>= pr) where- once = getStdRandom (unMeasure measure conds)- pr = print . finish+sample :: (Ord a, Show a) => 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 logit :: Floating a => a -> a logit !x = 1 / (1 + exp (- x))
+ Language/Hakaru/Lambda.hs view
@@ -0,0 +1,22 @@+-- The lambda-calculus part of the language, which can be shared+module Language.Hakaru.Lambda(lit, dbl, lam, app, fix, ifThenElse) where++lit :: (Eq a) => a -> a+lit = id++-- raw lit is a pain to use. These are nicer+dbl :: Double -> Double+dbl = lit++lam :: (a -> b) -> (a -> b)+lam f = f++app :: (a -> b) -> a -> b+app f x = f x++fix :: ((a -> b) -> (a -> b)) -> (a -> b)+fix g = f where f = g f++ifThenElse :: Bool -> a -> a -> a+ifThenElse True t _ = t+ifThenElse False _ f = f
Language/Hakaru/Metropolis.hs view
@@ -1,21 +1,17 @@ {-# LANGUAGE RankNTypes, NoMonomorphismRestriction, BangPatterns, DeriveDataTypeable, GADTs, ScopedTypeVariables, ExistentialQuantification, StandaloneDeriving #-}+{-# OPTIONS -Wall #-} module Language.Hakaru.Metropolis where import System.Random (RandomGen, StdGen, randomR, getStdGen)-import System.IO -import Control.Monad import Data.Dynamic-import Data.Function (on) import Data.Maybe import qualified Data.Map.Strict as M--import RandomChoice-import Visual+import Language.Hakaru.Types {- @@ -28,254 +24,160 @@ -} type DistVal = Dynamic--data Dist a = Dist {logDensity :: a -> Likelihood,- sample :: forall g. RandomGen g => g -> (a, g)}-deriving instance Typeable1 Dist -data XRP = forall e. Typeable e => XRP (e, Dist e)+-- and what does XRP stand for?+data XRP where+ XRP :: Typeable e => (Density e, Dist e) -> XRP -unXRP :: Typeable a => XRP -> Maybe (a, Dist a)+unXRP :: Typeable a => XRP -> Maybe (Density a, Dist a) unXRP (XRP (e,f)) = cast (e,f) -type Var a = Int--type Likelihood = Double type Visited = Bool type Observed = Bool-type Cond = Maybe DistVal+type LL = LogLikelihood type Subloc = Int type Name = [Subloc]-type Database = M.Map Name (XRP, Likelihood, Visited, Observed)-newtype Measure a = Measure {unMeasure :: (RandomGen g) =>- (Name- ,Database- ,(Likelihood, Likelihood)- ,[Cond]- ,g- ) -> (a- ,Database- ,(Likelihood, Likelihood)- ,[Cond]- ,g)}- deriving (Typeable)+data DBEntry = DBEntry {+ xrp :: XRP, + llhd :: LL, + vis :: Visited,+ observed :: Observed }+type Database = M.Map Name DBEntry --- n is structural_name--- d is database--- ll is likelihood of expression--- conds is the observed data--- g is the random seed+data SamplerState g 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 +type Sampler a = forall g. (RandomGen g) => SamplerState g -> (a, SamplerState g) -lit :: (Eq a, Typeable a) => a -> a-lit = id+sreturn :: a -> Sampler a+sreturn x s = (x, s) +sbind :: Sampler a -> (a -> Sampler b) -> Sampler b+sbind s k = \ st -> let (v, s') = s st in k v s'++smap :: (a -> b) -> Sampler a -> Sampler b+smap f s = sbind s (\a -> sreturn (f a))++newtype Measure a = Measure {unMeasure :: Name -> Sampler a }+ deriving (Typeable)+ return_ :: a -> Measure a-return_ x = Measure (\ (n, d, l, conds, g) -> (x, d, l, conds, g))+return_ x = Measure $ \ _ -> sreturn x -makeXRP :: (Typeable a, RandomGen g) => Cond -> Dist a- -> Name -> Database -> g- -> (a, Database, Likelihood, Likelihood, g)-makeXRP obs dist' n db g =+updateXRP :: Typeable a => Name -> Cond -> Dist a -> Sampler a+updateXRP n obs dist' s@(S {ldb = db, seed = g}) = case M.lookup n db of- Just (xd, lb, b, ob) ->+ Just (DBEntry xd lb _ ob) -> let Just (xb, dist) = unXRP xd- (x,l) = case obs of- Just xd ->- let Just x = fromDynamic xd- in (x, logDensity dist x)+ (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 (XRP (x,dist),- l',- True,- ob) db- in (x, d1, l', 0, g)+ 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}) Nothing ->- let (xnew, l, g1) = case obs of+ let (xnew2, l, g2) = case obs of Just xdnew -> let Just xnew = fromDynamic xdnew in (xnew, logDensity dist' xnew, g) Nothing ->- (xnew, logDensity dist' xnew, g1)- where (xnew, g1) = sample dist' g- d1 = M.insert n (XRP (xnew, dist'),- l,- True,- isJust obs) db- in (xnew, d1, l, l, g1)--updateLikelihood :: (Typeable a, RandomGen g) => - Likelihood -> Likelihood ->- (a, Database, Likelihood, Likelihood, g) ->- [Cond] ->- (a, Database, (Likelihood, Likelihood), [Cond], g)-updateLikelihood llTotal llFresh (x,d,l,lf,g) conds =- (x, d, (llTotal+l, llFresh+lf), conds, g)--dirac :: (Eq a, Typeable a) => a -> Cond -> Measure a-dirac theta obs = Measure $ \(n, d, (llTotal,llFresh), conds, g) ->- let dist' = Dist {logDensity = (\ x -> if x == theta then 0 else log 0),- sample = (\ g -> (theta,g))}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds--bern :: Double -> Cond -> Measure Bool-bern p obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- let dist' = Dist {logDensity = (\ x -> log (if x then p else 1 - p)),- sample = (\ g -> case randomR (0, 1) g of- (t, g') -> (t <= p, g'))}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds--poisson :: Double -> Cond -> Measure Int-poisson l obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- let poissonLogDensity l x | l > 0 && x> 0 = (fromIntegral x)*(log l) - lnFact x - l- poissonLogDensity l x | x==0 = -l- poissonLogDensity _ _ = log 0- dist' = Dist {logDensity = poissonLogDensity l,- sample = poisson_rng l}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds--gamma :: Double -> Double -> Cond -> Measure Double-gamma shape scale obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- let dist' = Dist {logDensity = gammaLogDensity shape scale,- sample = gamma_rng shape scale}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds--beta :: Double -> Double -> Cond -> Measure Double-beta a b obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- let dist' = Dist {logDensity = betaLogDensity a b,- sample = beta_rng a b}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds--uniform :: Double -> Double -> Cond -> Measure Double-uniform lo hi obs = Measure $ \(n, d, (llTotal,llFresh), conds, g) ->- let uniformLogDensity lo hi x | lo <= x && x <= hi = log (recip (hi - lo))- uniformLogDensity _ _ x = log 0- dist' = Dist {logDensity = uniformLogDensity lo hi,- sample = (\ g -> randomR (lo, hi) g)}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds--normal :: Double -> Double -> Cond -> Measure Double-normal mu sd obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- let dist' = Dist {logDensity = normalLogDensity mu sd,- sample = normal_rng mu sd}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds--laplace :: Double -> Double -> Cond -> Measure Double-laplace mu sd obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- let dist' = Dist {logDensity = laplaceLogDensity mu sd,- sample = laplace_rng mu sd}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds+ 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}) -categorical :: (Eq a, Typeable a) => [(a,Double)] - -> Cond -> Measure a-categorical list obs = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- let categoricalLogDensity list x = log $ fromMaybe 0 (lookup x list)- categoricalSample list g = (elem, g1)- where- (p, g1) = randomR (0, total) g- elem = fst $ head $ filter (\(_,p0) -> p <= p0) sumList- sumList = scanl1 (\acc (a, b) -> (a, b + snd(acc))) list- total = sum $ map snd list- dist' = Dist {logDensity = categoricalLogDensity list,- sample = categoricalSample list}- xrp = makeXRP obs dist' n d g- in updateLikelihood llTotal llFresh xrp conds+updateLogLikelihood :: RandomGen g => + LL -> LL -> SamplerState g ->+ (LL, LL)+updateLogLikelihood llTotal llFresh s =+ let (l,lf) = llh2 s in (llTotal+l, llFresh+lf) -factor :: Likelihood -> Measure ()-factor l = Measure $ \(n, d, (llTotal, llFresh), conds, g) ->- ((), d, (llTotal + l, llFresh), conds, g)+factor :: LL -> Measure ()+factor l = Measure $ \ _ -> \ s ->+ let (llTotal, llFresh) = llh2 s+ in ((), s {llh2 = (llTotal + l, llFresh)}) -resample :: RandomGen g => XRP -> g ->- (XRP, Likelihood, Likelihood, Likelihood, g)-resample (XRP (x, dist)) g =- let (x', g1) = sample dist g- fwd = logDensity dist x'- rvs = logDensity dist x- l' = fwd- in (XRP (x', dist), l', fwd, rvs, g1)+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)) bind :: Measure a -> (a -> Measure b) -> Measure b-bind (Measure m) cont = Measure $ \ (n,d,ll,conds,g) ->- let (v, d1, ll1, conds1, g1) = m (0:n, d, ll, conds, g)- in unMeasure (cont v) (1:n, d1, ll1, conds1, g1)+bind (Measure m) cont = Measure $ \ n ->+ sbind (m (0:n)) (\ a -> unMeasure (cont a) (1:n)) -conditioned :: (Cond -> Measure a) -> Measure a-conditioned f = Measure $ \ (n,d,ll,cond:conds,g) ->- unMeasure (f cond) (n, d, ll, conds, g)+conditioned :: Typeable a => Dist a -> Measure a+conditioned dist = Measure $ \ n -> + \s@(S {cnds = cond:conds }) ->+ updateXRP n cond dist s{cnds = conds} -unconditioned :: (Cond -> Measure a) -> Measure a-unconditioned f = f Nothing+unconditioned :: Typeable a => Dist a -> Measure a+unconditioned dist = Measure $ \ n ->+ updateXRP n Nothing dist instance Monad Measure where return = return_- (>>=) = bind--lam :: (a -> b) -> (a -> b)-lam f = f--app :: (a -> b) -> a -> b-app f x = f x--fix :: ((a -> b) -> (a -> b)) -> (a -> b)-fix g = f where f = g f--ifThenElse :: Bool -> a -> a -> a-ifThenElse True t _ = t-ifThenElse False _ f = f+ (>>=) = bind -run :: Measure a -> [Cond] -> IO (a, Database, Likelihood)-run (Measure prog) conds = do+run :: Measure a -> [Cond] -> IO (a, Database, LL)+run (Measure prog) cds = do g <- getStdGen- let (v, d, ll, conds1, g') =- prog ([0], M.empty, (0,0), conds, g)+ let (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, Likelihood, Likelihood, Likelihood, g)-traceUpdate (Measure prog) d conds g = do- let d1 = M.map (\ (x, l, _, ob) -> (x, l, False, ob)) d- let (v, d2, (llTotal, llFresh), conds1, g1) =- prog ([0], d1, (0,0), conds, g)- let (d3, stale_d) = M.partition (\ (_, _, v, _) -> v) d2- let llStale = M.foldl' (\ llStale (_,l,_,_) -> llStale + l)- 0 stale_d+ -> (a, Database, LL, LL, LL, g)+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)+ 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) initialStep :: Measure a -> [Cond] -> IO (a, Database,- Likelihood, Likelihood, Likelihood, StdGen)-initialStep prog conds = do+ LL, LL, LL, StdGen)+initialStep prog cds = do g <- getStdGen- return $ traceUpdate prog M.empty conds g+ return $ traceUpdate prog M.empty cds g -- TODO: Make a way of passing user-provided proposal distributions-updateDB :: (RandomGen g) => - Name -> Database -> Observed -> XRP -> g- -> (Database, Likelihood, Likelihood, Likelihood, g)-updateDB name db ob xd g = (db', l', fwd, rvs, g)- where db' = M.insert name (x', l', True, ob) db- (x', l', fwd, rvs, g1) = resample xd g+resample :: RandomGen g => Name -> Database -> Observed -> XRP -> g ->+ (Database, LL, LL, LL, g)+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) transition :: (Typeable a, RandomGen g) => Measure a -> [Cond]- -> a -> Database -> Likelihood -> g -> [a]-transition prog conds v db ll g =+ -> a -> Database -> LL -> g -> [a]+transition prog cds v db ll g = let dbSize = M.size db -- choose an unconditioned choice- (condDb, uncondDb) = M.partition (\ (_, _, _, ob) -> ob) db+ (_, uncondDb) = M.partition observed db (choice, g1) = randomR (0, (M.size uncondDb) -1) g- (name, (xd, l, _, ob)) = M.elemAt choice uncondDb- (db', l', fwd, rvs, g2) = updateDB name db ob xd g1- (v', db2, llTotal, llFresh, llStale, g3) = traceUpdate prog db' conds g2+ (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)@@ -283,12 +185,16 @@ (u, g4) = randomR (0 :: Double, 1) g3 in if (log u < a) then- v' : (transition prog conds v' db2 llTotal g4)+ v' : (transition prog cds v' db2 llTotal g4) else- v : (transition prog conds v db ll g4)+ v : (transition prog cds v db ll g4) mcmc :: Typeable a => Measure a -> [Cond] -> IO [a]-mcmc prog conds = do- (v, d, llTotal, llFresh, llStale, g) <- initialStep prog conds- return $ transition prog conds v d llTotal g+mcmc prog cds = do+ (v, d, llTotal, _, _, g) <- initialStep prog cds+ return $ 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)
+ Language/Hakaru/Mixture.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE RankNTypes, BangPatterns #-}+{-# OPTIONS -W #-}++module Language.Hakaru.Mixture (Prob, point, empty, scale,+ Mixture(..), toList, mnull, mmap, cross, mode) where++import Data.Monoid+import Data.Ord (comparing)+import Data.List (maximumBy)+import qualified Data.Map.Strict as M+import Data.Number.LogFloat hiding (isInfinite)+import Text.Show (showListWith)+import Numeric (showFFloat)++type Prob = LogFloat++-- Mixtures (the results of importance sampling)++newtype Mixture k = Mixture { unMixture :: M.Map k Prob }++instance (Show k) => Show (Mixture k) where+ showsPrec d (Mixture m) = showParen (d > 0) $+ showString "Mixture $ fromList " . showListWith s (M.toList m)+ where s (k,p) = showChar '('+ . shows k+ . showChar ','+ . (if isInfinite l || -42 < l && l < 42+ then showFFloat Nothing (fromLogFloat p :: Double)+ else showString "logToLogFloat " . showsPrec 11 l)+ . showChar ')'+ where l = logFromLogFloat p :: Double++instance (Ord k) => Monoid (Mixture k) where+ mempty = empty+ mappend m1 m2 = Mixture (M.unionWith (+) (unMixture m1) (unMixture m2))+ mconcat ms = Mixture (M.unionsWith (+) (map unMixture ms))++empty :: Mixture k+empty = Mixture M.empty++toList :: Mixture k -> [(k, Prob)]+toList = M.toList . unMixture++mnull :: Mixture k -> Bool+mnull = all (0>=) . M.elems . unMixture++point :: k -> Prob -> Mixture k+point k !v = Mixture (M.singleton k v)++scale :: Prob -> Mixture k -> Mixture k+scale !v = Mixture . M.map (v *) . unMixture++mmap :: (Ord k2) => (k1 -> k2) -> Mixture k1 -> Mixture k2+mmap f = Mixture . M.mapKeysWith (+) f . unMixture++cross :: (Ord k) => (k1 -> k2 -> k) -> Mixture k1 -> Mixture k2 -> Mixture k+cross f m1 m2 = mconcat [ mmap (`f` k) (scale v m1)+ | (k,v) <- M.toList (unMixture m2) ]++mode :: Mixture k -> (k, Prob)+mode (Mixture m) = maximumBy (comparing snd) (M.toList m)
+ Language/Hakaru/Sampler.hs view
@@ -0,0 +1,26 @@+{-# LANGUAGE RankNTypes #-}+{-# OPTIONS -W #-}++module Language.Hakaru.Sampler (Sampler, deterministic, sbind, smap) where++import Language.Hakaru.Mixture (Mixture, mnull, empty, scale, point)+import Language.Hakaru.Distribution (choose)+import System.Random (RandomGen)++-- Sampling procedures that produce one sample++type Sampler a = forall g. (RandomGen g) => g -> (Mixture a, g)++deterministic :: Mixture a -> Sampler a+deterministic m g = (m, g)++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) } } }++smap :: (a -> b) -> Sampler a -> Sampler b+smap f s = sbind s (\a -> deterministic (point (f a) 1))
Language/Hakaru/Symbolic.hs view
@@ -1,78 +1,185 @@-{-# LANGUAGE GADTs, TypeFamilies #-}-{-# OPTIONS -W #-}+{-# LANGUAGE GADTs, TypeFamilies, ScopedTypeVariables, FlexibleContexts #-}+{-# OPTIONS -Wall #-} module Language.Hakaru.Symbolic where +import Prelude hiding (Real)++data Real data Prob data Measure a data Dist a+data Exact --- Symbolic AST (from Syntax.hs)-class Symbolic repr where- real :: Double -> repr Double- bool :: Bool -> repr Bool- add, minus, mul, exp :: repr Double -> repr Double -> repr Double- sqrt, cos, sin :: repr Double -> repr Double- bind :: repr (Measure a) -> (repr a -> repr (Measure a)) - -> repr (Measure a)- ret :: repr a -> repr (Measure a)- uniformD, uniformC, normal :: repr Double -> repr Double -> repr (Dist Double)+class IntComp repr where+ int :: Integer -> repr Integer++class BoolComp repr where+ bool :: Bool -> repr Bool++class RealComp repr where+ real :: Rational -> repr Real+ exp :: repr Real -> repr Real -> repr Real+ sqrt, cos, sin :: repr Real -> repr Real++-- polymorphic operations and integer powering+-- probably should restrict 'a' to be some kind of numeric type?+class SymbComp repr where+ add, minus, mul :: repr a -> repr a -> repr a+ pow :: repr Real -> repr Integer -> repr Real+ scale :: repr Integer -> repr Real -> repr Real++class MeasMonad repr where+ bind :: repr (Measure a) -> (repr a -> repr (Measure b)) + -> repr (Measure b)+ ret :: repr a -> repr (Measure a)++class Distrib repr where+ uniform, normal :: repr Real -> repr Real -> repr (Dist Real)+ uniformD :: repr Integer -> repr Integer -> repr (Dist Integer)++class Conditioning repr where conditioned, unconditioned :: repr (Dist a) -> repr (Measure a) -- Printer (to Maple)+data Pos = Front | Back type VarCounter = Int-newtype Maple a = Maple { unMaple :: Bool -> VarCounter -> String }+newtype Maple a = Maple { unMaple :: Pos -> VarCounter -> String } -instance Symbolic Maple where- real x = Maple $ \_ _ -> show x- bool x = Maple $ \_ _ -> show x- add e1 e2 = Maple $ \f h -> unMaple e1 f h ++ "+" ++ unMaple e2 f h- minus e1 e2 = Maple $ \f h -> unMaple e1 f h ++ "-" ++ unMaple e2 f h - mul e1 e2 = Maple $ \f h -> unMaple e1 f h ++ "*" ++ unMaple e2 f h- exp e1 e2 = Maple $ \f h -> unMaple e1 f h ++ "^" ++ unMaple e2 f h- sqrt e = Maple $ \f h -> "sqrt(" ++ unMaple e f h ++ ")"- cos e = Maple $ \f h -> "cos(" ++ unMaple e f h ++ ")"- sin e = Maple $ \f h -> "sin(" ++ unMaple e f h ++ ")"- bind m c = Maple $ \f h -> unMaple m True h ++ - unMaple (c (Maple $ \_ _ -> ("x" ++ show h))) (f) (succ h)- ++ unMaple m False h - uniformC e1 e2 = Maple $ \f h -> if f == True then - show (1/((read (unMaple e2 f h) :: Double) - - (read (unMaple e1 f h) :: Double))) ++ " * Int (" - else ", x" ++ show h ++ "=" ++ unMaple e1 f h ++ ".." ++ - unMaple e2 f h ++ ")"- uniformD e1 e2 = Maple $ \f h -> if f == True then - show (1/((read (unMaple e2 f h) :: Double) - - (read (unMaple e1 f h) :: Double))) ++ " * Sum (" - else ", x" ++ show h ++ "=" ++ unMaple e1 f h ++ ".." ++ - unMaple e2 f h ++ ")" - normal e1 e2 = Maple $ \f h -> if f == True then - "Int (PDF (Normal (" ++ unMaple e1 f h ++ ", " ++- unMaple e2 f h ++ ", x" ++ show h ++ ") * " - else ", x" ++ show h ++ "=" ++ unMaple e1 f h ++ ".." ++ - unMaple e2 f h ++ ")" - unconditioned e = Maple $ \f h -> unMaple e f h- conditioned e = Maple $ \f h -> unMaple e f h - ret e = Maple $ \f h -> "g(" ++ unMaple e f h ++ ")"+type family MPL a+type instance MPL Real = Rational+type instance MPL Integer = Integer+type instance MPL Bool = Bool -view e = unMaple e True 0+-- if it weren't for the constraints, we could/should use Applicative+pure :: Show (MPL a) => MPL a -> Maple a+pure x = Maple $ \_ _ -> show x +liftA1 :: (String -> String) -> Maple a -> Maple a+liftA1 pr x = Maple $ \f h -> pr (unMaple x f h)++liftA2 :: (String -> String -> String) -> Maple a -> Maple a -> Maple a+liftA2 pr e1 e2 = Maple $ \f h -> pr (unMaple e1 f h) (unMaple e2 f h)++-- variant for ret+liftA1M :: (String -> String) -> Maple a -> Maple (Measure a)+liftA1M pr x = Maple $ \f h -> pr (unMaple x f h)++-- variant for pow+liftA2aba :: (String -> String -> String) -> Maple a -> Maple b -> Maple a+liftA2aba pr e1 e2 = Maple $ \f h -> pr (unMaple e1 f h) (unMaple e2 f h)++-- variant for scale+liftA2baa :: (String -> String -> String) -> Maple b -> Maple a -> Maple a+liftA2baa pr e1 e2 = Maple $ \f h -> pr (unMaple e1 f h) (unMaple e2 f h)++mkPr :: String -> (String -> String)+mkPr s t = s ++ "(" ++ t ++ ")"++d2m :: Maple (Dist a) -> Maple (Measure a)+d2m e = Maple $ unMaple e++infixPr :: String -> (String -> String -> String)+infixPr s a b = a ++ s ++ b++-- This is quite scary. Probably a mistake.+reify :: forall a. Read a => Pos -> VarCounter -> Maple a -> a+reify f h e = (read (unMaple e f h) :: a)++name :: String -> VarCounter -> String+name s h = s ++ show h++var :: String -> VarCounter -> Maple a+var s h = Maple $ \_ _ -> name s h++binder :: (String -> String -> Maybe String) -> + (String -> String -> VarCounter -> Maybe String) ->+ String -> + Maple a -> Maple a -> Maple (Dist a)+binder pr1 pr2 oper e1 e2 = Maple $ pr_+ where+ pr_ Front h = let x1 = unMaple e1 Front h+ x2 = unMaple e2 Front h in+ case (pr1 x1 x2, pr2 x1 x2 h) of+ (Just z, Just w) -> z ++ " * " ++ oper ++ " (" ++ w ++ " * "+ (Nothing, Just w) -> oper ++ " (" ++ w ++ " * "+ (Just z, Nothing) -> z ++ " * " ++ oper ++ " ("+ (Nothing, Nothing) -> oper ++ " ("+ pr_ Back h = ", " ++ (name "x" h) ++ "=" ++ unMaple e1 Back h +++ ".." ++ unMaple e2 Back h ++ ")"++instance RealComp Maple where+ -- serious problem here: all exact numbers will be printed as+ -- floats, which will really hamper the use of Maple in any + -- serious way. This needs a rethink.+ real = pure+ exp = liftA2 $ infixPr "^"+ sqrt = liftA1 $ mkPr "sqrt"+ cos = liftA1 $ mkPr "cos"+ sin = liftA1 $ mkPr "sin"++instance SymbComp Maple where+ add = liftA2 $ infixPr "+"+ minus = liftA2 $ infixPr "-"+ mul = liftA2 $ infixPr "*"+ pow = liftA2aba $ infixPr "^"+ scale = liftA2baa $ infixPr "*"++instance IntComp Maple where+ int = pure++instance BoolComp Maple where+ bool = pure++instance MeasMonad Maple where+ ret = liftA1M $ mkPr "g"+ bind m c = Maple $ \f h -> unMaple m Front h ++ + unMaple (c $ var "x" h) f (succ h) +++ unMaple m Back h ++instance Distrib Maple where+ uniform = binder (\e1 e2 -> Just $ + show (1/((read e2 :: Rational) - (read e1 :: Rational)))) + (\_ _ _ -> Nothing) "Int"+ uniformD = binder (\e1 e2 -> + let d = (read e2 :: Integer) - (read e1) in+ if d == 1 then Nothing+ else Just $ "(1 / " ++ show d ++ ")")+ (\_ _ _ -> Nothing) "Sum"+ normal = binder (\_ _ -> Nothing) + (\e1 e2 h -> Just $ "PDF (Normal (" ++ e1 ++ ", " ++ e2 ++ "), " ++ (name "x" h) ++ ")")+ "Int"++instance Conditioning Maple where+ unconditioned = d2m+ conditioned = d2m++view :: Maple a -> String+view e = unMaple e Front 0++lift :: Maple Integer -> Maple Real+lift x = Maple $ \f h -> unMaple x f h+ -- TEST CASES-exp1 = unconditioned (uniformC (real 1) (real 3)) `bind` \s ->- ret s+exp1, exp2, exp3, exp4 :: Maple (Measure Real)+exp1 = unconditioned (uniform (real 1) (real 3)) `bind` ret -- Borel's Paradox Simplified-exp2 = unconditioned (uniformD (real 1) (real 3)) `bind` \s ->- unconditioned (uniformC (real (-1)) (real 1)) `bind` \x ->- let y = s `mul` x in ret y+exp2 = unconditioned (uniformD (int 1) (int 3)) `bind` \s ->+ unconditioned (uniform (real (-1)) (real 1)) `bind` \x ->+ let y = s `scale` x in ret y -- Borel's Paradox-exp3 = unconditioned (uniformD (real 1) (real 2)) `bind` \s ->- unconditioned (uniformC (real (-1)) (real 1)) `bind` \x ->+exp3 = unconditioned (uniformD (int 1) (int 2)) `bind` \s ->+ unconditioned (uniform (real (-1)) (real 1)) `bind` \x -> let y = (Language.Hakaru.Symbolic.sqrt ((real 1 ) `minus` - (Language.Hakaru.Symbolic.exp s (real 2)))) `mul`- (Language.Hakaru.Symbolic.sin x) in ret y + (Language.Hakaru.Symbolic.pow (lift s) (int 2)))) `mul`+ (Language.Hakaru.Symbolic.sin x) in ret y +exp4 = unconditioned (normal (real 1) (real 4)) `bind` ret++test, test2, test3, test4 :: String test = view exp1 test2 = view exp2 test3 = view exp3+test4 = view exp4
− Language/Hakaru/Syntax.hs
@@ -1,185 +0,0 @@-{-# LANGUAGE TypeFamilies, ConstraintKinds, GADTs, FlexibleContexts #-}--module Language.Hakaru.Syntax where---- The syntax--import GHC.Exts (Constraint)---- TODO: The pretty-printing semantics--import qualified Text.PrettyPrint as PP---- The importance-sampling semantics--import Types (Cond, CSampler)-import Data.Dynamic (Typeable)-import qualified Data.Number.LogFloat as LF-import qualified Language.Hakaru.ImportanceSampler as IS---- The Metropolis-Hastings semantics--import qualified Language.Hakaru.Metropolis as MH---- The syntax--data Prob-data Measure a-data Dist a--class Mochastic repr where- type Type repr a :: Constraint- real :: Double -> repr Double- bool :: Bool -> repr Bool- add, mul :: repr Double -> repr Double -> repr Double- neg :: repr Double -> repr Double- neg = mul (real (-1))- logFloat, logToLogFloat- :: repr Double -> repr Prob- unbool :: repr Bool -> repr c -> repr c- -> repr c- pair :: repr a -> repr b -> repr (a, b)- unpair :: repr (a, b) -> (repr a -> repr b -> repr c)- -> repr c- inl :: repr a -> repr (Either a b)- inr :: repr b -> repr (Either a b)- uneither :: repr (Either a b) -> (repr a -> repr c) -> (repr b -> repr c)- -> repr c- nil :: repr [a]- cons :: repr a -> repr [a] -> repr [a]- unlist :: repr [a] -> repr c -> (repr a -> repr [a] -> repr c)- -> repr c- ret :: repr a -> repr (Measure a)- bind :: repr (Measure a) -> (repr a -> repr (Measure b))- -> repr (Measure b)- conditioned, unconditioned :: repr (Dist a) -> repr (Measure a)- factor :: repr Prob -> repr (Measure ())- dirac :: (Type repr a) => repr a -> repr (Dist a)- categorical :: (Type repr a) => repr [(a, Prob)] -> repr (Dist a)- bern :: (Type repr Bool) => repr Double -> repr (Dist Bool)- bern p = categorical $- cons (pair (bool True) (logFloat p)) $- cons (pair (bool False) (logFloat (add (real 1) (neg p)))) $- nil- normal, uniform- :: repr Double -> repr Double -> repr (Dist Double)- poisson :: repr Double -> repr (Dist Int)---- TODO: The initial (AST) "semantics"--- (Hey Oleg, is there any better way to deal with the Type constraint, so that--- the AST constructor doesn't have to take a repr constructor argument?)--data AST repr a where- Real :: Double -> AST repr Double- Unbool :: AST repr Bool -> AST repr c -> AST repr c -> AST repr c- Categorical :: (Type repr a) => AST repr [(a, Prob)] -> AST repr (Dist a)- -- ...--instance (Mochastic repr) => Mochastic (AST repr) where- type Type (AST repr) a = Type repr a- real = Real- unbool = Unbool- categorical = Categorical- -- ...--eval :: (Mochastic repr) => AST repr a -> repr a-eval (Real x) = real x-eval (Unbool b x y) = unbool (eval b) (eval x) (eval y)-eval (Categorical xps) = categorical (eval xps)--- ...---- TODO: The pretty-printing semantics--newtype PP a = PP (Int -> PP.Doc)---- The importance-sampling semantics--newtype IS a = IS (IS' a)-type family IS' a-type instance IS' (Measure a) = IS.Measure (IS' a)-type instance IS' (Dist a) = CSampler (IS' a)-type instance IS' [a] = [IS' a]-type instance IS' (a, b) = (IS' a, IS' b)-type instance IS' (Either a b) = Either (IS' a) (IS' b)-type instance IS' () = ()-type instance IS' Bool = Bool-type instance IS' Double = Double-type instance IS' Prob = LF.LogFloat-type instance IS' Int = Int--instance Mochastic IS where- type Type IS a = (Eq (IS' a), Typeable (IS' a))- real = IS- bool = IS- add (IS x) (IS y) = IS (x + y)- mul (IS x) (IS y) = IS (x * y)- neg (IS x) = IS (-x)- logFloat (IS x) = IS (LF.logFloat x)- logToLogFloat (IS x) = IS (LF.logToLogFloat x)- unbool (IS b) x y = if b then x else y- pair (IS x) (IS y) = IS (x, y)- unpair (IS (x, y)) c = c (IS x) (IS y)- inl (IS x) = IS (Left x)- inr (IS x) = IS (Right x)- uneither (IS e) c d = either (c . IS) (d . IS) e- nil = IS []- cons (IS x) (IS xs) = IS (x:xs)- unlist (IS []) n c = n- unlist (IS (x:xs)) n c = c (IS x) (IS xs)- ret (IS x) = IS (return x)- bind (IS m) k = IS (m >>= \x -> case k (IS x) of IS n -> n)- conditioned (IS dist) = IS (IS.conditioned dist)- unconditioned (IS dist) = IS (IS.unconditioned dist)- factor (IS p) = IS (IS.factor p)- dirac (IS x) = IS (IS.dirac x)- categorical (IS xps) = IS (IS.categorical xps)- bern (IS p) = IS (IS.bern p)- normal (IS m) (IS s) = IS (IS.normal m s)- uniform (IS lo) (IS hi) = IS (IS.uniformC lo hi)- poisson (IS l) = IS (IS.poisson l)---- The Metropolis-Hastings semantics--newtype MH a = MH (MH' a)-type family MH' a-type instance MH' (Measure a) = MH.Measure (MH' a)-type instance MH' (Dist a) = MH.Cond -> MH.Measure (MH' a)-type instance MH' [a] = [MH' a]-type instance MH' (a, b) = (MH' a, MH' b)-type instance MH' (Either a b) = Either (MH' a) (MH' b)-type instance MH' () = ()-type instance MH' Bool = Bool-type instance MH' Double = Double-type instance MH' Prob = MH.Likelihood-type instance MH' Int = Int--instance Mochastic MH where- type Type MH a = (Eq (MH' a), Typeable (MH' a), Show (MH' a))- real = MH- bool = MH- add (MH x) (MH y) = MH (x + y)- mul (MH x) (MH y) = MH (x * y)- neg (MH x) = MH (-x)- logFloat (MH x) = MH (LF.logFromLogFloat (LF.logFloat x))- logToLogFloat (MH x) = MH (LF.logFromLogFloat (LF.logToLogFloat x))- unbool (MH b) x y = if b then x else y- pair (MH x) (MH y) = MH (x, y)- unpair (MH (x, y)) c = c (MH x) (MH y)- inl (MH x) = MH (Left x)- inr (MH x) = MH (Right x)- uneither (MH e) c d = either (c . MH) (d . MH) e- nil = MH []- cons (MH x) (MH xs) = MH (x:xs)- unlist (MH []) n c = n- unlist (MH (x:xs)) n c = c (MH x) (MH xs)- ret (MH x) = MH (return x)- bind (MH m) k = MH (m >>= \x -> case k (MH x) of MH n -> n)- conditioned (MH dist) = MH (MH.conditioned dist)- unconditioned (MH dist) = MH (MH.unconditioned dist)- factor (MH p) = MH (MH.factor p)- dirac (MH x) = MH (MH.dirac x)- categorical (MH xps) = MH (MH.categorical xps)- bern (MH p) = MH (MH.bern p)- normal (MH m) (MH s) = MH (MH.normal m s)- uniform (MH lo) (MH hi) = MH (MH.uniform lo hi)- poisson = error "poisson: not implemented for MH" -- TODO
+ Language/Hakaru/Types.hs view
@@ -0,0 +1,29 @@+{-# LANGUAGE RankNTypes, BangPatterns, DeriveDataTypeable, StandaloneDeriving #-}+{-# OPTIONS -W #-}++module Language.Hakaru.Types where++import Data.Dynamic+import System.Random++-- Basic types for conditioning and conditioned sampler+data Density a = Lebesgue !a | Discrete !a deriving Typeable+type Cond = Maybe Dynamic++fromDiscrete :: Density t -> t+fromDiscrete (Discrete a) = a+fromDiscrete _ = error "got a non-discrete sampler"++fromLebesgue :: Density t -> t+fromLebesgue (Lebesgue a) = a+fromLebesgue _ = error "got a discrete sampler"++fromDensity :: Density t -> t+fromDensity (Discrete a) = a+fromDensity (Lebesgue a) = a++type LogLikelihood = Double+data Dist a = Dist {logDensity :: Density a -> LogLikelihood,+ distSample :: forall g.+ RandomGen g => g -> (Density a, g)}+deriving instance Typeable1 Dist
+ Language/Hakaru/Util/Coda.hs view
@@ -0,0 +1,12 @@+module Language.Hakaru.Util.Coda where++import Statistics.Autocorrelation+import qualified Data.Packed.Vector as V+import qualified Data.Vector.Generic as G++effectiveSampleSize :: [Double] -> Double+effectiveSampleSize samples = n / (1 + 2*(G.sum rho))+ where n = fromIntegral (V.dim vec)+ vec = V.fromList samples+ cov = autocovariance vec+ rho = G.map (/ G.head cov) cov
+ Language/Hakaru/Util/Csv.hs view
@@ -0,0 +1,40 @@+{-# LANGUAGE TypeOperators #-}++module Language.Hakaru.Util.Csv ((:::)((:::)), decodeFile, decodeGZipFile,+ decodeFileStream, decodeGZipFileStream) where++import Data.Csv ( HasHeader(..), FromRecord(..), FromField(..)+ , ToRecord(..), ToField(..), decode)+import qualified Data.Csv.Streaming as CS (decode)+import Codec.Compression.GZip (decompress)+import qualified Data.Foldable as F+import qualified Data.ByteString.Lazy as B+import qualified Data.Vector as V+import Control.Applicative ((<*>), (<$>))++data a ::: b = a ::: b deriving (Eq, Ord, Read, Show)+infixr 5 :::++instance (FromField a, FromRecord b) => FromRecord (a ::: b) where+ parseRecord v | V.null v = fail "too few fields in input record"+ | otherwise = (:::) <$> parseField (V.head v) <*> parseRecord (V.tail v)++instance (ToField a, ToRecord b) => ToRecord (a ::: b) where+ toRecord (a ::: b) = V.cons (toField a) (toRecord b)++decodeBytes :: FromRecord a => B.ByteString -> [a]+decodeBytes bs = case decode HasHeader bs of+ Left _ -> []+ Right v -> V.toList v++decodeFile :: FromRecord a => FilePath -> IO [a]+decodeFile = fmap decodeBytes . B.readFile++decodeGZipFile :: FromRecord a => FilePath -> IO [a]+decodeGZipFile = fmap (decodeBytes . decompress) . B.readFile++decodeFileStream :: FromRecord a => FilePath -> IO [a]+decodeFileStream = fmap (F.toList . CS.decode HasHeader) . B.readFile++decodeGZipFileStream :: FromRecord a => FilePath -> IO [a]+decodeGZipFileStream = fmap (F.toList . CS.decode HasHeader . decompress) . B.readFile
+ Language/Hakaru/Util/Extras.hs view
@@ -0,0 +1,78 @@+{-|+ 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.+-}++module Language.Hakaru.Util.Extras where++import qualified Data.Sequence as S+import System.Random+import Data.Maybe+import qualified Data.Foldable as F++extract :: S.Seq a -> Int -> Maybe (S.Seq a, a)+extract s i | S.null r = Nothing+ | otherwise = Just (a S.>< c, b)+ 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+ return $ extract s i++{-| + Given a sequence, return a *sorted* sequence of+ n randomly selected elements from *distinct positions* in the sequence+-}++randomElems :: Ord a => S.Seq a -> Int -> IO (S.Seq a)+randomElems = randomElemsTR S.empty++randomElemsTR :: Ord a => S.Seq a -> S.Seq a -> Int -> IO (S.Seq a)+randomElemsTR ixs s n+ | n == S.length s = return $ S.unstableSort s+ | n == 1 = do (_,i) <- fmap fromJust (randomExtract s)+ return.S.unstableSort $ i S.<| ixs+ | otherwise = do (s',i) <- fmap fromJust (randomExtract s)+ (randomElemsTR $! (i S.<| ixs)) s' (n-1)++{-|+ Chop a sequence at the given indices. + Assume number of indices given < length of sequence to be chopped+-}++pieces :: S.Seq a -> S.Seq Int -> [S.Seq a]+pieces s ixs = let f (ps,r,x) y = let (p,r') = S.splitAt (y-x) r+ in (p:ps,r',y)+ g (a,b,_) = b:a+ in g $ F.foldl f ([],s,0) ixs++{-|+ Given n, chop a sequence at m random points+ where m = min (length-1, n-1)+-}++randomPieces :: Int -> S.Seq a -> IO [S.Seq a]+randomPieces n s+ | n >= l = return $ F.toList $ fmap S.singleton s+ | otherwise = do ixs <- randomElems (S.fromList [1..l-1]) (n-1)+ return $ pieces s ixs+ where l = S.length s++{-|+ > pairs [1,2,3,4]+ [(1,2),(1,3),(1,4),(2,3),(2,4),(3,4)]+ > pairs [1,2,4,4]+ [(1,2),(1,4),(1,4),(2,4),(2,4),(4,4)]+-}++pairs :: [a] -> [(a,a)]+pairs [] = []+pairs (x:xs) = (zip (repeat x) xs) ++ pairs xs++l2Norm :: Floating a => [a] -> a+l2Norm l = sqrt.sum $ zipWith (*) l l
+ Language/Hakaru/Util/Visual.hs view
@@ -0,0 +1,43 @@+{-# LANGUAGE OverloadedStrings #-}++module Language.Hakaru.Util.Visual where++import System.IO+import Control.Monad++import Data.Aeson+import Data.List+import qualified Data.Text as T+import qualified Data.ByteString.Lazy.Char8 as B+import qualified Data.ByteString.Char8 as BS++plot :: Show a => [a] -> String -> IO ()+plot samples filename = do+ h <- openFile filename WriteMode+ hPrint h samples+ hClose h++batchPrint :: Show a => Int -> [a] -> IO ()+batchPrint n l = do+ let batch = take n l+ print batch+ when (length batch == n) $ batchPrint n (drop n l)++viz :: ToJSON a => Int -> [String] -> [[a]] -> IO ()+viz n name samples = viz' n 50 0 name samples++viz' :: ToJSON a => Int -> Int -> Int -> [String] -> [[a]] -> IO ()+viz' n c cur name samples = do+ putStrLn batch+ when (c+cur < n) $+ viz' n c (cur+c) name (drop c samples)+ where+ total = "total_samples" .= n+ current_sample = "current_sample" .= cur+ chunk = object (zipWith (\ name s -> T.pack name .= s)+ name+ (transpose $ take c samples))+ batch = B.unpack $ encode+ (object ["rvars" .= chunk,+ total,+ current_sample])
− Mixture.hs
@@ -1,61 +0,0 @@-{-# LANGUAGE RankNTypes, BangPatterns #-}-{-# OPTIONS -W #-}--module Mixture (Prob, point, empty, scale,- Mixture(..), toList, mnull, mmap, cross, mode) where--import Data.Monoid-import Data.Ord (comparing)-import Data.List (maximumBy)-import qualified Data.Map.Strict as M-import Data.Number.LogFloat hiding (isInfinite)-import Text.Show (showListWith)-import Numeric (showFFloat)--type Prob = LogFloat---- Mixtures (the results of importance sampling)--newtype Mixture k = Mixture { unMixture :: M.Map k Prob }--instance (Show k) => Show (Mixture k) where- showsPrec d (Mixture m) = showParen (d > 0) $- showString "Mixture $ fromList " . showListWith s (M.toList m)- where s (k,p) = showChar '('- . shows k- . showChar ','- . (if isInfinite l || -42 < l && l < 42- then showFFloat Nothing (fromLogFloat p :: Double)- else showString "logToLogFloat " . showsPrec 11 l)- . showChar ')'- where l = logFromLogFloat p :: Double--instance (Ord k) => Monoid (Mixture k) where- mempty = empty- mappend m1 m2 = Mixture (M.unionWith (+) (unMixture m1) (unMixture m2))- mconcat ms = Mixture (M.unionsWith (+) (map unMixture ms))--empty :: Mixture k-empty = Mixture M.empty--toList :: Mixture k -> [(k, Prob)]-toList = M.toList . unMixture--mnull :: Mixture k -> Bool-mnull = all (0>=) . M.elems . unMixture--point :: k -> Prob -> Mixture k-point k !v = Mixture (M.singleton k v)--scale :: Prob -> Mixture k -> Mixture k-scale !v = Mixture . M.map (v *) . unMixture--mmap :: (Ord k2) => (k1 -> k2) -> Mixture k1 -> Mixture k2-mmap f = Mixture . M.mapKeysWith (+) f . unMixture--cross :: (Ord k) => (k1 -> k2 -> k) -> Mixture k1 -> Mixture k2 -> Mixture k-cross f m1 m2 = mconcat [ mmap (`f` k) (scale v m1)- | (k,v) <- M.toList (unMixture m2) ]--mode :: Mixture k -> (k, Prob)-mode (Mixture m) = maximumBy (comparing snd) (M.toList m)
− RandomChoice.hs
@@ -1,145 +0,0 @@-{-# LANGUAGE BangPatterns #-}-module RandomChoice where--import System.Random-import Mixture-import Data.Maybe (fromMaybe)-import Data.List (findIndex, foldl')-import Numeric.SpecFunctions-import qualified Data.Vector.Unboxed as U-import qualified Data.Map.Strict as M-import qualified Data.Number.LogFloat as LF--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--choose :: (RandomGen g) => Mixture k -> g -> (k, Prob, g)-choose (Mixture m) g0 =- 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- 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)--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)--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 _ _ _ = error "normal: invalid parameters"--normalLogDensity mu sd x = (-tau * square (x - mu)- + log (tau / pi / 2)) / 2- where square y = y * y- tau = 1 / square sd--lnFact = logFactorial---- Makes use of Atkinson's algorithm as described in:--- Monte Carlo Statistical Methods pg. 55------ Further discussion at:--- http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/-poisson_rng :: (RandomGen g) => Double -> g -> (Int, g)-poisson_rng lambda g0 = make_poisson g0- where smu = sqrt lambda- b = 0.931 + 2.53*smu- a = -0.059 + 0.02483*b- vr = 0.9277 - 3.6224/(b - 2)- arep = 1.1239 + 1.1368/(b-3.4)- lnlam = log lambda-- make_poisson :: (RandomGen g) => g -> (Int,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 (us, v, k) of- (us,v,k) | us >= 0.07 && v <= vr -> (k, g2)- (_,_, k) | k < 0 -> make_poisson g2- (us,v,k) | us <= 0.013 && v > us -> make_poisson g2- (us,v,k) | accept_region us v k -> (k, g2)- _ -> make_poisson g2-- accept_region us v k = log (v * arep / (a/(us*us)+b)) <=- -lambda + (fromIntegral k)*lnlam - lnFact k---- 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 shape scale g | shape <= 0.0 = error "gamma: got a negative shape paramater"-gamma_rng shape scale g | scale <= 0.0 = error "gamma: got a negative scale paramater"-gamma_rng shape scale g | shape < 1.0 = (gvar2, g2)- where (gvar1, g1) = gamma_rng (shape + 1) scale g- (w, g2) = randomR (0,1) g1- gvar2 = scale * gvar1 * (w ** recip shape) -gamma_rng shape scale g = - let d = shape - 1/3- c = recip $ sqrt $ 9*d- -- Algorithm recommends inlining normal generator- n g = normal_rng 1 c g- (v, g2) = until (\x -> fst x > 0.0) (\ (_, g) -> normal_rng 1 c g) (n g)- 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 -> (scale*d*v3, g3)- False -> gamma_rng shape scale g3--gammaLogDensity shape scale x | x>= 0 && shape > 0 && scale > 0 =- scale * log shape - scale * x + (shape - 1) * log x - logGamma shape-gammaLogDensity _ _ _ = log 0--beta_rng :: (RandomGen g) => Double -> Double -> g -> (Double, g)--- Consider adding case for a <= 1 && b <= 1-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)--betaLogDensity a b x | x < 0 || x > 1 = error "beta: value must be between 0 and 1"-betaLogDensity a b x | a <= 0 || b <= 0 = error "beta: parameters must be positve" -betaLogDensity a b x = (logGamma (a + b)- - logGamma a- - logGamma b- + x * log (a - 1)- + (1 - x) * log (b - 1))--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)--laplaceLogDensity mu sd x = - log (2 * sd) - abs (x - mu) / sd---- 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 (a, total, g) = (map (/ total) a, g)--dirichletLogDensity a x | all (> 0) x = sum (zipWith logTerm a x) + logGamma (sum a)- where sum a = foldl' (+) 0 a- logTerm a x = (a-1) * log x - logGamma a-dirichletLogDensity _ _ = error "dirichlet: all values must be between 0 and 1"
− Sampler.hs
@@ -1,26 +0,0 @@-{-# LANGUAGE RankNTypes #-}-{-# OPTIONS -W #-}--module Sampler (Sampler, deterministic, sbind, smap) where--import Mixture (Mixture, mnull, empty, scale, point)-import RandomChoice (choose)-import System.Random (RandomGen)---- Sampling procedures that produce one sample--type Sampler a = forall g. (RandomGen g) => g -> (Mixture a, g)--deterministic :: Mixture a -> Sampler a-deterministic m g = (m, g)--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) } } }--smap :: (a -> b) -> Sampler a -> Sampler b-smap f s = sbind s (\a -> deterministic (point (f a) 1))
+ Tests/Distribution.hs view
@@ -0,0 +1,106 @@+{-# LANGUAGE BangPatterns #-}++module Tests.Distribution where++import Control.Monad+import qualified System.Random.MWC as MWC++import Language.Hakaru.Types+import Language.Hakaru.Util.Coda+import Language.Hakaru.Distribution hiding (choose)++import Test.QuickCheck+import Test.QuickCheck.Monadic as QM+import Test.Framework.Providers.QuickCheck2 (testProperty)++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++quickArg :: IO ()+quickArg = quickCheckWith stdArgs {maxSuccess = 2000} (\ x -> almostEqual tol x x)+ where tol :: Double+ 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]++betaTest = do+ Positive a <- QM.pick arbitrary+ Positive b <- QM.pick arbitrary+ g <- QM.run $ MWC.create+ samples <- QM.run $ replicateM 1000 $ distSample (beta 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 / (a + b)+ var a b = a*b / ((sq $ a + b) * (a + b + 1))++bernTest = do+ p <- QM.pick $ choose (0, 1)+ g <- QM.run $ MWC.create+ samples <- QM.run $ replicateM 1000 $ distSample (bern p) g+ let (mean, variance) = meanVariance (map fromDiscreteToNum samples)+ QM.assert $ (almostEqual tol mean (mu p)) && + (almostEqual tol variance (var p))+ 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
@@ -0,0 +1,118 @@+{-# LANGUAGE BangPatterns #-}++module Tests.ImportanceSampler where++import Data.Dynamic+import Language.Hakaru.Types+import Language.Hakaru.Lambda+import Language.Hakaru.Distribution+import Language.Hakaru.ImportanceSampler++-- import Test.QuickCheck.Monadic+import Tests.Models++-- Some test programs in our language++test_mixture :: IO ()+test_mixture = sample prog_mixture conds >>=+ print . take 10 >>+ putChar '\n' >>+ empiricalMeasure 1000 prog_mixture conds >>=+ print+ where conds = [Just (toDyn (Lebesgue 2 :: Density Double))]++prog_dup :: Measure (Bool, Bool)+prog_dup = do+ let c = unconditioned (bern 0.5)+ x <- c+ y <- c+ return (x,y)++prog_dbn :: Measure Bool+prog_dbn = do+ s0 <- unconditioned (bern 0.75)+ s1 <- unconditioned (if s0 then bern 0.75 else bern 0.25)+ _ <- conditioned (if s1 then bern 0.90 else bern 0.10)+ s2 <- unconditioned (if s1 then bern 0.75 else bern 0.25)+ _ <- conditioned (if s2 then bern 0.90 else bern 0.10)+ return s2++test_dbn :: IO ()+test_dbn = sample prog_dbn conds >>=+ print . take 10 >>+ putChar '\n' >>+ empiricalMeasure 1000 prog_dbn conds >>=+ print + where conds = [Just (toDyn (Discrete True)),+ Just (toDyn (Discrete True))]++prog_hmm :: Integer -> Measure Bool+prog_hmm n = do+ s <- unconditioned (bern 0.75) + loop_hmm n s++loop_hmm :: Integer -> (Bool -> Measure Bool)+loop_hmm !numLoops s = do+ _ <- conditioned (if s then bern 0.90 else bern 0.10)+ u <- unconditioned (if s then bern 0.75 else bern 0.25)+ if (numLoops > 1) then loop_hmm (numLoops - 1) u + else return s++test_hmm :: IO ()+test_hmm = sample (prog_hmm 2) conds >>=+ print . take 10 >>+ putChar '\n' >>+ empiricalMeasure 1000 (prog_hmm 2) conds >>=+ print + where conds = [Just (toDyn (Discrete True)),+ Just (toDyn (Discrete True))]++prog_carRoadModel :: Measure (Double, Double)+prog_carRoadModel = do+ speed <- unconditioned (uniform 5 15)+ let z0 = lit 0 + _ <- conditioned (normal z0 1)+ z1 <- unconditioned (normal (z0 + speed) 1)+ _ <- conditioned (normal z1 1)+ z2 <- unconditioned (normal (z1 + speed) 1) + _ <- conditioned (normal z2 1)+ z3 <- unconditioned (normal (z2 + speed) 1) + _ <- conditioned (normal z3 1)+ z4 <- unconditioned (normal (z3 + speed) 1) + return (z4, z3)++test_carRoadModel :: IO ()+test_carRoadModel = sample prog_carRoadModel conds >>=+ print . take 10 >>+ putChar '\n' >>+ empiricalMeasure 1000 prog_carRoadModel conds >>=+ print + where conds = [Just (toDyn (Lebesgue 0 :: Density Double)),+ Just (toDyn (Lebesgue 11 :: Density Double)), + Just (toDyn (Lebesgue 19 :: Density Double)),+ Just (toDyn (Lebesgue 33 :: Density Double))]++prog_categorical :: Measure Bool+prog_categorical = do + rain <- unconditioned (categorical [(True, 0.2), (False, 0.8)]) + sprinkler <- unconditioned (if rain+ then bern 0.01 else bern 0.4)+ _ <- conditioned (if rain+ then (if sprinkler then bern 0.99 else bern 0.8)+ else (if sprinkler then bern 0.90 else bern 0.1))+ return rain++test_categorical :: IO ()+test_categorical = sample prog_categorical conds >>=+ print . take 10 >>+ putChar '\n' >>+ empiricalMeasure 1000 prog_categorical conds >>=+ print + where conds = [Just (toDyn (Discrete True))]++prog_multiple_conditions :: Measure Double+prog_multiple_conditions = do+ b <- unconditioned (beta 1 1)+ _ <- conditioned (bern b)+ _ <- conditioned (bern b)+ return b
+ Tests/Metropolis.hs view
@@ -0,0 +1,45 @@+{-# LANGUAGE BangPatterns #-}++module Tests.Metropolis where++import Data.Dynamic++import Language.Hakaru.Types+import Language.Hakaru.Lambda+import Language.Hakaru.Metropolis+import Language.Hakaru.Distribution (bern, normal, uniform, beta)++-- import Test.QuickCheck.Monadic+-- import Distribution.TestSuite.QuickCheck+import Tests.Models++test_mixture :: IO [Bool]+test_mixture = mcmc prog_mixture [Just (toDyn (Lebesgue (-2) :: Density Double))]++test_multiple_conditions :: IO [Double]+test_multiple_conditions = do+ mcmc prog_multiple_conditions [Just (toDyn (Discrete True)),+ Just (toDyn (Discrete False))]++prog_two_normals :: Measure Bool+prog_two_normals = unconditioned (bern 0.5) `bind` \coin ->+ ifThenElse coin (conditioned (normal 0 1))+ (conditioned (normal 100 1)) `bind` \_ ->+ return coin++test_two_normals :: IO [Bool]+test_two_normals = mcmc prog_two_normals [Just (toDyn (Lebesgue 1 :: Density Double))]++test_normal :: IO [Double]+test_normal = mcmc (unconditioned (normal 1 3)) []++prog_joint = do+ bias <- unconditioned $ beta 1 1+ coin <- unconditioned $ bern bias+ return (bias, coin)++prog_condition = condition prog_joint True++test_condition :: IO [Double]+test_condition = mcmc prog_condition []+
+ Tests/Models.hs view
@@ -0,0 +1,26 @@+{-# LANGUAGE NoMonomorphismRestriction, GADTs #-}+module Tests.Models where++import Prelude hiding (Real)++import Language.Hakaru.Syntax2+import Language.Hakaru.Lambda+import qualified Language.Hakaru.Types as T++-- if we want to forgo the (D m) constraint, need to decorate the+-- program a little more.+prog_mixture :: (Meas m, D m ~ T.Dist ) => m Bool+prog_mixture = do+ c <- unconditioned (bern 0.5)+ _ <- conditioned (ifThenElse c (normal 1 1)+ (uniform 0 3))+ return c++prog_multiple_conditions :: (Meas m, D m ~ T.Dist) => m Double+prog_multiple_conditions = do+ b <- unconditioned (beta 1 1)+ _ <- conditioned (bern b)+ _ <- conditioned (bern b)+ return b++
+ Tests/Tests.hs view
@@ -0,0 +1,24 @@+module Main where++import Prelude hiding (Real)++import qualified Tests.ImportanceSampler as IS+import qualified Tests.Metropolis as MH+import qualified Tests.Distribution as D++import Language.Hakaru.Syntax2+import Language.Hakaru.Lambda+import qualified Language.Hakaru.Types as T++import Test.Framework (defaultMain, testGroup)+import Test.Framework.Providers.HUnit++import Test.HUnit++tests = [+ testGroup "Distribution checks" D.qtest, + testCase "alwaysPass" (1 @?= 1)+ --testCase "alwaysFail" (error "Fail!")+ ]++main = defaultMain tests
− Types.hs
@@ -1,13 +0,0 @@-{-# LANGUAGE RankNTypes, BangPatterns #-}-{-# OPTIONS -W #-}--module Types where--import Sampler (Sampler)--import Data.Dynamic---- Basic types for conditioning and conditioned sampler-data Cond = Unconditioned | Lebesgue !Dynamic | Discrete !Dynamic- deriving (Show)-newtype CSampler a = CSampler (Cond -> Sampler a)
− Util/Coda.hs
@@ -1,12 +0,0 @@-module Util.Coda where--import Statistics.Autocorrelation-import qualified Data.Packed.Vector as V-import qualified Data.Vector.Generic as G--effectiveSampleSize :: [Double] -> Double-effectiveSampleSize samples = n / (1 + 2*(G.sum rho))- where n = fromIntegral (V.dim vec)- vec = V.fromList samples- cov = autocovariance vec- rho = G.map (/ G.head cov) cov
− Util/Extras.hs
@@ -1,84 +0,0 @@-{-|- 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.--}--module Util.Extras where--import qualified Data.Sequence as S-import System.Random-import Data.Maybe-import qualified Data.Foldable as F--import Data.Dynamic-import Types--extract :: S.Seq a -> Int -> Maybe (S.Seq a, a)-extract s i | S.null r = Nothing- | otherwise = Just (a S.>< c, b)- 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- return $ extract s i--{-| - Given a sequence, return a *sorted* sequence of- n randomly selected elements from *distinct positions* in the sequence--}--randomElems :: Ord a => S.Seq a -> Int -> IO (S.Seq a)-randomElems = randomElemsTR S.empty--randomElemsTR :: Ord a => S.Seq a -> S.Seq a -> Int -> IO (S.Seq a)-randomElemsTR ixs s n- | n == S.length s = return $ S.unstableSort s- | n == 1 = do (_,i) <- fmap fromJust (randomExtract s)- return.S.unstableSort $ i S.<| ixs- | otherwise = do (s',i) <- fmap fromJust (randomExtract s)- (randomElemsTR $! (i S.<| ixs)) s' (n-1)--{-|- Chop a sequence at the given indices. - Assume number of indices given < length of sequence to be chopped--}--pieces :: S.Seq a -> S.Seq Int -> [S.Seq a]-pieces s ixs = let f (ps,r,x) y = let (p,r') = S.splitAt (y-x) r- in (p:ps,r',y)- g (a,b,_) = b:a- in g $ F.foldl f ([],s,0) ixs--{-|- Given n, chop a sequence at m random points- where m = min (length-1, n-1)--}--randomPieces :: Int -> S.Seq a -> IO [S.Seq a]-randomPieces n s- | n >= l = return $ F.toList $ fmap S.singleton s- | otherwise = do ixs <- randomElems (S.fromList [1..l-1]) (n-1)- return $ pieces s ixs- where l = S.length s--{-|- > pairs [1,2,3,4]- [(1,2),(1,3),(1,4),(2,3),(2,4),(3,4)]- > pairs [1,2,4,4]- [(1,2),(1,4),(1,4),(2,4),(2,4),(4,4)]--}--pairs :: [a] -> [(a,a)]-pairs [] = []-pairs (x:xs) = (zip (repeat x) xs) ++ pairs xs--l2Norm :: Floating a => [a] -> a-l2Norm l = sqrt.sum $ zipWith (*) l l--dataLoad [] = []-dataLoad (x:xs) = Lebesgue (toDyn (x :: Double)) : dataLoad xs
− Visual.hs
@@ -1,43 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}--module Visual where--import System.IO-import Control.Monad--import Data.Aeson-import Data.List-import qualified Data.Text as T-import qualified Data.ByteString.Lazy.Char8 as B-import qualified Data.ByteString.Char8 as BS--plot :: Show a => [a] -> String -> IO ()-plot samples filename = do- h <- openFile filename WriteMode- hPrint h samples- hClose h--batchPrint :: Show a => Int -> [a] -> IO ()-batchPrint n l = do- let batch = take n l- print batch- when (length batch == n) $ batchPrint n (drop n l)--viz :: ToJSON a => Int -> [String] -> [[a]] -> IO ()-viz n name samples = viz' n 50 0 name samples--viz' :: ToJSON a => Int -> Int -> Int -> [String] -> [[a]] -> IO ()-viz' n c cur name samples = do- putStrLn batch- when (c+cur < n) $- viz' n c (cur+c) name (drop c samples)- where- total = "total_samples" .= n- current_sample = "current_sample" .= cur- chunk = object (zipWith (\ name s -> T.pack name .= s)- name- (transpose $ take c samples))- batch = B.unpack $ encode- (object ["rvars" .= chunk,- total,- current_sample])
hakaru.cabal view
@@ -2,10 +2,10 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: hakaru-version: 0.1.2-synopsis: A probabilistic programming embedded DSL-description: Hakaru is an embedded DSL for performing probabilistic inference. It supports multiple inference backends.-homepage: http://www.indiana.edu/~ppaml+version: 0.1.3+synopsis: A probabilistic programming embedded DSL +-- description: +homepage: http://indiana.edu/~ppaml/ license: BSD3 license-file: LICENSE author: The Hakaru Team@@ -17,9 +17,83 @@ cabal-version: >=1.10 library- exposed-modules: Sampler, Types, Visual, Language.Hakaru.Syntax, Mixture, Language.Hakaru.Symbolic, RandomChoice, Util.Extras, Util.Coda, Examples.Tests, Language.Hakaru.ImportanceSampler, Language.Hakaru.Metropolis- -- other-modules: - other-extensions: RankNTypes, BangPatterns, OverloadedStrings, TypeFamilies, ConstraintKinds, GADTs, FlexibleContexts, TypeOperators, DataKinds, NoMonomorphismRestriction, DeriveDataTypeable, ScopedTypeVariables, ExistentialQuantification, StandaloneDeriving- build-depends: base >=4.6 && <4.7, aeson >=0.7 && <0.8, text >=1.1 && <1.2, bytestring >=0.10 && <0.11, pretty >=1.1 && <1.2, logfloat >=0.12 && <0.13, containers >=0.5 && <0.6, random >=1.0 && <1.1, math-functions >=0.1 && <0.2, vector >=0.10 && <0.11, cassava >=0.4 && <0.5, zlib >=0.5 && <0.6, statistics >=0.11 && <0.12, hmatrix >=0.16 && <0.17, parsec >=3.1 && <3.2+ exposed-modules: Language.Hakaru.Types,+ Language.Hakaru.Symbolic,+ Language.Hakaru.Arrow,+ Language.Hakaru.Mixture,+ Language.Hakaru.Sampler,+ Language.Hakaru.ImportanceSampler,+ Language.Hakaru.Metropolis,+ Language.Hakaru.Lambda,+ Language.Hakaru.Distribution,+ Language.Hakaru.Util.Csv,+ Language.Hakaru.Util.Extras,+ Language.Hakaru.Util.Visual,+ Language.Hakaru.Util.Coda+ other-extensions: RankNTypes, BangPatterns, GADTs, TypeFamilies, TypeOperators,+ ConstraintKinds, FlexibleContexts, NoMonomorphismRestriction,+ DeriveDataTypeable, ScopedTypeVariables, ExistentialQuantification,+ StandaloneDeriving, OverloadedStrings,+ FlexibleInstances, RebindableSyntax+ build-depends: base >=4.6 && <5.0,+ random >=1.0 && <1.1,+ transformers >=0.3 && <0.4,+ containers >=0.5 && <0.6,+ pretty >=1.1 && <1.2,+ logfloat >=0.12 && <0.13,+ hmatrix >=0.16 && <0.17,+ math-functions >=0.1 && <0.2,+ vector >=0.10 && <0.11,+ cassava >=0.4 && <0.5,+ zlib >=0.5 && <0.6,+ bytestring >=0.10 && <0.11,+ aeson >=0.7 && <0.8,+ text >=1.1 && <1.2,+ statistics >=0.11 && <0.14,+ parsec >=3.1 && <3.2,+ array >=0.4,+ mwc-random >=0.13 && <0.14,+ directory >=1.2 && <1.3,+ integration >= 0.2.0 && < 0.3.0,+ primitive >= 0.5 && < 0.6,+ parallel >=3.2 && <3.3,+ monad-loops >= 0.3.0.2 -- hs-source-dirs: default-language: Haskell2010+ ghc-options: -Wall++test-suite hakaru-test+ type: exitcode-stdio-1.0+ main-is: Tests/Tests.hs+ other-modules: Tests.ImportanceSampler,+ Tests.Models,+ Tests.Metropolis,+ Tests.Distribution+ build-depends: base, Cabal >= 1.10, + QuickCheck, + HUnit,+ test-framework,+ test-framework-quickcheck2,+ test-framework-hunit,+ random >=1.0 && <1.1,+ pretty >=1.1 && <1.2,+ containers >=0.5 && <0.6,+ logfloat >=0.12 && <0.13,+ math-functions >=0.1 && <0.2,+ statistics >=0.11 && <0.14,+ hmatrix >=0.16 && <0.17,+ vector >=0.10 && <0.11,+ hakaru >= 0.1.3,+ mwc-random >=0.13 && <0.14,+ primitive >= 0.5 && < 0.6,+ monad-loops >= 0.3.0.2+ default-language: Haskell2010++benchmark bench-all+ type: exitcode-stdio-1.0+ hs-source-dirs: Bench+ main-is: Bench.hs+ build-depends: base, deepseq, ghc-prim,+ criterion, hakaru >= 0.1.3+ ghc-options: -O2+ default-language: Haskell2010