monad-bayes 1.1.1 → 1.2.0

raw patch · 48 files changed

+674/−542 lines, 48 filesdep ~vectorPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: vector

API changes (from Hackage documentation)

Files

• CHANGELOG.md +5/−0
• benchmark/SSM.hs +7/−25
• benchmark/Single.hs +1/−55
• benchmark/Speed.hs +4/−4
• models/BetaBin.hs +5/−5
• models/ConjugatePriors.hs +4/−4
• models/Dice.hs +4/−4
• models/HMM.hs +7/−7
• models/Helper.hs +70/−0
• models/LDA.hs +4/−4
• models/LogReg.hs +2/−2
• models/NonlinearSSM.hs +2/−2
• models/NonlinearSSM/Algorithms.hs +56/−0
• models/Sprinkler.hs +2/−2
• monad-bayes.cabal +14/−3
• src/Control/Monad/Bayes/Class.hs +31/−31
• src/Control/Monad/Bayes/Density/Free.hs +20/−20
• src/Control/Monad/Bayes/Density/State.hs +13/−13
• src/Control/Monad/Bayes/Enumerator.hs +7/−7
• src/Control/Monad/Bayes/Inference/Lazy/MH.hs +7/−7
• src/Control/Monad/Bayes/Inference/Lazy/WIS.hs +5/−5
• src/Control/Monad/Bayes/Inference/MCMC.hs +8/−8
• src/Control/Monad/Bayes/Inference/PMMH.hs +12/−12
• src/Control/Monad/Bayes/Inference/RMSMC.hs +11/−11
• src/Control/Monad/Bayes/Inference/SMC.hs +6/−6
• src/Control/Monad/Bayes/Inference/SMC2.hs +12/−12
• src/Control/Monad/Bayes/Inference/TUI.hs +3/−3
• src/Control/Monad/Bayes/Integrator.hs +26/−24
• src/Control/Monad/Bayes/Population.hs +50/−54
• src/Control/Monad/Bayes/Sampler/Lazy.hs +16/−16
• src/Control/Monad/Bayes/Sampler/Strict.hs +17/−17
• src/Control/Monad/Bayes/Sequential/Coroutine.hs +21/−21
• src/Control/Monad/Bayes/Traced/Basic.hs +25/−25
• src/Control/Monad/Bayes/Traced/Common.hs +12/−12
• src/Control/Monad/Bayes/Traced/Dynamic.hs +29/−29
• src/Control/Monad/Bayes/Traced/Static.hs +27/−27
• src/Control/Monad/Bayes/Weighted.hs +21/−23
• src/Math/Integrators/StormerVerlet.hs +1/−1
• test/Spec.hs +6/−1
• test/TestAdvanced.hs +6/−6
• test/TestBenchmarks.hs +33/−0
• test/TestEnumerator.hs +2/−2
• test/TestInference.hs +14/−14
• test/TestIntegrator.hs +5/−5
• test/TestPopulation.hs +5/−5
• test/TestSSMFixtures.hs +28/−0
• test/TestSequential.hs +4/−4
• test/TestWeighted.hs +4/−4

@@ -1,3 +1,8 @@+# 1.2.0++- Renamed monad transformers idiomatically+  (https://github.com/tweag/monad-bayes/pull/295)+ # 1.1.1  - add fixture tests for benchmark models

@@ -1,43 +1,25 @@ module Main where +import Control.Monad (forM_) import Control.Monad.Bayes.Inference.MCMC import Control.Monad.Bayes.Inference.PMMH as PMMH (pmmh) import Control.Monad.Bayes.Inference.RMSMC (rmsmcDynamic) import Control.Monad.Bayes.Inference.SMC import Control.Monad.Bayes.Inference.SMC2 as SMC2 (smc2) import Control.Monad.Bayes.Population-import Control.Monad.Bayes.Population (population, resampleMultinomial)+import Control.Monad.Bayes.Population (resampleMultinomial, runPopulationT) import Control.Monad.Bayes.Sampler.Strict (sampleIO, sampleIOfixed, sampleWith) import Control.Monad.Bayes.Weighted (unweighted) import Control.Monad.IO.Class (MonadIO (liftIO)) import NonlinearSSM (generateData, model, param)+import NonlinearSSM.Algorithms import System.Random.Stateful (mkStdGen, newIOGenM)  main :: IO () main = sampleIOfixed $ do-  let t = 5   dat <- generateData t   let ys = map snd dat-  liftIO $ print "SMC"-  smcRes <- population $ smc SMCConfig {numSteps = t, numParticles = 10, resampler = resampleMultinomial} (param >>= model ys)-  liftIO $ print $ show smcRes-  liftIO $ print "RM-SMC"-  smcrmRes <--    population $-      rmsmcDynamic-        MCMCConfig {numMCMCSteps = 10, numBurnIn = 0, proposal = SingleSiteMH}-        SMCConfig {numSteps = t, numParticles = 10, resampler = resampleSystematic}-        (param >>= model ys)-  liftIO $ print $ show smcrmRes-  liftIO $ print "PMMH"-  pmmhRes <--    unweighted $-      pmmh-        MCMCConfig {numMCMCSteps = 2, numBurnIn = 0, proposal = SingleSiteMH}-        SMCConfig {numSteps = t, numParticles = 3, resampler = resampleSystematic}-        param-        (model ys)-  liftIO $ print $ show pmmhRes-  liftIO $ print "SMC2"-  smc2Res <- population $ smc2 t 3 2 1 param (model ys)-  liftIO $ print $ show smc2Res+  forM_ [SMC, RMSMCDynamic, PMMH, SMC2] $ \alg -> do+    liftIO $ print alg+    result <- runAlgFixed ys alg+    liftIO $ putStrLn result

@@ -1,22 +1,9 @@ {-# LANGUAGE DerivingStrategies #-} {-# LANGUAGE ImportQualifiedPost #-} -import Control.Monad.Bayes.Class (MonadMeasure)-import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (..), Proposal (SingleSiteMH))-import Control.Monad.Bayes.Inference.RMSMC (rmsmcBasic)-import Control.Monad.Bayes.Inference.SMC-  ( SMCConfig (SMCConfig, numParticles, numSteps, resampler),-    smc,-  )-import Control.Monad.Bayes.Population import Control.Monad.Bayes.Sampler.Strict-import Control.Monad.Bayes.Traced hiding (model)-import Control.Monad.Bayes.Weighted-import Control.Monad.ST (runST) import Data.Time (diffUTCTime, getCurrentTime)-import HMM qualified-import LDA qualified-import LogReg qualified+import Helper import Options.Applicative   ( Applicative (liftA2),     ParserInfo,@@ -30,47 +17,6 @@     option,     short,   )--data Model = LR Int | HMM Int | LDA (Int, Int)-  deriving stock (Show, Read)--parseModel :: String -> Maybe Model-parseModel s =-  case s of-    'L' : 'R' : n -> Just $ LR (read n)-    'H' : 'M' : 'M' : n -> Just $ HMM (read n)-    'L' : 'D' : 'A' : n -> Just $ LDA (5, read n)-    _ -> Nothing--getModel :: MonadMeasure m => Model -> (Int, m String)-getModel model = (size model, program model)-  where-    size (LR n) = n-    size (HMM n) = n-    size (LDA (d, w)) = d * w-    program (LR n) = show <$> (LogReg.logisticRegression (runST $ sampleSTfixed (LogReg.syntheticData n)))-    program (HMM n) = show <$> (HMM.hmm (runST $ sampleSTfixed (HMM.syntheticData n)))-    program (LDA (d, w)) = show <$> (LDA.lda (runST $ sampleSTfixed (LDA.syntheticData d w)))--data Alg = SMC | MH | RMSMC-  deriving stock (Read, Show)--runAlg :: Model -> Alg -> SamplerIO String-runAlg model alg =-  case alg of-    SMC ->-      let n = 100-          (k, m) = getModel model-       in show <$> population (smc SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic} m)-    MH ->-      let t = 100-          (_, m) = getModel model-       in show <$> unweighted (mh t m)-    RMSMC ->-      let n = 10-          t = 1-          (k, m) = getModel model-       in show <$> population (rmsmcBasic MCMCConfig {numMCMCSteps = t, numBurnIn = 0, proposal = SingleSiteMH} (SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic}) m)  infer :: Model -> Alg -> IO () infer model alg = do

@@ -8,7 +8,7 @@ import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (MCMCConfig, numBurnIn, numMCMCSteps, proposal), Proposal (SingleSiteMH)) import Control.Monad.Bayes.Inference.RMSMC (rmsmcDynamic) import Control.Monad.Bayes.Inference.SMC (SMCConfig (SMCConfig, numParticles, numSteps, resampler), smc)-import Control.Monad.Bayes.Population (population, resampleSystematic)+import Control.Monad.Bayes.Population (resampleSystematic, runPopulationT) import Control.Monad.Bayes.Sampler.Strict (SamplerIO, sampleIOfixed) import Control.Monad.Bayes.Traced (mh) import Control.Monad.Bayes.Weighted (unweighted)@@ -40,7 +40,7 @@   show (HMM xs) = "HMM" ++ show (length xs)   show (LDA xs) = "LDA" ++ show (length $ head xs) -buildModel :: MonadMeasure m => Model -> m String+buildModel :: (MonadMeasure m) => Model -> m String buildModel (LR dataset) = show <$> LogReg.logisticRegression dataset buildModel (HMM dataset) = show <$> HMM.hmm dataset buildModel (LDA dataset) = show <$> LDA.lda dataset@@ -59,10 +59,10 @@  runAlg :: Model -> Alg -> SamplerIO String runAlg model (MH n) = show <$> unweighted (mh n (buildModel model))-runAlg model (SMC n) = show <$> population (smc SMCConfig {numSteps = (modelLength model), numParticles = n, resampler = resampleSystematic} (buildModel model))+runAlg model (SMC n) = show <$> runPopulationT (smc SMCConfig {numSteps = (modelLength model), numParticles = n, resampler = resampleSystematic} (buildModel model)) runAlg model (RMSMC n t) =   show-    <$> population+    <$> runPopulationT       ( rmsmcDynamic           MCMCConfig {numMCMCSteps = t, numBurnIn = 0, proposal = SingleSiteMH}           SMCConfig {numSteps = modelLength model, numParticles = n, resampler = resampleSystematic}

@@ -17,14 +17,14 @@ import Pipes.Prelude qualified as P hiding (show)  -- | Beta-binomial model as an i.i.d. sequence conditionally on weight.-latent :: MonadDistribution m => Int -> m [Bool]+latent :: (MonadDistribution m) => Int -> m [Bool] latent n = do   weight <- uniform 0 1   replicateM n (bernoulli weight)  -- | Beta-binomial as a random process. -- Equivalent to the above by De Finetti's theorem.-urn :: MonadDistribution m => Int -> m [Bool]+urn :: (MonadDistribution m) => Int -> m [Bool] urn n = flip evalStateT (1, 1) $ do   replicateM n do     (a, b) <- get@@ -36,7 +36,7 @@  -- | Beta-binomial as a random process. -- This time using the Pipes library, for a more pure functional style-urnP :: MonadDistribution m => Int -> m [Bool]+urnP :: (MonadDistribution m) => Int -> m [Bool] urnP n = P.toListM $ P.take n <-< P.unfoldr toss (1, 1)   where     toss (a, b) = do@@ -47,7 +47,7 @@  -- | A beta-binomial model where the first three states are True,True,False. -- The resulting distribution is on the remaining outcomes.-cond :: MonadMeasure m => m [Bool] -> m [Bool]+cond :: (MonadMeasure m) => m [Bool] -> m [Bool] cond d = do   ~(first : second : third : rest) <- d   condition first@@ -56,7 +56,7 @@   return rest  -- | The final conditional model, abstracting the representation.-model :: MonadMeasure m => (Int -> m [Bool]) -> Int -> m Int+model :: (MonadMeasure m) => (Int -> m [Bool]) -> Int -> m Int model repr n = fmap count $ cond $ repr (n + 3)   where     -- Post-processing by counting the number of True values.

@@ -41,16 +41,16 @@     a' = a + s     b' = b + fromIntegral n - s -bernoulliPdf :: Floating a => a -> Bool -> Log a+bernoulliPdf :: (Floating a) => a -> Bool -> Log a bernoulliPdf p x = let numBool = if x then 1.0 else 0 in Exp $ log (p ** numBool * (1 - p) ** (1 - numBool)) -betaBernoulli' :: MonadMeasure m => (Double, Double) -> Bayesian m Double Bool+betaBernoulli' :: (MonadMeasure m) => (Double, Double) -> Bayesian m Double Bool betaBernoulli' (a, b) = Bayesian (beta a b) bernoulli bernoulliPdf -normalNormal' :: MonadMeasure m => Double -> (Double, Double) -> Bayesian m Double Double+normalNormal' :: (MonadMeasure m) => Double -> (Double, Double) -> Bayesian m Double Double normalNormal' var (mu0, var0) = Bayesian (normal mu0 (sqrt var0)) (`normal` (sqrt var)) (`normalPdf` (sqrt var)) -gammaNormal' :: MonadMeasure m => (Double, Double) -> Bayesian m Double Double+gammaNormal' :: (MonadMeasure m) => (Double, Double) -> Bayesian m Double Double gammaNormal' (a, b) = Bayesian (gamma a (recip b)) (normal 0 . sqrt . recip) (normalPdf 0 . sqrt . recip)  normalNormalAnalytic ::

@@ -12,23 +12,23 @@   )  -- | A toss of a six-sided die.-die :: MonadDistribution m => m Int+die :: (MonadDistribution m) => m Int die = uniformD [1 .. 6]  -- | A sum of outcomes of n independent tosses of six-sided dice.-dice :: MonadDistribution m => Int -> m Int+dice :: (MonadDistribution m) => Int -> m Int dice 1 = die dice n = liftA2 (+) die (dice (n - 1))  -- | Toss of two dice where the output is greater than 4.-diceHard :: MonadMeasure m => m Int+diceHard :: (MonadMeasure m) => m Int diceHard = do   result <- dice 2   condition (result > 4)   return result  -- | Toss of two dice with an artificial soft constraint.-diceSoft :: MonadMeasure m => m Int+diceSoft :: (MonadMeasure m) => m Int diceSoft = do   result <- dice 2   score (1 / fromIntegral result)

@@ -39,7 +39,7 @@   ]  -- | The transition model.-trans :: MonadDistribution m => Int -> m Int+trans :: (MonadDistribution m) => Int -> m Int trans 0 = categorical $ fromList [0.1, 0.4, 0.5] trans 1 = categorical $ fromList [0.2, 0.6, 0.2] trans 2 = categorical $ fromList [0.15, 0.7, 0.15]@@ -53,7 +53,7 @@ emissionMean _ = error "unreachable"  -- | Initial state distribution-start :: MonadDistribution m => m Int+start :: (MonadDistribution m) => m Int start = uniformD [0, 1, 2]  -- | Example HMM from http://dl.acm.org/citation.cfm?id=2804317@@ -67,7 +67,7 @@     f [] k = start >>= k []     f (y : ys) k = f ys (\xs x -> expand x y >>= k (x : xs)) -syntheticData :: MonadDistribution m => Int -> m [Double]+syntheticData :: (MonadDistribution m) => Int -> m [Double] syntheticData n = replicateM n syntheticPoint   where     syntheticPoint = uniformD [0, 1, 2]@@ -75,14 +75,14 @@ -- | Equivalent model, but using pipes for simplicity  -- | Prior expressed as a stream-hmmPrior :: MonadDistribution m => Producer Int m b+hmmPrior :: (MonadDistribution m) => Producer Int m b hmmPrior = do   x <- lift start   yield x   Pipes.unfoldr (fmap (Right . (\k -> (k, k))) . trans) x  -- | Observations expressed as a stream-hmmObservations :: Functor m => [a] -> Producer (Maybe a) m ()+hmmObservations :: (Functor m) => [a] -> Producer (Maybe a) m () hmmObservations dataset = each (Nothing : (Just <$> reverse dataset))  -- | Posterior expressed as a stream@@ -93,12 +93,12 @@     hmmPrior     (hmmObservations dataset)   where-    hmmLikelihood :: MonadFactor f => (Int, Maybe Double) -> f ()+    hmmLikelihood :: (MonadFactor f) => (Int, Maybe Double) -> f ()     hmmLikelihood (l, o) = when (isJust o) (factor $ normalPdf (emissionMean l) 1 (fromJust o))      zipWithM f p1 p2 = Pipes.zip p1 p2 >-> Pipes.chain f >-> Pipes.map fst -hmmPosteriorPredictive :: MonadDistribution m => [Double] -> Producer Double m ()+hmmPosteriorPredictive :: (MonadDistribution m) => [Double] -> Producer Double m () hmmPosteriorPredictive dataset =   Pipes.hoist enumerateToDistribution (hmmPosterior dataset)     >-> Pipes.mapM (\x -> normal (emissionMean x) 1)

@@ -0,0 +1,70 @@+{-# LANGUAGE DerivingStrategies #-}+{-# LANGUAGE ImportQualifiedPost #-}++module Helper where++import Control.Monad.Bayes.Class (MonadMeasure)+import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (..), Proposal (SingleSiteMH))+import Control.Monad.Bayes.Inference.RMSMC (rmsmcBasic)+import Control.Monad.Bayes.Inference.SMC+  ( SMCConfig (SMCConfig, numParticles, numSteps, resampler),+    smc,+  )+import Control.Monad.Bayes.Population+import Control.Monad.Bayes.Sampler.Strict+import Control.Monad.Bayes.Traced hiding (model)+import Control.Monad.Bayes.Weighted+import Control.Monad.ST (runST)+import HMM qualified+import LDA qualified+import LogReg qualified++data Model = LR Int | HMM Int | LDA (Int, Int)+  deriving stock (Show, Read)++parseModel :: String -> Maybe Model+parseModel s =+  case s of+    'L' : 'R' : n -> Just $ LR (read n)+    'H' : 'M' : 'M' : n -> Just $ HMM (read n)+    'L' : 'D' : 'A' : n -> Just $ LDA (5, read n)+    _ -> Nothing++serializeModel :: Model -> Maybe String+serializeModel (LR n) = Just $ "LR" ++ show n+serializeModel (HMM n) = Just $ "HMM" ++ show n+serializeModel (LDA (5, n)) = Just $ "LDA" ++ show n+serializeModel (LDA _) = Nothing++data Alg = SMC | MH | RMSMC+  deriving stock (Read, Show, Eq, Ord, Enum, Bounded)++getModel :: (MonadMeasure m) => Model -> (Int, m String)+getModel model = (size model, program model)+  where+    size (LR n) = n+    size (HMM n) = n+    size (LDA (d, w)) = d * w+    program (LR n) = show <$> (LogReg.logisticRegression (runST $ sampleSTfixed (LogReg.syntheticData n)))+    program (HMM n) = show <$> (HMM.hmm (runST $ sampleSTfixed (HMM.syntheticData n)))+    program (LDA (d, w)) = show <$> (LDA.lda (runST $ sampleSTfixed (LDA.syntheticData d w)))++runAlg :: Model -> Alg -> SamplerIO String+runAlg model alg =+  case alg of+    SMC ->+      let n = 100+          (k, m) = getModel model+       in show <$> runPopulationT (smc SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic} m)+    MH ->+      let t = 100+          (_, m) = getModel model+       in show <$> unweighted (mh t m)+    RMSMC ->+      let n = 10+          t = 1+          (k, m) = getModel model+       in show <$> runPopulationT (rmsmcBasic MCMCConfig {numMCMCSteps = t, numBurnIn = 0, proposal = SingleSiteMH} (SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic}) m)++runAlgFixed :: Model -> Alg -> IO String+runAlgFixed model alg = sampleIOfixed $ runAlg model alg

@@ -43,17 +43,17 @@     words "bear wolf bear python bear wolf bear wolf bear wolf"   ] -wordDistPrior :: MonadDistribution m => m (V.Vector Double)+wordDistPrior :: (MonadDistribution m) => m (V.Vector Double) wordDistPrior = dirichlet $ V.replicate (length vocabulary) 1 -topicDistPrior :: MonadDistribution m => m (V.Vector Double)+topicDistPrior :: (MonadDistribution m) => m (V.Vector Double) topicDistPrior = dirichlet $ V.replicate (length topics) 1  wordIndex :: Map.Map Text Int wordIndex = Map.fromList $ zip vocabulary [0 ..]  lda ::-  MonadMeasure m =>+  (MonadMeasure m) =>   Documents ->   m (Map.Map Text (V.Vector (Text, Double)), [(Text, V.Vector (Text, Double))]) lda docs = do@@ -73,7 +73,7 @@       zip (fmap (foldr1 (\x y -> x <> " " <> y)) docs) (fmap (V.zip $ V.fromList ["topic1", "topic2"]) td)     ) -syntheticData :: MonadDistribution m => Int -> Int -> m [[Text]]+syntheticData :: (MonadDistribution m) => Int -> Int -> m [[Text]] syntheticData d w = List.replicateM d (List.replicateM w syntheticWord)   where     syntheticWord = uniformD vocabulary

@@ -13,7 +13,7 @@   ) import Numeric.Log (Log (Exp)) -logisticRegression :: MonadMeasure m => [(Double, Bool)] -> m Double+logisticRegression :: (MonadMeasure m) => [(Double, Bool)] -> m Double logisticRegression dat = do   m <- normal 0 1   b <- normal 0 1@@ -27,7 +27,7 @@   sigmoid 8  -- make a synthetic dataset by randomly choosing input-label pairs-syntheticData :: MonadDistribution m => Int -> m [(Double, Bool)]+syntheticData :: (MonadDistribution m) => Int -> m [(Double, Bool)] syntheticData n = replicateM n do   x <- uniform (-1) 1   label <- bernoulli 0.5

@@ -7,7 +7,7 @@     normalPdf,   ) -param :: MonadDistribution m => m (Double, Double)+param :: (MonadDistribution m) => m (Double, Double) param = do   let a = 0.01   let b = 0.01@@ -43,7 +43,7 @@   return $ reverse xs  generateData ::-  MonadDistribution m =>+  (MonadDistribution m) =>   -- | T   Int ->   -- | list of latent and observable states from t=1

@@ -0,0 +1,56 @@+module NonlinearSSM.Algorithms where++import Control.Monad.Bayes.Class (MonadDistribution)+import Control.Monad.Bayes.Inference.MCMC+import Control.Monad.Bayes.Inference.PMMH as PMMH (pmmh)+import Control.Monad.Bayes.Inference.RMSMC (rmsmc, rmsmcBasic, rmsmcDynamic)+import Control.Monad.Bayes.Inference.SMC+import Control.Monad.Bayes.Inference.SMC2 as SMC2 (smc2)+import Control.Monad.Bayes.Population+import Control.Monad.Bayes.Weighted (unweighted)+import NonlinearSSM++data Alg = SMC | RMSMC | RMSMCDynamic | RMSMCBasic | PMMH | SMC2+  deriving (Show, Read, Eq, Ord, Enum, Bounded)++algs :: [Alg]+algs = [minBound .. maxBound]++type SSMData = [Double]++t :: Int+t = 5++-- FIXME refactor such that it can be reused in ssm benchmark+runAlgFixed :: (MonadDistribution m) => SSMData -> Alg -> m String+runAlgFixed ys SMC = fmap show $ runPopulationT $ smc SMCConfig {numSteps = t, numParticles = 10, resampler = resampleMultinomial} (param >>= model ys)+runAlgFixed ys RMSMC =+  fmap show $+    runPopulationT $+      rmsmc+        MCMCConfig {numMCMCSteps = 10, numBurnIn = 0, proposal = SingleSiteMH}+        SMCConfig {numSteps = t, numParticles = 10, resampler = resampleSystematic}+        (param >>= model ys)+runAlgFixed ys RMSMCBasic =+  fmap show $+    runPopulationT $+      rmsmcBasic+        MCMCConfig {numMCMCSteps = 10, numBurnIn = 0, proposal = SingleSiteMH}+        SMCConfig {numSteps = t, numParticles = 10, resampler = resampleSystematic}+        (param >>= model ys)+runAlgFixed ys RMSMCDynamic =+  fmap show $+    runPopulationT $+      rmsmcDynamic+        MCMCConfig {numMCMCSteps = 10, numBurnIn = 0, proposal = SingleSiteMH}+        SMCConfig {numSteps = t, numParticles = 10, resampler = resampleSystematic}+        (param >>= model ys)+runAlgFixed ys PMMH =+  fmap show $+    unweighted $+      pmmh+        MCMCConfig {numMCMCSteps = 2, numBurnIn = 0, proposal = SingleSiteMH}+        SMCConfig {numSteps = t, numParticles = 3, resampler = resampleSystematic}+        param+        (model ys)+runAlgFixed ys SMC2 = fmap show $ runPopulationT $ smc2 t 3 2 1 param (model ys)

@@ -3,7 +3,7 @@ import Control.Monad (when) import Control.Monad.Bayes.Class -hard :: MonadMeasure m => m Bool+hard :: (MonadMeasure m) => m Bool hard = do   rain <- bernoulli 0.3   sprinkler <- bernoulli $ if rain then 0.1 else 0.4@@ -15,7 +15,7 @@   condition (not wet)   return rain -soft :: MonadMeasure m => m Bool+soft :: (MonadMeasure m) => m Bool soft = do   rain <- bernoulli 0.3   when rain (factor 0.2)

@@ -1,6 +1,6 @@ cabal-version:      2.2 name:               monad-bayes-version:            1.1.1+version:            1.2.0 license:            MIT license-file:       LICENSE.md copyright:          2015-2020 Adam Scibior@@ -61,7 +61,7 @@     , scientific       ^>=0.3     , statistics       >=0.14.0  && <0.17     , text             >=1.2     && <2.1-    , vector           ^>=0.12.0+    , vector           >=0.12.0  && <0.14     , vty              ^>=5.38  common test-deps@@ -131,6 +131,7 @@   hs-source-dirs:     benchmark models   other-modules:     Dice+    Helper     HMM     LDA     LogReg@@ -161,9 +162,15 @@   other-modules:     BetaBin     ConjugatePriors+    Helper     HMM+    LDA+    LogReg+    NonlinearSSM+    NonlinearSSM.Algorithms     Sprinkler     TestAdvanced+    TestBenchmarks     TestDistribution     TestEnumerator     TestInference@@ -172,6 +179,7 @@     TestPopulation     TestSampler     TestSequential+    TestSSMFixtures     TestStormerVerlet     TestWeighted @@ -198,7 +206,10 @@   type:             exitcode-stdio-1.0   main-is:          SSM.hs   hs-source-dirs:   models benchmark-  other-modules:    NonlinearSSM+  other-modules:+    NonlinearSSM+    NonlinearSSM.Algorithms+   default-language: Haskell2010   build-depends:     , base

@@ -107,7 +107,7 @@ import Statistics.Distribution.Uniform (uniformDistr)  -- | Monads that can draw random variables.-class Monad m => MonadDistribution m where+class (Monad m) => MonadDistribution m where   -- | Draw from a uniform distribution.   random ::     -- | \(\sim \mathcal{U}(0, 1)\)@@ -163,7 +163,7 @@    -- | Draw from a categorical distribution.   categorical ::-    Vector v Double =>+    (Vector v Double) =>     -- | event probabilities     v Double ->     -- | outcome category@@ -208,7 +208,7 @@    -- | Draw from a Dirichlet distribution.   dirichlet ::-    Vector v Double =>+    (Vector v Double) =>     -- | concentration parameters @as@     v Double ->     -- | \(\sim \mathrm{Dir}(\mathrm{as})\)@@ -226,7 +226,7 @@  -- | Draw from a discrete distribution using a sequence of draws from -- Bernoulli.-fromPMF :: MonadDistribution m => (Int -> Double) -> m Int+fromPMF :: (MonadDistribution m) => (Int -> Double) -> m Int fromPMF p = f 0 1   where     f i r = do@@ -241,7 +241,7 @@ discrete = fromPMF . probability  -- | Monads that can score different execution paths.-class Monad m => MonadFactor m where+class (Monad m) => MonadFactor m where   -- | Record a likelihood.   score ::     -- | likelihood of the execution path@@ -250,7 +250,7 @@  -- | Synonym for 'score'. factor ::-  MonadFactor m =>+  (MonadFactor m) =>   -- | likelihood of the execution path   Log Double ->   m ()@@ -258,17 +258,17 @@  -- | synonym for pretty type signatures, but note that (A -> Distribution B) won't work as intended: for that, use Kernel -- Also note that the use of RankNTypes means performance may take a hit: really the main point of these signatures is didactic-type Distribution a = forall m. MonadDistribution m => m a+type Distribution a = forall m. (MonadDistribution m) => m a -type Measure a = forall m. MonadMeasure m => m a+type Measure a = forall m. (MonadMeasure m) => m a -type Kernel a b = forall m. MonadMeasure m => a -> m b+type Kernel a b = forall m. (MonadMeasure m) => a -> m b  -- | Hard conditioning.-condition :: MonadFactor m => Bool -> m ()+condition :: (MonadFactor m) => Bool -> m () condition b = score $ if b then 1 else 0 -independent :: Applicative m => Int -> m a -> m [a]+independent :: (Applicative m) => Int -> m a -> m [a] independent = replicateM  -- | Monads that support both sampling and scoring.@@ -290,7 +290,7 @@ poissonPdf rate n = Exp $ logProbability (Poisson.poisson rate) (fromIntegral n)  -- | multivariate normal-mvNormal :: MonadDistribution m => V.Vector Double -> Matrix Double -> m (V.Vector Double)+mvNormal :: (MonadDistribution m) => V.Vector Double -> Matrix Double -> m (V.Vector Double) mvNormal mu bigSigma = do   let n = length mu   ss <- replicateM n (normal 0 1)@@ -312,7 +312,7 @@   factor $ product $ fmap (likelihood z) os   return z -priorPredictive :: Monad m => Bayesian m a b -> m b+priorPredictive :: (Monad m) => Bayesian m a b -> m b priorPredictive bm = prior bm >>= generative bm  posteriorPredictive ::@@ -342,32 +342,32 @@ ---------------------------------------------------------------------------- -- Instances that lift probabilistic effects to standard tranformers. -instance MonadDistribution m => MonadDistribution (IdentityT m) where+instance (MonadDistribution m) => MonadDistribution (IdentityT m) where   random = lift random   bernoulli = lift . bernoulli -instance MonadFactor m => MonadFactor (IdentityT m) where+instance (MonadFactor m) => MonadFactor (IdentityT m) where   score = lift . score -instance MonadMeasure m => MonadMeasure (IdentityT m)+instance (MonadMeasure m) => MonadMeasure (IdentityT m) -instance MonadDistribution m => MonadDistribution (ExceptT e m) where+instance (MonadDistribution m) => MonadDistribution (ExceptT e m) where   random = lift random   uniformD = lift . uniformD -instance MonadFactor m => MonadFactor (ExceptT e m) where+instance (MonadFactor m) => MonadFactor (ExceptT e m) where   score = lift . score -instance MonadMeasure m => MonadMeasure (ExceptT e m)+instance (MonadMeasure m) => MonadMeasure (ExceptT e m) -instance MonadDistribution m => MonadDistribution (ReaderT r m) where+instance (MonadDistribution m) => MonadDistribution (ReaderT r m) where   random = lift random   bernoulli = lift . bernoulli -instance MonadFactor m => MonadFactor (ReaderT r m) where+instance (MonadFactor m) => MonadFactor (ReaderT r m) where   score = lift . score -instance MonadMeasure m => MonadMeasure (ReaderT r m)+instance (MonadMeasure m) => MonadMeasure (ReaderT r m)  instance (Monoid w, MonadDistribution m) => MonadDistribution (WriterT w m) where   random = lift random@@ -379,31 +379,31 @@  instance (Monoid w, MonadMeasure m) => MonadMeasure (WriterT w m) -instance MonadDistribution m => MonadDistribution (StateT s m) where+instance (MonadDistribution m) => MonadDistribution (StateT s m) where   random = lift random   bernoulli = lift . bernoulli   categorical = lift . categorical   uniformD = lift . uniformD -instance MonadFactor m => MonadFactor (StateT s m) where+instance (MonadFactor m) => MonadFactor (StateT s m) where   score = lift . score -instance MonadMeasure m => MonadMeasure (StateT s m)+instance (MonadMeasure m) => MonadMeasure (StateT s m) -instance MonadDistribution m => MonadDistribution (ListT m) where+instance (MonadDistribution m) => MonadDistribution (ListT m) where   random = lift random   bernoulli = lift . bernoulli   categorical = lift . categorical -instance MonadFactor m => MonadFactor (ListT m) where+instance (MonadFactor m) => MonadFactor (ListT m) where   score = lift . score -instance MonadMeasure m => MonadMeasure (ListT m)+instance (MonadMeasure m) => MonadMeasure (ListT m) -instance MonadDistribution m => MonadDistribution (ContT r m) where+instance (MonadDistribution m) => MonadDistribution (ContT r m) where   random = lift random -instance MonadFactor m => MonadFactor (ContT r m) where+instance (MonadFactor m) => MonadFactor (ContT r m) where   score = lift . score -instance MonadMeasure m => MonadMeasure (ContT r m)+instance (MonadMeasure m) => MonadMeasure (ContT r m)

@@ -12,13 +12,13 @@ -- Stability   : experimental -- Portability : GHC ----- 'Density' is a free monad transformer over random sampling.+-- 'DensityT' is a free monad transformer over random sampling. module Control.Monad.Bayes.Density.Free-  ( Density,+  ( DensityT (..),     hoist,     interpret,     withRandomness,-    density,+    runDensityT,     traced,   ) where@@ -36,40 +36,40 @@ -- | Free monad transformer over random sampling. -- -- Uses the Church-encoded version of the free monad for efficiency.-newtype Density m a = Density {runDensity :: FT SamF m a}+newtype DensityT m a = DensityT {getDensityT :: FT SamF m a}   deriving newtype (Functor, Applicative, Monad, MonadTrans) -instance MonadFree SamF (Density m) where-  wrap = Density . wrap . fmap runDensity+instance MonadFree SamF (DensityT m) where+  wrap = DensityT . wrap . fmap getDensityT -instance Monad m => MonadDistribution (Density m) where-  random = Density $ liftF (Random id)+instance (Monad m) => MonadDistribution (DensityT m) where+  random = DensityT $ liftF (Random id) --- | Hoist 'Density' through a monad transform.-hoist :: (Monad m, Monad n) => (forall x. m x -> n x) -> Density m a -> Density n a-hoist f (Density m) = Density (hoistFT f m)+-- | Hoist 'DensityT' through a monad transform.+hoist :: (Monad m, Monad n) => (forall x. m x -> n x) -> DensityT m a -> DensityT n a+hoist f (DensityT m) = DensityT (hoistFT f m)  -- | Execute random sampling in the transformed monad.-interpret :: MonadDistribution m => Density m a -> m a-interpret (Density m) = iterT f m+interpret :: (MonadDistribution m) => DensityT m a -> m a+interpret (DensityT m) = iterT f m   where     f (Random k) = random >>= k  -- | Execute computation with supplied values for random choices.-withRandomness :: Monad m => [Double] -> Density m a -> m a-withRandomness randomness (Density m) = evalStateT (iterTM f m) randomness+withRandomness :: (Monad m) => [Double] -> DensityT m a -> m a+withRandomness randomness (DensityT m) = evalStateT (iterTM f m) randomness   where     f (Random k) = do       xs <- get       case xs of-        [] -> error "Density: the list of randomness was too short"+        [] -> error "DensityT: the list of randomness was too short"         y : ys -> put ys >> k y  -- | Execute computation with supplied values for a subset of random choices. -- Return the output value and a record of all random choices used, whether -- taken as input or drawn using the transformed monad.-density :: MonadDistribution m => [Double] -> Density m a -> m (a, [Double])-density randomness (Density m) =+runDensityT :: (MonadDistribution m) => [Double] -> DensityT m a -> m (a, [Double])+runDensityT randomness (DensityT m) =   runWriterT $ evalStateT (iterTM f $ hoistFT lift m) randomness   where     f (Random k) = do@@ -84,5 +84,5 @@       k x  -- | Like 'density', but use an arbitrary sampling monad.-traced :: MonadDistribution m => [Double] -> Density Identity a -> m (a, [Double])-traced randomness m = density randomness $ hoist (return . runIdentity) m+traced :: (MonadDistribution m) => [Double] -> DensityT Identity a -> m (a, [Double])+traced randomness m = runDensityT randomness $ hoist (return . runIdentity) m

@@ -13,21 +13,21 @@ import Control.Monad.State (MonadState (get, put), StateT, evalStateT) import Control.Monad.Writer -newtype Density m a = Density {runDensity :: WriterT [Double] (StateT [Double] m) a} deriving newtype (Functor, Applicative, Monad)+newtype DensityT m a = DensityT {getDensityT :: WriterT [Double] (StateT [Double] m) a} deriving newtype (Functor, Applicative, Monad) -instance MonadTrans Density where-  lift = Density . lift . lift+instance MonadTrans DensityT where+  lift = DensityT . lift . lift -instance Monad m => MonadState [Double] (Density m) where-  get = Density $ lift $ get-  put = Density . lift . put+instance (Monad m) => MonadState [Double] (DensityT m) where+  get = DensityT $ lift $ get+  put = DensityT . lift . put -instance Monad m => MonadWriter [Double] (Density m) where-  tell = Density . tell-  listen = Density . listen . runDensity-  pass = Density . pass . runDensity+instance (Monad m) => MonadWriter [Double] (DensityT m) where+  tell = DensityT . tell+  listen = DensityT . listen . getDensityT+  pass = DensityT . pass . getDensityT -instance MonadDistribution m => MonadDistribution (Density m) where+instance (MonadDistribution m) => MonadDistribution (DensityT m) where   random = do     trace <- get     x <- case trace of@@ -36,5 +36,5 @@     tell [x]     pure x -density :: Monad m => Density m b -> [Double] -> m (b, [Double])-density (Density m) = evalStateT (runWriterT m)+runDensityT :: (Monad m) => DensityT m b -> [Double] -> m (b, [Double])+runDensityT (DensityT m) = evalStateT (runWriterT m)

@@ -80,7 +80,7 @@ evidence = Log.sum . map snd . logExplicit  -- | Normalized probability mass of a specific value.-mass :: Ord a => Enumerator a -> a -> Double+mass :: (Ord a) => Enumerator a -> a -> Double mass d = f   where     f a = fromMaybe 0 $ lookup a m@@ -95,7 +95,7 @@ -- The resulting list is sorted ascendingly according to values. -- -- > enumerator = compact . explicit-enumerator, enumerate :: Ord a => Enumerator a -> [(a, Double)]+enumerator, enumerate :: (Ord a) => Enumerator a -> [(a, Double)] enumerator d = filter ((/= 0) . snd) $ compact (zip xs ws)   where     (xs, ws) = second (map (exp . ln) . normalize) $ unzip (logExplicit d)@@ -107,20 +107,20 @@ expectation :: (a -> Double) -> Enumerator a -> Double expectation f = Prelude.sum . map (\(x, w) -> f x * (exp . ln) w) . normalizeWeights . logExplicit -normalize :: Fractional b => [b] -> [b]+normalize :: (Fractional b) => [b] -> [b] normalize xs = map (/ z) xs   where     z = Prelude.sum xs  -- | Divide all weights by their sum.-normalizeWeights :: Fractional b => [(a, b)] -> [(a, b)]+normalizeWeights :: (Fractional b) => [(a, b)] -> [(a, b)] normalizeWeights ls = zip xs ps   where     (xs, ws) = unzip ls     ps = normalize ws  -- | 'compact' followed by removing values with zero weight.-normalForm :: Ord a => Enumerator a -> [(a, Double)]+normalForm :: (Ord a) => Enumerator a -> [(a, Double)] normalForm = filter ((/= 0) . snd) . compact . explicit  toEmpirical :: (Fractional b, Ord a, Ord b) => [a] -> [(a, b)]@@ -139,10 +139,10 @@ removeZeros :: Enumerator a -> Enumerator a removeZeros (Enumerator (WriterT a)) = Enumerator $ WriterT $ filter ((\(Product x) -> x /= 0) . snd) a -instance Ord a => Eq (Enumerator a) where+instance (Ord a) => Eq (Enumerator a) where   p == q = normalForm p == normalForm q -instance Ord a => AEq (Enumerator a) where+instance (Ord a) => AEq (Enumerator a) where   p === q = xs == ys && ps === qs     where       (xs, ps) = unzip (normalForm p)

@@ -7,18 +7,18 @@  import Control.Monad.Bayes.Class (Log (ln)) import Control.Monad.Bayes.Sampler.Lazy-  ( Sampler (runSampler),+  ( SamplerT (runSamplerT),     Tree (..),     Trees (..),     randomTree,   )-import Control.Monad.Bayes.Weighted (Weighted, weighted)+import Control.Monad.Bayes.Weighted (WeightedT, runWeightedT) import Control.Monad.Extra (iterateM) import Control.Monad.State.Lazy (MonadState (get, put), runState) import System.Random (RandomGen (split), getStdGen, newStdGen) import System.Random qualified as R -mh :: forall a. Double -> Weighted Sampler a -> IO [(a, Log Double)]+mh :: forall a. Double -> WeightedT SamplerT a -> IO [(a, Log Double)] mh p m = do   -- Top level: produce a stream of samples.   -- Split the random number generator in two@@ -27,7 +27,7 @@   g <- newStdGen >> getStdGen   let (g1, g2) = split g   let t = randomTree g1-  let (x, w) = runSampler (weighted m) t+  let (x, w) = runSamplerT (runWeightedT m) t   -- Now run step over and over to get a stream of (tree,result,weight)s.   let (samples, _) = runState (iterateM step (t, x, w)) g2   -- The stream of seeds is used to produce a stream of result/weight pairs.@@ -50,7 +50,7 @@       let t' = mutateTree p g1 t       -- Rerun the model with the new tree, to get a new       -- weight w'.-      let (x', w') = runSampler (weighted m) t'+      let (x', w') = runSamplerT (runWeightedT m) t'       -- MH acceptance ratio. This is the probability of either       -- returning the new seed or the old one.       let ratio = w' / w@@ -61,7 +61,7 @@         else return (t, x, w)  -- Replace the labels of a tree randomly, with probability p-mutateTree :: forall g. RandomGen g => Double -> g -> Tree -> Tree+mutateTree :: forall g. (RandomGen g) => Double -> g -> Tree -> Tree mutateTree p g (Tree a ts) =   let (a', g') = (R.random g :: (Double, g))       (a'', g'') = R.random g'@@ -70,7 +70,7 @@           lazyUniforms = mutateTrees p g'' ts         } -mutateTrees :: RandomGen g => Double -> g -> Trees -> Trees+mutateTrees :: (RandomGen g) => Double -> g -> Trees -> Trees mutateTrees p g (Trees t ts) =   let (g1, g2) = split g    in Trees

@@ -1,15 +1,15 @@ module Control.Monad.Bayes.Inference.Lazy.WIS where -import Control.Monad.Bayes.Sampler.Lazy (Sampler, weightedsamples)-import Control.Monad.Bayes.Weighted (Weighted)+import Control.Monad.Bayes.Sampler.Lazy (SamplerT, weightedsamples)+import Control.Monad.Bayes.Weighted (WeightedT) import Numeric.Log (Log (Exp)) import System.Random (Random (randoms), getStdGen, newStdGen) --- | Weighted Importance Sampling+-- | WeightedT Importance Sampling  -- | Likelihood weighted importance sampling first draws n weighted samples, --    and then samples a stream of results from that regarded as an empirical distribution-lwis :: Int -> Weighted Sampler a -> IO [a]+lwis :: Int -> WeightedT SamplerT a -> IO [a] lwis n m = do   xws <- weightedsamples m   let xws' = take n $ accumulate xws 0@@ -18,6 +18,6 @@   rs <- randoms <$> getStdGen   return $ fmap (\r -> fst $ head $ filter ((>= Exp (log r) * max') . snd) xws') rs   where-    accumulate :: Num t => [(a, t)] -> t -> [(a, t)]+    accumulate :: (Num t) => [(a, t)] -> t -> [(a, t)]     accumulate ((x, w) : xws) a = (x, w + a) : (x, w + a) : accumulate xws (w + a)     accumulate [] _ = []

@@ -21,7 +21,7 @@   ) import Control.Monad.Bayes.Traced.Dynamic qualified as Dynamic import Control.Monad.Bayes.Traced.Static qualified as Static-import Control.Monad.Bayes.Weighted (Weighted, unweighted)+import Control.Monad.Bayes.Weighted (WeightedT, unweighted) import Pipes ((>->)) import Pipes qualified as P import Pipes.Prelude qualified as P@@ -33,25 +33,25 @@ defaultMCMCConfig :: MCMCConfig defaultMCMCConfig = MCMCConfig {proposal = SingleSiteMH, numMCMCSteps = 1, numBurnIn = 0} -mcmc :: MonadDistribution m => MCMCConfig -> Static.Traced (Weighted m) a -> m [a]+mcmc :: (MonadDistribution m) => MCMCConfig -> Static.TracedT (WeightedT m) a -> m [a] mcmc (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Static.mh numMCMCSteps m -mcmcBasic :: MonadDistribution m => MCMCConfig -> Basic.Traced (Weighted m) a -> m [a]+mcmcBasic :: (MonadDistribution m) => MCMCConfig -> Basic.TracedT (WeightedT m) a -> m [a] mcmcBasic (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Basic.mh numMCMCSteps m -mcmcDynamic :: MonadDistribution m => MCMCConfig -> Dynamic.Traced (Weighted m) a -> m [a]+mcmcDynamic :: (MonadDistribution m) => MCMCConfig -> Dynamic.TracedT (WeightedT m) a -> m [a] mcmcDynamic (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Dynamic.mh numMCMCSteps m  -- -- | draw iid samples until you get one that has non-zero likelihood-independentSamples :: Monad m => Static.Traced m a -> P.Producer (MHResult a) m (Trace a)-independentSamples (Static.Traced _w d) =+independentSamples :: (Monad m) => Static.TracedT m a -> P.Producer (MHResult a) m (Trace a)+independentSamples (Static.TracedT _w d) =   P.repeatM d     >-> P.takeWhile' ((== 0) . probDensity)     >-> P.map (MHResult False)  -- | convert a probabilistic program into a producer of samples-mcmcP :: MonadDistribution m => MCMCConfig -> Static.Traced m a -> P.Producer (MHResult a) m ()-mcmcP MCMCConfig {..} m@(Static.Traced w _) = do+mcmcP :: (MonadDistribution m) => MCMCConfig -> Static.TracedT m a -> P.Producer (MHResult a) m ()+mcmcP MCMCConfig {..} m@(Static.TracedT w _) = do   initialValue <- independentSamples m >-> P.drain   ( P.unfoldr (fmap (Right . (\k -> (k, trace k))) . mhTransWithBool w) initialValue       >-> P.drop numBurnIn

@@ -22,30 +22,30 @@ import Control.Monad.Bayes.Inference.MCMC (MCMCConfig, mcmc) import Control.Monad.Bayes.Inference.SMC (SMCConfig (), smc) import Control.Monad.Bayes.Population as Pop-  ( Population,+  ( PopulationT,     hoist,-    population,     pushEvidence,+    runPopulationT,   )-import Control.Monad.Bayes.Sequential.Coroutine (Sequential)-import Control.Monad.Bayes.Traced.Static (Traced)+import Control.Monad.Bayes.Sequential.Coroutine (SequentialT)+import Control.Monad.Bayes.Traced.Static (TracedT) import Control.Monad.Bayes.Weighted import Control.Monad.Trans (lift) import Numeric.Log (Log)  -- | Particle Marginal Metropolis-Hastings sampling. pmmh ::-  MonadDistribution m =>+  (MonadDistribution m) =>   MCMCConfig ->-  SMCConfig (Weighted m) ->-  Traced (Weighted m) a1 ->-  (a1 -> Sequential (Population (Weighted m)) a2) ->+  SMCConfig (WeightedT m) ->+  TracedT (WeightedT m) a1 ->+  (a1 -> SequentialT (PopulationT (WeightedT m)) a2) ->   m [[(a2, Log Double)]] pmmh mcmcConf smcConf param model =   mcmc     mcmcConf     ( param-        >>= population+        >>= runPopulationT           . pushEvidence           . Pop.hoist lift           . smc smcConf@@ -54,9 +54,9 @@  -- | Particle Marginal Metropolis-Hastings sampling from a Bayesian model pmmhBayesianModel ::-  MonadMeasure m =>+  (MonadMeasure m) =>   MCMCConfig ->-  SMCConfig (Weighted m) ->-  (forall m'. MonadMeasure m' => Bayesian m' a1 a2) ->+  SMCConfig (WeightedT m) ->+  (forall m'. (MonadMeasure m') => Bayesian m' a1 a2) ->   m [[(a2, Log Double)]] pmmhBayesianModel mcmcConf smcConf bm = pmmh mcmcConf smcConf (prior bm) (generative bm)

@@ -24,7 +24,7 @@ import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (..)) import Control.Monad.Bayes.Inference.SMC import Control.Monad.Bayes.Population-  ( Population,+  ( PopulationT,     spawn,     withParticles,   )@@ -33,7 +33,7 @@ import Control.Monad.Bayes.Traced.Basic qualified as TrBas import Control.Monad.Bayes.Traced.Dynamic qualified as TrDyn import Control.Monad.Bayes.Traced.Static as Tr-  ( Traced,+  ( TracedT,     marginal,     mhStep,   )@@ -42,12 +42,12 @@  -- | Resample-move Sequential Monte Carlo. rmsmc ::-  MonadDistribution m =>+  (MonadDistribution m) =>   MCMCConfig ->   SMCConfig m ->   -- | model-  Sequential (Traced (Population m)) a ->-  Population m a+  SequentialT (TracedT (PopulationT m)) a ->+  PopulationT m a rmsmc (MCMCConfig {..}) (SMCConfig {..}) =   marginal     . S.sequentially (composeCopies numMCMCSteps mhStep . TrStat.hoist resampler) numSteps@@ -56,12 +56,12 @@ -- | Resample-move Sequential Monte Carlo with a more efficient -- tracing representation. rmsmcBasic ::-  MonadDistribution m =>+  (MonadDistribution m) =>   MCMCConfig ->   SMCConfig m ->   -- | model-  Sequential (TrBas.Traced (Population m)) a ->-  Population m a+  SequentialT (TrBas.TracedT (PopulationT m)) a ->+  PopulationT m a rmsmcBasic (MCMCConfig {..}) (SMCConfig {..}) =   TrBas.marginal     . S.sequentially (composeCopies numMCMCSteps TrBas.mhStep . TrBas.hoist resampler) numSteps@@ -71,12 +71,12 @@ -- where only random variables since last resampling are considered -- for rejuvenation. rmsmcDynamic ::-  MonadDistribution m =>+  (MonadDistribution m) =>   MCMCConfig ->   SMCConfig m ->   -- | model-  Sequential (TrDyn.Traced (Population m)) a ->-  Population m a+  SequentialT (TrDyn.TracedT (PopulationT m)) a ->+  PopulationT m a rmsmcDynamic (MCMCConfig {..}) (SMCConfig {..}) =   TrDyn.marginal     . S.sequentially (TrDyn.freeze . composeCopies numMCMCSteps TrDyn.mhStep . TrDyn.hoist resampler) numSteps

@@ -22,14 +22,14 @@  import Control.Monad.Bayes.Class (MonadDistribution, MonadMeasure) import Control.Monad.Bayes.Population-  ( Population,+  ( PopulationT,     pushEvidence,     withParticles,   ) import Control.Monad.Bayes.Sequential.Coroutine as Coroutine  data SMCConfig m = SMCConfig-  { resampler :: forall x. Population m x -> Population m x,+  { resampler :: forall x. PopulationT m x -> PopulationT m x,     numSteps :: Int,     numParticles :: Int   }@@ -37,10 +37,10 @@ -- | Sequential importance resampling. -- Basically an SMC template that takes a custom resampler. smc ::-  MonadDistribution m =>+  (MonadDistribution m) =>   SMCConfig m ->-  Coroutine.Sequential (Population m) a ->-  Population m a+  Coroutine.SequentialT (PopulationT m) a ->+  PopulationT m a smc SMCConfig {..} =   Coroutine.sequentially resampler numSteps     . Coroutine.hoistFirst (withParticles numParticles)@@ -49,5 +49,5 @@ -- Weights are normalized at each timestep and the total weight is pushed -- as a score into the transformed monad. smcPush ::-  MonadMeasure m => SMCConfig m -> Coroutine.Sequential (Population m) a -> Population m a+  (MonadMeasure m) => SMCConfig m -> Coroutine.SequentialT (PopulationT m) a -> PopulationT m a smcPush config = smc config {resampler = (pushEvidence . resampler config)}

@@ -27,33 +27,33 @@ import Control.Monad.Bayes.Inference.MCMC import Control.Monad.Bayes.Inference.RMSMC (rmsmc) import Control.Monad.Bayes.Inference.SMC (SMCConfig (SMCConfig, numParticles, numSteps, resampler), smcPush)-import Control.Monad.Bayes.Population as Pop (Population, population, resampleMultinomial)-import Control.Monad.Bayes.Sequential.Coroutine (Sequential)+import Control.Monad.Bayes.Population as Pop (PopulationT, resampleMultinomial, runPopulationT)+import Control.Monad.Bayes.Sequential.Coroutine (SequentialT) import Control.Monad.Bayes.Traced import Control.Monad.Trans (MonadTrans (..)) import Numeric.Log (Log)  -- | Helper monad transformer for preprocessing the model for 'smc2'.-newtype SMC2 m a = SMC2 (Sequential (Traced (Population m)) a)+newtype SMC2 m a = SMC2 (SequentialT (TracedT (PopulationT m)) a)   deriving newtype (Functor, Applicative, Monad) -setup :: SMC2 m a -> Sequential (Traced (Population m)) a+setup :: SMC2 m a -> SequentialT (TracedT (PopulationT m)) a setup (SMC2 m) = m  instance MonadTrans SMC2 where   lift = SMC2 . lift . lift . lift -instance MonadDistribution m => MonadDistribution (SMC2 m) where+instance (MonadDistribution m) => MonadDistribution (SMC2 m) where   random = lift random -instance Monad m => MonadFactor (SMC2 m) where+instance (Monad m) => MonadFactor (SMC2 m) where   score = SMC2 . score -instance MonadDistribution m => MonadMeasure (SMC2 m)+instance (MonadDistribution m) => MonadMeasure (SMC2 m)  -- | Sequential Monte Carlo squared. smc2 ::-  MonadDistribution m =>+  (MonadDistribution m) =>   -- | number of time steps   Int ->   -- | number of inner particles@@ -63,12 +63,12 @@   -- | number of MH transitions   Int ->   -- | model parameters-  Sequential (Traced (Population m)) b ->+  SequentialT (TracedT (PopulationT m)) b ->   -- | model-  (b -> Sequential (Population (SMC2 m)) a) ->-  Population m [(a, Log Double)]+  (b -> SequentialT (PopulationT (SMC2 m)) a) ->+  PopulationT m [(a, Log Double)] smc2 k n p t param m =   rmsmc     MCMCConfig {numMCMCSteps = t, proposal = SingleSiteMH, numBurnIn = 0}     SMCConfig {numParticles = p, numSteps = k, resampler = resampleMultinomial}-    (param >>= setup . population . smcPush (SMCConfig {numSteps = k, numParticles = n, resampler = resampleMultinomial}) . m)+    (param >>= setup . runPopulationT . smcPush (SMCConfig {numSteps = k, numParticles = n, resampler = resampleMultinomial}) . m)

@@ -18,7 +18,7 @@ import Control.Monad.Bayes.Enumerator (toEmpirical) import Control.Monad.Bayes.Inference.MCMC import Control.Monad.Bayes.Sampler.Strict (SamplerIO, sampleIO)-import Control.Monad.Bayes.Traced (Traced)+import Control.Monad.Bayes.Traced (TracedT) import Control.Monad.Bayes.Traced.Common hiding (burnIn) import Control.Monad.Bayes.Weighted import Data.Scientific (FPFormat (Exponent), formatScientific, fromFloatDigits)@@ -106,7 +106,7 @@     . (fmap (second (formatScientific Exponent (Just 3) . fromFloatDigits)))     . toEmpirical -showVal :: Show a => [a] -> Widget n+showVal :: (Show a) => [a] -> Widget n showVal = txt . T.pack . (\case [] -> ""; a -> show $ head a)  -- | handler for events received by the TUI@@ -130,7 +130,7 @@       (attrName "highlight", fg yellow)     ] -tui :: Show a => Int -> Traced (Weighted SamplerIO) a -> ([a] -> Widget ()) -> IO ()+tui :: (Show a) => Int -> TracedT (WeightedT SamplerIO) a -> ([a] -> Widget ()) -> IO () tui burnIn distribution visualizer = void do   eventChan <- B.newBChan 10   initialVty <- buildVty

@@ -35,7 +35,7 @@ import Control.Foldl (Fold) import Control.Foldl qualified as Foldl import Control.Monad.Bayes.Class (MonadDistribution (bernoulli, random, uniformD))-import Control.Monad.Bayes.Weighted (Weighted, weighted)+import Control.Monad.Bayes.Weighted (WeightedT, runWeightedT) import Control.Monad.Cont   ( Cont,     ContT (ContT),@@ -49,13 +49,15 @@ import Statistics.Distribution qualified as Statistics import Statistics.Distribution.Uniform qualified as Statistics -newtype Integrator a = Integrator {getCont :: Cont Double a}+newtype Integrator a = Integrator {getIntegrator :: Cont Double a}   deriving newtype (Functor, Applicative, Monad) -integrator, runIntegrator :: (a -> Double) -> Integrator a -> Double-integrator f (Integrator a) = runCont a f-runIntegrator = integrator+runIntegrator :: (a -> Double) -> Integrator a -> Double+runIntegrator f (Integrator a) = runCont a f +integrator :: ((a -> Double) -> Double) -> Integrator a+integrator = Integrator . cont+ instance MonadDistribution Integrator where   random = fromDensityFunction $ Statistics.density $ Statistics.uniformDistr 0 1   bernoulli p = Integrator $ cont (\f -> p * f True + (1 - p) * f False)@@ -68,45 +70,45 @@   where     integralWithQuadrature = result . last . (\z -> trap z 0 1) -fromMassFunction :: Foldable f => (a -> Double) -> f a -> Integrator a+fromMassFunction :: (Foldable f) => (a -> Double) -> f a -> Integrator a fromMassFunction f support = Integrator $ cont \g ->   foldl' (\acc x -> acc + f x * g x) 0 support -empirical :: Foldable f => f a -> Integrator a+empirical :: (Foldable f) => f a -> Integrator a empirical = Integrator . cont . flip weightedAverage   where     weightedAverage :: (Foldable f, Fractional r) => (a -> r) -> f a -> r     weightedAverage f = Foldl.fold (weightedAverageFold f) -    weightedAverageFold :: Fractional r => (a -> r) -> Fold a r+    weightedAverageFold :: (Fractional r) => (a -> r) -> Fold a r     weightedAverageFold f = Foldl.premap f averageFold -    averageFold :: Fractional a => Fold a a+    averageFold :: (Fractional a) => Fold a a     averageFold = (/) <$> Foldl.sum <*> Foldl.genericLength  expectation :: Integrator Double -> Double-expectation = integrator id+expectation = runIntegrator id  variance :: Integrator Double -> Double-variance nu = integrator (^ 2) nu - expectation nu ^ 2+variance nu = runIntegrator (^ 2) nu - expectation nu ^ 2  momentGeneratingFunction :: Integrator Double -> Double -> Double-momentGeneratingFunction nu t = integrator (\x -> exp (t * x)) nu+momentGeneratingFunction nu t = runIntegrator (\x -> exp (t * x)) nu  cumulantGeneratingFunction :: Integrator Double -> Double -> Double cumulantGeneratingFunction nu = log . momentGeneratingFunction nu -normalize :: Weighted Integrator a -> Integrator a+normalize :: WeightedT Integrator a -> Integrator a normalize m =-  let m' = weighted m-      z = integrator (ln . exp . snd) m'+  let m' = runWeightedT m+      z = runIntegrator (ln . exp . snd) m'    in do-        (x, d) <- weighted m+        (x, d) <- runWeightedT m         Integrator $ cont $ \f -> (f () * (ln $ exp d)) / z         return x  cdf :: Integrator Double -> Double -> Double-cdf nu x = integrator (negativeInfinity `to` x) nu+cdf nu x = runIntegrator (negativeInfinity `to` x) nu   where     negativeInfinity :: Double     negativeInfinity = negate (1 / 0)@@ -117,14 +119,14 @@       | otherwise = 0  volume :: Integrator Double -> Double-volume = integrator (const 1)+volume = runIntegrator (const 1)  containing :: (Num a, Eq b) => [b] -> b -> a containing xs x   | x `elem` xs = 1   | otherwise = 0 -instance Num a => Num (Integrator a) where+instance (Num a) => Num (Integrator a) where   (+) = liftA2 (+)   (-) = liftA2 (-)   (*) = liftA2 (*)@@ -132,13 +134,13 @@   signum = fmap signum   fromInteger = pure . fromInteger -probability :: Ord a => (a, a) -> Integrator a -> Double-probability (lower, upper) = integrator (\x -> if x < upper && x >= lower then 1 else 0)+probability :: (Ord a) => (a, a) -> Integrator a -> Double+probability (lower, upper) = runIntegrator (\x -> if x < upper && x >= lower then 1 else 0) -enumeratorWith :: Ord a => Set a -> Integrator a -> [(a, Double)]+enumeratorWith :: (Ord a) => Set a -> Integrator a -> [(a, Double)] enumeratorWith ls meas =   [ ( val,-      integrator+      runIntegrator         (\x -> if x == val then 1 else 0)         meas     )@@ -149,7 +151,7 @@   (Enum a, Ord a, Fractional a) =>   Int ->   a ->-  Weighted Integrator a ->+  WeightedT Integrator a ->   [(a, Double)] histogram nBins binSize model = do   x <- take nBins [1 ..]

@@ -13,11 +13,10 @@ -- Stability   : experimental -- Portability : GHC ----- 'Population' turns a single sample into a collection of weighted samples.+-- 'PopulationT' turns a single sample into a collection of weighted samples. module Control.Monad.Bayes.Population-  ( Population,-    population,-    runPopulation,+  ( PopulationT (..),+    runPopulationT,     explicitPopulation,     fromWeightedList,     spawn,@@ -47,11 +46,11 @@     factor,   ) import Control.Monad.Bayes.Weighted-  ( Weighted,+  ( WeightedT,     applyWeight,     extractWeight,-    weighted,-    withWeight,+    runWeightedT,+    weightedT,   ) import Control.Monad.List (ListT (..), MonadIO, MonadTrans (..)) import Data.List (unfoldr)@@ -64,46 +63,43 @@ import Prelude hiding (all, sum)  -- | A collection of weighted samples, or particles.-newtype Population m a = Population (Weighted (ListT m) a)+newtype PopulationT m a = PopulationT {getPopulationT :: WeightedT (ListT m) a}   deriving newtype (Functor, Applicative, Monad, MonadIO, MonadDistribution, MonadFactor, MonadMeasure) -instance MonadTrans Population where-  lift = Population . lift . lift+instance MonadTrans PopulationT where+  lift = PopulationT . lift . lift  -- | Explicit representation of the weighted sample with weights in the log -- domain.-population, runPopulation :: Population m a -> m [(a, Log Double)]-population (Population m) = runListT $ weighted m---- | deprecated synonym-runPopulation = population+runPopulationT :: PopulationT m a -> m [(a, Log Double)]+runPopulationT = runListT . runWeightedT . getPopulationT  -- | Explicit representation of the weighted sample.-explicitPopulation :: Functor m => Population m a -> m [(a, Double)]-explicitPopulation = fmap (map (second (exp . ln))) . population+explicitPopulation :: (Functor m) => PopulationT m a -> m [(a, Double)]+explicitPopulation = fmap (map (second (exp . ln))) . runPopulationT --- | Initialize 'Population' with a concrete weighted sample.-fromWeightedList :: Monad m => m [(a, Log Double)] -> Population m a-fromWeightedList = Population . withWeight . ListT+-- | Initialize 'PopulationT' with a concrete weighted sample.+fromWeightedList :: (Monad m) => m [(a, Log Double)] -> PopulationT m a+fromWeightedList = PopulationT . weightedT . ListT  -- | Increase the sample size by a given factor. -- The weights are adjusted such that their sum is preserved. -- It is therefore safe to use 'spawn' in arbitrary places in the program -- without introducing bias.-spawn :: Monad m => Int -> Population m ()+spawn :: (Monad m) => Int -> PopulationT m () spawn n = fromWeightedList $ pure $ replicate n ((), 1 / fromIntegral n) -withParticles :: Monad m => Int -> Population m a -> Population m a+withParticles :: (Monad m) => Int -> PopulationT m a -> PopulationT m a withParticles n = (spawn n >>)  resampleGeneric ::-  MonadDistribution m =>+  (MonadDistribution m) =>   -- | resampler   (V.Vector Double -> m [Int]) ->-  Population m a ->-  Population m a+  PopulationT m a ->+  PopulationT m a resampleGeneric resampler m = fromWeightedList $ do-  pop <- population m+  pop <- runPopulationT m   let (xs, ps) = unzip pop   let n = length xs   let z = Log.sum ps@@ -150,8 +146,8 @@ -- The total weight is preserved. resampleSystematic ::   (MonadDistribution m) =>-  Population m a ->-  Population m a+  PopulationT m a ->+  PopulationT m a resampleSystematic = resampleGeneric (\ps -> (`systematic` ps) <$> random)  -- | Stratified sampler.@@ -171,7 +167,7 @@ -- and \(w^{(k)}\) are the weights. -- -- The conditional variance of stratified sampling is always smaller than that of multinomial sampling and it is also unbiased - see  [Comparison of Resampling Schemes for Particle Filtering](https://arxiv.org/abs/cs/0507025).-stratified :: MonadDistribution m => V.Vector Double -> m [Int]+stratified :: (MonadDistribution m) => V.Vector Double -> m [Int] stratified weights = do   let bigN = V.length weights   dithers <- V.replicateM bigN (uniform 0.0 1.0)@@ -192,32 +188,32 @@ -- The total weight is preserved. resampleStratified ::   (MonadDistribution m) =>-  Population m a ->-  Population m a+  PopulationT m a ->+  PopulationT m a resampleStratified = resampleGeneric stratified  -- | Multinomial sampler.  Sample from \(0, \ldots, n - 1\) \(n\) -- times drawn at random according to the weights where \(n\) is the -- length of vector of weights.-multinomial :: MonadDistribution m => V.Vector Double -> m [Int]+multinomial :: (MonadDistribution m) => V.Vector Double -> m [Int] multinomial ps = replicateM (V.length ps) (categorical ps)  -- | Resample the population using the underlying monad and a multinomial resampling scheme. -- The total weight is preserved. resampleMultinomial ::   (MonadDistribution m) =>-  Population m a ->-  Population m a+  PopulationT m a ->+  PopulationT m a resampleMultinomial = resampleGeneric multinomial --- | Separate the sum of weights into the 'Weighted' transformer.+-- | Separate the sum of weights into the 'WeightedT' transformer. -- Weights are normalized after this operation. extractEvidence ::-  Monad m =>-  Population m a ->-  Population (Weighted m) a+  (Monad m) =>+  PopulationT m a ->+  PopulationT (WeightedT m) a extractEvidence m = fromWeightedList $ do-  pop <- lift $ population m+  pop <- lift $ runPopulationT m   let (xs, ps) = unzip pop   let z = sum ps   let ws = map (if z > 0 then (/ z) else const (1 / fromIntegral (length ps))) ps@@ -227,27 +223,27 @@ -- | Push the evidence estimator as a score to the transformed monad. -- Weights are normalized after this operation. pushEvidence ::-  MonadFactor m =>-  Population m a ->-  Population m a+  (MonadFactor m) =>+  PopulationT m a ->+  PopulationT m a pushEvidence = hoist applyWeight . extractEvidence  -- | A properly weighted single sample, that is one picked at random according -- to the weights, with the sum of all weights. proper ::   (MonadDistribution m) =>-  Population m a ->-  Weighted m a+  PopulationT m a ->+  WeightedT m a proper m = do-  pop <- population $ extractEvidence m+  pop <- runPopulationT $ extractEvidence m   let (xs, ps) = unzip pop   index <- logCategorical $ V.fromList ps   let x = xs !! index   return x  -- | Model evidence estimator, also known as pseudo-marginal likelihood.-evidence :: (Monad m) => Population m a -> m (Log Double)-evidence = extractWeight . population . extractEvidence+evidence :: (Monad m) => PopulationT m a -> m (Log Double)+evidence = extractWeight . runPopulationT . extractEvidence  -- | Picks one point from the population and uses model evidence as a 'score' -- in the transformed monad.@@ -255,12 +251,12 @@ -- introducing bias. collapse ::   (MonadMeasure m) =>-  Population m a ->+  PopulationT m a ->   m a collapse = applyWeight . proper --- | Population average of a function, computed using unnormalized weights.-popAvg :: (Monad m) => (a -> Double) -> Population m a -> m Double+-- | PopulationT average of a function, computed using unnormalized weights.+popAvg :: (Monad m) => (a -> Double) -> PopulationT m a -> m Double popAvg f p = do   xs <- explicitPopulation p   let ys = map (\(x, w) -> f x * w) xs@@ -269,8 +265,8 @@  -- | Applies a transformation to the inner monad. hoist ::-  Monad n =>+  (Monad n) =>   (forall x. m x -> n x) ->-  Population m a ->-  Population n a-hoist f = fromWeightedList . f . population+  PopulationT m a ->+  PopulationT n a+hoist f = fromWeightedList . f . runPopulationT

@@ -8,7 +8,7 @@  import Control.Monad (ap) import Control.Monad.Bayes.Class (MonadDistribution (random))-import Control.Monad.Bayes.Weighted (Weighted, weighted)+import Control.Monad.Bayes.Weighted (WeightedT, runWeightedT) import Numeric.Log (Log (..)) import System.Random   ( RandomGen (split),@@ -35,7 +35,7 @@ -- | A probability distribution over a is -- | a function 'Tree -> a' -- | The idea is that it uses up bits of the tree as it runs-newtype Sampler a = Sampler {runSampler :: Tree -> a}+newtype SamplerT a = SamplerT {runSamplerT :: Tree -> a}   deriving (Functor)  -- | Two key things to do with trees:@@ -45,35 +45,35 @@ splitTree (Tree r (Trees t ts)) = (t, Tree r ts)  -- | Preliminaries for the simulation methods. Generate a tree with uniform random labels. This uses 'split' to split a random seed-randomTree :: RandomGen g => g -> Tree+randomTree :: (RandomGen g) => g -> Tree randomTree g = let (a, g') = R.random g in Tree a (randomTrees g') -randomTrees :: RandomGen g => g -> Trees+randomTrees :: (RandomGen g) => g -> Trees randomTrees g = let (g1, g2) = split g in Trees (randomTree g1) (randomTrees g2) -instance Applicative Sampler where-  pure = Sampler . const+instance Applicative SamplerT where+  pure = SamplerT . const   (<*>) = ap  -- | probabilities for a monad. -- | Sequencing is done by splitting the tree -- | and using different bits for different computations.-instance Monad Sampler where+instance Monad SamplerT where   return = pure-  (Sampler m) >>= f = Sampler \g ->+  (SamplerT m) >>= f = SamplerT \g ->     let (g1, g2) = splitTree g-        (Sampler m') = f (m g1)+        (SamplerT m') = f (m g1)      in m' g2 -instance MonadDistribution Sampler where-  random = Sampler \(Tree r _) -> r+instance MonadDistribution SamplerT where+  random = SamplerT \(Tree r _) -> r -sampler :: Sampler a -> IO a-sampler m = newStdGen *> (runSampler m . randomTree <$> getStdGen)+sampler :: SamplerT a -> IO a+sampler m = newStdGen *> (runSamplerT m . randomTree <$> getStdGen) -independent :: Monad m => m a -> m [a]+independent :: (Monad m) => m a -> m [a] independent = sequence . repeat  -- | 'weightedsamples' runs a probability measure and gets out a stream of (result,weight) pairs-weightedsamples :: Weighted Sampler a -> IO [(a, Log Double)]-weightedsamples = sampler . independent . weighted+weightedsamples :: WeightedT SamplerT a -> IO [(a, Log Double)]+weightedsamples = sampler . independent . runWeightedT

@@ -15,7 +15,7 @@ -- 'SamplerIO' and 'SamplerST' are instances of 'MonadDistribution'. Apply a 'MonadFactor' -- transformer to obtain a 'MonadMeasure' that can execute probabilistic models. module Control.Monad.Bayes.Sampler.Strict-  ( Sampler,+  ( SamplerT (..),     SamplerIO,     SamplerST,     sampleIO,@@ -48,27 +48,27 @@  -- | The sampling interpretation of a probabilistic program -- Here m is typically IO or ST-newtype Sampler g m a = Sampler (ReaderT g m a) deriving (Functor, Applicative, Monad, MonadIO)+newtype SamplerT g m a = SamplerT {runSamplerT :: ReaderT g m a} deriving (Functor, Applicative, Monad, MonadIO) --- | convenient type synonym to show specializations of Sampler+-- | convenient type synonym to show specializations of SamplerT -- to particular pairs of monad and RNG-type SamplerIO = Sampler (IOGenM StdGen) IO+type SamplerIO = SamplerT (IOGenM StdGen) IO --- | convenient type synonym to show specializations of Sampler+-- | convenient type synonym to show specializations of SamplerT -- to particular pairs of monad and RNG-type SamplerST s = Sampler (STGenM StdGen s) (ST s)+type SamplerST s = SamplerT (STGenM StdGen s) (ST s) -instance StatefulGen g m => MonadDistribution (Sampler g m) where-  random = Sampler (ReaderT uniformDouble01M)+instance (StatefulGen g m) => MonadDistribution (SamplerT g m) where+  random = SamplerT (ReaderT uniformDouble01M) -  uniform a b = Sampler (ReaderT $ uniformRM (a, b))-  normal m s = Sampler (ReaderT (MWC.normal m s))-  gamma shape scale = Sampler (ReaderT $ MWC.gamma shape scale)-  beta a b = Sampler (ReaderT $ MWC.beta a b)+  uniform a b = SamplerT (ReaderT $ uniformRM (a, b))+  normal m s = SamplerT (ReaderT (MWC.normal m s))+  gamma shape scale = SamplerT (ReaderT $ MWC.gamma shape scale)+  beta a b = SamplerT (ReaderT $ MWC.beta a b) -  bernoulli p = Sampler (ReaderT $ MWC.bernoulli p)-  categorical ps = Sampler (ReaderT $ MWC.categorical ps)-  geometric p = Sampler (ReaderT $ MWC.geometric0 p)+  bernoulli p = SamplerT (ReaderT $ MWC.bernoulli p)+  categorical ps = SamplerT (ReaderT $ MWC.categorical ps)+  geometric p = SamplerT (ReaderT $ MWC.geometric0 p)  -- | Sample with a random number generator of your choice e.g. the one -- from `System.Random`.@@ -77,8 +77,8 @@ -- >>> import System.Random.Stateful hiding (random) -- >>> newIOGenM (mkStdGen 1729) >>= sampleWith random -- 4.690861245089605e-2-sampleWith :: Sampler g m a -> g -> m a-sampleWith (Sampler m) = runReaderT m+sampleWith :: SamplerT g m a -> g -> m a+sampleWith (SamplerT m) = runReaderT m  -- | initialize random seed using system entropy, and sample sampleIO, sampler :: SamplerIO a -> IO a

@@ -11,9 +11,9 @@ -- Stability   : experimental -- Portability : GHC ----- 'Sequential' represents a computation that can be suspended.+-- 'SequentialT' represents a computation that can be suspended. module Control.Monad.Bayes.Sequential.Coroutine-  ( Sequential,+  ( SequentialT,     suspend,     finish,     advance,@@ -48,55 +48,55 @@ -- useful for implementation of Sequential Monte Carlo related methods. -- All the probabilistic effects are lifted from the transformed monad, but -- also `suspend` is inserted after each `factor`.-newtype Sequential m a = Sequential {runSequential :: Coroutine (Await ()) m a}+newtype SequentialT m a = SequentialT {runSequentialT :: Coroutine (Await ()) m a}   deriving newtype (Functor, Applicative, Monad, MonadTrans, MonadIO)  extract :: Await () a -> a extract (Await f) = f () -instance MonadDistribution m => MonadDistribution (Sequential m) where+instance (MonadDistribution m) => MonadDistribution (SequentialT m) where   random = lift random   bernoulli = lift . bernoulli   categorical = lift . categorical  -- | Execution is 'suspend'ed after each 'score'.-instance MonadFactor m => MonadFactor (Sequential m) where+instance (MonadFactor m) => MonadFactor (SequentialT m) where   score w = lift (score w) >> suspend -instance MonadMeasure m => MonadMeasure (Sequential m)+instance (MonadMeasure m) => MonadMeasure (SequentialT m)  -- | A point where the computation is paused.-suspend :: Monad m => Sequential m ()-suspend = Sequential await+suspend :: (Monad m) => SequentialT m ()+suspend = SequentialT await  -- | Remove the remaining suspension points.-finish :: Monad m => Sequential m a -> m a-finish = pogoStick extract . runSequential+finish :: (Monad m) => SequentialT m a -> m a+finish = pogoStick extract . runSequentialT  -- | Execute to the next suspension point. -- If the computation is finished, do nothing. -- -- > finish = finish . advance-advance :: Monad m => Sequential m a -> Sequential m a-advance = Sequential . bounce extract . runSequential+advance :: (Monad m) => SequentialT m a -> SequentialT m a+advance = SequentialT . bounce extract . runSequentialT  -- | Return True if no more suspension points remain.-finished :: Monad m => Sequential m a -> m Bool-finished = fmap isRight . resume . runSequential+finished :: (Monad m) => SequentialT m a -> m Bool+finished = fmap isRight . resume . runSequentialT  -- | Transform the inner monad. -- This operation only applies to computation up to the first suspension.-hoistFirst :: (forall x. m x -> m x) -> Sequential m a -> Sequential m a-hoistFirst f = Sequential . Coroutine . f . resume . runSequential+hoistFirst :: (forall x. m x -> m x) -> SequentialT m a -> SequentialT m a+hoistFirst f = SequentialT . Coroutine . f . resume . runSequentialT  -- | Transform the inner monad. -- The transformation is applied recursively through all the suspension points. hoist ::   (Monad m, Monad n) =>   (forall x. m x -> n x) ->-  Sequential m a ->-  Sequential n a-hoist f = Sequential . mapMonad f . runSequential+  SequentialT m a ->+  SequentialT n a+hoist f = SequentialT . mapMonad f . runSequentialT  -- | Apply a function a given number of times. composeCopies :: Int -> (a -> a) -> (a -> a)@@ -106,12 +106,12 @@ -- Applies a given transformation after each time step. sequentially,   sis ::-    Monad m =>+    (Monad m) =>     -- | transformation     (forall x. m x -> m x) ->     -- | number of time steps     Int ->-    Sequential m a ->+    SequentialT m a ->     m a sequentially f k = finish . composeCopies k (advance . hoistFirst f) 

@@ -9,7 +9,7 @@ -- Stability   : experimental -- Portability : GHC module Control.Monad.Bayes.Traced.Basic-  ( Traced,+  ( TracedT,     hoist,     marginal,     mhStep,@@ -23,7 +23,7 @@     MonadFactor (..),     MonadMeasure,   )-import Control.Monad.Bayes.Density.Free (Density)+import Control.Monad.Bayes.Density.Free (DensityT) import Control.Monad.Bayes.Traced.Common   ( Trace (..),     bind,@@ -31,56 +31,56 @@     scored,     singleton,   )-import Control.Monad.Bayes.Weighted (Weighted)+import Control.Monad.Bayes.Weighted (WeightedT) import Data.Functor.Identity (Identity) import Data.List.NonEmpty as NE (NonEmpty ((:|)), toList)  -- | Tracing monad that records random choices made in the program.-data Traced m a = Traced+data TracedT m a = TracedT   { -- | Run the program with a modified trace.-    model :: Weighted (Density Identity) a,+    model :: WeightedT (DensityT Identity) a,     -- | Record trace and output.     traceDist :: m (Trace a)   } -instance Monad m => Functor (Traced m) where-  fmap f (Traced m d) = Traced (fmap f m) (fmap (fmap f) d)+instance (Monad m) => Functor (TracedT m) where+  fmap f (TracedT m d) = TracedT (fmap f m) (fmap (fmap f) d) -instance Monad m => Applicative (Traced m) where-  pure x = Traced (pure x) (pure (pure x))-  (Traced mf df) <*> (Traced mx dx) = Traced (mf <*> mx) (liftA2 (<*>) df dx)+instance (Monad m) => Applicative (TracedT m) where+  pure x = TracedT (pure x) (pure (pure x))+  (TracedT mf df) <*> (TracedT mx dx) = TracedT (mf <*> mx) (liftA2 (<*>) df dx) -instance Monad m => Monad (Traced m) where-  (Traced mx dx) >>= f = Traced my dy+instance (Monad m) => Monad (TracedT m) where+  (TracedT mx dx) >>= f = TracedT my dy     where       my = mx >>= model . f       dy = dx `bind` (traceDist . f) -instance MonadDistribution m => MonadDistribution (Traced m) where-  random = Traced random (fmap singleton random)+instance (MonadDistribution m) => MonadDistribution (TracedT m) where+  random = TracedT random (fmap singleton random) -instance MonadFactor m => MonadFactor (Traced m) where-  score w = Traced (score w) (score w >> pure (scored w))+instance (MonadFactor m) => MonadFactor (TracedT m) where+  score w = TracedT (score w) (score w >> pure (scored w)) -instance MonadMeasure m => MonadMeasure (Traced m)+instance (MonadMeasure m) => MonadMeasure (TracedT m) -hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a-hoist f (Traced m d) = Traced m (f d)+hoist :: (forall x. m x -> m x) -> TracedT m a -> TracedT m a+hoist f (TracedT m d) = TracedT m (f d)  -- | Discard the trace and supporting infrastructure.-marginal :: Monad m => Traced m a -> m a-marginal (Traced _ d) = fmap output d+marginal :: (Monad m) => TracedT m a -> m a+marginal (TracedT _ d) = fmap output d  -- | A single step of the Trace Metropolis-Hastings algorithm.-mhStep :: MonadDistribution m => Traced m a -> Traced m a-mhStep (Traced m d) = Traced m d'+mhStep :: (MonadDistribution m) => TracedT m a -> TracedT m a+mhStep (TracedT m d) = TracedT m d'   where     d' = d >>= mhTrans' m  -- | Full run of the Trace Metropolis-Hastings algorithm with a specified -- number of steps.-mh :: MonadDistribution m => Int -> Traced m a -> m [a]-mh n (Traced m d) = fmap (map output . NE.toList) (f n)+mh :: (MonadDistribution m) => Int -> TracedT m a -> m [a]+mh n (TracedT m d) = fmap (map output . NE.toList) (f n)   where     f k       | k <= 0 = fmap (:| []) d

@@ -26,10 +26,10 @@   ) import Control.Monad.Bayes.Density.Free qualified as Free import Control.Monad.Bayes.Density.State qualified as State-import Control.Monad.Bayes.Weighted as Weighted-  ( Weighted,+import Control.Monad.Bayes.Weighted as WeightedT+  ( WeightedT,     hoist,-    weighted,+    runWeightedT,   ) import Control.Monad.Writer (WriterT (WriterT, runWriterT)) import Data.Functor.Identity (Identity (runIdentity))@@ -74,14 +74,14 @@ scored :: Log Double -> Trace () scored w = Trace {variables = [], output = (), probDensity = w} -bind :: Monad m => m (Trace a) -> (a -> m (Trace b)) -> m (Trace b)+bind :: (Monad m) => m (Trace a) -> (a -> m (Trace b)) -> m (Trace b) bind dx f = do   t1 <- dx   t2 <- f (output t1)   return $ t2 {variables = variables t1 ++ variables t2, probDensity = probDensity t1 * probDensity t2}  -- | A single Metropolis-corrected transition of single-site Trace MCMC.-mhTrans :: MonadDistribution m => (Weighted (State.Density m)) a -> Trace a -> m (Trace a)+mhTrans :: (MonadDistribution m) => (WeightedT (State.DensityT m)) a -> Trace a -> m (Trace a) mhTrans m t@Trace {variables = us, probDensity = p} = do   let n = length us   us' <- do@@ -90,16 +90,16 @@     case splitAt i us of       (xs, _ : ys) -> return $ xs ++ (u' : ys)       _ -> error "impossible"-  ((b, q), vs) <- State.density (weighted m) us'+  ((b, q), vs) <- State.runDensityT (runWeightedT m) us'   let ratio = (exp . ln) $ min 1 (q * fromIntegral n / (p * fromIntegral (length vs)))   accept <- bernoulli ratio   return $ if accept then Trace vs b q else t -mhTransFree :: MonadDistribution m => Weighted (Free.Density m) a -> Trace a -> m (Trace a)+mhTransFree :: (MonadDistribution m) => WeightedT (Free.DensityT m) a -> Trace a -> m (Trace a) mhTransFree m t = trace <$> mhTransWithBool m t  -- | A single Metropolis-corrected transition of single-site Trace MCMC.-mhTransWithBool :: MonadDistribution m => Weighted (Free.Density m) a -> Trace a -> m (MHResult a)+mhTransWithBool :: (MonadDistribution m) => WeightedT (Free.DensityT m) a -> Trace a -> m (MHResult a) mhTransWithBool m t@Trace {variables = us, probDensity = p} = do   let n = length us   us' <- do@@ -108,17 +108,17 @@     case splitAt i us of       (xs, _ : ys) -> return $ xs ++ (u' : ys)       _ -> error "impossible"-  ((b, q), vs) <- runWriterT $ weighted $ Weighted.hoist (WriterT . Free.density us') m+  ((b, q), vs) <- runWriterT $ runWeightedT $ WeightedT.hoist (WriterT . Free.runDensityT us') m   let ratio = (exp . ln) $ min 1 (q * fromIntegral n / (p * fromIntegral (length vs)))   accept <- bernoulli ratio   return if accept then MHResult True (Trace vs b q) else MHResult False t  -- | A variant of 'mhTrans' with an external sampling monad.-mhTrans' :: MonadDistribution m => Weighted (Free.Density Identity) a -> Trace a -> m (Trace a)-mhTrans' m = mhTransFree (Weighted.hoist (Free.hoist (return . runIdentity)) m)+mhTrans' :: (MonadDistribution m) => WeightedT (Free.DensityT Identity) a -> Trace a -> m (Trace a)+mhTrans' m = mhTransFree (WeightedT.hoist (Free.hoist (return . runIdentity)) m)  -- | burn in an MCMC chain for n steps (which amounts to dropping samples of the end of the list)-burnIn :: Functor m => Int -> m [a] -> m [a]+burnIn :: (Functor m) => Int -> m [a] -> m [a] burnIn n = fmap dropEnd   where     dropEnd ls = let len = length ls in take (len - n) ls

@@ -9,7 +9,7 @@ -- Stability   : experimental -- Portability : GHC module Control.Monad.Bayes.Traced.Dynamic-  ( Traced,+  ( TracedT,     hoist,     marginal,     freeze,@@ -24,7 +24,7 @@     MonadFactor (..),     MonadMeasure,   )-import Control.Monad.Bayes.Density.Free (Density)+import Control.Monad.Bayes.Density.Free (DensityT) import Control.Monad.Bayes.Traced.Common   ( Trace (..),     bind,@@ -32,75 +32,75 @@     scored,     singleton,   )-import Control.Monad.Bayes.Weighted (Weighted)+import Control.Monad.Bayes.Weighted (WeightedT) import Control.Monad.Trans (MonadTrans (..)) import Data.List.NonEmpty as NE (NonEmpty ((:|)), toList)  -- | A tracing monad where only a subset of random choices are traced and this -- subset can be adjusted dynamically.-newtype Traced m a = Traced {runTraced :: m (Weighted (Density m) a, Trace a)}+newtype TracedT m a = TracedT {runTraced :: m (WeightedT (DensityT m) a, Trace a)} -pushM :: Monad m => m (Weighted (Density m) a) -> Weighted (Density m) a+pushM :: (Monad m) => m (WeightedT (DensityT m) a) -> WeightedT (DensityT m) a pushM = join . lift . lift -instance Monad m => Functor (Traced m) where-  fmap f (Traced c) = Traced $ do+instance (Monad m) => Functor (TracedT m) where+  fmap f (TracedT c) = TracedT $ do     (m, t) <- c     let m' = fmap f m     let t' = fmap f t     return (m', t') -instance Monad m => Applicative (Traced m) where-  pure x = Traced $ pure (pure x, pure x)-  (Traced cf) <*> (Traced cx) = Traced $ do+instance (Monad m) => Applicative (TracedT m) where+  pure x = TracedT $ pure (pure x, pure x)+  (TracedT cf) <*> (TracedT cx) = TracedT $ do     (mf, tf) <- cf     (mx, tx) <- cx     return (mf <*> mx, tf <*> tx) -instance Monad m => Monad (Traced m) where-  (Traced cx) >>= f = Traced $ do+instance (Monad m) => Monad (TracedT m) where+  (TracedT cx) >>= f = TracedT $ do     (mx, tx) <- cx     let m = mx >>= pushM . fmap fst . runTraced . f     t <- return tx `bind` (fmap snd . runTraced . f)     return (m, t) -instance MonadTrans Traced where-  lift m = Traced $ fmap ((,) (lift $ lift m) . pure) m+instance MonadTrans TracedT where+  lift m = TracedT $ fmap ((,) (lift $ lift m) . pure) m -instance MonadDistribution m => MonadDistribution (Traced m) where-  random = Traced $ fmap ((,) random . singleton) random+instance (MonadDistribution m) => MonadDistribution (TracedT m) where+  random = TracedT $ fmap ((,) random . singleton) random -instance MonadFactor m => MonadFactor (Traced m) where-  score w = Traced $ fmap (score w,) (score w >> pure (scored w))+instance (MonadFactor m) => MonadFactor (TracedT m) where+  score w = TracedT $ fmap (score w,) (score w >> pure (scored w)) -instance MonadMeasure m => MonadMeasure (Traced m)+instance (MonadMeasure m) => MonadMeasure (TracedT m) -hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a-hoist f (Traced c) = Traced (f c)+hoist :: (forall x. m x -> m x) -> TracedT m a -> TracedT m a+hoist f (TracedT c) = TracedT (f c)  -- | Discard the trace and supporting infrastructure.-marginal :: Monad m => Traced m a -> m a-marginal (Traced c) = fmap (output . snd) c+marginal :: (Monad m) => TracedT m a -> m a+marginal (TracedT c) = fmap (output . snd) c  -- | Freeze all traced random choices to their current values and stop tracing -- them.-freeze :: Monad m => Traced m a -> Traced m a-freeze (Traced c) = Traced $ do+freeze :: (Monad m) => TracedT m a -> TracedT m a+freeze (TracedT c) = TracedT $ do   (_, t) <- c   let x = output t   return (return x, pure x)  -- | A single step of the Trace Metropolis-Hastings algorithm.-mhStep :: MonadDistribution m => Traced m a -> Traced m a-mhStep (Traced c) = Traced $ do+mhStep :: (MonadDistribution m) => TracedT m a -> TracedT m a+mhStep (TracedT c) = TracedT $ do   (m, t) <- c   t' <- mhTransFree m t   return (m, t')  -- | Full run of the Trace Metropolis-Hastings algorithm with a specified -- number of steps.-mh :: MonadDistribution m => Int -> Traced m a -> m [a]-mh n (Traced c) = do+mh :: (MonadDistribution m) => Int -> TracedT m a -> m [a]+mh n (TracedT c) = do   (m, t) <- c   let f k         | k <= 0 = return (t :| [])

@@ -10,7 +10,7 @@ -- Stability   : experimental -- Portability : GHC module Control.Monad.Bayes.Traced.Static-  ( Traced (..),+  ( TracedT (..),     hoist,     marginal,     mhStep,@@ -24,7 +24,7 @@     MonadFactor (..),     MonadMeasure,   )-import Control.Monad.Bayes.Density.Free (Density)+import Control.Monad.Bayes.Density.Free (DensityT) import Control.Monad.Bayes.Traced.Common   ( Trace (..),     bind,@@ -32,7 +32,7 @@     scored,     singleton,   )-import Control.Monad.Bayes.Weighted (Weighted)+import Control.Monad.Bayes.Weighted (WeightedT) import Control.Monad.Trans (MonadTrans (..)) import Data.List.NonEmpty as NE (NonEmpty ((:|)), toList) @@ -40,45 +40,45 @@ -- -- The random choices that are not to be traced should be lifted from the -- transformed monad.-data Traced m a = Traced-  { model :: Weighted (Density m) a,+data TracedT m a = TracedT+  { model :: WeightedT (DensityT m) a,     traceDist :: m (Trace a)   } -instance Monad m => Functor (Traced m) where-  fmap f (Traced m d) = Traced (fmap f m) (fmap (fmap f) d)+instance (Monad m) => Functor (TracedT m) where+  fmap f (TracedT m d) = TracedT (fmap f m) (fmap (fmap f) d) -instance Monad m => Applicative (Traced m) where-  pure x = Traced (pure x) (pure (pure x))-  (Traced mf df) <*> (Traced mx dx) = Traced (mf <*> mx) (liftA2 (<*>) df dx)+instance (Monad m) => Applicative (TracedT m) where+  pure x = TracedT (pure x) (pure (pure x))+  (TracedT mf df) <*> (TracedT mx dx) = TracedT (mf <*> mx) (liftA2 (<*>) df dx) -instance Monad m => Monad (Traced m) where-  (Traced mx dx) >>= f = Traced my dy+instance (Monad m) => Monad (TracedT m) where+  (TracedT mx dx) >>= f = TracedT my dy     where       my = mx >>= model . f       dy = dx `bind` (traceDist . f) -instance MonadTrans Traced where-  lift m = Traced (lift $ lift m) (fmap pure m)+instance MonadTrans TracedT where+  lift m = TracedT (lift $ lift m) (fmap pure m) -instance MonadDistribution m => MonadDistribution (Traced m) where-  random = Traced random (fmap singleton random)+instance (MonadDistribution m) => MonadDistribution (TracedT m) where+  random = TracedT random (fmap singleton random) -instance MonadFactor m => MonadFactor (Traced m) where-  score w = Traced (score w) (score w >> pure (scored w))+instance (MonadFactor m) => MonadFactor (TracedT m) where+  score w = TracedT (score w) (score w >> pure (scored w)) -instance MonadMeasure m => MonadMeasure (Traced m)+instance (MonadMeasure m) => MonadMeasure (TracedT m) -hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a-hoist f (Traced m d) = Traced m (f d)+hoist :: (forall x. m x -> m x) -> TracedT m a -> TracedT m a+hoist f (TracedT m d) = TracedT m (f d)  -- | Discard the trace and supporting infrastructure.-marginal :: Monad m => Traced m a -> m a-marginal (Traced _ d) = fmap output d+marginal :: (Monad m) => TracedT m a -> m a+marginal (TracedT _ d) = fmap output d  -- | A single step of the Trace Metropolis-Hastings algorithm.-mhStep :: MonadDistribution m => Traced m a -> Traced m a-mhStep (Traced m d) = Traced m d'+mhStep :: (MonadDistribution m) => TracedT m a -> TracedT m a+mhStep (TracedT m d) = TracedT m d'   where     d' = d >>= mhTransFree m @@ -111,8 +111,8 @@ -- [True,True,True] -- -- Of course, it will need to be run more than twice to get a reasonable estimate.-mh :: MonadDistribution m => Int -> Traced m a -> m [a]-mh n (Traced m d) = fmap (map output . NE.toList) (f n)+mh :: (MonadDistribution m) => Int -> TracedT m a -> m [a]+mh n (TracedT m d) = fmap (map output . NE.toList) (f n)   where     f k       | k <= 0 = fmap (:| []) d

@@ -11,17 +11,16 @@ -- Stability   : experimental -- Portability : GHC ----- 'Weighted' is an instance of 'MonadFactor'. Apply a 'MonadDistribution' transformer to+-- 'WeightedT' is an instance of 'MonadFactor'. Apply a 'MonadDistribution' transformer to -- obtain a 'MonadMeasure' that can execute probabilistic models. module Control.Monad.Bayes.Weighted-  ( Weighted,-    withWeight,-    weighted,+  ( WeightedT,+    weightedT,     extractWeight,     unweighted,     applyWeight,     hoist,-    runWeighted,+    runWeightedT,   ) where @@ -35,46 +34,45 @@ import Numeric.Log (Log)  -- | Execute the program using the prior distribution, while accumulating likelihood.-newtype Weighted m a = Weighted (StateT (Log Double) m a)+newtype WeightedT m a = WeightedT (StateT (Log Double) m a)   -- StateT is more efficient than WriterT   deriving newtype (Functor, Applicative, Monad, MonadIO, MonadTrans, MonadDistribution) -instance Monad m => MonadFactor (Weighted m) where-  score w = Weighted (modify (* w))+instance (Monad m) => MonadFactor (WeightedT m) where+  score w = WeightedT (modify (* w)) -instance MonadDistribution m => MonadMeasure (Weighted m)+instance (MonadDistribution m) => MonadMeasure (WeightedT m)  -- | Obtain an explicit value of the likelihood for a given value.-weighted, runWeighted :: Weighted m a -> m (a, Log Double)-weighted (Weighted m) = runStateT m 1-runWeighted = weighted+runWeightedT :: WeightedT m a -> m (a, Log Double)+runWeightedT (WeightedT m) = runStateT m 1  -- | Compute the sample and discard the weight. -- -- This operation introduces bias.-unweighted :: Functor m => Weighted m a -> m a-unweighted = fmap fst . weighted+unweighted :: (Functor m) => WeightedT m a -> m a+unweighted = fmap fst . runWeightedT  -- | Compute the weight and discard the sample.-extractWeight :: Functor m => Weighted m a -> m (Log Double)-extractWeight = fmap snd . weighted+extractWeight :: (Functor m) => WeightedT m a -> m (Log Double)+extractWeight = fmap snd . runWeightedT  -- | Embed a random variable with explicitly given likelihood. ----- > weighted . withWeight = id-withWeight :: (Monad m) => m (a, Log Double) -> Weighted m a-withWeight m = Weighted $ do+-- > runWeightedT . weightedT = id+weightedT :: (Monad m) => m (a, Log Double) -> WeightedT m a+weightedT m = WeightedT $ do   (x, w) <- lift m   modify (* w)   return x  -- | Use the weight as a factor in the transformed monad.-applyWeight :: MonadFactor m => Weighted m a -> m a+applyWeight :: (MonadFactor m) => WeightedT m a -> m a applyWeight m = do-  (x, w) <- weighted m+  (x, w) <- runWeightedT m   factor w   return x  -- | Apply a transformation to the transformed monad.-hoist :: (forall x. m x -> n x) -> Weighted m a -> Weighted n a-hoist t (Weighted m) = Weighted $ mapStateT t m+hoist :: (forall x. m x -> n x) -> WeightedT m a -> WeightedT n a+hoist t (WeightedT m) = WeightedT $ mapStateT t m

@@ -53,7 +53,7 @@ -- of solutions corrensdonded to times that was requested. -- It takes Vector of time points as a parameter and returns a vector of results integrateV ::-  PrimMonad m =>+  (PrimMonad m) =>   -- | Internal integrator   Integrator a ->   -- | initial  value

@@ -5,6 +5,7 @@ import Test.Hspec.QuickCheck (prop) import Test.QuickCheck (ioProperty, property, (==>)) import TestAdvanced qualified+import TestBenchmarks qualified import TestDistribution qualified import TestEnumerator qualified import TestInference qualified@@ -12,6 +13,7 @@ import TestPipes (hmms) import TestPipes qualified import TestPopulation qualified+import TestSSMFixtures qualified import TestSampler qualified import TestSequential qualified import TestStormerVerlet qualified@@ -55,7 +57,7 @@         -- Because of rounding issues, require the variance to be a bit bigger than 0         -- See https://github.com/tweag/monad-bayes/issues/275         var > 0.1 ==> property $ TestIntegrator.normalVariance mean (sqrt var) ~== var-  describe "Sampler mean and variance" do+  describe "SamplerT mean and variance" do     it "gets right mean and variance" $       TestSampler.testMeanAndVariance `shouldBe` True   describe "Integrator Volume" do@@ -166,3 +168,6 @@       passed6 `shouldBe` True       passed7 <- TestAdvanced.passed7       passed7 `shouldBe` True++  TestBenchmarks.test+  TestSSMFixtures.test

@@ -15,7 +15,7 @@ mcmcConfig :: MCMCConfig mcmcConfig = MCMCConfig {numMCMCSteps = 0, numBurnIn = 0, proposal = SingleSiteMH} -smcConfig :: MonadDistribution m => SMCConfig m+smcConfig :: (MonadDistribution m) => SMCConfig m smcConfig = SMCConfig {numSteps = 0, numParticles = 1000, resampler = resampleMultinomial}  passed1, passed2, passed3, passed4, passed5, passed6, passed7 :: IO Bool@@ -23,16 +23,16 @@   sample <- sampleIOfixed $ mcmc MCMCConfig {numMCMCSteps = 10000, numBurnIn = 5000, proposal = SingleSiteMH} random   return $ abs (0.5 - (expectation id $ fromList $ toEmpirical sample)) < 0.01 passed2 = do-  sample <- sampleIOfixed $ population $ smc (SMCConfig {numSteps = 0, numParticles = 10000, resampler = resampleMultinomial}) random+  sample <- sampleIOfixed $ runPopulationT $ smc (SMCConfig {numSteps = 0, numParticles = 10000, resampler = resampleMultinomial}) random   return $ close 0.5 sample passed3 = do-  sample <- sampleIOfixed $ population $ rmsmcDynamic mcmcConfig smcConfig random+  sample <- sampleIOfixed $ runPopulationT $ rmsmcDynamic mcmcConfig smcConfig random   return $ close 0.5 sample passed4 = do-  sample <- sampleIOfixed $ population $ rmsmcBasic mcmcConfig smcConfig random+  sample <- sampleIOfixed $ runPopulationT $ rmsmcBasic mcmcConfig smcConfig random   return $ close 0.5 sample passed5 = do-  sample <- sampleIOfixed $ population $ rmsmc mcmcConfig smcConfig random+  sample <- sampleIOfixed $ runPopulationT $ rmsmc mcmcConfig smcConfig random   return $ close 0.5 sample passed6 = do   sample <-@@ -48,7 +48,7 @@ close :: Double -> [(Double, Log Double)] -> Bool  passed7 = do-  sample <- fmap join $ sampleIOfixed $ fmap (fmap (\(x, y) -> fmap (second (* y)) x)) $ population $ smc2 0 100 100 100 random (normal 0)+  sample <- fmap join $ sampleIOfixed $ fmap (fmap (\(x, y) -> fmap (second (* y)) x)) $ runPopulationT $ smc2 0 100 100 100 random (normal 0)   return $ close 0.0 sample  close n sample = abs (n - (expectation id $ fromList $ toEmpiricalWeighted sample)) < 0.01

@@ -0,0 +1,33 @@+module TestBenchmarks where++import Control.Monad (forM_)+import Data.Maybe (fromJust)+import Helper+import System.IO (readFile')+import System.IO.Error (catchIOError, isDoesNotExistError)+import Test.Hspec++fixtureToFilename :: Model -> Alg -> String+fixtureToFilename model alg = fromJust (serializeModel model) ++ "-" ++ show alg ++ ".txt"++models :: [Model]+models = [LR 10, HMM 10, LDA (5, 10)]++algs :: [Alg]+algs = [minBound .. maxBound]++test :: SpecWith ()+test = describe "Benchmarks" $ forM_ models $ \model -> forM_ algs $ testFixture model++testFixture :: Model -> Alg -> SpecWith ()+testFixture model alg = do+  let filename = "test/fixtures/" ++ fixtureToFilename model alg+  it ("should agree with the fixture " ++ filename) $ do+    fixture <- catchIOError (readFile' filename) $ \e ->+      if isDoesNotExistError e+        then return ""+        else ioError e+    sampled <- runAlgFixed model alg+    -- Reset in case of fixture update or creation+    writeFile filename sampled+    fixture `shouldBe` sampled

@@ -15,13 +15,13 @@ import Numeric.Log (Log (ln)) import Sprinkler (hard, soft) -unnorm :: MonadDistribution m => m Int+unnorm :: (MonadDistribution m) => m Int unnorm = categorical $ V.fromList [0.5, 0.8]  passed1 :: Bool passed1 = (exp . ln) (evidence unnorm) ~== 1 -agg :: MonadDistribution m => m Int+agg :: (MonadDistribution m) => m Int agg = do   x <- uniformD [0, 1]   y <- uniformD [2, 1]

@@ -20,33 +20,33 @@ import Control.Monad.Bayes.Integrator (normalize) import Control.Monad.Bayes.Integrator qualified as Integrator import Control.Monad.Bayes.Population-import Control.Monad.Bayes.Sampler.Strict (Sampler, sampleIOfixed)+import Control.Monad.Bayes.Sampler.Strict (SamplerT, sampleIOfixed) import Control.Monad.Bayes.Sampler.Strict qualified as Sampler-import Control.Monad.Bayes.Weighted (Weighted)-import Control.Monad.Bayes.Weighted qualified as Weighted+import Control.Monad.Bayes.Weighted (WeightedT)+import Control.Monad.Bayes.Weighted qualified as WeightedT import Data.AEq (AEq ((~==))) import Numeric.Log (Log) import Sprinkler (soft) import System.Random.Stateful (IOGenM, StdGen) -sprinkler :: MonadMeasure m => m Bool+sprinkler :: (MonadMeasure m) => m Bool sprinkler = Sprinkler.soft  -- | Count the number of particles produced by SMC checkParticles :: Int -> Int -> IO Int checkParticles observations particles =-  sampleIOfixed (fmap length (population $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleMultinomial} Sprinkler.soft))+  sampleIOfixed (fmap length (runPopulationT $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleMultinomial} Sprinkler.soft))  checkParticlesSystematic :: Int -> Int -> IO Int checkParticlesSystematic observations particles =-  sampleIOfixed (fmap length (population $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleSystematic} Sprinkler.soft))+  sampleIOfixed (fmap length (runPopulationT $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleSystematic} Sprinkler.soft))  checkParticlesStratified :: Int -> Int -> IO Int checkParticlesStratified observations particles =-  sampleIOfixed (fmap length (population $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleStratified} Sprinkler.soft))+  sampleIOfixed (fmap length (runPopulationT $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleStratified} Sprinkler.soft))  checkTerminateSMC :: IO [(Bool, Log Double)]-checkTerminateSMC = sampleIOfixed (population $ smc SMCConfig {numSteps = 2, numParticles = 5, resampler = resampleMultinomial} sprinkler)+checkTerminateSMC = sampleIOfixed (runPopulationT $ smc SMCConfig {numSteps = 2, numParticles = 5, resampler = resampleMultinomial} sprinkler)  checkPreserveSMC :: Bool checkPreserveSMC =@@ -54,8 +54,8 @@     ~== enumerator sprinkler  expectationNearNumeric ::-  Weighted Integrator.Integrator Double ->-  Weighted Integrator.Integrator Double ->+  WeightedT Integrator.Integrator Double ->+  WeightedT Integrator.Integrator Double ->   Double expectationNearNumeric x y =   let e1 = Integrator.expectation $ normalize x@@ -63,12 +63,12 @@    in (abs (e1 - e2))  expectationNearSampling ::-  Weighted (Sampler (IOGenM StdGen) IO) Double ->-  Weighted (Sampler (IOGenM StdGen) IO) Double ->+  WeightedT (SamplerT (IOGenM StdGen) IO) Double ->+  WeightedT (SamplerT (IOGenM StdGen) IO) Double ->   IO Double expectationNearSampling x y = do-  e1 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ Weighted.weighted x-  e2 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ Weighted.weighted y+  e1 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ WeightedT.runWeightedT x+  e2 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ WeightedT.runWeightedT y   return (abs (e1 - e2))  testNormalNormal :: [Double] -> IO Bool

@@ -13,7 +13,7 @@   ) import Control.Monad.Bayes.Integrator import Control.Monad.Bayes.Sampler.Strict-import Control.Monad.Bayes.Weighted (weighted)+import Control.Monad.Bayes.Weighted (runWeightedT) import Control.Monad.ST (runST) import Data.AEq (AEq ((~==))) import Data.List (sortOn)@@ -32,7 +32,7 @@ volumeIsOne :: [Double] -> Bool volumeIsOne = (~== 1.0) . volume . uniformD -agg :: MonadDistribution m => m Int+agg :: (MonadDistribution m) => m Int agg = do   x <- uniformD [0, 1]   y <- uniformD [2, 1]@@ -98,7 +98,7 @@ passed8 =   1     == ( volume $-           fmap (ln . exp . snd) $ weighted do+           fmap (ln . exp . snd) $ runWeightedT do              x <- bernoulli 0.5              factor $ if x then 0.2 else 0.1              return x@@ -137,11 +137,11 @@     ~== 1 -- sampler and integrator agree on a non-trivial model passed14 =-  let sample = runST $ sampleSTfixed $ fmap sampleMean $ replicateM 10000 $ weighted $ model1+  let sample = runST $ sampleSTfixed $ fmap sampleMean $ replicateM 10000 $ runWeightedT $ model1       quadrature = expectation $ normalize $ model1    in abs (sample - quadrature) < 0.01 -model1 :: MonadMeasure m => m Double+model1 :: (MonadMeasure m) => m Double model1 = do   x <- random   y <- random

@@ -3,20 +3,20 @@ import Control.Monad.Bayes.Class (MonadDistribution, MonadMeasure) import Control.Monad.Bayes.Enumerator (enumerator, expectation) import Control.Monad.Bayes.Population as Population-  ( Population,+  ( PopulationT,     collapse,     popAvg,-    population,     pushEvidence,     resampleMultinomial,+    runPopulationT,     spawn,   ) import Control.Monad.Bayes.Sampler.Strict (sampleIOfixed) import Data.AEq (AEq ((~==))) import Sprinkler (soft) -weightedSampleSize :: MonadDistribution m => Population m a -> m Int-weightedSampleSize = fmap length . population+weightedSampleSize :: (MonadDistribution m) => PopulationT m a -> m Int+weightedSampleSize = fmap length . runPopulationT  popSize :: IO Int popSize =@@ -26,7 +26,7 @@ manySize =   sampleIOfixed (weightedSampleSize $ spawn 5 >> sprinkler >> spawn 3) -sprinkler :: MonadMeasure m => m Bool+sprinkler :: (MonadMeasure m) => m Bool sprinkler = Sprinkler.soft  sprinklerExact :: [(Bool, Double)]

@@ -0,0 +1,28 @@+module TestSSMFixtures where++import Control.Monad.Bayes.Sampler.Strict (sampleIOfixed)+import NonlinearSSM+import NonlinearSSM.Algorithms+import System.IO (readFile')+import System.IO.Error (catchIOError, isDoesNotExistError)+import Test.Hspec++fixtureToFilename :: Alg -> FilePath+fixtureToFilename alg = "test/fixtures/SSM-" ++ show alg ++ ".txt"++testFixture :: Alg -> SpecWith ()+testFixture alg = do+  let filename = fixtureToFilename alg+  it ("should agree with the fixture " ++ filename) $ do+    ys <- sampleIOfixed $ generateData t+    fixture <- catchIOError (readFile' filename) $ \e ->+      if isDoesNotExistError e+        then return ""+        else ioError e+    sampled <- sampleIOfixed $ runAlgFixed (map fst ys) alg+    -- Reset in case of fixture update or creation+    writeFile filename sampled+    fixture `shouldBe` sampled++test :: SpecWith ()+test = describe "TestSSMFixtures" $ mapM_ testFixture algs

@@ -10,7 +10,7 @@ import Data.AEq (AEq ((~==))) import Sprinkler (soft) -twoSync :: MonadMeasure m => m Int+twoSync :: (MonadMeasure m) => m Int twoSync = do   x <- uniformD [0, 1]   factor (fromIntegral x)@@ -18,7 +18,7 @@   factor (fromIntegral y)   return (x + y) -finishedTwoSync :: MonadMeasure m => Int -> m Bool+finishedTwoSync :: (MonadMeasure m) => Int -> m Bool finishedTwoSync n = finished (run n twoSync)   where     run 0 d = d@@ -30,7 +30,7 @@ checkTwoSync 2 = mass (finishedTwoSync 2) True ~== 1 checkTwoSync _ = error "Unexpected argument" -sprinkler :: MonadMeasure m => m Bool+sprinkler :: (MonadMeasure m) => m Bool sprinkler = Sprinkler.soft  checkPreserve :: Bool@@ -42,7 +42,7 @@ pFinished 2 = 1 pFinished _ = error "Unexpected argument" -isFinished :: MonadMeasure m => Int -> m Bool+isFinished :: (MonadMeasure m) => Int -> m Bool isFinished n = finished (run n sprinkler)   where     run 0 d = d

@@ -8,13 +8,13 @@     factor,   ) import Control.Monad.Bayes.Sampler.Strict (sampleIOfixed)-import Control.Monad.Bayes.Weighted (weighted)+import Control.Monad.Bayes.Weighted (runWeightedT) import Control.Monad.State (unless, when) import Data.AEq (AEq ((~==))) import Data.Bifunctor (second) import Numeric.Log (Log (Exp, ln)) -model :: MonadMeasure m => m (Int, Double)+model :: (MonadMeasure m) => m (Int, Double) model = do   n <- uniformD [0, 1, 2]   unless (n == 0) (factor 0.5)@@ -22,8 +22,8 @@   when (n == 2) (factor $ (Exp . log) (x * x))   return (n, x) -result :: MonadDistribution m => m ((Int, Double), Double)-result = second (exp . ln) <$> weighted model+result :: (MonadDistribution m) => m ((Int, Double), Double)+result = second (exp . ln) <$> runWeightedT model  passed :: IO Bool passed = fmap check (sampleIOfixed result)