diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,11 +1,52 @@
-# Unreleased
+# 1.0.0 (2022-09-10)
 
-- No unreleased changes so far.
+- host website from repo
+- host notebooks from repo
+- use histogram-fill
 
-# 0.1.1.0 (2020-04-08)
+# 0.2.0 (2022-07-26)
 
+- rename various functions to match the names of the corresponding types (e.g. `Enumerator` goes with `enumerator`)
+- add configs as arguments to inference methods `smc` and `mcmc`
+- add rudimentary tests for all inference methods
+- put `mcmc` as inference method in new module `Control.Monad.Bayes.Inference.MCMC`
+- update history of changelog in line with semantic versioning conventions
+- bumped to GHC 9.2.3
+
+# 0.1.5 (2022-07-26)
+
+- Refactor of sampler to be parametric in the choice of a pair of IO monad and RNG
+
+# 0.1.4 (2022-06-15)
+
+Addition of new helper functions, plotting tools, tests, and Integrator monad.
+
+- helpers include: `toEmpirical` (list of samples to empirical distribution) and `toBins` (simple histogramming)
+- `Integrator` is an instance of `MonadSample` for numerical integration
+- `notebooks` now contains working notebook-based tutorials and examples
+- new tests, including with conjugate distributions to compare analytic solution against inferred posterior
+- `models` directory is cleaned up. New sequential models using `pipes` package to represent monadic streams
+
+# 0.1.3 (2022-06-08)
+
+Clean up of unused functions and broken code
+
+- remove unused functions in `Weighted` and `Population`
+- remove broken models in `models`
+- explicit imports
+- added some global language extensions
+
+# 0.1.2 (2022-06-08)
+
+Add documentation
+
+- docs written in markdown
+- docs built by sphinx
+
+# 0.1.1 (2020-04-08)
+
 - New exported function: `Control.Monad.Bayes.Class` now exports `discrete`.
 
-# 0.1.0.0 (2020-02-17)
+# 0.1.0 (2020-02-17)
 
 Initial release.
diff --git a/benchmark/SSM.hs b/benchmark/SSM.hs
--- a/benchmark/SSM.hs
+++ b/benchmark/SSM.hs
@@ -1,29 +1,43 @@
 module Main where
 
-import Control.Monad.Bayes.Inference.PMMH as PMMH
-import Control.Monad.Bayes.Inference.RMSMC
+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
+import Control.Monad.Bayes.Inference.SMC2 as SMC2 (smc2)
 import Control.Monad.Bayes.Population
-import Control.Monad.Bayes.Sampler
-import Control.Monad.Bayes.Weighted
-import Control.Monad.IO.Class
-import NonlinearSSM
+import Control.Monad.Bayes.Population (population, resampleMultinomial)
+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 System.Random.Stateful (mkStdGen, newIOGenM)
 
 main :: IO ()
-main = sampleIO $ do
+main = sampleIOfixed $ do
   let t = 5
   dat <- generateData t
   let ys = map snd dat
   liftIO $ print "SMC"
-  smcRes <- runPopulation $ smcMultinomial t 10 (param >>= model ys)
+  smcRes <- population $ smc SMCConfig {numSteps = t, numParticles = 10, resampler = resampleMultinomial} (param >>= model ys)
   liftIO $ print $ show smcRes
   liftIO $ print "RM-SMC"
-  smcrmRes <- runPopulation $ rmsmcLocal t 10 10 (param >>= model ys)
+  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 <- prior $ pmmh 2 t 3 param (model ys)
+  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 <- runPopulation $ smc2 t 3 2 1 param (model ys)
+  smc2Res <- population $ smc2 t 3 2 1 param (model ys)
   liftIO $ print $ show smc2Res
diff --git a/benchmark/Single.hs b/benchmark/Single.hs
--- a/benchmark/Single.hs
+++ b/benchmark/Single.hs
@@ -1,19 +1,38 @@
-import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Inference.RMSMC
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE ImportQualifiedPost #-}
+
+import Control.Monad.Bayes.Class (MonadInfer)
+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
-import Control.Monad.Bayes.Traced
+import Control.Monad.Bayes.Sampler.Strict
+import Control.Monad.Bayes.Traced hiding (model)
 import Control.Monad.Bayes.Weighted
-import Data.Time
-import qualified HMM
-import qualified LDA
-import qualified LogReg
+import Control.Monad.ST (runST)
+import Data.Time (diffUTCTime, getCurrentTime)
+import HMM qualified
+import LDA qualified
+import LogReg qualified
 import Options.Applicative
-import System.Random.MWC (createSystemRandom)
+  ( Applicative (liftA2),
+    ParserInfo,
+    auto,
+    execParser,
+    fullDesc,
+    help,
+    info,
+    long,
+    maybeReader,
+    option,
+    short,
+  )
 
 data Model = LR Int | HMM Int | LDA (Int, Int)
-  deriving (Show, Read)
+  deriving stock (Show, Read)
 
 parseModel :: String -> Maybe Model
 parseModel s =
@@ -29,14 +48,12 @@
     size (LR n) = n
     size (HMM n) = n
     size (LDA (d, w)) = d * w
-    synthesize :: SamplerST a -> (a -> b) -> b
-    synthesize dataGen prog = prog (sampleSTfixed dataGen)
-    program (LR n) = show <$> synthesize (LogReg.syntheticData n) LogReg.logisticRegression
-    program (HMM n) = show <$> synthesize (HMM.syntheticData n) HMM.hmm
-    program (LDA (d, w)) = show <$> synthesize (LDA.syntheticData d w) LDA.lda
+    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 (Read, Show)
+  deriving stock (Read, Show)
 
 runAlg :: Model -> Alg -> SamplerIO String
 runAlg model alg =
@@ -44,21 +61,20 @@
     SMC ->
       let n = 100
           (k, m) = getModel model
-       in show <$> runPopulation (smcSystematic k n m)
+       in show <$> population (smc SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic} m)
     MH ->
       let t = 100
           (_, m) = getModel model
-       in show <$> prior (mh t m)
+       in show <$> unweighted (mh t m)
     RMSMC ->
       let n = 10
           t = 1
           (k, m) = getModel model
-       in show <$> runPopulation (rmsmcBasic k n t m)
+       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
-  g <- createSystemRandom
-  x <- sampleIOwith (runAlg model alg) g
+  x <- sampleIOfixed (runAlg model alg)
   print x
 
 opts :: ParserInfo (Model, Alg)
diff --git a/benchmark/Speed.hs b/benchmark/Speed.hs
--- a/benchmark/Speed.hs
+++ b/benchmark/Speed.hs
@@ -1,123 +1,38 @@
-import Control.Monad (unless)
-import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Inference.RMSMC
-import Control.Monad.Bayes.Inference.SMC
-import Control.Monad.Bayes.Population
-import Control.Monad.Bayes.Sampler
-import Control.Monad.Bayes.Traced
-import Control.Monad.Bayes.Weighted
-import Criterion.Main
-import Criterion.Types
-import GHC.IO.Handle
--- import NonlinearSSM
-import qualified HMM
-import qualified LDA
-import qualified LogReg
-import System.Exit
-import System.Process hiding (env)
-import System.Random.MWC (GenIO, createSystemRandom)
-
--- | Path to the Anglican project with benchmarks.
-anglicanPath :: String
-anglicanPath = "/scratch/ams240/repos/anglican-white-paper/experiments"
-
--- | Running Leiningen repl process.
--- Leiningen takes a lot of time to start so we keep a repl
--- running to speed up benchmarking.
--- Contains in order input handle, output handle, and process handle.
-data LeinProc = LeinProc Handle Handle ProcessHandle
-
-anglicanModelName :: Model -> String
-anglicanModelName (LR _) = "logisticRegression"
-anglicanModelName (HMM _) = "hmm"
-anglicanModelName (LDA _) = "lda"
-
-clojureBool :: Bool -> String
-clojureBool False = "false"
-clojureBool True = "true"
-
--- | Format Haskell [Bool] as a Clojure Boolean vector.
-clojureBoolVector :: [Bool] -> String
-clojureBoolVector = clojureVector . map clojureBool
-
--- | Format a Haskell list as a Clojure vector.
-clojureVector :: [String] -> String
-clojureVector xs = "[" ++ unwords xs ++ "]"
-
--- | Format a Haskell list as a Clojure vector.
-clojureShowVector :: Show a => [a] -> String
-clojureShowVector xs = clojureVector $ map show xs
-
--- | Insert data into an Anglican model.
-anglicanData :: LeinProc -> Model -> IO ()
-anglicanData leinProc model = do
-  anglican leinProc ["(use 'nstools.ns)\n"]
-  anglican leinProc ["(ns+ " ++ anglicanModelName model ++ ")\n"]
-  case model of
-    LR dataset -> do
-      let (xs, labels) = unzip dataset
-      anglican leinProc ["(ns-unmap *ns* 'xs)\n"]
-      anglican leinProc ["(def xs " ++ clojureShowVector xs ++ ")\n", "nil\n"]
-      anglican leinProc ["(ns-unmap *ns* 'labels)\n"]
-      anglican leinProc ["(def labels " ++ clojureBoolVector labels ++ ")\n", "nil\n"]
-    HMM observations -> do
-      anglican leinProc ["(ns-unmap *ns* 'observations)\n"]
-      anglican leinProc ["(def observations " ++ clojureShowVector observations ++ ")\n", "nil\n"]
-    LDA docs -> do
-      anglican leinProc ["(ns-unmap *ns* 'docs)\n"]
-      anglican leinProc ["(def docs " ++ clojureVector (map clojureShowVector docs) ++ ")\n", "nil\n"]
-  anglican leinProc ["(use 'nstools.ns)\n"]
-  anglican leinProc ["(ns+ anglican.core)\n"]
-
-anglican :: LeinProc -> [String] -> IO ()
-anglican (LeinProc input output _) cmds = do
-  -- execute command
-  mapM_ (hPutStr input) cmds
-  hFlush input
-  -- wait until all samples are produced
-  waitForAnglican output
-
--- | Wait until an Anglican program finishes
-waitForAnglican :: Handle -> IO ()
-waitForAnglican handle = run
-  where
-    run = do
-      l <- hGetLine handle
-      -- traceIO l
-      -- We recognize that Anglican is finished when the repl emits a line "nil"
-      unless (l == "nil") run
-
--- | Start a Leiningen process in a directory that contains benchmarks.
-startLein :: IO LeinProc
-startLein = do
-  let setup = (shell "lein repl") {cwd = Just anglicanPath, std_in = CreatePipe, std_out = CreatePipe}
-  (Just input, Just output, _, process) <- createProcess setup
-  -- wait until Leiningen starts producing output to make sure it's ready
-  _ <- hWaitForInput output (-1) -- wait for output indefinitely
-  let leinProc = LeinProc input output process
-  anglican leinProc ["(use 'nstools.ns)\n"]
-  return leinProc
-
--- | Path to the WebPPL project with benchmarks.
-webpplPath :: String
-webpplPath = "/scratch/ams240/repos/anglican-white-paper/experiments/WebPPL"
-
--- | Format Haskell list as a Javascript list.
-javascriptList :: Show a => [a] -> String
-javascriptList xs = unwords (map show xs)
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# OPTIONS_GHC -Wall #-}
 
-webpplModelName :: Model -> String
-webpplModelName (LR _) = "logisticRegression"
-webpplModelName (HMM _) = "hmm"
-webpplModelName (LDA _) = "lda"
+module Main (main) where
 
--- | Environment to execute benchmarks in.
-data Env = Env {rng :: GenIO, lein :: LeinProc}
+import Control.Monad.Bayes.Class (MonadInfer, MonadSample)
+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.Sampler.Strict (SamplerIO, sampleIOfixed)
+import Control.Monad.Bayes.Traced (mh)
+import Control.Monad.Bayes.Weighted (unweighted)
+import Criterion.Main
+  ( Benchmark,
+    Benchmarkable,
+    bench,
+    defaultConfig,
+    defaultMainWith,
+    nfIO,
+  )
+import Criterion.Types (Config (csvFile, rawDataFile))
+import Data.Functor (void)
+import Data.Text qualified as T
+import HMM qualified
+import LDA qualified
+import LogReg qualified
+import System.Process.Typed (runProcess)
+import System.Random.Stateful (IOGenM, StatefulGen, StdGen, mkStdGen, newIOGenM)
 
-data ProbProgSys = MonadBayes | Anglican | WebPPL
-  deriving (Show)
+data ProbProgSys = MonadBayes
+  deriving stock (Show)
 
-data Model = LR [(Double, Bool)] | HMM [Double] | LDA [[String]]
+data Model = LR [(Double, Bool)] | HMM [Double] | LDA [[T.Text]]
 
 instance Show Model where
   show (LR xs) = "LR" ++ show (length xs)
@@ -142,67 +57,43 @@
   show (RMSMC n t) = "RMSMC" ++ show n ++ "-" ++ show t
 
 runAlg :: Model -> Alg -> SamplerIO String
-runAlg model (MH n) = show <$> prior (mh n (buildModel model))
-runAlg model (SMC n) = show <$> runPopulation (smcSystematic (modelLength model) n (buildModel model))
-runAlg model (RMSMC n t) = show <$> runPopulation (rmsmcLocal (modelLength model) n t (buildModel model))
+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 (RMSMC n t) =
+  show
+    <$> population
+      ( rmsmcDynamic
+          MCMCConfig {numMCMCSteps = t, numBurnIn = 0, proposal = SingleSiteMH}
+          SMCConfig {numSteps = modelLength model, numParticles = n, resampler = resampleSystematic}
+          (buildModel model)
+      )
 
-prepareBenchmarkable :: GenIO -> ProbProgSys -> Model -> Alg -> Benchmarkable
-prepareBenchmarkable g MonadBayes model alg = nfIO $ sampleIOwith (runAlg model alg) g
-prepareBenchmarkable _ Anglican _ _ = error "Anglican benchmarks not available"
-prepareBenchmarkable _ WebPPL _ _ = error "WebPPL benchmarks not available"
+prepareBenchmarkable :: ProbProgSys -> Model -> Alg -> Benchmarkable
+prepareBenchmarkable MonadBayes model alg = nfIO $ sampleIOfixed (runAlg model alg)
 
-prepareBenchmark :: Env -> ProbProgSys -> Model -> Alg -> Benchmark
-prepareBenchmark e MonadBayes model alg =
+prepareBenchmark :: ProbProgSys -> Model -> Alg -> Benchmark
+prepareBenchmark MonadBayes model alg =
   bench (show MonadBayes ++ sep ++ show model ++ sep ++ show alg) $
-    prepareBenchmarkable (rng e) MonadBayes model alg
-  where
-    sep = "_"
-prepareBenchmark e Anglican model alg = env prepareData (const $ bench name $ whnfIO collect)
-  where
-    name = show Anglican ++ sep ++ show model ++ sep ++ show alg
-    sep = "_"
-    algString (MH n) = "-a lmh -n " ++ show n
-    algString (SMC n) = "-a smc -n " ++ show n
-    algString (RMSMC _ _) = error "Anglican does not support resample-move SMC"
-    prepareData = anglicanData (lein e) model
-    collect =
-      anglican (lein e) ["(time (m! " ++ anglicanModelName model ++ " " ++ algString alg ++ "))\n"]
-prepareBenchmark _ WebPPL model alg = bench name $ whnfIO run
+    prepareBenchmarkable MonadBayes model alg
   where
-    name = show WebPPL ++ sep ++ show model ++ sep ++ show alg
-    sep = "_"
-    algString (MH n) = "--alg MCMC --samples " ++ show n ++ " --rejuv 0 "
-    algString (SMC n) = "--alg SMC --samples " ++ show n ++ " --rejuv 0 "
-    algString (RMSMC n t) = "--alg SMC --samples " ++ show n ++ " --rejuv " ++ show t ++ " "
-    dataString (LR dataset) = let (xs, labels) = unzip dataset in "--xs='" ++ javascriptList xs ++ "' --labels='" ++ javascriptList (map (\b -> if b then (1 :: Int) else 0) labels) ++ "'"
-    dataString (HMM obs) = "--obs='" ++ javascriptList obs ++ "'"
-    dataString (LDA docs) = unwords $ zipWith (\i doc -> "--doc" ++ show (i :: Int) ++ "='" ++ unwords doc ++ "'") [1 .. 5] docs
-    run = do
-      let command = "node " ++ webpplModelName model ++ ".js " ++ algString alg ++ dataString model
-      (_, _, _, process) <- createProcess $ (shell command) {cwd = Just webpplPath, std_out = NoStream, std_err = NoStream}
-      exitCode <- waitForProcess process
-      case exitCode of
-        ExitSuccess -> return ()
-        ExitFailure i -> error $ "WebPPL terminated with exit code " ++ show i
+    sep = "_" :: String
 
 -- | Checks if the requested benchmark is implemented.
 supported :: (ProbProgSys, Model, Alg) -> Bool
-supported (Anglican, _, RMSMC _ _) = False
+supported (_, _, RMSMC _ _) = True
 supported _ = True
 
 systems :: [ProbProgSys]
 systems =
   [ MonadBayes
-    -- Anglican,
-    -- WebPPL
   ]
 
-lengthBenchmarks :: Env -> [(Double, Bool)] -> [Double] -> [[String]] -> [Benchmark]
-lengthBenchmarks e lrData hmmData ldaData = benchmarks
+lengthBenchmarks :: [(Double, Bool)] -> [Double] -> [[T.Text]] -> [Benchmark]
+lengthBenchmarks lrData hmmData ldaData = benchmarks
   where
-    lrLengths = 10 : map (* 100) [1 .. 10]
-    hmmLengths = 10 : map (* 100) [1 .. 10]
-    ldaLengths = 5 : map (* 50) [1 .. 10]
+    lrLengths = 10 : map (* 100) [1 :: Int .. 10]
+    hmmLengths = 10 : map (* 100) [1 :: Int .. 10]
+    ldaLengths = 5 : map (* 50) [1 :: Int .. 10]
     models =
       map (LR . (`take` lrData)) lrLengths
         ++ map (HMM . (`take` hmmData)) hmmLengths
@@ -212,7 +103,7 @@
         SMC 100,
         RMSMC 10 1
       ]
-    benchmarks = map (uncurry3 (prepareBenchmark e)) $ filter supported xs
+    benchmarks = map (uncurry3 (prepareBenchmark)) $ filter supported xs
       where
         uncurry3 f (x, y, z) = f x y z
         xs = do
@@ -221,12 +112,12 @@
           a <- algs
           return (s, m, a)
 
-samplesBenchmarks :: Env -> [(Double, Bool)] -> [Double] -> [[String]] -> [Benchmark]
-samplesBenchmarks e lrData hmmData ldaData = benchmarks
+samplesBenchmarks :: [(Double, Bool)] -> [Double] -> [[T.Text]] -> [Benchmark]
+samplesBenchmarks lrData hmmData ldaData = benchmarks
   where
-    lrLengths = [50]
-    hmmLengths = [20]
-    ldaLengths = [10]
+    lrLengths = [50 :: Int]
+    hmmLengths = [20 :: Int]
+    ldaLengths = [10 :: Int]
     models =
       map (LR . (`take` lrData)) lrLengths
         ++ map (HMM . (`take` hmmData)) hmmLengths
@@ -234,7 +125,7 @@
     algs =
       map (\x -> MH (100 * x)) [1 .. 10] ++ map (\x -> SMC (100 * x)) [1 .. 10]
         ++ map (\x -> RMSMC 10 (10 * x)) [1 .. 10]
-    benchmarks = map (uncurry3 (prepareBenchmark e)) $ filter supported xs
+    benchmarks = map (uncurry3 (prepareBenchmark)) $ filter supported xs
       where
         uncurry3 f (x, y, z) = f x y z
         xs = do
@@ -245,13 +136,11 @@
 
 main :: IO ()
 main = do
-  g <- createSystemRandom
-  l <- startLein
-  let e = Env g l
-  lrData <- sampleIOwith (LogReg.syntheticData 1000) g
-  hmmData <- sampleIOwith (HMM.syntheticData 1000) g
-  ldaData <- sampleIOwith (LDA.syntheticData 5 1000) g
+  lrData <- sampleIOfixed (LogReg.syntheticData 1000)
+  hmmData <- sampleIOfixed (HMM.syntheticData 1000)
+  ldaData <- sampleIOfixed (LDA.syntheticData 5 1000)
   let configLength = defaultConfig {csvFile = Just "speed-length.csv", rawDataFile = Just "raw.dat"}
-  defaultMainWith configLength (lengthBenchmarks e lrData hmmData ldaData)
+  defaultMainWith configLength (lengthBenchmarks lrData hmmData ldaData)
   let configSamples = defaultConfig {csvFile = Just "speed-samples.csv", rawDataFile = Just "raw.dat"}
-  defaultMainWith configSamples (samplesBenchmarks e lrData hmmData ldaData)
+  defaultMainWith configSamples (samplesBenchmarks lrData hmmData ldaData)
+  void $ runProcess "python plots.py"
diff --git a/models/BetaBin.hs b/models/BetaBin.hs
new file mode 100644
--- /dev/null
+++ b/models/BetaBin.hs
@@ -0,0 +1,64 @@
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}
+{-# OPTIONS_GHC -Wno-missing-export-lists #-}
+
+module BetaBin where
+
+-- The beta-binomial model in latent variable and urn model representations.
+-- The two formulations should be exactly equivalent, but only urn works with Dist.
+import Control.Monad (replicateM)
+import Control.Monad.Bayes.Class
+  ( MonadInfer,
+    MonadSample (bernoulli, uniform),
+    condition,
+  )
+import Control.Monad.State.Lazy (evalStateT, get, put)
+import Pipes ((<-<))
+import Pipes.Prelude qualified as P hiding (show)
+
+-- | Beta-binomial model as an i.i.d. sequence conditionally on weight.
+latent :: MonadSample 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 :: MonadSample m => Int -> m [Bool]
+urn n = flip evalStateT (1, 1) $ do
+  replicateM n do
+    (a, b) <- get
+    let weight = a / (a + b)
+    outcome <- bernoulli weight
+    let (a', b') = if outcome then (a + 1, b) else (a, b + 1)
+    put (a', b')
+    return outcome
+
+-- | Beta-binomial as a random process.
+-- This time using the Pipes library, for a more pure functional style
+urnP :: MonadSample m => Int -> m [Bool]
+urnP n = P.toListM $ P.take n <-< P.unfoldr toss (1, 1)
+  where
+    toss (a, b) = do
+      let weight = a / (a + b)
+      outcome <- bernoulli weight
+      let (a', b') = if outcome then (a + 1, b) else (a, b + 1)
+      return $ Right (outcome, (a', b'))
+
+-- | A beta-binomial model where the first three states are True,True,False.
+-- The resulting distribution is on the remaining outcomes.
+cond :: MonadInfer m => m [Bool] -> m [Bool]
+cond d = do
+  ~(first : second : third : rest) <- d
+  condition first
+  condition second
+  condition (not third)
+  return rest
+
+-- | The final conditional model, abstracting the representation.
+model :: MonadInfer 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.
+    count :: [Bool] -> Int
+    count = length . filter id
diff --git a/models/ConjugatePriors.hs b/models/ConjugatePriors.hs
new file mode 100644
--- /dev/null
+++ b/models/ConjugatePriors.hs
@@ -0,0 +1,66 @@
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE TypeFamilies #-}
+
+module ConjugatePriors where
+
+import Control.Applicative (Applicative (liftA2))
+import Control.Foldl (fold)
+import Control.Foldl qualified as F
+import Control.Monad.Bayes.Class (Bayesian (..), MonadInfer, MonadSample (bernoulli, beta, gamma, normal), normalPdf)
+import Numeric.Log (Log (Exp))
+import Prelude
+
+type GammaParams = (Double, Double)
+
+type BetaParams = (Double, Double)
+
+type NormalParams = (Double, Double)
+
+-- | Posterior on the precision of the normal after the points are observed
+gammaNormalAnalytic ::
+  (MonadInfer m, Foldable t, Functor t) =>
+  GammaParams ->
+  t Double ->
+  m Double
+
+-- | Exact posterior for the model.
+-- For derivation see Kevin Murphy's
+-- "Conjugate Bayesian analysis of the Gaussian distribution"
+-- section 4.
+gammaNormalAnalytic (a, b) points = gamma a' (recip b')
+  where
+    a' = a + fromIntegral (length points) / 2
+    b' = b + sum (fmap (** 2) points) / 2
+
+-- | Posterior on beta after the bernoulli sample
+betaBernoulliAnalytic :: (MonadInfer m, Foldable t) => BetaParams -> t Bool -> m Double
+betaBernoulliAnalytic (a, b) points = beta a' b'
+  where
+    (n, s) = fold (liftA2 (,) F.length (F.premap (\case True -> 1; False -> 0) F.sum)) points
+    a' = a + s
+    b' = b + fromIntegral n - s
+
+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' :: MonadInfer m => (Double, Double) -> Bayesian m Double Bool
+betaBernoulli' (a, b) = Bayesian (beta a b) bernoulli bernoulliPdf
+
+normalNormal' :: MonadInfer m => Double -> (Double, Double) -> Bayesian m Double Double
+normalNormal' var (mu0, var0) = Bayesian (normal mu0 (sqrt var0)) (`normal` (sqrt var)) (`normalPdf` (sqrt var))
+
+gammaNormal' :: MonadInfer m => (Double, Double) -> Bayesian m Double Double
+gammaNormal' (a, b) = Bayesian (gamma a (recip b)) (normal 0 . sqrt . recip) (normalPdf 0 . sqrt . recip)
+
+normalNormalAnalytic ::
+  (MonadInfer m, Foldable t) =>
+  Double ->
+  NormalParams ->
+  t Double ->
+  m Double
+normalNormalAnalytic sigma_2 (mu0, sigma0_2) points = normal mu' (sqrt sigma_2')
+  where
+    (n, s) = fold (liftA2 (,) F.length F.sum) points
+    mu' = sigma_2' * (mu0 / sigma0_2 + s / sigma_2)
+    sigma_2' = recip (recip sigma0_2 + fromIntegral n / sigma_2)
diff --git a/models/Dice.hs b/models/Dice.hs
--- a/models/Dice.hs
+++ b/models/Dice.hs
@@ -1,10 +1,15 @@
-module Dice where
+module Dice (diceHard, diceSoft) where
 
 -- A toy model for dice rolling from http://dl.acm.org/citation.cfm?id=2804317
 -- Exact results can be obtained using Dist monad
 
-import Control.Monad (liftM2)
+import Control.Applicative (liftA2)
 import Control.Monad.Bayes.Class
+  ( MonadCond (score),
+    MonadInfer,
+    MonadSample (uniformD),
+    condition,
+  )
 
 -- | A toss of a six-sided die.
 die :: MonadSample m => m Int
@@ -13,7 +18,7 @@
 -- | A sum of outcomes of n independent tosses of six-sided dice.
 dice :: MonadSample m => Int -> m Int
 dice 1 = die
-dice n = liftM2 (+) die (dice (n -1))
+dice n = liftA2 (+) die (dice (n - 1))
 
 -- | Toss of two dice where the output is greater than 4.
 diceHard :: MonadInfer m => m Int
diff --git a/models/HMM.hs b/models/HMM.hs
--- a/models/HMM.hs
+++ b/models/HMM.hs
@@ -1,20 +1,21 @@
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE TypeFamilies #-}
-
 -- HMM from Anglican (https://bitbucket.org/probprog/anglican-white-paper)
 
-module HMM
-  ( values,
-    hmm,
-    syntheticData,
-  )
-where
-
---Hidden Markov Models
+module HMM where
 
-import Control.Monad (replicateM)
+import Control.Monad (replicateM, when)
 import Control.Monad.Bayes.Class
+  ( MonadCond,
+    MonadInfer,
+    MonadSample (categorical, normal, uniformD),
+    factor,
+    normalPdf,
+  )
+import Control.Monad.Bayes.Enumerator (enumerateToDistribution)
+import Data.Maybe (fromJust, isJust)
 import Data.Vector (fromList)
+import Pipes (MFunctor (hoist), MonadTrans (lift), each, yield, (>->))
+import Pipes.Core (Producer)
+import qualified Pipes.Prelude as Pipes
 
 -- | Observed values
 values :: [Double]
@@ -70,3 +71,34 @@
 syntheticData n = replicateM n syntheticPoint
   where
     syntheticPoint = uniformD [0, 1, 2]
+
+-- | Equivalent model, but using pipes for simplicity
+
+-- | Prior expressed as a stream
+hmmPrior :: MonadSample 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 dataset = each (Nothing : (Just <$> reverse dataset))
+
+-- | Posterior expressed as a stream
+hmmPosterior :: (MonadInfer m) => [Double] -> Producer Int m ()
+hmmPosterior dataset =
+  zipWithM
+    hmmLikelihood
+    hmmPrior
+    (hmmObservations dataset)
+  where
+    hmmLikelihood :: MonadCond 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 :: MonadSample m => [Double] -> Producer Double m ()
+hmmPosteriorPredictive dataset =
+  Pipes.hoist enumerateToDistribution (hmmPosterior dataset)
+    >-> Pipes.mapM (\x -> normal (emissionMean x) 1)
diff --git a/models/LDA.hs b/models/LDA.hs
--- a/models/LDA.hs
+++ b/models/LDA.hs
@@ -1,21 +1,40 @@
+{-# LANGUAGE ImportQualifiedPost #-}
+
 -- LDA model from Anglican
 -- (https://bitbucket.org/probprog/anglican-white-paper)
 
+-- This model is just a toy/reference implementation.
+-- A more serious one would not store documents as lists of words.
+-- The point is just to showcase the model
+
 module LDA where
 
-import qualified Control.Monad as List (replicateM)
+import Control.Monad qualified as List (replicateM)
 import Control.Monad.Bayes.Class
-import qualified Data.Map as Map
-import Data.Vector as V hiding (length, mapM, mapM_, zip)
-import Numeric.Log
+  ( MonadInfer,
+    MonadSample (categorical, dirichlet, uniformD),
+    factor,
+  )
+import Control.Monad.Bayes.Sampler.Strict (sampleIO, sampleIOfixed)
+import Control.Monad.Bayes.Traced (mh)
+import Control.Monad.Bayes.Weighted (unweighted)
+import Data.Map qualified as Map
+import Data.Text (Text, words)
+import Data.Vector as V (Vector, replicate, (!))
+import Data.Vector qualified as V hiding (length, mapM, mapM_)
+import Numeric.Log (Log (Exp))
+import Text.Pretty.Simple (pPrint)
+import Prelude hiding (words)
 
-vocabulary :: [String]
+vocabulary :: [Text]
 vocabulary = ["bear", "wolf", "python", "prolog"]
 
-topics :: [String]
+topics :: [Text]
 topics = ["topic1", "topic2"]
 
-documents :: [[String]]
+type Documents = [[Text]]
+
+documents :: Documents
 documents =
   [ words "bear wolf bear wolf bear wolf python wolf bear wolf",
     words "python prolog python prolog python prolog python prolog python prolog",
@@ -24,31 +43,42 @@
     words "bear wolf bear python bear wolf bear wolf bear wolf"
   ]
 
-wordDistPrior :: MonadSample m => m (Vector Double)
+wordDistPrior :: MonadSample m => m (V.Vector Double)
 wordDistPrior = dirichlet $ V.replicate (length vocabulary) 1
 
-topicDistPrior :: MonadSample m => m (Vector Double)
+topicDistPrior :: MonadSample m => m (V.Vector Double)
 topicDistPrior = dirichlet $ V.replicate (length topics) 1
 
-wordIndex :: Map.Map String Int
+wordIndex :: Map.Map Text Int
 wordIndex = Map.fromList $ zip vocabulary [0 ..]
 
-lda :: MonadInfer m => [[String]] -> m [Int]
+lda ::
+  MonadInfer m =>
+  Documents ->
+  m (Map.Map Text (V.Vector (Text, Double)), [(Text, V.Vector (Text, Double))])
 lda docs = do
   word_dist_for_topic <- do
-    ts <- mapM (const wordDistPrior) [0 .. length topics]
-    return $ Map.fromList $ zip [0 .. length topics] ts
+    ts <- List.replicateM (length topics) wordDistPrior
+    return $ Map.fromList $ zip topics ts
   let obs doc = do
-        topic_dist <- fmap categorical topicDistPrior
+        topic_dist <- topicDistPrior
         let f word = do
-              topic <- topic_dist
+              topic <- (fmap (topics !!) . categorical) topic_dist
               factor $ (Exp . log) $ (word_dist_for_topic Map.! topic) V.! (wordIndex Map.! word)
         mapM_ f doc
-  mapM_ obs docs
-  -- return samples since Discrete is not NFData
-  mapM (categorical . snd) $ Map.toList word_dist_for_topic
+        return topic_dist
+  td <- mapM obs docs
+  return
+    ( fmap (V.zip (V.fromList vocabulary)) word_dist_for_topic,
+      zip (fmap (foldr1 (\x y -> x <> " " <> y)) docs) (fmap (V.zip $ V.fromList ["topic1", "topic2"]) td)
+    )
 
-syntheticData :: MonadSample m => Int -> Int -> m [[String]]
+syntheticData :: MonadSample m => Int -> Int -> m [[Text]]
 syntheticData d w = List.replicateM d (List.replicateM w syntheticWord)
   where
     syntheticWord = uniformD vocabulary
+
+runLDA :: IO ()
+runLDA = do
+  s <- sampleIOfixed $ unweighted $ mh 1000 $ lda documents
+  pPrint (head s)
diff --git a/models/LogReg.hs b/models/LogReg.hs
--- a/models/LogReg.hs
+++ b/models/LogReg.hs
@@ -1,35 +1,41 @@
+{-# LANGUAGE BlockArguments #-}
+
 -- Logistic regression model from Anglican
 -- (https://bitbucket.org/probprog/anglican-white-paper)
 
-module LogReg where
+module LogReg (logisticRegression, syntheticData, xs, labels) where
 
 import Control.Monad (replicateM)
 import Control.Monad.Bayes.Class
-import Numeric.Log
-
-xs :: [Double]
-xs = [-10, -5, 2, 6, 10]
-
-labels :: [Bool]
-labels = [False, False, True, True, True]
+  ( MonadInfer,
+    MonadSample (bernoulli, gamma, normal, uniform),
+    factor,
+  )
+import Numeric.Log (Log (Exp))
 
-logisticRegression :: (MonadInfer m) => [(Double, Bool)] -> m Double
+logisticRegression :: MonadInfer m => [(Double, Bool)] -> m Double
 logisticRegression dat = do
   m <- normal 0 1
   b <- normal 0 1
   sigma <- gamma 1 1
   let y x = normal (m * x + b) sigma
-      sigmoid x = y x >>= \t -> return $ 1 / (1 + exp (- t))
+      sigmoid x = y x >>= \t -> return $ 1 / (1 + exp (-t))
       obs x label = do
         p <- sigmoid x
         factor $ (Exp . log) $ if label then p else 1 - p
   mapM_ (uncurry obs) dat
   sigmoid 8
 
+-- make a synthetic dataset by randomly choosing input-label pairs
 syntheticData :: MonadSample m => Int -> m [(Double, Bool)]
-syntheticData n = replicateM n syntheticPoint
-  where
-    syntheticPoint = do
-      x <- uniform (-1) 1
-      label <- bernoulli 0.5
-      return (x, label)
+syntheticData n = replicateM n do
+  x <- uniform (-1) 1
+  label <- bernoulli 0.5
+  return (x, label)
+
+-- a tiny test dataset, for sanity-checking
+xs :: [Double]
+xs = [-10, -5, 2, 6, 10]
+
+labels :: [Bool]
+labels = [False, False, True, True, True]
diff --git a/models/NonlinearSSM.hs b/models/NonlinearSSM.hs
--- a/models/NonlinearSSM.hs
+++ b/models/NonlinearSSM.hs
@@ -1,6 +1,11 @@
 module NonlinearSSM where
 
 import Control.Monad.Bayes.Class
+  ( MonadInfer,
+    MonadSample (gamma, normal),
+    factor,
+    normalPdf,
+  )
 
 param :: MonadSample m => m (Double, Double)
 param = do
@@ -51,7 +56,7 @@
         let n = length acc
         x' <- normal (mean x n) sigmaX
         y' <- normal (sq x' / 20) sigmaY
-        simulate (k -1) x' ((x', y') : acc)
+        simulate (k - 1) x' ((x', y') : acc)
   x0 <- normal 0 (sqrt 5)
   xys <- simulate t x0 []
   return $ reverse xys
diff --git a/models/Sprinkler.hs b/models/Sprinkler.hs
--- a/models/Sprinkler.hs
+++ b/models/Sprinkler.hs
@@ -1,4 +1,4 @@
-module Sprinkler where
+module Sprinkler (hard, soft) where
 
 import Control.Monad (when)
 import Control.Monad.Bayes.Class
diff --git a/monad-bayes.cabal b/monad-bayes.cabal
--- a/monad-bayes.cabal
+++ b/monad-bayes.cabal
@@ -1,183 +1,262 @@
 cabal-version:      2.0
 name:               monad-bayes
-version:            0.1.1.0
+version:            1.0.0
 license:            MIT
 license-file:       LICENSE.md
 copyright:          2015-2020 Adam Scibior
-maintainer:         leonhard.markert@tweag.io
+maintainer:         dominic.steinitz@tweag.io
 author:             Adam Scibior <adscib@gmail.com>
 stability:          experimental
-tested-with:        ghc ==8.4.4 ghc ==8.6.5 ghc ==8.8.1
+tested-with:        GHC ==9.2.2
 homepage:           http://github.com/tweag/monad-bayes#readme
 bug-reports:        https://github.com/tweag/monad-bayes/issues
 synopsis:           A library for probabilistic programming.
 description:
-    A library for probabilistic programming using probability monads. The
-    emphasis is on composition of inference algorithms implemented in
-    terms of monad transformers.
+  A library for probabilistic programming using probability monads. The
+  emphasis is on composition of inference algorithms implemented in
+  terms of monad transformers.
 
 category:           Statistics
 build-type:         Simple
 extra-source-files: CHANGELOG.md
 
 source-repository head
-    type:     git
-    location: https://github.com/tweag/monad-bayes.git
+  type:     git
+  location: https://github.com/tweag/monad-bayes.git
 
 flag dev
-    description: Turn on development settings.
-    default:     False
-    manual:      True
+  description: Turn on development settings.
+  default:     False
+  manual:      True
 
 library
-    exposed-modules:
-        Control.Monad.Bayes.Class
-        Control.Monad.Bayes.Enumerator
-        Control.Monad.Bayes.Free
-        Control.Monad.Bayes.Helpers
-        Control.Monad.Bayes.Inference.PMMH
-        Control.Monad.Bayes.Inference.RMSMC
-        Control.Monad.Bayes.Inference.SMC
-        Control.Monad.Bayes.Inference.SMC2
-        Control.Monad.Bayes.Population
-        Control.Monad.Bayes.Sampler
-        Control.Monad.Bayes.Sequential
-        Control.Monad.Bayes.Traced
-        Control.Monad.Bayes.Traced.Basic
-        Control.Monad.Bayes.Traced.Dynamic
-        Control.Monad.Bayes.Traced.Static
-        Control.Monad.Bayes.Weighted
+  exposed-modules:
+    Control.Monad.Bayes.Class
+    Control.Monad.Bayes.Density.Free
+    Control.Monad.Bayes.Density.State
+    Control.Monad.Bayes.Enumerator
+    Control.Monad.Bayes.Inference.Lazy.MH
+    Control.Monad.Bayes.Inference.Lazy.WIS
+    Control.Monad.Bayes.Inference.MCMC
+    Control.Monad.Bayes.Inference.PMMH
+    Control.Monad.Bayes.Inference.RMSMC
+    Control.Monad.Bayes.Inference.SMC
+    Control.Monad.Bayes.Inference.SMC2
+    Control.Monad.Bayes.Inference.TUI
+    Control.Monad.Bayes.Integrator
+    Control.Monad.Bayes.Population
+    Control.Monad.Bayes.Sampler.Lazy
+    Control.Monad.Bayes.Sampler.Strict
+    Control.Monad.Bayes.Sequential.Coroutine
+    Control.Monad.Bayes.Traced
+    Control.Monad.Bayes.Traced.Basic
+    Control.Monad.Bayes.Traced.Dynamic
+    Control.Monad.Bayes.Traced.Static
+    Control.Monad.Bayes.Weighted
+    Math.Integrators.StormerVerlet
 
-    hs-source-dirs:     src
-    other-modules:      Control.Monad.Bayes.Traced.Common
-    default-language:   Haskell2010
-    default-extensions:
-        MultiParamTypeClasses RankNTypes FlexibleContexts FlexibleInstances
-        GeneralizedNewtypeDeriving TypeFamilies StandaloneDeriving GADTs
-        TupleSections
+  hs-source-dirs:     src
+  other-modules:      Control.Monad.Bayes.Traced.Common
+  default-language:   Haskell2010
+  build-depends:
+      base             >=4.11   && <4.17
+    , brick
+    , containers       >=0.5.10 && <0.7
+    , foldl
+    , free             >=5.0.2  && <5.2
+    , histogram-fill
+    , ieee754          ^>=0.8.0
+    , integration
+    , lens
+    , linear
+    , log-domain       >=0.12   && <0.14
+    , math-functions   >=0.2.1  && <0.4
+    , matrix
+    , monad-coroutine  ^>=0.9.0
+    , monad-extras
+    , mtl              ^>=2.2.2
+    , mwc-random       >=0.13.6 && <0.16
+    , pipes
+    , pretty-simple
+    , primitive
+    , random
+    , safe             ^>=0.3.17
+    , scientific
+    , statistics       >=0.14.0 && <0.17
+    , text
+    , vector           ^>=0.12.0
+    , vty
 
-    other-extensions:   ScopedTypeVariables DeriveFunctor
-    build-depends:
-        base >=4.11 && <4.14,
-        containers >=0.5.10 && <0.7,
-        free >=5.0.2 && <5.2,
-        ieee754 >=0.8.0 && <0.9,
-        log-domain >=0.12 && <0.14,
-        math-functions >=0.2.1 && <0.4,
-        monad-coroutine >=0.9.0 && <0.10,
-        mtl >=2.2.2 && <2.3,
-        mwc-random >=0.13.6 && <0.15,
-        safe >=0.3.17 && <0.4,
-        statistics >=0.14.0 && <0.16,
-        transformers >=0.5.2 && <0.6,
-        vector >=0.12.0 && <0.13
+  default-extensions:
+    BlockArguments
+    FlexibleContexts
+    ImportQualifiedPost
+    LambdaCase
+    OverloadedStrings
+    TupleSections
 
-    if flag(dev)
-        ghc-options:
-            -Wall -Wcompat -Wincomplete-record-updates
-            -Wincomplete-uni-patterns -Wnoncanonical-monad-instances
+  if flag(dev)
+    ghc-options:
+      -O2 -Wall -Wno-missing-local-signatures -Wno-trustworthy-safe
+      -Wno-missing-import-lists -Wno-implicit-prelude
+      -Wno-monomorphism-restriction
 
-    else
-        ghc-options: -Wall
+  else
+    ghc-options: -Wall
 
 executable example
-    main-is:          Single.hs
-    hs-source-dirs:   benchmark models
-    other-modules:
-        Dice
-        HMM
-        LDA
-        LogReg
+  main-is:            Single.hs
+  hs-source-dirs:     benchmark models
+  other-modules:
+    Dice
+    HMM
+    LDA
+    LogReg
 
-    default-language: Haskell2010
-    build-depends:
-        base -any,
-        containers -any,
-        log-domain -any,
-        monad-bayes -any,
-        mwc-random -any,
-        optparse-applicative -any,
-        time -any,
-        vector -any
+  default-language:   Haskell2010
+  build-depends:
+      base
+    , containers
+    , log-domain
+    , math-functions
+    , monad-bayes
+    , mwc-random
+    , optparse-applicative
+    , pipes
+    , pretty-simple
+    , random
+    , text
+    , time
+    , vector
 
-    if flag(dev)
-        ghc-options:
-            -Wall -Wcompat -Wincomplete-record-updates
-            -Wincomplete-uni-patterns -Wnoncanonical-monad-instances
+  if flag(dev)
+    ghc-options:
+      -O2 -Wall -Wcompat -Wincomplete-record-updates
+      -Wincomplete-uni-patterns -Wnoncanonical-monad-instances
 
-    else
-        ghc-options: -Wall
+  else
+    ghc-options: -Wall
 
+  default-extensions:
+    BlockArguments
+    FlexibleContexts
+    ImportQualifiedPost
+    LambdaCase
+    OverloadedStrings
+    TupleSections
+
 test-suite monad-bayes-test
-    type:             exitcode-stdio-1.0
-    main-is:          Spec.hs
-    hs-source-dirs:   test models
-    other-modules:
-        Sprinkler
-        TestEnumerator
-        TestInference
-        TestPopulation
-        TestSequential
-        TestWeighted
+  type:               exitcode-stdio-1.0
+  main-is:            Spec.hs
+  hs-source-dirs:     test models
+  other-modules:
+    BetaBin
+    ConjugatePriors
+    HMM
+    Sprinkler
+    TestAdvanced
+    TestDistribution
+    TestEnumerator
+    TestInference
+    TestIntegrator
+    TestPipes
+    TestPopulation
+    TestSampler
+    TestSequential
+    TestStormerVerlet
+    TestWeighted
 
-    default-language: Haskell2010
-    build-depends:
-        base -any,
-        QuickCheck -any,
-        hspec -any,
-        ieee754 -any,
-        log-domain -any,
-        math-functions -any,
-        monad-bayes -any,
-        mtl -any,
-        transformers -any,
-        vector -any
+  default-language:   Haskell2010
+  build-depends:
+      base
+    , containers
+    , foldl
+    , hspec
+    , ieee754
+    , lens
+    , linear
+    , log-domain
+    , math-functions
+    , matrix
+    , monad-bayes
+    , mtl
+    , mwc-random
+    , pipes
+    , pretty-simple
+    , profunctors
+    , QuickCheck
+    , random
+    , statistics
+    , text
+    , transformers
+    , vector
 
-    if flag(dev)
-        ghc-options:
-            -Wall -Wcompat -Wincomplete-record-updates
-            -Wincomplete-uni-patterns -Wnoncanonical-monad-instances
+  if flag(dev)
+    ghc-options:
+      -Wall -Wno-missing-local-signatures -Wno-unsafe
+      -Wno-missing-import-lists -Wno-implicit-prelude
 
-    else
-        ghc-options: -Wall
+  else
+    ghc-options: -Wall
 
+  default-extensions:
+    BlockArguments
+    FlexibleContexts
+    ImportQualifiedPost
+    LambdaCase
+    OverloadedStrings
+    TupleSections
+
 benchmark ssm-bench
-    type:               exitcode-stdio-1.0
-    main-is:            SSM.hs
-    hs-source-dirs:     models benchmark
-    other-modules:      NonlinearSSM
-    default-language:   Haskell2010
-    default-extensions: RankNTypes
-    build-depends:
-        base -any,
-        monad-bayes -any
+  type:             exitcode-stdio-1.0
+  main-is:          SSM.hs
+  hs-source-dirs:   models benchmark
+  other-modules:    NonlinearSSM
+  default-language: Haskell2010
+  build-depends:
+      base
+    , monad-bayes
+    , pretty-simple
+    , random
 
 benchmark speed-bench
-    type:               exitcode-stdio-1.0
-    main-is:            Speed.hs
-    hs-source-dirs:     models benchmark
-    other-modules:
-        HMM
-        LDA
-        LogReg
+  type:               exitcode-stdio-1.0
+  main-is:            Speed.hs
+  hs-source-dirs:     models benchmark
+  other-modules:
+    HMM
+    LDA
+    LogReg
 
-    default-language:   Haskell2010
-    default-extensions: RankNTypes
-    build-depends:
-        base -any,
-        abstract-par -any,
-        containers -any,
-        criterion -any,
-        log-domain -any,
-        monad-bayes -any,
-        mwc-random -any,
-        process -any,
-        vector -any
+  default-language:   Haskell2010
+  build-depends:
+      abstract-par
+    , base
+    , containers
+    , criterion
+    , log-domain
+    , monad-bayes
+    , mwc-random
+    , pipes
+    , pretty-simple
+    , process
+    , random
+    , text
+    , typed-process
+    , vector
 
-    if flag(dev)
-        ghc-options:
-            -Wall -Wcompat -Wincomplete-record-updates
-            -Wincomplete-uni-patterns -Wnoncanonical-monad-instances
+  if flag(dev)
+    ghc-options:
+      -Wall -Wno-missing-local-signatures -Wno-unsafe
+      -Wno-missing-import-lists -Wno-implicit-prelude
 
-    else
-        ghc-options: -Wall
+  else
+    ghc-options: -Wall
+
+  default-extensions:
+    BlockArguments
+    FlexibleContexts
+    ImportQualifiedPost
+    LambdaCase
+    OverloadedStrings
+    TupleSections
diff --git a/src/Control/Monad/Bayes/Class.hs b/src/Control/Monad/Bayes/Class.hs
--- a/src/Control/Monad/Bayes/Class.hs
+++ b/src/Control/Monad/Bayes/Class.hs
@@ -1,3 +1,8 @@
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE RecordWildCards #-}
+{-# OPTIONS_GHC -Wno-deprecations #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Class
 -- Description : Types for probabilistic modelling
@@ -7,8 +12,8 @@
 -- Stability   : experimental
 -- Portability : GHC
 --
--- This module defines 'MonadInfer', which can be used to represent a simple model
--- like the following:
+-- This module defines 'MonadInfer', which can be used to represent any probabilistic program,
+-- such as the following:
 --
 -- @
 -- import Control.Monad (when)
@@ -52,28 +57,52 @@
     MonadInfer,
     discrete,
     normalPdf,
+    Bayesian (..),
+    posterior,
+    priorPredictive,
+    posteriorPredictive,
+    independent,
+    mvNormal,
+    Histogram,
+    histogram,
+    histogramToList,
+    Distribution,
+    Measure,
+    Kernel,
+    Log (ln, Exp),
   )
 where
 
-import Control.Monad (when)
-import Control.Monad.Trans.Class
-import Control.Monad.Trans.Cont
-import Control.Monad.Trans.Identity
-import Control.Monad.Trans.List
-import Control.Monad.Trans.Maybe
-import Control.Monad.Trans.RWS hiding (tell)
-import Control.Monad.Trans.Reader
-import Control.Monad.Trans.State
-import Control.Monad.Trans.Writer
-import qualified Data.Vector as V
-import Data.Vector.Generic as VG
-import Numeric.Log
+import Control.Arrow (Arrow (second))
+import Control.Monad (replicateM, when)
+import Control.Monad.Cont (ContT)
+import Control.Monad.Except (ExceptT, lift)
+import Control.Monad.Identity (IdentityT)
+import Control.Monad.List (ListT)
+import Control.Monad.Reader (ReaderT)
+import Control.Monad.State (StateT)
+import Control.Monad.Writer (WriterT)
+import Data.Histogram qualified as H
+import Data.Histogram.Fill qualified as H
+import Data.Matrix
+  ( Matrix,
+    cholDecomp,
+    colVector,
+    getCol,
+    multStd,
+  )
+import Data.Vector qualified as V
+import Data.Vector.Generic as VG (Vector, map, mapM, null, sum, (!))
+import Numeric.Log (Log (..))
 import Statistics.Distribution
+  ( ContDistr (logDensity, quantile),
+    DiscreteDistr (probability),
+  )
 import Statistics.Distribution.Beta (betaDistr)
 import Statistics.Distribution.Gamma (gammaDistr)
 import Statistics.Distribution.Geometric (geometric0)
 import Statistics.Distribution.Normal (normalDistr)
-import qualified Statistics.Distribution.Poisson as Poisson
+import Statistics.Distribution.Poisson qualified as Poisson
 import Statistics.Distribution.Uniform (uniformDistr)
 
 -- | Monads that can draw random variables.
@@ -138,7 +167,7 @@
     v Double ->
     -- | outcome category
     m Int
-  categorical ps = fromPMF (ps !)
+  categorical ps = if VG.null ps then error "empty input list" else fromPMF (ps !)
 
   -- | Draw from a categorical distribution in the log domain.
   logCategorical ::
@@ -226,10 +255,21 @@
   m ()
 factor = score
 
+-- | 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. MonadSample m => m a
+
+type Measure a = forall m. MonadInfer m => m a
+
+type Kernel a b = forall m. MonadInfer m => a -> m b
+
 -- | Hard conditioning.
 condition :: MonadCond m => Bool -> m ()
 condition b = score $ if b then 1 else 0
 
+independent :: Applicative m => Int -> m a -> m [a]
+independent = replicateM
+
 -- | Monads that support both sampling and scoring.
 class (MonadSample m, MonadCond m) => MonadInfer m
 
@@ -245,6 +285,55 @@
   Log Double
 normalPdf mu sigma x = Exp $ logDensity (normalDistr mu sigma) x
 
+-- | multivariate normal
+mvNormal :: MonadSample m => V.Vector Double -> Matrix Double -> m (V.Vector Double)
+mvNormal mu bigSigma = do
+  let n = length mu
+  ss <- replicateM n (normal 0 1)
+  let bigL = cholDecomp bigSigma
+  let ts = (colVector mu) + bigL `multStd` (colVector $ V.fromList ss)
+  return $ getCol 1 ts
+
+-- | a useful datatype for expressing bayesian models
+data Bayesian m z o = Bayesian
+  { prior :: m z, -- prior over latent variable Z
+    generative :: z -> m o, -- distribution over observations given Z=z
+    likelihood :: z -> o -> Log Double -- p(o|z)
+  }
+
+posterior :: (MonadInfer m, Foldable f, Functor f) => Bayesian m z o -> f o -> m z
+posterior Bayesian {..} os = do
+  z <- prior
+  factor $ product $ fmap (likelihood z) os
+  return z
+
+priorPredictive :: Monad m => Bayesian m a b -> m b
+priorPredictive bm = prior bm >>= generative bm
+
+posteriorPredictive ::
+  (MonadInfer m, Foldable f, Functor f) =>
+  Bayesian m a b ->
+  f b ->
+  m b
+posteriorPredictive bm os = posterior bm os >>= generative bm
+
+-- helper funcs
+--------------------
+
+type Histogram = H.Histogram H.BinD Double
+
+histogram :: Int -> [(Double, Log Double)] -> Histogram
+histogram n v = H.fillBuilder buildr $ fmap (second (ln . exp)) v
+  where
+    v1 = fmap fst v
+    mi = Prelude.minimum v1
+    ma = Prelude.maximum v1
+    bins = H.binD mi n ma
+    buildr = H.mkWeighted bins
+
+histogramToList :: Histogram -> [(Double, Double)]
+histogramToList = H.asList
+
 ----------------------------------------------------------------------------
 -- Instances that lift probabilistic effects to standard tranformers.
 
@@ -257,13 +346,14 @@
 
 instance MonadInfer m => MonadInfer (IdentityT m)
 
-instance MonadSample m => MonadSample (MaybeT m) where
+instance MonadSample m => MonadSample (ExceptT e m) where
   random = lift random
+  uniformD = lift . uniformD
 
-instance MonadCond m => MonadCond (MaybeT m) where
+instance MonadCond m => MonadCond (ExceptT e m) where
   score = lift . score
 
-instance MonadInfer m => MonadInfer (MaybeT m)
+instance MonadInfer m => MonadInfer (ExceptT e m)
 
 instance MonadSample m => MonadSample (ReaderT r m) where
   random = lift random
@@ -288,19 +378,12 @@
   random = lift random
   bernoulli = lift . bernoulli
   categorical = lift . categorical
+  uniformD = lift . uniformD
 
 instance MonadCond m => MonadCond (StateT s m) where
   score = lift . score
 
 instance MonadInfer m => MonadInfer (StateT s m)
-
-instance (MonadSample m, Monoid w) => MonadSample (RWST r w s m) where
-  random = lift random
-
-instance (MonadCond m, Monoid w) => MonadCond (RWST r w s m) where
-  score = lift . score
-
-instance (MonadInfer m, Monoid w) => MonadInfer (RWST r w s m)
 
 instance MonadSample m => MonadSample (ListT m) where
   random = lift random
diff --git a/src/Control/Monad/Bayes/Density/Free.hs b/src/Control/Monad/Bayes/Density/Free.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Density/Free.hs
@@ -0,0 +1,88 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE RankNTypes #-}
+
+-- |
+-- Module      : Control.Monad.Bayes.Density.Free
+-- Description : Free monad transformer over random sampling
+-- Copyright   : (c) Adam Scibior, 2015-2020
+-- License     : MIT
+-- Maintainer  : leonhard.markert@tweag.io
+-- Stability   : experimental
+-- Portability : GHC
+--
+-- 'Density' is a free monad transformer over random sampling.
+module Control.Monad.Bayes.Density.Free
+  ( Density,
+    hoist,
+    interpret,
+    withRandomness,
+    density,
+    traced,
+  )
+where
+
+import Control.Monad.Bayes.Class (MonadSample (random))
+import Control.Monad.RWS
+import Control.Monad.State (evalStateT)
+import Control.Monad.Trans.Free.Church (FT, MonadFree (..), hoistFT, iterT, iterTM, liftF)
+import Control.Monad.Writer (WriterT (..))
+import Data.Functor.Identity (Identity, runIdentity)
+
+-- | Random sampling functor.
+newtype SamF a = Random (Double -> a) deriving (Functor)
+
+-- | 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}
+  deriving newtype (Functor, Applicative, Monad, MonadTrans)
+
+instance MonadFree SamF (Density m) where
+  wrap = Density . wrap . fmap runDensity
+
+instance Monad m => MonadSample (Density m) where
+  random = Density $ 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)
+
+-- | Execute random sampling in the transformed monad.
+interpret :: MonadSample m => Density m a -> m a
+interpret (Density 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
+  where
+    f (Random k) = do
+      xs <- get
+      case xs of
+        [] -> error "Density: 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 :: MonadSample m => [Double] -> Density m a -> m (a, [Double])
+density randomness (Density m) =
+  runWriterT $ evalStateT (iterTM f $ hoistFT lift m) randomness
+  where
+    f (Random k) = do
+      -- This block runs in StateT [Double] (WriterT [Double]) m.
+      -- StateT propagates consumed randomness while WriterT records
+      -- randomness used, whether old or new.
+      xs <- get
+      x <- case xs of
+        [] -> random
+        y : ys -> put ys >> return y
+      tell [x]
+      k x
+
+-- | Like 'density', but use an arbitrary sampling monad.
+traced :: MonadSample m => [Double] -> Density Identity a -> m (a, [Double])
+traced randomness m = density randomness $ hoist (return . runIdentity) m
diff --git a/src/Control/Monad/Bayes/Density/State.hs b/src/Control/Monad/Bayes/Density/State.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Density/State.hs
@@ -0,0 +1,40 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+
+-- |
+-- Slower than Control.Monad.Bayes.Density.Free, so not used by default,
+-- but more elementary to understand. Just uses standard
+-- monad transformer techniques.
+module Control.Monad.Bayes.Density.State where
+
+import Control.Monad.Bayes.Class (MonadSample (random))
+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)
+
+instance MonadTrans Density where
+  lift = Density . lift . lift
+
+instance Monad m => MonadState [Double] (Density m) where
+  get = Density $ lift $ get
+  put = Density . lift . put
+
+instance Monad m => MonadWriter [Double] (Density m) where
+  tell = Density . tell
+  listen = Density . listen . runDensity
+  pass = Density . pass . runDensity
+
+instance MonadSample m => MonadSample (Density m) where
+  random = do
+    trace <- get
+    x <- case trace of
+      [] -> random
+      r : xs -> put xs >> pure r
+    tell [x]
+    pure x
+
+density :: Monad m => Density m b -> [Double] -> m (b, [Double])
+density (Density m) = evalStateT (runWriterT m)
diff --git a/src/Control/Monad/Bayes/Enumerator.hs b/src/Control/Monad/Bayes/Enumerator.hs
--- a/src/Control/Monad/Bayes/Enumerator.hs
+++ b/src/Control/Monad/Bayes/Enumerator.hs
@@ -1,3 +1,7 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE ImportQualifiedPost #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Enumerator
 -- Description : Exhaustive enumeration of discrete random variables
@@ -13,28 +17,40 @@
     evidence,
     mass,
     compact,
+    enumerator,
     enumerate,
     expectation,
     normalForm,
+    toEmpirical,
+    toEmpiricalWeighted,
+    normalizeWeights,
+    enumerateToDistribution,
+    removeZeros,
+    fromList,
   )
 where
 
 import Control.Applicative (Alternative)
 import Control.Arrow (second)
-import Control.Monad (MonadPlus)
 import Control.Monad.Bayes.Class
-import Control.Monad.Trans.Writer
-import Data.AEq ((===), AEq, (~==))
-import qualified Data.Map as Map
-import Data.Maybe
-import Data.Monoid
-import qualified Data.Vector.Generic as V
-import Numeric.Log as Log
+  ( MonadCond (..),
+    MonadInfer,
+    MonadSample (bernoulli, categorical, logCategorical, random),
+  )
+import Control.Monad.Writer
+import Data.AEq (AEq, (===), (~==))
+import Data.List (sortOn)
+import Data.Map qualified as Map
+import Data.Maybe (fromMaybe)
+import Data.Ord (Down (Down))
+import Data.Vector qualified as VV
+import Data.Vector.Generic qualified as V
+import Numeric.Log as Log (Log (..), sum)
 
 -- | An exact inference transformer that integrates
 -- discrete random variables by enumerating all execution paths.
 newtype Enumerator a = Enumerator (WriterT (Product (Log Double)) [] a)
-  deriving (Functor, Applicative, Monad, Alternative, MonadPlus)
+  deriving newtype (Functor, Applicative, Monad, Alternative, MonadPlus)
 
 instance MonadSample Enumerator where
   random = error "Infinitely supported random variables not supported in Enumerator"
@@ -68,33 +84,36 @@
 mass d = f
   where
     f a = fromMaybe 0 $ lookup a m
-    m = enumerate d
+    m = enumerator d
 
 -- | Aggregate weights of equal values.
 -- The resulting list is sorted ascendingly according to values.
-compact :: (Num r, Ord a) => [(a, r)] -> [(a, r)]
-compact = Map.toAscList . Map.fromListWith (+)
+compact :: (Num r, Ord a, Ord r) => [(a, r)] -> [(a, r)]
+compact = sortOn (Down . snd) . Map.toAscList . Map.fromListWith (+)
 
 -- | Aggregate and normalize of weights.
 -- The resulting list is sorted ascendingly according to values.
 --
--- > enumerate = compact . explicit
-enumerate :: Ord a => Enumerator a -> [(a, Double)]
-enumerate d = compact (zip xs ws)
+-- > enumerator = compact . explicit
+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)
 
+-- | deprecated synonym
+enumerate = enumerator
+
 -- | Expectation of a given function computed using normalized weights.
 expectation :: (a -> Double) -> Enumerator a -> Double
 expectation f = Prelude.sum . map (\(x, w) -> f x * (exp . ln) w) . normalizeWeights . logExplicit
 
-normalize :: [Log Double] -> [Log Double]
+normalize :: Fractional b => [b] -> [b]
 normalize xs = map (/ z) xs
   where
-    z = Log.sum xs
+    z = Prelude.sum xs
 
 -- | Divide all weights by their sum.
-normalizeWeights :: [(a, Log Double)] -> [(a, Log Double)]
+normalizeWeights :: Fractional b => [(a, b)] -> [(a, b)]
 normalizeWeights ls = zip xs ps
   where
     (xs, ws) = unzip ls
@@ -103,6 +122,22 @@
 -- | 'compact' followed by removing values with zero weight.
 normalForm :: Ord a => Enumerator a -> [(a, Double)]
 normalForm = filter ((/= 0) . snd) . compact . explicit
+
+toEmpirical :: (Fractional b, Ord a, Ord b) => [a] -> [(a, b)]
+toEmpirical ls = normalizeWeights $ compact (zip ls (repeat 1))
+
+toEmpiricalWeighted :: (Fractional b, Ord a, Ord b) => [(a, b)] -> [(a, b)]
+toEmpiricalWeighted = normalizeWeights . compact
+
+enumerateToDistribution :: (MonadSample n) => Enumerator a -> n a
+enumerateToDistribution model = do
+  let samples = logExplicit model
+  let (support, logprobs) = unzip samples
+  i <- logCategorical $ VV.fromList logprobs
+  return $ support !! i
+
+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
   p == q = normalForm p == normalForm q
diff --git a/src/Control/Monad/Bayes/Free.hs b/src/Control/Monad/Bayes/Free.hs
deleted file mode 100644
--- a/src/Control/Monad/Bayes/Free.hs
+++ /dev/null
@@ -1,86 +0,0 @@
--- |
--- Module      : Control.Monad.Bayes.Free
--- Description : Free monad transformer over random sampling
--- Copyright   : (c) Adam Scibior, 2015-2020
--- License     : MIT
--- Maintainer  : leonhard.markert@tweag.io
--- Stability   : experimental
--- Portability : GHC
---
--- 'FreeSampler' is a free monad transformer over random sampling.
-module Control.Monad.Bayes.Free
-  ( FreeSampler,
-    hoist,
-    interpret,
-    withRandomness,
-    withPartialRandomness,
-    runWith,
-  )
-where
-
-import Control.Monad.Bayes.Class
-import Control.Monad.State (evalStateT, get, put)
-import Control.Monad.Trans (MonadTrans (..))
-import Control.Monad.Trans.Free.Church (FT, MonadFree (..), hoistFT, iterT, iterTM, liftF)
-import Control.Monad.Writer (WriterT (..), tell)
-import Data.Functor.Identity (Identity, runIdentity)
-
--- | Random sampling functor.
-newtype SamF a = Random (Double -> a)
-
-instance Functor SamF where
-  fmap f (Random k) = Random (f . k)
-
--- | Free monad transformer over random sampling.
---
--- Uses the Church-encoded version of the free monad for efficiency.
-newtype FreeSampler m a = FreeSampler {runFreeSampler :: FT SamF m a}
-  deriving (Functor, Applicative, Monad, MonadTrans)
-
-instance Monad m => MonadFree SamF (FreeSampler m) where
-  wrap = FreeSampler . wrap . fmap runFreeSampler
-
-instance Monad m => MonadSample (FreeSampler m) where
-  random = FreeSampler $ liftF (Random id)
-
--- | Hoist 'FreeSampler' through a monad transform.
-hoist :: (Monad m, Monad n) => (forall x. m x -> n x) -> FreeSampler m a -> FreeSampler n a
-hoist f (FreeSampler m) = FreeSampler (hoistFT f m)
-
--- | Execute random sampling in the transformed monad.
-interpret :: MonadSample m => FreeSampler m a -> m a
-interpret (FreeSampler m) = iterT f m
-  where
-    f (Random k) = random >>= k
-
--- | Execute computation with supplied values for random choices.
-withRandomness :: Monad m => [Double] -> FreeSampler m a -> m a
-withRandomness randomness (FreeSampler m) = evalStateT (iterTM f m) randomness
-  where
-    f (Random k) = do
-      xs <- get
-      case xs of
-        [] -> error "FreeSampler: 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.
-withPartialRandomness :: MonadSample m => [Double] -> FreeSampler m a -> m (a, [Double])
-withPartialRandomness randomness (FreeSampler m) =
-  runWriterT $ evalStateT (iterTM f $ hoistFT lift m) randomness
-  where
-    f (Random k) = do
-      -- This block runs in StateT [Double] (WriterT [Double]) m.
-      -- StateT propagates consumed randomness while WriterT records
-      -- randomness used, whether old or new.
-      xs <- get
-      x <- case xs of
-        [] -> random
-        y : ys -> put ys >> return y
-      tell [x]
-      k x
-
--- | Like 'withPartialRandomness', but use an arbitrary sampling monad.
-runWith :: MonadSample m => [Double] -> FreeSampler Identity a -> m (a, [Double])
-runWith randomness m = withPartialRandomness randomness $ hoist (return . runIdentity) m
diff --git a/src/Control/Monad/Bayes/Helpers.hs b/src/Control/Monad/Bayes/Helpers.hs
deleted file mode 100644
--- a/src/Control/Monad/Bayes/Helpers.hs
+++ /dev/null
@@ -1,77 +0,0 @@
--- |
--- Module      : Control.Monad.Bayes.Helpers
--- Description : Helper functions for working with inference monads
--- Copyright   : (c) Adam Scibior, 2015-2020
--- License     : MIT
--- Maintainer  : leonhard.markert@tweag.io
--- Stability   : experimental
--- Portability : GHC
-module Control.Monad.Bayes.Helpers
-  ( W,
-    hoistW,
-    P,
-    hoistP,
-    S,
-    hoistS,
-    F,
-    hoistF,
-    T,
-    hoistT,
-    hoistWF,
-    hoistSP,
-    hoistSTP,
-  )
-where
-
-import Control.Monad.Bayes.Free as Free
-import Control.Monad.Bayes.Population as Pop
-import Control.Monad.Bayes.Sequential as Seq
-import Control.Monad.Bayes.Traced as Tr
-import Control.Monad.Bayes.Weighted as Weighted
-
-type W = Weighted
-
-type P = Population
-
-type S = Sequential
-
-type F = FreeSampler
-
-type T = Traced
-
-hoistW :: (forall x. m x -> n x) -> W m a -> W n a
-hoistW = Weighted.hoist
-
-hoistP ::
-  (Monad m, Monad n) =>
-  (forall x. m x -> n x) ->
-  P m a ->
-  P n a
-hoistP = Pop.hoist
-
-hoistS :: (forall x. m x -> m x) -> S m a -> S m a
-hoistS = Seq.hoistFirst
-
-hoistF :: (Monad m, Monad n) => (forall x. m x -> n x) -> F m a -> F n a
-hoistF = Free.hoist
-
-hoistWF ::
-  (Monad m, Monad n) =>
-  (forall x. m x -> n x) ->
-  W (F m) a ->
-  W (F n) a
-hoistWF m = hoistW $ hoistF m
-
-hoistSP ::
-  Monad m =>
-  (forall x. m x -> m x) ->
-  S (P m) a ->
-  S (P m) a
-hoistSP m = hoistS $ hoistP m
-
-hoistSTP ::
-  Monad m =>
-  (forall x. m x -> m x) ->
-  S (T (P m)) a ->
-  S (T (P m)) a
-hoistSTP m = hoistS $ hoistT $ hoistP m
diff --git a/src/Control/Monad/Bayes/Inference/Lazy/MH.hs b/src/Control/Monad/Bayes/Inference/Lazy/MH.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/Lazy/MH.hs
@@ -0,0 +1,79 @@
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# OPTIONS_GHC -Wno-name-shadowing #-}
+
+module Control.Monad.Bayes.Inference.Lazy.MH where
+
+import Control.Monad.Bayes.Class (Log (ln))
+import Control.Monad.Bayes.Sampler.Lazy
+  ( Sampler (runSampler),
+    Tree (..),
+    Trees (..),
+    randomTree,
+  )
+import Control.Monad.Bayes.Weighted (Weighted, weighted)
+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 p m = do
+  -- Top level: produce a stream of samples.
+  -- Split the random number generator in two
+  -- One part is used as the first seed for the simulation,
+  -- and one part is used for the randomness in the MH algorithm.
+  g <- newStdGen >> getStdGen
+  let (g1, g2) = split g
+  let t = randomTree g1
+  let (x, w) = runSampler (weighted 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.
+  return $ map (\(_, x, w) -> (x, w)) samples
+  where
+    --   where
+    {- NB There are three kinds of randomness in the step function.
+    1. The start tree 't', which is the source of randomness for simulating the
+    program m to start with. This is sort-of the point in the "state space".
+    2. The randomness needed to propose a new tree ('g1')
+    3. The randomness needed to decide whether to accept or reject that ('g2')
+    The tree t is an argument and result,
+    but we use a state monad ('get'/'put') to deal with the other randomness '(g,g1,g2)' -}
+
+    -- step :: RandomGen g => (Tree, a, Log Double) -> State g (Tree, a, Log Double)
+    step (t, x, w) = do
+      -- Randomly change some sites
+      g <- get
+      let (g1, g2) = split g
+      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'
+      -- MH acceptance ratio. This is the probability of either
+      -- returning the new seed or the old one.
+      let ratio = w' / w
+      let (r, g2') = R.random g2
+      put g2'
+      if r < min 1 (exp $ ln ratio)
+        then return (t', x', w')
+        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 p g (Tree a ts) =
+  let (a', g') = (R.random g :: (Double, g))
+      (a'', g'') = R.random g'
+   in Tree
+        { currentUniform = if a' < p then a'' else a,
+          lazyUniforms = mutateTrees p g'' ts
+        }
+
+mutateTrees :: RandomGen g => Double -> g -> Trees -> Trees
+mutateTrees p g (Trees t ts) =
+  let (g1, g2) = split g
+   in Trees
+        { headTree = mutateTree p g1 t,
+          tailTrees = mutateTrees p g2 ts
+        }
diff --git a/src/Control/Monad/Bayes/Inference/Lazy/WIS.hs b/src/Control/Monad/Bayes/Inference/Lazy/WIS.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/Lazy/WIS.hs
@@ -0,0 +1,23 @@
+module Control.Monad.Bayes.Inference.Lazy.WIS where
+
+import Control.Monad.Bayes.Sampler.Lazy (Sampler, weightedsamples)
+import Control.Monad.Bayes.Weighted (Weighted)
+import Numeric.Log (Log (Exp))
+import System.Random (Random (randoms), getStdGen, newStdGen)
+
+-- | Weighted 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 n m = do
+  xws <- weightedsamples m
+  let xws' = take n $ accumulate xws 0
+  let max' = snd $ last xws'
+  _ <- newStdGen
+  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 ((x, w) : xws) a = (x, w + a) : (x, w + a) : accumulate xws (w + a)
+    accumulate [] _ = []
diff --git a/src/Control/Monad/Bayes/Inference/MCMC.hs b/src/Control/Monad/Bayes/Inference/MCMC.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/MCMC.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE RecordWildCards #-}
+
+-- |
+-- Module      : Control.Monad.Bayes.Inference.MCMC
+-- Description : Markov Chain Monte Carlo (MCMC)
+-- Copyright   : (c) Adam Scibior, 2015-2020
+-- License     : MIT
+-- Maintainer  : tweag.io
+-- Stability   : experimental
+-- Portability : GHC
+module Control.Monad.Bayes.Inference.MCMC where
+
+import Control.Monad.Bayes.Class
+import qualified Control.Monad.Bayes.Traced.Basic as Basic
+import Control.Monad.Bayes.Traced.Common
+import qualified Control.Monad.Bayes.Traced.Dynamic as Dynamic
+import qualified Control.Monad.Bayes.Traced.Static as Static
+import Control.Monad.Bayes.Weighted
+import Pipes ((>->))
+import qualified Pipes as P
+import qualified Pipes.Prelude as P
+
+data Proposal = SingleSiteMH
+
+data MCMCConfig = MCMCConfig {proposal :: Proposal, numMCMCSteps :: Int, numBurnIn :: Int}
+
+defaultMCMCConfig :: MCMCConfig
+defaultMCMCConfig = MCMCConfig {proposal = SingleSiteMH, numMCMCSteps = 1, numBurnIn = 0}
+
+mcmc :: MonadSample m => MCMCConfig -> Static.Traced (Weighted m) a -> m [a]
+mcmc (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Static.mh numMCMCSteps m
+
+mcmcBasic :: MonadSample m => MCMCConfig -> Basic.Traced (Weighted m) a -> m [a]
+mcmcBasic (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Basic.mh numMCMCSteps m
+
+mcmcDynamic :: MonadSample m => MCMCConfig -> Dynamic.Traced (Weighted 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) =
+  P.repeatM d
+    >-> P.takeWhile' ((== 0) . probDensity)
+    >-> P.map (MHResult False)
+
+-- | convert a probabilistic program into a producer of samples
+mcmcP :: MonadSample m => MCMCConfig -> Static.Traced m a -> P.Producer (MHResult a) m ()
+mcmcP MCMCConfig {..} m@(Static.Traced w _) = do
+  initialValue <- independentSamples m >-> P.drain
+  ( P.unfoldr (fmap (Right . (\k -> (k, trace k))) . mhTransWithBool w) initialValue
+      >-> P.drop numBurnIn
+    )
diff --git a/src/Control/Monad/Bayes/Inference/PMMH.hs b/src/Control/Monad/Bayes/Inference/PMMH.hs
--- a/src/Control/Monad/Bayes/Inference/PMMH.hs
+++ b/src/Control/Monad/Bayes/Inference/PMMH.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE RankNTypes #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Inference.PMMH
 -- Description : Particle Marginal Metropolis-Hastings (PMMH)
@@ -12,30 +14,49 @@
 -- Christophe Andrieu, Arnaud Doucet, and Roman Holenstein. 2010. Particle Markov chain Monte Carlo Methods. /Journal of the Royal Statistical Society/ 72 (2010), 269-342. <http://www.stats.ox.ac.uk/~doucet/andrieu_doucet_holenstein_PMCMC.pdf>
 module Control.Monad.Bayes.Inference.PMMH
   ( pmmh,
+    pmmhBayesianModel,
   )
 where
 
-import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Inference.SMC
+import Control.Monad.Bayes.Class (Bayesian (generative), MonadInfer, MonadSample, prior)
+import Control.Monad.Bayes.Inference.MCMC (MCMCConfig, mcmc)
+import Control.Monad.Bayes.Inference.SMC (SMCConfig (), smc)
 import Control.Monad.Bayes.Population as Pop
-import Control.Monad.Bayes.Sequential
-import Control.Monad.Bayes.Traced
+  ( Population,
+    hoist,
+    population,
+    pushEvidence,
+  )
+import Control.Monad.Bayes.Sequential.Coroutine (Sequential)
+import Control.Monad.Bayes.Traced.Static (Traced)
+import Control.Monad.Bayes.Weighted
 import Control.Monad.Trans (lift)
-import Numeric.Log
+import Numeric.Log (Log)
 
 -- | Particle Marginal Metropolis-Hastings sampling.
 pmmh ::
+  MonadSample m =>
+  MCMCConfig ->
+  SMCConfig (Weighted m) ->
+  Traced (Weighted m) a1 ->
+  (a1 -> Sequential (Population (Weighted m)) a2) ->
+  m [[(a2, Log Double)]]
+pmmh mcmcConf smcConf param model =
+  mcmc
+    mcmcConf
+    ( param
+        >>= population
+          . pushEvidence
+          . Pop.hoist lift
+          . smc smcConf
+          . model
+    )
+
+-- | Particle Marginal Metropolis-Hastings sampling from a Bayesian model
+pmmhBayesianModel ::
   MonadInfer m =>
-  -- | number of Metropolis-Hastings steps
-  Int ->
-  -- | number of time steps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | model parameters prior
-  Traced m b ->
-  -- | model
-  (b -> Sequential (Population m) a) ->
-  m [[(a, Log Double)]]
-pmmh t k n param model =
-  mh t (param >>= runPopulation . pushEvidence . Pop.hoist lift . smcSystematic k n . model)
+  MCMCConfig ->
+  SMCConfig (Weighted m) ->
+  (forall m'. MonadInfer m' => Bayesian m' a1 a2) ->
+  m [[(a2, Log Double)]]
+pmmhBayesianModel mcmcConf smcConf bm = pmmh mcmcConf smcConf (prior bm) (generative bm)
diff --git a/src/Control/Monad/Bayes/Inference/RMSMC.hs b/src/Control/Monad/Bayes/Inference/RMSMC.hs
--- a/src/Control/Monad/Bayes/Inference/RMSMC.hs
+++ b/src/Control/Monad/Bayes/Inference/RMSMC.hs
@@ -1,3 +1,6 @@
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# LANGUAGE RecordWildCards #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Inference.RMSMC
 -- Description : Resample-Move Sequential Monte Carlo (RM-SMC)
@@ -12,73 +15,76 @@
 -- Walter Gilks and Carlo Berzuini. 2001. Following a moving target - Monte Carlo inference for dynamic Bayesian models. /Journal of the Royal Statistical Society/ 63 (2001), 127-146. <http://www.mathcs.emory.edu/~whalen/Papers/BNs/MonteCarlo-DBNs.pdf>
 module Control.Monad.Bayes.Inference.RMSMC
   ( rmsmc,
-    rmsmcLocal,
+    rmsmcDynamic,
     rmsmcBasic,
   )
 where
 
-import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Helpers
+import Control.Monad.Bayes.Class (MonadSample)
+import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (..))
+import Control.Monad.Bayes.Inference.SMC
 import Control.Monad.Bayes.Population
-import Control.Monad.Bayes.Sequential as Seq
-import Control.Monad.Bayes.Traced as Tr
-import qualified Control.Monad.Bayes.Traced.Basic as TrBas
-import qualified Control.Monad.Bayes.Traced.Dynamic as TrDyn
+  ( Population,
+    spawn,
+    withParticles,
+  )
+import Control.Monad.Bayes.Sequential.Coroutine as Seq
+import Control.Monad.Bayes.Sequential.Coroutine qualified as S
+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,
+    marginal,
+    mhStep,
+  )
+import Control.Monad.Bayes.Traced.Static qualified as TrStat
+import Data.Monoid (Endo (..))
 
 -- | Resample-move Sequential Monte Carlo.
 rmsmc ::
   MonadSample m =>
-  -- | number of timesteps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | number of Metropolis-Hastings transitions after each resampling
-  Int ->
+  MCMCConfig ->
+  SMCConfig m ->
   -- | model
   Sequential (Traced (Population m)) a ->
   Population m a
-rmsmc k n t =
+rmsmc (MCMCConfig {..}) (SMCConfig {..}) =
   marginal
-    . sis (composeCopies t mhStep . hoistT resampleSystematic) k
-    . hoistS (hoistT (spawn n >>))
+    . S.sequentially (composeCopies numMCMCSteps mhStep . TrStat.hoist resampler) numSteps
+    . S.hoistFirst (TrStat.hoist (spawn numParticles >>))
 
 -- | Resample-move Sequential Monte Carlo with a more efficient
 -- tracing representation.
 rmsmcBasic ::
   MonadSample m =>
-  -- | number of timesteps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | number of Metropolis-Hastings transitions after each resampling
-  Int ->
+  MCMCConfig ->
+  SMCConfig m ->
   -- | model
   Sequential (TrBas.Traced (Population m)) a ->
   Population m a
-rmsmcBasic k n t =
+rmsmcBasic (MCMCConfig {..}) (SMCConfig {..}) =
   TrBas.marginal
-    . sis (composeCopies t TrBas.mhStep . TrBas.hoistT resampleSystematic) k
-    . hoistS (TrBas.hoistT (spawn n >>))
+    . S.sequentially (composeCopies numMCMCSteps TrBas.mhStep . TrBas.hoist resampler) numSteps
+    . S.hoistFirst (TrBas.hoist (withParticles numParticles))
 
 -- | A variant of resample-move Sequential Monte Carlo
 -- where only random variables since last resampling are considered
 -- for rejuvenation.
-rmsmcLocal ::
+rmsmcDynamic ::
   MonadSample m =>
-  -- | number of timesteps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | number of Metropolis-Hastings transitions after each resampling
-  Int ->
+  MCMCConfig ->
+  SMCConfig m ->
   -- | model
   Sequential (TrDyn.Traced (Population m)) a ->
   Population m a
-rmsmcLocal k n t =
+rmsmcDynamic (MCMCConfig {..}) (SMCConfig {..}) =
   TrDyn.marginal
-    . sis (TrDyn.freeze . composeCopies t TrDyn.mhStep . TrDyn.hoistT resampleSystematic) k
-    . hoistS (TrDyn.hoistT (spawn n >>))
+    . S.sequentially (TrDyn.freeze . composeCopies numMCMCSteps TrDyn.mhStep . TrDyn.hoist resampler) numSteps
+    . S.hoistFirst (TrDyn.hoist (withParticles numParticles))
 
 -- | Apply a function a given number of times.
 composeCopies :: Int -> (a -> a) -> (a -> a)
-composeCopies k f = foldr (.) id (replicate k f)
+composeCopies k = withEndo (mconcat . replicate k)
+
+withEndo :: (Endo a -> Endo b) -> (a -> a) -> b -> b
+withEndo f = appEndo . f . Endo
diff --git a/src/Control/Monad/Bayes/Inference/SMC.hs b/src/Control/Monad/Bayes/Inference/SMC.hs
--- a/src/Control/Monad/Bayes/Inference/SMC.hs
+++ b/src/Control/Monad/Bayes/Inference/SMC.hs
@@ -1,3 +1,6 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE RecordWildCards #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Inference.SMC
 -- Description : Sequential Monte Carlo (SMC)
@@ -11,83 +14,40 @@
 --
 -- Arnaud Doucet and Adam M. Johansen. 2011. A tutorial on particle filtering and smoothing: fifteen years later. In /The Oxford Handbook of Nonlinear Filtering/, Dan Crisan and Boris Rozovskii (Eds.). Oxford University Press, Chapter 8.
 module Control.Monad.Bayes.Inference.SMC
-  ( sir,
-    smcMultinomial,
-    smcSystematic,
-    smcMultinomialPush,
-    smcSystematicPush,
+  ( smc,
+    smcPush,
+    SMCConfig (..),
   )
 where
 
-import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Class (MonadInfer, MonadSample)
 import Control.Monad.Bayes.Population
-import Control.Monad.Bayes.Sequential as Seq
+  ( Population,
+    pushEvidence,
+    withParticles,
+  )
+import Control.Monad.Bayes.Sequential.Coroutine as Coroutine
 
+data SMCConfig m = SMCConfig
+  { resampler :: forall x. Population m x -> Population m x,
+    numSteps :: Int,
+    numParticles :: Int
+  }
+
 -- | Sequential importance resampling.
 -- Basically an SMC template that takes a custom resampler.
-sir ::
-  Monad m =>
-  -- | resampler
-  (forall x. Population m x -> Population m x) ->
-  -- | number of timesteps
-  Int ->
-  -- | population size
-  Int ->
-  -- | model
-  Sequential (Population m) a ->
-  Population m a
-sir resampler k n = sis resampler k . Seq.hoistFirst (spawn n >>)
-
--- | Sequential Monte Carlo with multinomial resampling at each timestep.
--- Weights are not normalized.
-smcMultinomial ::
-  MonadSample m =>
-  -- | number of timesteps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | model
-  Sequential (Population m) a ->
-  Population m a
-smcMultinomial = sir resampleMultinomial
-
--- | Sequential Monte Carlo with systematic resampling at each timestep.
--- Weights are not normalized.
-smcSystematic ::
+smc ::
   MonadSample m =>
-  -- | number of timesteps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | model
-  Sequential (Population m) a ->
+  SMCConfig m ->
+  Coroutine.Sequential (Population m) a ->
   Population m a
-smcSystematic = sir resampleSystematic
+smc SMCConfig {..} =
+  Coroutine.sequentially resampler numSteps
+    . Coroutine.hoistFirst (withParticles numParticles)
 
 -- | Sequential Monte Carlo with multinomial resampling at each timestep.
 -- Weights are normalized at each timestep and the total weight is pushed
 -- as a score into the transformed monad.
-smcMultinomialPush ::
-  MonadInfer m =>
-  -- | number of timesteps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | model
-  Sequential (Population m) a ->
-  Population m a
-smcMultinomialPush = sir (pushEvidence . resampleMultinomial)
-
--- | Sequential Monte Carlo with systematic resampling at each timestep.
--- Weights are normalized at each timestep and the total weight is pushed
--- as a score into the transformed monad.
-smcSystematicPush ::
-  MonadInfer m =>
-  -- | number of timesteps
-  Int ->
-  -- | number of particles
-  Int ->
-  -- | model
-  Sequential (Population m) a ->
-  Population m a
-smcSystematicPush = sir (pushEvidence . resampleSystematic)
+smcPush ::
+  MonadInfer m => SMCConfig m -> Coroutine.Sequential (Population m) a -> Population m a
+smcPush config = smc config {resampler = (pushEvidence . resampler config)}
diff --git a/src/Control/Monad/Bayes/Inference/SMC2.hs b/src/Control/Monad/Bayes/Inference/SMC2.hs
--- a/src/Control/Monad/Bayes/Inference/SMC2.hs
+++ b/src/Control/Monad/Bayes/Inference/SMC2.hs
@@ -1,3 +1,6 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Inference.SMC2
 -- Description : Sequential Monte Carlo squared (SMC²)
@@ -12,22 +15,29 @@
 -- Nicolas Chopin, Pierre E. Jacob, and Omiros Papaspiliopoulos. 2013. SMC²: an efficient algorithm for sequential analysis of state space models. /Journal of the Royal Statistical Society Series B: Statistical Methodology/ 75 (2013), 397-426. Issue 3. <https://doi.org/10.1111/j.1467-9868.2012.01046.x>
 module Control.Monad.Bayes.Inference.SMC2
   ( smc2,
+    SMC2,
   )
 where
 
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Helpers
-import Control.Monad.Bayes.Inference.RMSMC
-import Control.Monad.Bayes.Inference.SMC
-import Control.Monad.Bayes.Population as Pop
-import Control.Monad.Trans
-import Numeric.Log
+  ( MonadCond (..),
+    MonadInfer,
+    MonadSample (random),
+  )
+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.Traced
+import Control.Monad.Trans (MonadTrans (..))
+import Numeric.Log (Log)
 
 -- | Helper monad transformer for preprocessing the model for 'smc2'.
-newtype SMC2 m a = SMC2 (S (T (P m)) a)
-  deriving (Functor, Applicative, Monad)
+newtype SMC2 m a = SMC2 (Sequential (Traced (Population m)) a)
+  deriving newtype (Functor, Applicative, Monad)
 
-setup :: SMC2 m a -> S (T (P m)) a
+setup :: SMC2 m a -> Sequential (Traced (Population m)) a
 setup (SMC2 m) = m
 
 instance MonadTrans SMC2 where
@@ -53,9 +63,12 @@
   -- | number of MH transitions
   Int ->
   -- | model parameters
-  S (T (P m)) b ->
+  Sequential (Traced (Population m)) b ->
   -- | model
-  (b -> S (P (SMC2 m)) a) ->
-  P m [(a, Log Double)]
-smc2 k n p t param model =
-  rmsmc k p t (param >>= setup . runPopulation . smcSystematicPush k n . model)
+  (b -> Sequential (Population (SMC2 m)) a) ->
+  Population 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)
diff --git a/src/Control/Monad/Bayes/Inference/TUI.hs b/src/Control/Monad/Bayes/Inference/TUI.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/TUI.hs
@@ -0,0 +1,176 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# OPTIONS_GHC -Wno-type-defaults #-}
+
+module Control.Monad.Bayes.Inference.TUI where
+
+import Brick
+import Brick qualified as B
+import Brick.BChan qualified as B
+import Brick.Widgets.Border
+import Brick.Widgets.Border.Style
+import Brick.Widgets.Center
+import Brick.Widgets.ProgressBar qualified as B
+import Control.Arrow (Arrow (..))
+import Control.Concurrent (forkIO)
+import Control.Foldl qualified as Fold
+import Control.Monad (void)
+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.Common
+import Control.Monad.Bayes.Weighted
+import Data.Scientific (FPFormat (Exponent), formatScientific, fromFloatDigits)
+import Data.Text qualified as T
+import Data.Text.Lazy qualified as TL
+import Data.Text.Lazy.IO qualified as TL
+import GHC.Float (double2Float)
+import Graphics.Vty
+import Graphics.Vty qualified as V
+import Numeric.Log (Log (ln))
+import Pipes (runEffect, (>->))
+import Pipes qualified as P
+import Pipes.Prelude qualified as P
+import Text.Pretty.Simple (pShow, pShowNoColor)
+
+data MCMCData a = MCMCData
+  { numSteps :: Int,
+    numSuccesses :: Int,
+    samples :: [a],
+    lk :: [Double],
+    totalSteps :: Int
+  }
+  deriving stock (Show)
+
+-- | Brick is a terminal user interface (TUI)
+-- which we use to display inference algorithms in progress
+
+-- | draw the brick app
+drawUI :: ([a] -> Widget n) -> MCMCData a -> [Widget n]
+drawUI handleSamples state = [ui]
+  where
+    completionBar =
+      updateAttrMap
+        ( B.mapAttrNames
+            [ (doneAttr, B.progressCompleteAttr),
+              (toDoAttr, B.progressIncompleteAttr)
+            ]
+        )
+        $ toBar $ fromIntegral $ numSteps state
+
+    likelihoodBar =
+      updateAttrMap
+        ( B.mapAttrNames
+            [ (doneAttr, B.progressCompleteAttr),
+              (toDoAttr, B.progressIncompleteAttr)
+            ]
+        )
+        $ B.progressBar
+          (Just $ "Mean likelihood for last 1000 samples: " <> take 10 (show (head $ lk state <> [0])))
+          (double2Float (Fold.fold Fold.mean $ take 1000 $ lk state) / double2Float (maximum $ 0 : lk state))
+
+    displayStep c = Just $ "Step " <> show c
+    numFailures = numSteps state - numSuccesses state
+    toBar v = B.progressBar (displayStep v) (v / fromIntegral (totalSteps state))
+    displaySuccessesAndFailures =
+      withBorderStyle unicode $
+        borderWithLabel (str "Successes and failures") $
+          center (str (show $ numSuccesses state))
+            <+> vBorder
+            <+> center (str (show numFailures))
+    warning =
+      if numSteps state > 1000 && (fromIntegral (numSuccesses state) / fromIntegral (numSteps state)) < 0.1
+        then withAttr (attrName "highlight") $ str "Warning: acceptance rate is rather low.\nThis probably means that your proposal isn't good."
+        else str ""
+
+    ui =
+      (str "Progress: " <+> completionBar)
+        <=> (str "Likelihood: " <+> likelihoodBar)
+        <=> str "\n"
+        <=> displaySuccessesAndFailures
+        <=> warning
+        <=> handleSamples (samples state)
+
+noVisual :: b -> Widget n
+noVisual = const emptyWidget
+
+showEmpirical :: (Show a, Ord a) => [a] -> Widget n
+showEmpirical =
+  txt
+    . T.pack
+    . TL.unpack
+    . pShow
+    . (fmap (second (formatScientific Exponent (Just 3) . fromFloatDigits)))
+    . toEmpirical
+
+showVal :: Show a => [a] -> Widget n
+showVal = txt . T.pack . (\case [] -> ""; a -> show $ head a)
+
+-- | handler for events received by the TUI
+appEvent :: s -> B.BrickEvent n1 s -> B.EventM n2 (B.Next s)
+appEvent p (B.VtyEvent e) =
+  case e of
+    V.EvKey (V.KChar 'q') [] -> do
+      B.halt p
+    _ -> B.continue p
+appEvent _ (B.AppEvent d) = B.continue d
+appEvent _ _ = error "unknown event"
+
+doneAttr, toDoAttr :: B.AttrName
+doneAttr = B.attrName "theBase" <> B.attrName "done"
+toDoAttr = B.attrName "theBase" <> B.attrName "remaining"
+
+theMap :: B.AttrMap
+theMap =
+  B.attrMap
+    V.defAttr
+    [ (B.attrName "theBase", bg V.brightBlack),
+      (doneAttr, V.black `on` V.white),
+      (toDoAttr, V.white `on` V.black),
+      (attrName "highlight", fg yellow)
+    ]
+
+tui :: Show a => Int -> Traced (Weighted SamplerIO) a -> ([a] -> Widget ()) -> IO ()
+tui burnIn distribution visualizer = void do
+  eventChan <- B.newBChan 10
+  initialVty <- buildVty
+  _ <- forkIO $ run (mcmcP MCMCConfig {numBurnIn = burnIn, proposal = SingleSiteMH, numMCMCSteps = -1} distribution) eventChan n
+  samples <-
+    B.customMain
+      initialVty
+      buildVty
+      (Just eventChan)
+      ( ( B.App
+            { B.appDraw = drawUI visualizer,
+              B.appChooseCursor = B.showFirstCursor,
+              B.appHandleEvent = appEvent,
+              B.appStartEvent = return,
+              B.appAttrMap = const theMap
+            }
+        )
+      )
+      (initialState n)
+  TL.writeFile "data/tui_output.txt" (pShowNoColor samples)
+  return samples
+  where
+    buildVty = V.mkVty V.defaultConfig
+    n = 100000
+    initialState n = MCMCData {numSteps = 0, samples = [], lk = [], numSuccesses = 0, totalSteps = n}
+
+    run prod chan i =
+      runEffect $
+        P.hoist (sampleIO . unweighted) prod
+          >-> P.scan
+            ( \mcmcdata@(MCMCData ns nsc smples lk _) a ->
+                mcmcdata
+                  { numSteps = ns + 1,
+                    numSuccesses = nsc + if success a then 1 else 0,
+                    samples = output (trace a) : smples,
+                    lk = exp (ln (probDensity (trace a))) : lk
+                  }
+            )
+            (initialState i)
+            id
+          >-> P.take i
+          >-> P.mapM_ (B.writeBChan chan)
diff --git a/src/Control/Monad/Bayes/Integrator.hs b/src/Control/Monad/Bayes/Integrator.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Integrator.hs
@@ -0,0 +1,167 @@
+{-# LANGUAGE ApplicativeDo #-}
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# OPTIONS_GHC -Wno-type-defaults #-}
+{-# OPTIONS_GHC -Wno-unused-top-binds #-}
+
+-- |
+-- This is adapted from https://jtobin.io/giry-monad-implementation
+-- but brought into the monad-bayes framework (i.e. Integrator is an instance of MonadInfer)
+-- It's largely for debugging other inference methods and didactic use,
+-- because brute force integration of measures is
+-- only practical for small programs
+module Control.Monad.Bayes.Integrator
+  ( probability,
+    variance,
+    expectation,
+    cdf,
+    empirical,
+    enumeratorWith,
+    histogram,
+    plotCdf,
+    volume,
+    normalize,
+    Integrator,
+    momentGeneratingFunction,
+    cumulantGeneratingFunction,
+    integrator,
+    runIntegrator,
+  )
+where
+
+import Control.Applicative (Applicative (..))
+import Control.Foldl (Fold)
+import Control.Foldl qualified as Foldl
+import Control.Monad.Bayes.Class (MonadSample (bernoulli, random, uniformD))
+import Control.Monad.Bayes.Weighted (Weighted, weighted)
+import Control.Monad.Cont
+  ( Cont,
+    ContT (ContT),
+    cont,
+    runCont,
+  )
+import Data.Foldable (Foldable (foldl'))
+import Data.Set (Set, elems)
+import Numeric.Integration.TanhSinh (Result (result), trap)
+import Numeric.Log (Log (ln))
+import Statistics.Distribution qualified as Statistics
+import Statistics.Distribution.Uniform qualified as Statistics
+
+newtype Integrator a = Integrator {getCont :: Cont Double a}
+  deriving newtype (Functor, Applicative, Monad)
+
+integrator, runIntegrator :: (a -> Double) -> Integrator a -> Double
+integrator f (Integrator a) = runCont a f
+runIntegrator = integrator
+
+instance MonadSample Integrator where
+  random = fromDensityFunction $ Statistics.density $ Statistics.uniformDistr 0 1
+  bernoulli p = Integrator $ cont (\f -> p * f True + (1 - p) * f False)
+  uniformD ls = fromMassFunction (const (1 / fromIntegral (length ls))) ls
+
+fromDensityFunction :: (Double -> Double) -> Integrator Double
+fromDensityFunction d = Integrator $
+  cont $ \f ->
+    integralWithQuadrature (\x -> f x * d x)
+  where
+    integralWithQuadrature = result . last . (\z -> trap z 0 1)
+
+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 = 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 f = Foldl.premap f averageFold
+
+    averageFold :: Fractional a => Fold a a
+    averageFold = (/) <$> Foldl.sum <*> Foldl.genericLength
+
+expectation :: Integrator Double -> Double
+expectation = integrator id
+
+variance :: Integrator Double -> Double
+variance nu = integrator (^ 2) nu - expectation nu ^ 2
+
+momentGeneratingFunction :: Integrator Double -> Double -> Double
+momentGeneratingFunction nu t = integrator (\x -> exp (t * x)) nu
+
+cumulantGeneratingFunction :: Integrator Double -> Double -> Double
+cumulantGeneratingFunction nu = log . momentGeneratingFunction nu
+
+normalize :: Weighted Integrator a -> Integrator a
+normalize m =
+  let m' = weighted m
+      z = integrator (ln . exp . snd) m'
+   in do
+        (x, d) <- weighted m
+        Integrator $ cont $ \f -> (f () * (ln $ exp d)) / z
+        return x
+
+cdf :: Integrator Double -> Double -> Double
+cdf nu x = integrator (negativeInfinity `to` x) nu
+  where
+    negativeInfinity :: Double
+    negativeInfinity = negate (1 / 0)
+
+    to :: (Num a, Ord a) => a -> a -> a -> a
+    to a b k
+      | k >= a && k <= b = 1
+      | otherwise = 0
+
+volume :: Integrator Double -> Double
+volume = integrator (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
+  (+) = liftA2 (+)
+  (-) = liftA2 (-)
+  (*) = liftA2 (*)
+  abs = fmap abs
+  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)
+
+enumeratorWith :: Ord a => Set a -> Integrator a -> [(a, Double)]
+enumeratorWith ls meas =
+  [ ( val,
+      integrator
+        (\x -> if x == val then 1 else 0)
+        meas
+    )
+    | val <- elems ls
+  ]
+
+histogram ::
+  (Enum a, Ord a, Fractional a) =>
+  Int ->
+  a ->
+  Weighted Integrator a ->
+  [(a, Double)]
+histogram nBins binSize model = do
+  x <- take nBins [1 ..]
+  let transform k = (k - (fromIntegral nBins / 2)) * binSize
+  return
+    ( (fst)
+        (transform x, transform (x + 1)),
+      probability (transform x, transform (x + 1)) $ normalize model
+    )
+
+plotCdf :: Int -> Double -> Double -> Integrator Double -> [(Double, Double)]
+plotCdf nBins binSize middlePoint model = do
+  x <- take nBins [1 ..]
+  let transform k = (k - (fromIntegral nBins / 2)) * binSize + middlePoint
+  return (transform x, cdf model (transform x))
diff --git a/src/Control/Monad/Bayes/Population.hs b/src/Control/Monad/Bayes/Population.hs
--- a/src/Control/Monad/Bayes/Population.hs
+++ b/src/Control/Monad/Bayes/Population.hs
@@ -1,3 +1,9 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# LANGUAGE RankNTypes #-}
+{-# OPTIONS_GHC -Wno-deprecations #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Population
 -- Description : Representation of distributions using multiple samples
@@ -10,51 +16,71 @@
 -- 'Population' turns a single sample into a collection of weighted samples.
 module Control.Monad.Bayes.Population
   ( Population,
+    population,
     runPopulation,
     explicitPopulation,
     fromWeightedList,
     spawn,
+    multinomial,
     resampleMultinomial,
+    systematic,
     resampleSystematic,
+    stratified,
+    resampleStratified,
     extractEvidence,
     pushEvidence,
     proper,
     evidence,
+    hoist,
     collapse,
-    mapPopulation,
-    normalize,
     popAvg,
-    flatten,
-    hoist,
+    withParticles,
   )
 where
 
 import Control.Arrow (second)
 import Control.Monad (replicateM)
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Weighted hiding (flatten, hoist)
-import Control.Monad.Trans
-import Control.Monad.Trans.List
-import qualified Data.List
-import qualified Data.Vector as V
-import Numeric.Log
+  ( MonadCond,
+    MonadInfer,
+    MonadSample (categorical, logCategorical, random, uniform),
+    factor,
+  )
+import Control.Monad.Bayes.Weighted
+  ( Weighted,
+    applyWeight,
+    extractWeight,
+    weighted,
+    withWeight,
+  )
+import Control.Monad.List (ListT (..), MonadIO, MonadTrans (..))
+import Data.List (unfoldr)
+import Data.List qualified
+import Data.Maybe (catMaybes)
+import Data.Vector ((!))
+import Data.Vector qualified as V
+import Numeric.Log (Log, ln, sum)
+import Numeric.Log qualified as Log
 import Prelude hiding (all, sum)
 
 -- | A collection of weighted samples, or particles.
 newtype Population m a = Population (Weighted (ListT m) a)
-  deriving (Functor, Applicative, Monad, MonadIO, MonadSample, MonadCond, MonadInfer)
+  deriving newtype (Functor, Applicative, Monad, MonadIO, MonadSample, MonadCond, MonadInfer)
 
 instance MonadTrans Population where
   lift = Population . lift . lift
 
 -- | Explicit representation of the weighted sample with weights in the log
 -- domain.
-runPopulation :: Functor m => Population m a -> m [(a, Log Double)]
-runPopulation (Population m) = runListT $ runWeighted m
+population, runPopulation :: Population m a -> m [(a, Log Double)]
+population (Population m) = runListT $ weighted m
 
+-- | deprecated synonym
+runPopulation = population
+
 -- | Explicit representation of the weighted sample.
 explicitPopulation :: Functor m => Population m a -> m [(a, Double)]
-explicitPopulation = fmap (map (second (exp . ln))) . runPopulation
+explicitPopulation = fmap (map (second (exp . ln))) . population
 
 -- | Initialize 'Population' with a concrete weighted sample.
 fromWeightedList :: Monad m => m [(a, Log Double)] -> Population m a
@@ -67,6 +93,9 @@
 spawn :: Monad m => Int -> Population m ()
 spawn n = fromWeightedList $ pure $ replicate n ((), 1 / fromIntegral n)
 
+withParticles :: Monad m => Int -> Population m a -> Population m a
+withParticles n = (spawn n >>)
+
 resampleGeneric ::
   MonadSample m =>
   -- | resampler
@@ -74,10 +103,10 @@
   Population m a ->
   Population m a
 resampleGeneric resampler m = fromWeightedList $ do
-  pop <- runPopulation m
+  pop <- population m
   let (xs, ps) = unzip pop
   let n = length xs
-  let z = sum ps
+  let z = Log.sum ps
   if z > 0
     then do
       let weights = V.fromList (map (exp . ln . (/ z)) ps)
@@ -85,10 +114,26 @@
       let xvec = V.fromList xs
       let offsprings = map (xvec V.!) ancestors
       return $ map (,z / fromIntegral n) offsprings
-    else-- if all weights are zero do not resample
+    else -- if all weights are zero do not resample
       return pop
 
--- | Systematic resampling helper.
+-- | Systematic sampler.
+-- Sample \(n\) values from \((0,1]\) as follows
+-- \[
+-- \begin{aligned}
+-- u^{(1)} &\sim U\left(0, \frac{1}{n}\right] \\
+-- u^{(i)} &=u^{(1)}+\frac{i-1}{n}, \quad i=2,3, \ldots, n
+-- \end{aligned}
+-- \]
+-- and then pick integers \(m\) according to
+-- \[
+-- Q^{(m-1)}<u^{(n)} \leq Q^{(m)}
+-- \]
+-- where
+-- \[
+-- Q^{(m)}=\sum_{k=1}^{m} w^{(k)}
+-- \]
+-- and \(w^{(k)}\) are the weights. See also [Comparison of Resampling Schemes for Particle Filtering](https://arxiv.org/abs/cs/0507025).
 systematic :: Double -> V.Vector Double -> [Int]
 systematic u ps = f 0 (u / fromIntegral n) 0 0 []
   where
@@ -98,7 +143,7 @@
     f i _ _ _ acc | i == n = acc
     f i v j q acc =
       if v < q
-        then f (i + 1) (v + inc) j q (j -1 : acc)
+        then f (i + 1) (v + inc) j q (j - 1 : acc)
         else f i v (j + 1) (q + prob j) acc
 
 -- | Resample the population using the underlying monad and a systematic resampling scheme.
@@ -109,7 +154,51 @@
   Population m a
 resampleSystematic = resampleGeneric (\ps -> (`systematic` ps) <$> random)
 
--- | Multinomial resampler.
+-- | Stratified sampler.
+--
+-- Sample \(n\) values from \((0,1]\) as follows
+-- \[
+-- u^{(i)} \sim U\left(\frac{i-1}{n}, \frac{i}{n}\right], \quad i=1,2, \ldots, n
+-- \]
+-- and then pick integers \(m\) according to
+-- \[
+-- Q^{(m-1)}<u^{(n)} \leq Q^{(m)}
+-- \]
+-- where
+-- \[
+-- Q^{(m)}=\sum_{k=1}^{m} w^{(k)}
+-- \]
+-- 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 :: MonadSample m => V.Vector Double -> m [Int]
+stratified weights = do
+  let bigN = V.length weights
+  dithers <- V.replicateM bigN (uniform 0.0 1.0)
+  let positions =
+        V.map (/ fromIntegral bigN) $
+          V.zipWith (+) dithers (V.map fromIntegral $ V.fromList [0 .. bigN - 1])
+      cumulativeSum = V.scanl (+) 0.0 weights
+      coalg (i, j)
+        | i < bigN =
+          if (positions ! i) < (cumulativeSum ! j)
+            then Just (Just j, (i + 1, j))
+            else Just (Nothing, (i, j + 1))
+        | otherwise =
+          Nothing
+  return $ map (\i -> i - 1) $ catMaybes $ unfoldr coalg (0, 0)
+
+-- | Resample the population using the underlying monad and a stratified resampling scheme.
+-- The total weight is preserved.
+resampleStratified ::
+  (MonadSample m) =>
+  Population m a ->
+  Population 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 :: MonadSample m => V.Vector Double -> m [Int]
 multinomial ps = replicateM (V.length ps) (categorical ps)
 
@@ -128,7 +217,7 @@
   Population m a ->
   Population (Weighted m) a
 extractEvidence m = fromWeightedList $ do
-  pop <- lift $ runPopulation m
+  pop <- lift $ population 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
@@ -150,7 +239,7 @@
   Population m a ->
   Weighted m a
 proper m = do
-  pop <- runPopulation $ extractEvidence m
+  pop <- population $ extractEvidence m
   let (xs, ps) = unzip pop
   index <- logCategorical $ V.fromList ps
   let x = xs !! index
@@ -158,7 +247,7 @@
 
 -- | Model evidence estimator, also known as pseudo-marginal likelihood.
 evidence :: (Monad m) => Population m a -> m (Log Double)
-evidence = extractWeight . runPopulation . extractEvidence
+evidence = extractWeight . population . extractEvidence
 
 -- | Picks one point from the population and uses model evidence as a 'score'
 -- in the transformed monad.
@@ -170,19 +259,6 @@
   m a
 collapse = applyWeight . proper
 
--- | Applies a random transformation to a population.
-mapPopulation ::
-  (Monad m) =>
-  ([(a, Log Double)] -> m [(a, Log Double)]) ->
-  Population m a ->
-  Population m a
-mapPopulation f m = fromWeightedList $ runPopulation m >>= f
-
--- | Normalizes the weights in the population so that their sum is 1.
--- This transformation introduces bias.
-normalize :: (Monad m) => Population m a -> Population m a
-normalize = hoist prior . extractEvidence
-
 -- | Population average of a function, computed using unnormalized weights.
 popAvg :: (Monad m) => (a -> Double) -> Population m a -> m Double
 popAvg f p = do
@@ -191,20 +267,10 @@
   let t = Data.List.sum ys
   return t
 
--- | Combine a population of populations into a single population.
-flatten :: Monad m => Population (Population m) a -> Population m a
-flatten m = Population $ withWeight $ ListT t
-  where
-    t = f <$> (runPopulation . runPopulation) m
-    f d = do
-      (x, p) <- d
-      (y, q) <- x
-      return (y, p * q)
-
 -- | Applies a transformation to the inner monad.
 hoist ::
-  (Monad m, Monad n) =>
+  Monad n =>
   (forall x. m x -> n x) ->
   Population m a ->
   Population n a
-hoist f = fromWeightedList . f . runPopulation
+hoist f = fromWeightedList . f . population
diff --git a/src/Control/Monad/Bayes/Sampler.hs b/src/Control/Monad/Bayes/Sampler.hs
deleted file mode 100644
--- a/src/Control/Monad/Bayes/Sampler.hs
+++ /dev/null
@@ -1,105 +0,0 @@
--- |
--- Module      : Control.Monad.Bayes.Sampler
--- Description : Pseudo-random sampling monads
--- Copyright   : (c) Adam Scibior, 2015-2020
--- License     : MIT
--- Maintainer  : leonhard.markert@tweag.io
--- Stability   : experimental
--- Portability : GHC
---
--- 'SamplerIO' and 'SamplerST' are instances of 'MonadSample'. Apply a 'MonadCond'
--- transformer to obtain a 'MonadInfer' that can execute probabilistic models.
-module Control.Monad.Bayes.Sampler
-  ( SamplerIO,
-    sampleIO,
-    sampleIOfixed,
-    sampleIOwith,
-    Seed,
-    SamplerST (SamplerST),
-    runSamplerST,
-    sampleST,
-    sampleSTfixed,
-  )
-where
-
-import Control.Monad.Bayes.Class
-import Control.Monad.ST (ST, runST, stToIO)
-import Control.Monad.State (State, state)
-import Control.Monad.Trans (MonadIO, lift)
-import Control.Monad.Trans.Reader (ReaderT, ask, mapReaderT, runReaderT)
-import System.Random.MWC
-import qualified System.Random.MWC.Distributions as MWC
-
--- | An 'IO' based random sampler using the MWC-Random package.
-newtype SamplerIO a = SamplerIO (ReaderT GenIO IO a)
-  deriving (Functor, Applicative, Monad, MonadIO)
-
--- | Initialize a pseudo-random number generator using randomness supplied by
--- the operating system.
--- For efficiency this operation should be applied at the very end, ideally
--- once per program.
-sampleIO :: SamplerIO a -> IO a
-sampleIO (SamplerIO m) = createSystemRandom >>= runReaderT m
-
--- | Like 'sampleIO', but with a fixed random seed.
--- Useful for reproducibility.
-sampleIOfixed :: SamplerIO a -> IO a
-sampleIOfixed (SamplerIO m) = create >>= runReaderT m
-
--- | Like 'sampleIO' but with a custom pseudo-random number generator.
-sampleIOwith :: SamplerIO a -> GenIO -> IO a
-sampleIOwith (SamplerIO m) = runReaderT m
-
-fromSamplerST :: SamplerST a -> SamplerIO a
-fromSamplerST (SamplerST m) = SamplerIO $ mapReaderT stToIO m
-
-instance MonadSample SamplerIO where
-  random = fromSamplerST random
-
--- | An 'ST' based random sampler using the @mwc-random@ package.
-newtype SamplerST a = SamplerST (forall s. ReaderT (GenST s) (ST s) a)
-
-runSamplerST :: SamplerST a -> ReaderT (GenST s) (ST s) a
-runSamplerST (SamplerST s) = s
-
-instance Functor SamplerST where
-  fmap f (SamplerST s) = SamplerST $ fmap f s
-
-instance Applicative SamplerST where
-  pure x = SamplerST $ pure x
-  (SamplerST f) <*> (SamplerST x) = SamplerST $ f <*> x
-
-instance Monad SamplerST where
-  (SamplerST x) >>= f = SamplerST $ x >>= runSamplerST . f
-
--- | Run the sampler with a supplied seed.
--- Note that 'State Seed' is much less efficient than 'SamplerST' for composing computation.
-sampleST :: SamplerST a -> State Seed a
-sampleST (SamplerST s) =
-  state $ \seed -> runST $ do
-    gen <- restore seed
-    y <- runReaderT s gen
-    finalSeed <- save gen
-    return (y, finalSeed)
-
--- | Run the sampler with a fixed random seed.
-sampleSTfixed :: SamplerST a -> a
-sampleSTfixed (SamplerST s) = runST $ do
-  gen <- create
-  runReaderT s gen
-
--- | Convert a distribution supplied by @mwc-random@.
-fromMWC :: (forall s. GenST s -> ST s a) -> SamplerST a
-fromMWC s = SamplerST $ ask >>= lift . s
-
-instance MonadSample SamplerST where
-  random = fromMWC System.Random.MWC.uniform
-
-  uniform a b = fromMWC $ uniformR (a, b)
-  normal m s = fromMWC $ MWC.normal m s
-  gamma shape scale = fromMWC $ MWC.gamma shape scale
-  beta a b = fromMWC $ MWC.beta a b
-
-  bernoulli p = fromMWC $ MWC.bernoulli p
-  categorical ps = fromMWC $ MWC.categorical ps
-  geometric p = fromMWC $ MWC.geometric0 p
diff --git a/src/Control/Monad/Bayes/Sampler/Lazy.hs b/src/Control/Monad/Bayes/Sampler/Lazy.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Sampler/Lazy.hs
@@ -0,0 +1,79 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- | This is a port of the implementation of LazyPPL: https://lazyppl.bitbucket.io/
+module Control.Monad.Bayes.Sampler.Lazy where
+
+import Control.Monad (ap, liftM)
+import Control.Monad.Bayes.Class (MonadSample (random))
+import Control.Monad.Bayes.Weighted (Weighted, weighted)
+import Numeric.Log (Log (..))
+import System.Random
+  ( RandomGen (split),
+    getStdGen,
+    newStdGen,
+  )
+import qualified System.Random as R
+
+-- | A 'Tree' is a lazy, infinitely wide and infinitely deep tree, labelled by Doubles
+-- | Our source of randomness will be a Tree, populated by uniform [0,1] choices for each label.
+-- | Often people just use a list or stream instead of a tree.
+-- | But a tree allows us to be lazy about how far we are going all the time.
+data Tree = Tree
+  { currentUniform :: Double,
+    lazyUniforms :: Trees
+  }
+
+-- | An infinite stream of 'Tree's.
+data Trees = Trees
+  { headTree :: Tree,
+    tailTrees :: Trees
+  }
+
+-- | 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}
+  deriving (Functor)
+
+-- | Two key things to do with trees:
+-- | Split tree splits a tree in two (bijectively)
+-- | Get the label at the head of the tree and discard the rest
+splitTree :: Tree -> (Tree, Tree)
+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 g = let (a, g') = R.random g in Tree a (randomTrees g')
+
+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
+  (<*>) = ap
+
+-- | probabilities for a monad.
+-- | Sequencing is done by splitting the tree
+-- | and using different bits for different computations.
+instance Monad Sampler where
+  return = pure
+  (Sampler m) >>= f = Sampler \g ->
+    let (g1, g2) = splitTree g
+        (Sampler m') = f (m g1)
+     in m' g2
+
+instance MonadSample Sampler where
+  random = Sampler \(Tree r _) -> r
+
+sampler :: Sampler a -> IO a
+sampler m = newStdGen *> (runSampler m . randomTree <$> getStdGen)
+
+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
diff --git a/src/Control/Monad/Bayes/Sampler/Strict.hs b/src/Control/Monad/Bayes/Sampler/Strict.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Sampler/Strict.hs
@@ -0,0 +1,101 @@
+{-# LANGUAGE ApplicativeDo #-}
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE ImportQualifiedPost #-}
+
+-- |
+-- Module      : Control.Monad.Bayes.Sampler
+-- Description : Pseudo-random sampling monads
+-- Copyright   : (c) Adam Scibior, 2015-2020
+-- License     : MIT
+-- Maintainer  : leonhard.markert@tweag.io
+-- Stability   : experimental
+-- Portability : GHC
+--
+-- 'SamplerIO' and 'SamplerST' are instances of 'MonadSample'. Apply a 'MonadCond'
+-- transformer to obtain a 'MonadInfer' that can execute probabilistic models.
+module Control.Monad.Bayes.Sampler.Strict
+  ( Sampler,
+    SamplerIO,
+    SamplerST,
+    sampleIO,
+    sampleIOfixed,
+    sampleWith,
+    sampleSTfixed,
+    sampleMean,
+    sampler,
+  )
+where
+
+import Control.Foldl qualified as F hiding (random)
+import Control.Monad.Bayes.Class
+  ( MonadSample
+      ( bernoulli,
+        beta,
+        categorical,
+        gamma,
+        geometric,
+        normal,
+        random,
+        uniform
+      ),
+  )
+import Control.Monad.Reader (MonadIO, ReaderT (..))
+import Control.Monad.ST (ST)
+import Numeric.Log (Log (ln))
+import System.Random.MWC.Distributions qualified as MWC
+import System.Random.Stateful (IOGenM (..), STGenM, StatefulGen, StdGen, initStdGen, mkStdGen, newIOGenM, newSTGenM, uniformDouble01M, uniformRM)
+
+-- | 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)
+
+-- | convenient type synonym to show specializations of Sampler
+-- to particular pairs of monad and RNG
+type SamplerIO = Sampler (IOGenM StdGen) IO
+
+-- | convenient type synonym to show specializations of Sampler
+-- to particular pairs of monad and RNG
+type SamplerST s = Sampler (STGenM StdGen s) (ST s)
+
+instance StatefulGen g m => MonadSample (Sampler g m) where
+  random = Sampler (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)
+
+  bernoulli p = Sampler (ReaderT $ MWC.bernoulli p)
+  categorical ps = Sampler (ReaderT $ MWC.categorical ps)
+  geometric p = Sampler (ReaderT $ MWC.geometric0 p)
+
+-- | Sample with a random number generator of your choice e.g. the one
+-- from `System.Random`.
+--
+-- >>> import Control.Monad.Bayes.Class
+-- >>> import System.Random.Stateful hiding (random)
+-- >>> newIOGenM (mkStdGen 1729) >>= sampleWith random
+-- 4.690861245089605e-2
+sampleWith :: StatefulGen g m => Sampler g m a -> g -> m a
+sampleWith (Sampler m) = runReaderT m
+
+-- | initialize random seed using system entropy, and sample
+sampleIO, sampler :: SamplerIO a -> IO a
+sampleIO x = initStdGen >>= newIOGenM >>= sampleWith x
+sampler = sampleIO
+
+-- | Run the sampler with a fixed random seed
+sampleIOfixed :: SamplerIO a -> IO a
+sampleIOfixed x = newIOGenM (mkStdGen 1729) >>= sampleWith x
+
+-- | Run the sampler with a fixed random seed
+sampleSTfixed :: SamplerST s b -> ST s b
+sampleSTfixed x = newSTGenM (mkStdGen 1729) >>= sampleWith x
+
+sampleMean :: [(Double, Log Double)] -> Double
+sampleMean samples =
+  let z = F.premap (ln . exp . snd) F.sum
+      w = (F.premap (\(x, y) -> x * ln (exp y)) F.sum)
+      s = (/) <$> w <*> z
+   in F.fold s samples
diff --git a/src/Control/Monad/Bayes/Sequential.hs b/src/Control/Monad/Bayes/Sequential.hs
deleted file mode 100644
--- a/src/Control/Monad/Bayes/Sequential.hs
+++ /dev/null
@@ -1,98 +0,0 @@
--- |
--- Module      : Control.Monad.Bayes.Sequential
--- Description : Suspendable probabilistic computation
--- Copyright   : (c) Adam Scibior, 2015-2020
--- License     : MIT
--- Maintainer  : leonhard.markert@tweag.io
--- Stability   : experimental
--- Portability : GHC
---
--- 'Sequential' represents a computation that can be suspended.
-module Control.Monad.Bayes.Sequential
-  ( Sequential,
-    suspend,
-    finish,
-    advance,
-    finished,
-    hoistFirst,
-    hoist,
-    sis,
-  )
-where
-
-import Control.Monad.Bayes.Class
-import Control.Monad.Coroutine hiding (suspend)
-import Control.Monad.Coroutine.SuspensionFunctors
-import Control.Monad.Trans
-import Data.Either
-
--- | Represents a computation that can be suspended at certain points.
--- The intermediate monadic effects can be extracted, which is particularly
--- 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}
-  deriving (Functor, Applicative, Monad, MonadTrans, MonadIO)
-
-extract :: Await () a -> a
-extract (Await f) = f ()
-
-instance MonadSample m => MonadSample (Sequential m) where
-  random = lift random
-  bernoulli = lift . bernoulli
-  categorical = lift . categorical
-
--- | Execution is 'suspend'ed after each 'score'.
-instance MonadCond m => MonadCond (Sequential m) where
-  score w = lift (score w) >> suspend
-
-instance MonadInfer m => MonadInfer (Sequential m)
-
--- | A point where the computation is paused.
-suspend :: Monad m => Sequential m ()
-suspend = Sequential await
-
--- | Remove the remaining suspension points.
-finish :: Monad m => Sequential m a -> m a
-finish = pogoStick extract . runSequential
-
--- | 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
-
--- | Return True if no more suspension points remain.
-finished :: Monad m => Sequential m a -> m Bool
-finished = fmap isRight . resume . runSequential
-
--- | 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
-
--- | 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
-
--- | Apply a function a given number of times.
-composeCopies :: Int -> (a -> a) -> (a -> a)
-composeCopies k f = foldr (.) id (replicate k f)
-
--- | Sequential importance sampling.
--- Applies a given transformation after each time step.
-sis ::
-  Monad m =>
-  -- | transformation
-  (forall x. m x -> m x) ->
-  -- | number of time steps
-  Int ->
-  Sequential m a ->
-  m a
-sis f k = finish . composeCopies k (advance . hoistFirst f)
diff --git a/src/Control/Monad/Bayes/Sequential/Coroutine.hs b/src/Control/Monad/Bayes/Sequential/Coroutine.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Sequential/Coroutine.hs
@@ -0,0 +1,119 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE RankNTypes #-}
+
+-- |
+-- Module      : Control.Monad.Bayes.Sequential
+-- Description : Suspendable probabilistic computation
+-- Copyright   : (c) Adam Scibior, 2015-2020
+-- License     : MIT
+-- Maintainer  : leonhard.markert@tweag.io
+-- Stability   : experimental
+-- Portability : GHC
+--
+-- 'Sequential' represents a computation that can be suspended.
+module Control.Monad.Bayes.Sequential.Coroutine
+  ( Sequential,
+    suspend,
+    finish,
+    advance,
+    finished,
+    hoistFirst,
+    hoist,
+    sequentially,
+    sis,
+  )
+where
+
+import Control.Monad.Bayes.Class
+  ( MonadCond (..),
+    MonadInfer,
+    MonadSample (bernoulli, categorical, random),
+  )
+import Control.Monad.Coroutine
+  ( Coroutine (..),
+    bounce,
+    mapMonad,
+    pogoStick,
+  )
+import Control.Monad.Coroutine.SuspensionFunctors
+  ( Await (..),
+    await,
+  )
+import Control.Monad.Trans (MonadIO, MonadTrans (..))
+import Data.Either (isRight)
+
+-- | Represents a computation that can be suspended at certain points.
+-- The intermediate monadic effects can be extracted, which is particularly
+-- 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}
+  deriving newtype (Functor, Applicative, Monad, MonadTrans, MonadIO)
+
+extract :: Await () a -> a
+extract (Await f) = f ()
+
+instance MonadSample m => MonadSample (Sequential m) where
+  random = lift random
+  bernoulli = lift . bernoulli
+  categorical = lift . categorical
+
+-- | Execution is 'suspend'ed after each 'score'.
+instance MonadCond m => MonadCond (Sequential m) where
+  score w = lift (score w) >> suspend
+
+instance MonadInfer m => MonadInfer (Sequential m)
+
+-- | A point where the computation is paused.
+suspend :: Monad m => Sequential m ()
+suspend = Sequential await
+
+-- | Remove the remaining suspension points.
+finish :: Monad m => Sequential m a -> m a
+finish = pogoStick extract . runSequential
+
+-- | 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
+
+-- | Return True if no more suspension points remain.
+finished :: Monad m => Sequential m a -> m Bool
+finished = fmap isRight . resume . runSequential
+
+-- | 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
+
+-- | 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
+
+-- | Apply a function a given number of times.
+composeCopies :: Int -> (a -> a) -> (a -> a)
+composeCopies k f = foldr (.) id (replicate k f)
+
+-- | Sequential importance sampling.
+-- Applies a given transformation after each time step.
+sequentially,
+  sis ::
+    Monad m =>
+    -- | transformation
+    (forall x. m x -> m x) ->
+    -- | number of time steps
+    Int ->
+    Sequential m a ->
+    m a
+sequentially f k = finish . composeCopies k (advance . hoistFirst f)
+
+-- | synonym
+sis = sequentially
diff --git a/src/Control/Monad/Bayes/Traced/Basic.hs b/src/Control/Monad/Bayes/Traced/Basic.hs
--- a/src/Control/Monad/Bayes/Traced/Basic.hs
+++ b/src/Control/Monad/Bayes/Traced/Basic.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE RankNTypes #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Traced.Basic
 -- Description : Distributions on full execution traces of full programs
@@ -8,7 +10,7 @@
 -- Portability : GHC
 module Control.Monad.Bayes.Traced.Basic
   ( Traced,
-    hoistT,
+    hoist,
     marginal,
     mhStep,
     mh,
@@ -17,19 +19,29 @@
 
 import Control.Applicative (liftA2)
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Free (FreeSampler)
+  ( MonadCond (..),
+    MonadInfer,
+    MonadSample (random),
+  )
+import Control.Monad.Bayes.Density.Free (Density)
 import Control.Monad.Bayes.Traced.Common
+  ( Trace (..),
+    bind,
+    mhTrans',
+    scored,
+    singleton,
+  )
 import Control.Monad.Bayes.Weighted (Weighted)
 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
-      { -- | Run the program with a modified trace.
-        model :: Weighted (FreeSampler Identity) a,
-        -- | Record trace and output.
-        traceDist :: m (Trace a)
-      }
+data Traced m a = Traced
+  { -- | Run the program with a modified trace.
+    model :: Weighted (Density 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)
@@ -52,8 +64,8 @@
 
 instance MonadInfer m => MonadInfer (Traced m)
 
-hoistT :: (forall x. m x -> m x) -> Traced m a -> Traced m a
-hoistT f (Traced m d) = Traced m (f d)
+hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a
+hoist f (Traced m d) = Traced m (f d)
 
 -- | Discard the trace and supporting infrastructure.
 marginal :: Monad m => Traced m a -> m a
@@ -68,10 +80,11 @@
 -- | Full run of the Trace Metropolis-Hastings algorithm with a specified
 -- number of steps.
 mh :: MonadSample m => Int -> Traced m a -> m [a]
-mh n (Traced m d) = fmap (map output) (f n)
+mh n (Traced m d) = fmap (map output . NE.toList) (f n)
   where
-    f 0 = fmap (: []) d
-    f k = do
-      ~(x : xs) <- f (k -1)
-      y <- mhTrans' m x
-      return (y : x : xs)
+    f k
+      | k <= 0 = fmap (:| []) d
+      | otherwise = do
+        (x :| xs) <- f (k - 1)
+        y <- mhTrans' m x
+        return (y :| x : xs)
diff --git a/src/Control/Monad/Bayes/Traced/Common.hs b/src/Control/Monad/Bayes/Traced/Common.hs
--- a/src/Control/Monad/Bayes/Traced/Common.hs
+++ b/src/Control/Monad/Bayes/Traced/Common.hs
@@ -7,78 +7,118 @@
 -- Stability   : experimental
 -- Portability : GHC
 module Control.Monad.Bayes.Traced.Common
-  ( Trace,
+  ( Trace (..),
     singleton,
-    output,
     scored,
     bind,
     mhTrans,
+    mhTransWithBool,
+    mhTransFree,
     mhTrans',
+    burnIn,
+    MHResult (..),
   )
 where
 
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Free as FreeSampler
+  ( MonadSample (bernoulli, random),
+    discrete,
+  )
+import qualified Control.Monad.Bayes.Density.Free as Free
+import qualified Control.Monad.Bayes.Density.State as State
 import Control.Monad.Bayes.Weighted as Weighted
-import Control.Monad.Trans.Writer
-import Data.Functor.Identity
+  ( Weighted,
+    hoist,
+    weighted,
+  )
+import Control.Monad.Writer (WriterT (WriterT, runWriterT))
+import Data.Functor.Identity (Identity (runIdentity))
 import Numeric.Log (Log, ln)
 import Statistics.Distribution.DiscreteUniform (discreteUniformAB)
 
--- | Collection of random variables sampled during the program's execution.
-data Trace a
-  = Trace
-      { -- | Sequence of random variables sampled during the program's execution.
-        variables :: [Double],
-        -- |
-        output :: a,
-        -- | The probability of observing this particular sequence.
-        density :: Log Double
-      }
+data MHResult a = MHResult
+  { success :: Bool,
+    trace :: Trace a
+  }
 
+-- | Collection of random variables sampler during the program's execution.
+data Trace a = Trace
+  { -- | Sequence of random variables sampler during the program's execution.
+    variables :: [Double],
+    --
+    output :: a,
+    -- | The probability of observing this particular sequence.
+    probDensity :: Log Double
+  }
+
 instance Functor Trace where
   fmap f t = t {output = f (output t)}
 
 instance Applicative Trace where
-  pure x = Trace {variables = [], output = x, density = 1}
+  pure x = Trace {variables = [], output = x, probDensity = 1}
   tf <*> tx =
     Trace
       { variables = variables tf ++ variables tx,
         output = output tf (output tx),
-        density = density tf * density tx
+        probDensity = probDensity tf * probDensity tx
       }
 
 instance Monad Trace where
   t >>= f =
     let t' = f (output t)
-     in t' {variables = variables t ++ variables t', density = density t * density t'}
+     in t' {variables = variables t ++ variables t', probDensity = probDensity t * probDensity t'}
 
 singleton :: Double -> Trace Double
-singleton u = Trace {variables = [u], output = u, density = 1}
+singleton u = Trace {variables = [u], output = u, probDensity = 1}
 
 scored :: Log Double -> Trace ()
-scored w = Trace {variables = [], output = (), density = w}
+scored w = Trace {variables = [], output = (), probDensity = w}
 
 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, density = density t1 * density t2}
+  return $ t2 {variables = variables t1 ++ variables t2, probDensity = probDensity t1 * probDensity t2}
 
 -- | A single Metropolis-corrected transition of single-site Trace MCMC.
-mhTrans :: MonadSample m => Weighted (FreeSampler m) a -> Trace a -> m (Trace a)
-mhTrans m t@Trace {variables = us, density = p} = do
+mhTrans :: MonadSample m => (Weighted (State.Density m)) a -> Trace a -> m (Trace a)
+mhTrans m t@Trace {variables = us, probDensity = p} = do
   let n = length us
   us' <- do
-    i <- discrete $ discreteUniformAB 0 (n -1)
+    i <- discrete $ discreteUniformAB 0 (n - 1)
     u' <- random
-    let (xs, _ : ys) = splitAt i us
-    return $ xs ++ (u' : ys)
-  ((b, q), vs) <- runWriterT $ runWeighted $ Weighted.hoist (WriterT . withPartialRandomness us') m
+    case splitAt i us of
+      (xs, _ : ys) -> return $ xs ++ (u' : ys)
+      _ -> error "impossible"
+  ((b, q), vs) <- State.density (weighted 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 :: MonadSample m => Weighted (Free.Density 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 :: MonadSample m => Weighted (Free.Density m) a -> Trace a -> m (MHResult a)
+mhTransWithBool m t@Trace {variables = us, probDensity = p} = do
+  let n = length us
+  us' <- do
+    i <- discrete $ discreteUniformAB 0 (n - 1)
+    u' <- random
+    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
+  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' :: MonadSample m => Weighted (FreeSampler Identity) a -> Trace a -> m (Trace a)
-mhTrans' m = mhTrans (Weighted.hoist (FreeSampler.hoist (return . runIdentity)) m)
+mhTrans' :: MonadSample m => Weighted (Free.Density Identity) a -> Trace a -> m (Trace a)
+mhTrans' m = mhTransFree (Weighted.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 n = fmap dropEnd
+  where
+    dropEnd ls = let len = length ls in take (len - n) ls
diff --git a/src/Control/Monad/Bayes/Traced/Dynamic.hs b/src/Control/Monad/Bayes/Traced/Dynamic.hs
--- a/src/Control/Monad/Bayes/Traced/Dynamic.hs
+++ b/src/Control/Monad/Bayes/Traced/Dynamic.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE RankNTypes #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Traced.Dynamic
 -- Description : Distributions on execution traces that can be dynamically frozen
@@ -8,7 +10,7 @@
 -- Portability : GHC
 module Control.Monad.Bayes.Traced.Dynamic
   ( Traced,
-    hoistT,
+    hoist,
     marginal,
     freeze,
     mhStep,
@@ -18,16 +20,27 @@
 
 import Control.Monad (join)
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Free (FreeSampler)
+  ( MonadCond (..),
+    MonadInfer,
+    MonadSample (random),
+  )
+import Control.Monad.Bayes.Density.Free (Density)
 import Control.Monad.Bayes.Traced.Common
+  ( Trace (..),
+    bind,
+    mhTransFree,
+    scored,
+    singleton,
+  )
 import Control.Monad.Bayes.Weighted (Weighted)
 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 (FreeSampler m) a, Trace a)}
+newtype Traced m a = Traced {runTraced :: m (Weighted (Density m) a, Trace a)}
 
-pushM :: Monad m => m (Weighted (FreeSampler m) a) -> Weighted (FreeSampler m) a
+pushM :: Monad m => m (Weighted (Density m) a) -> Weighted (Density m) a
 pushM = join . lift . lift
 
 instance Monad m => Functor (Traced m) where
@@ -62,8 +75,8 @@
 
 instance MonadInfer m => MonadInfer (Traced m)
 
-hoistT :: (forall x. m x -> m x) -> Traced m a -> Traced m a
-hoistT f (Traced c) = Traced (f c)
+hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a
+hoist f (Traced c) = Traced (f c)
 
 -- | Discard the trace and supporting infrastructure.
 marginal :: Monad m => Traced m a -> m a
@@ -81,7 +94,7 @@
 mhStep :: MonadSample m => Traced m a -> Traced m a
 mhStep (Traced c) = Traced $ do
   (m, t) <- c
-  t' <- mhTrans m t
+  t' <- mhTransFree m t
   return (m, t')
 
 -- | Full run of the Trace Metropolis-Hastings algorithm with a specified
@@ -89,11 +102,10 @@
 mh :: MonadSample m => Int -> Traced m a -> m [a]
 mh n (Traced c) = do
   (m, t) <- c
-  let f 0 = return [t]
-      f k = do
-        ~(x : xs) <- f (k -1)
-        y <- mhTrans m x
-        return (y : x : xs)
-  ts <- f n
-  let xs = map output ts
-  return xs
+  let f k
+        | k <= 0 = return (t :| [])
+        | otherwise = do
+          (x :| xs) <- f (k - 1)
+          y <- mhTransFree m x
+          return (y :| x : xs)
+  fmap (map output . NE.toList) (f n)
diff --git a/src/Control/Monad/Bayes/Traced/Static.hs b/src/Control/Monad/Bayes/Traced/Static.hs
--- a/src/Control/Monad/Bayes/Traced/Static.hs
+++ b/src/Control/Monad/Bayes/Traced/Static.hs
@@ -1,3 +1,6 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE RecordWildCards #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Traced.Static
 -- Description : Distributions on execution traces of full programs
@@ -7,8 +10,8 @@
 -- Stability   : experimental
 -- Portability : GHC
 module Control.Monad.Bayes.Traced.Static
-  ( Traced,
-    hoistT,
+  ( Traced (..),
+    hoist,
     marginal,
     mhStep,
     mh,
@@ -17,20 +20,30 @@
 
 import Control.Applicative (liftA2)
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Free (FreeSampler)
+  ( MonadCond (..),
+    MonadInfer,
+    MonadSample (random),
+  )
+import Control.Monad.Bayes.Density.Free (Density)
 import Control.Monad.Bayes.Traced.Common
+  ( Trace (..),
+    bind,
+    mhTransFree,
+    scored,
+    singleton,
+  )
 import Control.Monad.Bayes.Weighted (Weighted)
 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.
 --
 -- The random choices that are not to be traced should be lifted from the
 -- transformed monad.
-data Traced m a
-  = Traced
-      { model :: Weighted (FreeSampler m) a,
-        traceDist :: m (Trace a)
-      }
+data Traced m a = Traced
+  { model :: Weighted (Density 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)
@@ -56,8 +69,8 @@
 
 instance MonadInfer m => MonadInfer (Traced m)
 
-hoistT :: (forall x. m x -> m x) -> Traced m a -> Traced m a
-hoistT f (Traced m d) = Traced m (f d)
+hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a
+hoist f (Traced m d) = Traced m (f d)
 
 -- | Discard the trace and supporting infrastructure.
 marginal :: Monad m => Traced m a -> m a
@@ -67,15 +80,16 @@
 mhStep :: MonadSample m => Traced m a -> Traced m a
 mhStep (Traced m d) = Traced m d'
   where
-    d' = d >>= mhTrans m
+    d' = d >>= mhTransFree m
 
 -- | Full run of the Trace Metropolis-Hastings algorithm with a specified
--- number of steps.
+-- number of steps. Newest samples are at the head of the list.
 mh :: MonadSample m => Int -> Traced m a -> m [a]
-mh n (Traced m d) = fmap (map output) (f n)
+mh n (Traced m d) = fmap (map output . NE.toList) (f n)
   where
-    f 0 = fmap (: []) d
-    f k = do
-      ~(x : xs) <- f (k -1)
-      y <- mhTrans m x
-      return (y : x : xs)
+    f k
+      | k <= 0 = fmap (:| []) d
+      | otherwise = do
+        (x :| xs) <- f (k - 1)
+        y <- mhTransFree m x
+        return (y :| x : xs)
diff --git a/src/Control/Monad/Bayes/Weighted.hs b/src/Control/Monad/Bayes/Weighted.hs
--- a/src/Control/Monad/Bayes/Weighted.hs
+++ b/src/Control/Monad/Bayes/Weighted.hs
@@ -1,3 +1,7 @@
+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+{-# LANGUAGE RankNTypes #-}
+
 -- |
 -- Module      : Control.Monad.Bayes.Weighted
 -- Description : Probability monad accumulating the likelihood
@@ -12,24 +16,28 @@
 module Control.Monad.Bayes.Weighted
   ( Weighted,
     withWeight,
-    runWeighted,
+    weighted,
     extractWeight,
-    prior,
-    flatten,
+    unweighted,
     applyWeight,
     hoist,
+    runWeighted,
   )
 where
 
 import Control.Monad.Bayes.Class
-import Control.Monad.Trans (MonadIO, MonadTrans (..))
-import Control.Monad.Trans.State (StateT (..), mapStateT, modify)
+  ( MonadCond (..),
+    MonadInfer,
+    MonadSample,
+    factor,
+  )
+import Control.Monad.State (MonadIO, MonadTrans, StateT (..), lift, mapStateT, modify)
 import Numeric.Log (Log)
 
 -- | Execute the program using the prior distribution, while accumulating likelihood.
 newtype Weighted m a = Weighted (StateT (Log Double) m a)
   -- StateT is more efficient than WriterT
-  deriving (Functor, Applicative, Monad, MonadIO, MonadTrans, MonadSample)
+  deriving newtype (Functor, Applicative, Monad, MonadIO, MonadTrans, MonadSample)
 
 instance Monad m => MonadCond (Weighted m) where
   score w = Weighted (modify (* w))
@@ -37,36 +45,33 @@
 instance MonadSample m => MonadInfer (Weighted m)
 
 -- | Obtain an explicit value of the likelihood for a given value.
-runWeighted :: (Functor m) => Weighted m a -> m (a, Log Double)
-runWeighted (Weighted m) = runStateT m 1
+weighted, runWeighted :: Weighted m a -> m (a, Log Double)
+weighted (Weighted m) = runStateT m 1
+runWeighted = weighted
 
 -- | Compute the sample and discard the weight.
 --
 -- This operation introduces bias.
-prior :: Functor m => Weighted m a -> m a
-prior = fmap fst . runWeighted
+unweighted :: Functor m => Weighted m a -> m a
+unweighted = fmap fst . weighted
 
 -- | Compute the weight and discard the sample.
 extractWeight :: Functor m => Weighted m a -> m (Log Double)
-extractWeight = fmap snd . runWeighted
+extractWeight = fmap snd . weighted
 
 -- | Embed a random variable with explicitly given likelihood.
 --
--- > runWeighted . withWeight = id
+-- > weighted . withWeight = id
 withWeight :: (Monad m) => m (a, Log Double) -> Weighted m a
 withWeight m = Weighted $ do
   (x, w) <- lift m
   modify (* w)
   return x
 
--- | Combine weights from two different levels.
-flatten :: Monad m => Weighted (Weighted m) a -> Weighted m a
-flatten m = withWeight $ (\((x, p), q) -> (x, p * q)) <$> runWeighted (runWeighted m)
-
 -- | Use the weight as a factor in the transformed monad.
 applyWeight :: MonadCond m => Weighted m a -> m a
 applyWeight m = do
-  (x, w) <- runWeighted m
+  (x, w) <- weighted m
   factor w
   return x
 
diff --git a/src/Math/Integrators/StormerVerlet.hs b/src/Math/Integrators/StormerVerlet.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Integrators/StormerVerlet.hs
@@ -0,0 +1,77 @@
+module Math.Integrators.StormerVerlet
+  ( integrateV,
+    stormerVerlet2H,
+    Integrator,
+  )
+where
+
+import Control.Lens
+import Control.Monad.Primitive
+import Data.Vector (Vector, (!))
+import qualified Data.Vector as V
+import Data.Vector.Mutable
+import Linear (V2 (..))
+
+-- | Integrator function
+-- -   \Phi [h] |->  y_0 -> y_1
+type Integrator a =
+  -- | Step
+  Double ->
+  -- | Initial value
+  a ->
+  -- | Next value
+  a
+
+-- | Störmer-Verlet integration scheme for systems of the form
+-- \(\mathbb{H}(p,q) = T(p) + V(q)\)
+stormerVerlet2H ::
+  (Applicative f, Num (f a), Show (f a), Fractional a) =>
+  -- | Step size
+  a ->
+  -- | \(\frac{\partial H}{\partial q}\)
+  (f a -> f a) ->
+  -- | \(\frac{\partial H}{\partial p}\)
+  (f a -> f a) ->
+  -- | Current \((p, q)\) as a 2-dimensional vector
+  V2 (f a) ->
+  -- | New \((p, q)\) as a 2-dimensional vector
+  V2 (f a)
+stormerVerlet2H hh nablaQ nablaP prev =
+  V2 qNew pNew
+  where
+    h2 = hh / 2
+    hhs = pure hh
+    hh2s = pure h2
+    qsPrev = prev ^. _1
+    psPrev = prev ^. _2
+    pp2 = psPrev - hh2s * nablaQ qsPrev
+    qNew = qsPrev + hhs * nablaP pp2
+    pNew = pp2 - hh2s * nablaQ qNew
+
+-- |
+-- Integrate ODE equation using fixed steps set by a vector, and returns a vector
+-- 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 =>
+  -- | Internal integrator
+  Integrator a ->
+  -- | initial  value
+  a ->
+  -- | vector of time points
+  Vector Double ->
+  -- | vector of solution
+  m (Vector a)
+integrateV integrator initial times = do
+  out <- new (V.length times)
+  write out 0 initial
+  compute initial 1 out
+  V.unsafeFreeze out
+  where
+    compute y i out
+      | i == V.length times = return ()
+      | otherwise = do
+        let h = (times ! i) - (times ! (i - 1))
+            y' = integrator h y
+        write out i y'
+        compute y' (i + 1) out
diff --git a/test/Spec.hs b/test/Spec.hs
--- a/test/Spec.hs
+++ b/test/Spec.hs
@@ -1,67 +1,166 @@
-import Test.Hspec
-import Test.Hspec.QuickCheck
-import Test.QuickCheck
-import qualified TestEnumerator
-import qualified TestInference
-import qualified TestPopulation
-import qualified TestSequential
-import qualified TestWeighted
+{-# LANGUAGE ImportQualifiedPost #-}
 
+import Data.AEq (AEq ((~==)))
+import Test.Hspec (context, describe, hspec, it, shouldBe)
+import Test.Hspec.QuickCheck (prop)
+import Test.QuickCheck (ioProperty, property, (==>))
+import TestAdvanced qualified
+import TestDistribution qualified
+import TestEnumerator qualified
+import TestInference qualified
+import TestIntegrator qualified
+import TestPipes (hmms)
+import TestPipes qualified
+import TestPopulation qualified
+import TestSampler qualified
+import TestSequential qualified
+import TestStormerVerlet qualified
+import TestWeighted qualified
+
 main :: IO ()
-main = hspec $ do
-  describe "Weighted"
-    $ it "accumulates likelihood correctly"
-    $ do
-      passed <- TestWeighted.passed
-      passed `shouldBe` True
-  describe "Dist" $ do
+main = hspec do
+  describe "Stormer Verlet" $
+    it "conserves energy" $
+      do
+        p1 <- TestStormerVerlet.passed1
+        p1 `shouldBe` True
+  describe "Distribution" $
+    it "gives correct mean, variance and covariance" $
+      do
+        p1 <- TestDistribution.passed1
+        p1 `shouldBe` True
+        p2 <- TestDistribution.passed2
+        p2 `shouldBe` True
+        p3 <- TestDistribution.passed3
+        p3 `shouldBe` True
+  describe "Weighted" $
+    it "accumulates likelihood correctly" $
+      do
+        passed <- TestWeighted.passed
+        passed `shouldBe` True
+  describe "Enumerator" do
     it "sorts samples and aggregates weights" $
       TestEnumerator.passed2 `shouldBe` True
     it "gives correct answer for the sprinkler model" $
       TestEnumerator.passed3 `shouldBe` True
     it "computes expectation correctly" $
       TestEnumerator.passed4 `shouldBe` True
-  describe "Population" $ do
-    context "controlling population" $ do
-      it "preserves the population when not explicitly altered" $ do
+  describe "Integrator Expectation" do
+    prop "expectation numerically" $
+      \mean var ->
+        var > 0 ==> property $ TestIntegrator.normalExpectation mean (sqrt var) ~== mean
+  describe "Integrator Variance" do
+    prop "variance numerically" $
+      \mean var ->
+        var > 0 ==> property $ TestIntegrator.normalVariance mean (sqrt var) ~== var
+  describe "Sampler mean and variance" do
+    it "gets right mean and variance" $
+      TestSampler.testMeanAndVariance `shouldBe` True
+  describe "Integrator Volume" do
+    prop "volume sums to 1" $
+      property $ \case
+        [] -> True
+        ls -> (TestIntegrator.volumeIsOne ls)
+
+  describe "Integrator" do
+    it "" $
+      all
+        (== True)
+        [ TestIntegrator.passed1,
+          TestIntegrator.passed2,
+          TestIntegrator.passed3,
+          TestIntegrator.passed4,
+          TestIntegrator.passed5,
+          TestIntegrator.passed6,
+          TestIntegrator.passed7,
+          TestIntegrator.passed8,
+          TestIntegrator.passed9,
+          TestIntegrator.passed10,
+          TestIntegrator.passed11,
+          TestIntegrator.passed12,
+          TestIntegrator.passed13,
+          TestIntegrator.passed14
+        ]
+        `shouldBe` True
+
+  describe "Population" do
+    context "controlling population" do
+      it "preserves the population when not explicitly altered" do
         popSize <- TestPopulation.popSize
         popSize `shouldBe` 5
-      it "multiplies the number of samples when spawn invoked twice" $ do
+      it "multiplies the number of samples when spawn invoked twice" do
         manySize <- TestPopulation.manySize
         manySize `shouldBe` 15
       it "correctly computes population average" $
         TestPopulation.popAvgCheck `shouldBe` True
-    context "distribution-preserving transformations" $ do
-      it "collapse preserves the distribution" $ do
+    context "distribution-preserving transformations" do
+      it "collapse preserves the distribution" do
         TestPopulation.transCheck1 `shouldBe` True
         TestPopulation.transCheck2 `shouldBe` True
-      it "resample preserves the distribution" $ do
+      it "resample preserves the distribution" do
         TestPopulation.resampleCheck 1 `shouldBe` True
         TestPopulation.resampleCheck 2 `shouldBe` True
-  describe "Sequential" $ do
-    it "stops at every factor" $ do
+  describe "Sequential" do
+    it "stops at every factor" do
       TestSequential.checkTwoSync 0 `shouldBe` True
       TestSequential.checkTwoSync 1 `shouldBe` True
       TestSequential.checkTwoSync 2 `shouldBe` True
     it "preserves the distribution" $
       TestSequential.checkPreserve `shouldBe` True
-    it "produces correct intermediate weights" $ do
+    it "produces correct intermediate weights" do
       TestSequential.checkSync 0 `shouldBe` True
       TestSequential.checkSync 1 `shouldBe` True
       TestSequential.checkSync 2 `shouldBe` True
-  describe "SMC" $ do
+  describe "SMC" do
     it "terminates" $
       seq TestInference.checkTerminateSMC () `shouldBe` ()
     it "preserves the distribution on the sprinkler model" $
       TestInference.checkPreserveSMC `shouldBe` True
     prop "number of particles is equal to its second parameter" $
       \observations particles ->
-        observations >= 0 && particles >= 1 ==> ioProperty $ do
+        observations >= 0 && particles >= 1 ==> ioProperty do
           checkParticles <- TestInference.checkParticles observations particles
           return $ checkParticles == particles
-  describe "SMC with systematic resampling"
-    $ prop "number of particles is equal to its second parameter"
-    $ \observations particles ->
-      observations >= 0 && particles >= 1 ==> ioProperty $ do
-        checkParticles <- TestInference.checkParticlesSystematic observations particles
-        return $ checkParticles == particles
+  describe "SMC with systematic resampling" $
+    prop "number of particles is equal to its second parameter" $
+      \observations particles ->
+        observations >= 0 && particles >= 1 ==> ioProperty do
+          checkParticles <- TestInference.checkParticlesSystematic observations particles
+          return $ checkParticles == particles
+  describe "Equivalent Expectations" do
+    prop "Gamma Normal" $
+      ioProperty . TestInference.testGammaNormal
+    prop "Normal Normal" $
+      \n -> ioProperty (TestInference.testNormalNormal [max (-3) $ min 3 n])
+    prop "Beta Bernoulli" $
+      ioProperty . TestInference.testBetaBernoulli
+  describe "Pipes: Urn" do
+    it "Distributions are equivalent" do
+      TestPipes.urns 10 `shouldBe` True
+  describe "Pipes: HMM" do
+    prop "pipe model is equivalent to standard model" $
+      \num -> property $ hmms $ take 5 num
+
+  describe "SMC with stratified resampling" $
+    prop "number of particles is equal to its second parameter" $
+      \observations particles ->
+        observations >= 0 && particles >= 1 ==> ioProperty do
+          checkParticles <- TestInference.checkParticlesStratified observations particles
+          return $ checkParticles == particles
+
+  describe "Expectation from all inference methods" $
+    it "gives correct answer for the sprinkler model" do
+      passed1 <- TestAdvanced.passed1
+      passed1 `shouldBe` True
+      passed2 <- TestAdvanced.passed2
+      passed2 `shouldBe` True
+      passed3 <- TestAdvanced.passed3
+      passed3 `shouldBe` True
+      passed4 <- TestAdvanced.passed4
+      passed4 `shouldBe` True
+      passed5 <- TestAdvanced.passed5
+      passed5 `shouldBe` True
+      passed6 <- TestAdvanced.passed6
+      passed6 `shouldBe` True
+      passed7 <- TestAdvanced.passed7
+      passed7 `shouldBe` True
diff --git a/test/TestAdvanced.hs b/test/TestAdvanced.hs
new file mode 100644
--- /dev/null
+++ b/test/TestAdvanced.hs
@@ -0,0 +1,65 @@
+module TestAdvanced where
+
+import ConjugatePriors
+  ( betaBernoulli',
+    betaBernoulliAnalytic,
+    gammaNormal',
+    gammaNormalAnalytic,
+    normalNormal',
+    normalNormalAnalytic,
+  )
+import Control.Arrow
+import Control.Monad (join, replicateM)
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Enumerator
+import Control.Monad.Bayes.Inference.MCMC
+import Control.Monad.Bayes.Inference.PMMH
+import Control.Monad.Bayes.Inference.RMSMC
+import Control.Monad.Bayes.Inference.SMC
+import Control.Monad.Bayes.Inference.SMC2
+import Control.Monad.Bayes.Population
+import Control.Monad.Bayes.Sampler.Strict
+import Control.Monad.Bayes.Traced
+import Control.Monad.Bayes.Weighted
+import Numeric.Log (Log)
+
+mcmcConfig :: MCMCConfig
+mcmcConfig = MCMCConfig {numMCMCSteps = 0, numBurnIn = 0, proposal = SingleSiteMH}
+
+smcConfig :: MonadSample m => SMCConfig m
+smcConfig = SMCConfig {numSteps = 0, numParticles = 1000, resampler = resampleMultinomial}
+
+passed1, passed2, passed3, passed4, passed5, passed6, passed7 :: IO Bool
+passed1 = do
+  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
+  return $ close 0.5 sample
+passed3 = do
+  sample <- sampleIOfixed $ population $ rmsmcDynamic mcmcConfig smcConfig random
+  return $ close 0.5 sample
+passed4 = do
+  sample <- sampleIOfixed $ population $ rmsmcBasic mcmcConfig smcConfig random
+  return $ close 0.5 sample
+passed5 = do
+  sample <- sampleIOfixed $ population $ rmsmc mcmcConfig smcConfig random
+  return $ close 0.5 sample
+passed6 = do
+  sample <-
+    fmap join $
+      sampleIOfixed $
+        pmmh
+          mcmcConfig {numMCMCSteps = 100}
+          smcConfig {numSteps = 0, numParticles = 100}
+          random
+          (normal 0)
+  return $ close 0.0 sample
+
+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)
+  return $ close 0.0 sample
+
+close n sample = abs (n - (expectation id $ fromList $ toEmpiricalWeighted sample)) < 0.01
diff --git a/test/TestDistribution.hs b/test/TestDistribution.hs
new file mode 100644
--- /dev/null
+++ b/test/TestDistribution.hs
@@ -0,0 +1,69 @@
+{-# LANGUAGE ImportQualifiedPost #-}
+{-# LANGUAGE Trustworthy #-}
+
+module TestDistribution
+  ( passed1,
+    passed2,
+    passed3,
+  )
+where
+
+import Control.Monad (replicateM)
+import Control.Monad.Bayes.Class (MonadSample, mvNormal)
+import Control.Monad.Bayes.Sampler.Strict
+import Control.Monad.Identity (runIdentity)
+import Control.Monad.State (evalStateT)
+import Data.Matrix (fromList)
+import Data.Vector qualified as V
+import System.Random.MWC (toSeed)
+
+-- Test the sampled covariance is approximately the same as the
+-- specified covariance.
+passed1 :: IO Bool
+passed1 = sampleIOfixed $ do
+  let mu = (V.fromList [0.0, 0.0])
+      sigma11 = 2.0
+      sigma12 = 1.0
+      bigSigma = (fromList 2 2 [sigma11, sigma12, sigma12, sigma11])
+      nSamples = 200000
+      nSamples' = fromIntegral nSamples
+  ss <- replicateM nSamples $ (mvNormal mu bigSigma)
+  let xbar = (/ nSamples') $ sum $ fmap (V.! 0) ss
+      ybar = (/ nSamples') $ sum $ fmap (V.! 1) ss
+  let term1 = (/ nSamples') $ sum $ zipWith (*) (fmap (V.! 0) ss) (fmap (V.! 1) ss)
+  let term2 = xbar * ybar
+  return $ abs (sigma12 - (term1 - term2)) < 2e-2
+
+-- Test the sampled means are approximately the same as the specified
+-- means.
+passed2 :: IO Bool
+passed2 = sampleIOfixed $ do
+  let mu = (V.fromList [0.0, 0.0])
+      sigma11 = 2.0
+      sigma12 = 1.0
+      bigSigma = (fromList 2 2 [sigma11, sigma12, sigma12, sigma11])
+      nSamples = 100000
+      nSamples' = fromIntegral nSamples
+  ss <- replicateM nSamples $ (mvNormal mu bigSigma)
+  let xbar = (/ nSamples') $ sum $ fmap (V.! 0) ss
+      ybar = (/ nSamples') $ sum $ fmap (V.! 1) ss
+  return $ abs xbar < 1e-2 && abs ybar < 1e-2
+
+-- Test the sampled variances are approximately the same as the
+-- specified variances.
+passed3 :: IO Bool
+passed3 = sampleIOfixed $ do
+  let mu = (V.fromList [0.0, 0.0])
+      sigma11 = 2.0
+      sigma12 = 1.0
+      bigSigma = (fromList 2 2 [sigma11, sigma12, sigma12, sigma11])
+      nSamples = 200000
+      nSamples' = fromIntegral nSamples
+  ss <- replicateM nSamples $ (mvNormal mu bigSigma)
+  let xbar = (/ nSamples') $ sum $ fmap (V.! 0) ss
+      ybar = (/ nSamples') $ sum $ fmap (V.! 1) ss
+  let xbar2 = (/ nSamples') $ sum $ fmap (\x -> x * x) $ fmap (V.! 0) ss
+      ybar2 = (/ nSamples') $ sum $ fmap (\x -> x * x) $ fmap (V.! 1) ss
+  let xvar = xbar2 - xbar * xbar
+  let yvar = ybar2 - ybar * ybar
+  return $ abs (xvar - sigma11) < 1e-2 && abs (yvar - sigma11) < 2e-2
diff --git a/test/TestEnumerator.hs b/test/TestEnumerator.hs
--- a/test/TestEnumerator.hs
+++ b/test/TestEnumerator.hs
@@ -1,11 +1,19 @@
-module TestEnumerator where
+{-# LANGUAGE ImportQualifiedPost #-}
 
+module TestEnumerator (passed1, passed2, passed3, passed4) where
+
 import Control.Monad.Bayes.Class
+  ( MonadSample (categorical, uniformD),
+  )
 import Control.Monad.Bayes.Enumerator
-import Data.AEq
-import qualified Data.Vector as V
-import Numeric.Log
-import Sprinkler
+  ( enumerator,
+    evidence,
+    expectation,
+  )
+import Data.AEq (AEq ((~==)))
+import Data.Vector qualified as V
+import Numeric.Log (Log (ln))
+import Sprinkler (hard, soft)
 
 unnorm :: MonadSample m => m Int
 unnorm = categorical $ V.fromList [0.5, 0.8]
@@ -20,10 +28,10 @@
   return (x + y)
 
 passed2 :: Bool
-passed2 = enumerate agg ~== [(1, 0.25), (2, 0.5), (3, 0.25)]
+passed2 = enumerator agg ~== [(2, 0.5), (1, 0.25), (3, 0.25)]
 
 passed3 :: Bool
-passed3 = enumerate Sprinkler.hard ~== enumerate Sprinkler.soft
+passed3 = enumerator Sprinkler.hard ~== enumerator Sprinkler.soft
 
 passed4 :: Bool
 passed4 =
diff --git a/test/TestInference.hs b/test/TestInference.hs
--- a/test/TestInference.hs
+++ b/test/TestInference.hs
@@ -1,16 +1,34 @@
+{-# LANGUAGE ImportQualifiedPost #-}
 {-# LANGUAGE Rank2Types #-}
 {-# LANGUAGE TypeFamilies #-}
+{-# OPTIONS_GHC -Wno-missing-export-lists #-}
 
 module TestInference where
 
-import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Enumerator
+import ConjugatePriors
+  ( betaBernoulli',
+    betaBernoulliAnalytic,
+    gammaNormal',
+    gammaNormalAnalytic,
+    normalNormal',
+    normalNormalAnalytic,
+  )
+import Control.Monad (replicateM)
+import Control.Monad.Bayes.Class (MonadInfer, posterior)
+import Control.Monad.Bayes.Enumerator (enumerator)
 import Control.Monad.Bayes.Inference.SMC
+import Control.Monad.Bayes.Integrator (normalize)
+import Control.Monad.Bayes.Integrator qualified as Integrator
 import Control.Monad.Bayes.Population
-import Control.Monad.Bayes.Sampler
-import Data.AEq
-import Numeric.Log
-import Sprinkler
+import Control.Monad.Bayes.Population (collapse, runPopulation)
+import Control.Monad.Bayes.Sampler.Strict (Sampler, 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 Data.AEq (AEq ((~==)))
+import Numeric.Log (Log)
+import Sprinkler (soft)
+import System.Random.Stateful (IOGenM, StdGen, mkStdGen, newIOGenM)
 
 sprinkler :: MonadInfer m => m Bool
 sprinkler = Sprinkler.soft
@@ -18,16 +36,62 @@
 -- | Count the number of particles produced by SMC
 checkParticles :: Int -> Int -> IO Int
 checkParticles observations particles =
-  sampleIOfixed (fmap length (runPopulation $ smcMultinomial observations particles Sprinkler.soft))
+  sampleIOfixed (fmap length (population $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleMultinomial} Sprinkler.soft))
 
 checkParticlesSystematic :: Int -> Int -> IO Int
 checkParticlesSystematic observations particles =
-  sampleIOfixed (fmap length (runPopulation $ smcSystematic observations particles Sprinkler.soft))
+  sampleIOfixed (fmap length (population $ 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))
+
 checkTerminateSMC :: IO [(Bool, Log Double)]
-checkTerminateSMC = sampleIOfixed (runPopulation $ smcMultinomial 2 5 sprinkler)
+checkTerminateSMC = sampleIOfixed (population $ smc SMCConfig {numSteps = 2, numParticles = 5, resampler = resampleMultinomial} sprinkler)
 
 checkPreserveSMC :: Bool
 checkPreserveSMC =
-  (enumerate . collapse . smcMultinomial 2 2) sprinkler
-    ~== enumerate sprinkler
+  (enumerator . collapse . smc SMCConfig {numSteps = 2, numParticles = 2, resampler = resampleMultinomial}) sprinkler
+    ~== enumerator sprinkler
+
+expectationNearNumeric ::
+  Weighted Integrator.Integrator Double ->
+  Weighted Integrator.Integrator Double ->
+  Double
+expectationNearNumeric x y =
+  let e1 = Integrator.expectation $ normalize x
+      e2 = Integrator.expectation $ normalize y
+   in (abs (e1 - e2))
+
+expectationNearSampling ::
+  Weighted (Sampler (IOGenM StdGen) IO) Double ->
+  Weighted (Sampler (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
+  return (abs (e1 - e2))
+
+testNormalNormal :: [Double] -> IO Bool
+testNormalNormal n = do
+  let e =
+        expectationNearNumeric
+          (posterior (normalNormal' 1 (1, 10)) n)
+          (normalNormalAnalytic 1 (1, 10) n)
+  return (e < 1e-0)
+
+testGammaNormal :: [Double] -> IO Bool
+testGammaNormal n = do
+  let e =
+        expectationNearNumeric
+          (posterior (gammaNormal' (1, 1)) n)
+          (gammaNormalAnalytic (1, 1) n)
+  return (e < 1e-1)
+
+testBetaBernoulli :: [Bool] -> IO Bool
+testBetaBernoulli bs = do
+  let e =
+        expectationNearNumeric
+          (posterior (betaBernoulli' (1, 1)) bs)
+          (betaBernoulliAnalytic (1, 1) bs)
+  return (e < 1e-1)
diff --git a/test/TestIntegrator.hs b/test/TestIntegrator.hs
new file mode 100644
--- /dev/null
+++ b/test/TestIntegrator.hs
@@ -0,0 +1,151 @@
+{-# LANGUAGE BlockArguments #-}
+
+module TestIntegrator where
+
+import Control.Monad (replicateM)
+import Control.Monad.Bayes.Class
+  ( MonadCond (score),
+    MonadInfer,
+    MonadSample (bernoulli, gamma, normal, random, uniformD),
+    condition,
+    factor,
+    normalPdf,
+  )
+import Control.Monad.Bayes.Integrator
+import Control.Monad.Bayes.Sampler.Strict
+import Control.Monad.Bayes.Weighted (weighted)
+import Control.Monad.ST (ST, runST)
+import Data.AEq (AEq ((~==)))
+import Data.List (sortOn)
+import Data.Set (fromList)
+import Numeric.Log (Log (Exp, ln))
+import Sprinkler (hard, soft)
+import Statistics.Distribution (Distribution (cumulative))
+import Statistics.Distribution.Normal (normalDistr)
+
+normalExpectation :: Double -> Double -> Double
+normalExpectation mean std = expectation (normal mean std)
+
+normalVariance :: Double -> Double -> Double
+normalVariance mean std = variance (normal mean std)
+
+volumeIsOne :: [Double] -> Bool
+volumeIsOne = (~== 1.0) . volume . uniformD
+
+agg :: MonadSample m => m Int
+agg = do
+  x <- uniformD [0, 1]
+  y <- uniformD [2, 1]
+  return (x + y)
+
+within :: (Ord a, Num a) => a -> a -> a -> Bool
+within n x y = abs (x - y) < n
+
+passed1,
+  passed2,
+  passed3,
+  passed4,
+  passed5,
+  passed6,
+  passed7,
+  passed8,
+  passed9,
+  passed10,
+  passed11,
+  passed12,
+  passed13,
+  passed14 ::
+    Bool
+-- enumerator from Integrator works
+passed1 =
+  sortOn fst (enumeratorWith (fromList [3, 1, 2]) agg)
+    ~== sortOn fst [(2, 0.5), (1, 0.25), (3, 0.25)]
+-- hard and soft sprinkers are equivalent under enumerator from Integrator
+passed2 =
+  enumeratorWith (fromList [True, False]) (normalize (Sprinkler.hard))
+    ~== enumeratorWith (fromList [True, False]) (normalize (Sprinkler.soft))
+-- expectation is as expected
+passed3 =
+  expectation (fmap ((** 2) . (+ 1)) $ uniformD [0, 1]) == 2.5
+-- distribution is normalized
+passed4 = volume (uniformD [1, 2]) ~== 1.0
+-- enumerator is as expected
+passed5 =
+  sortOn fst (enumeratorWith (fromList [0, 1 :: Int]) (empirical [0 :: Int, 1, 1, 1]))
+    == sortOn fst [(1, 0.75), (0, 0.25)]
+-- normalization works right for enumerator, when there is conditioning
+passed6 =
+  sortOn fst [(2, 0.5), (3, 0.5), (1, 0.0)]
+    == sortOn
+      fst
+      ( enumeratorWith (fromList [1, 2, 3]) $
+          normalize $ do
+            x <- uniformD [1 :: Int, 2, 3]
+            condition (x > 1)
+            return x
+      )
+-- soft factor statements work with enumerator and normalization
+passed7 =
+  sortOn fst [(True, 0.75), (False, 0.25)]
+    ~== sortOn
+      fst
+      ( enumeratorWith (fromList [True, False]) $ normalize do
+          x <- bernoulli 0.5
+          factor $ if x then 0.3 else 0.1
+          return x
+      )
+-- volume of weight remains 1
+passed8 =
+  1
+    == ( volume $
+           fmap (ln . exp . snd) $ weighted do
+             x <- bernoulli 0.5
+             factor $ if x then 0.2 else 0.1
+             return x
+       )
+-- normal probability in positive region is half
+passed9 = probability (1, 1000) (normal 1 10) - 0.5 < 0.05
+-- cdf as expected
+passed10 = cdf (normal 5 5) 5 - 0.5 < 0.05
+-- cdf as expected
+passed11 =
+  (within 0.001)
+    ( cdf
+        ( do
+            x <- normal 0 1
+            return x
+        )
+        3
+    )
+    (cumulative (normalDistr 0 1) 3)
+-- volume as expected
+passed12 =
+  volume
+    ( do
+        x <- gamma 2 3
+        return x
+    )
+    ~== 1
+-- normalization preserves volume
+passed13 =
+  (volume . normalize)
+    ( do
+        x <- gamma 2 3
+        factor (normalPdf 0 1 x)
+        return x
+    )
+    ~== 1
+-- sampler and integrator agree on a non-trivial model
+passed14 =
+  let sample = runST $ sampleSTfixed $ fmap sampleMean $ replicateM 10000 $ weighted $ model1
+      quadrature = expectation $ normalize $ model1
+   in abs (sample - quadrature) < 0.01
+
+model1 :: MonadInfer m => m Double
+model1 = do
+  x <- random
+  y <- random
+  score (Exp $ log (f x + y))
+  return x
+  where
+    f x = cos (x ** 4) + x ** 3
diff --git a/test/TestPipes.hs b/test/TestPipes.hs
new file mode 100644
--- /dev/null
+++ b/test/TestPipes.hs
@@ -0,0 +1,21 @@
+{-# OPTIONS_GHC -Wno-monomorphism-restriction #-}
+
+module TestPipes where
+
+import BetaBin (urn, urnP)
+import Control.Monad.Bayes.Class ()
+import Control.Monad.Bayes.Enumerator (enumerator)
+import Data.AEq (AEq ((~==)))
+import HMM (hmm, hmmPosterior)
+import Pipes ((>->))
+import Pipes.Prelude (toListM)
+import qualified Pipes.Prelude as Pipes
+
+urns :: Int -> Bool
+urns n = enumerator (urn n) ~== enumerator (urnP n)
+
+hmms :: [Double] -> Bool
+hmms observations =
+  let hmmWithoutPipe = hmm observations
+      hmmWithPipe = reverse . init <$> toListM (hmmPosterior observations)
+   in enumerator hmmWithPipe ~== enumerator hmmWithoutPipe
diff --git a/test/TestPopulation.hs b/test/TestPopulation.hs
--- a/test/TestPopulation.hs
+++ b/test/TestPopulation.hs
@@ -1,42 +1,51 @@
-module TestPopulation where
+module TestPopulation (weightedSampleSize, popSize, manySize, sprinkler, sprinklerExact, transCheck1, transCheck2, resampleCheck, popAvgCheck) where
 
-import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Enumerator
+import Control.Monad.Bayes.Class (MonadInfer, MonadSample)
+import Control.Monad.Bayes.Enumerator (enumerator, expectation)
 import Control.Monad.Bayes.Population as Population
-import Control.Monad.Bayes.Sampler
-import Data.AEq
-import Sprinkler
+  ( Population,
+    collapse,
+    popAvg,
+    population,
+    pushEvidence,
+    resampleMultinomial,
+    spawn,
+  )
+import Control.Monad.Bayes.Sampler.Strict (sampleIOfixed)
+import Data.AEq (AEq ((~==)))
+import Sprinkler (soft)
+import System.Random.Stateful (mkStdGen, newIOGenM)
 
 weightedSampleSize :: MonadSample m => Population m a -> m Int
-weightedSampleSize = fmap length . runPopulation
+weightedSampleSize = fmap length . population
 
 popSize :: IO Int
-popSize = sampleIOfixed $ weightedSampleSize $ spawn 5 >> sprinkler
+popSize =
+  sampleIOfixed (weightedSampleSize $ spawn 5 >> sprinkler)
 
 manySize :: IO Int
-manySize = sampleIOfixed $ weightedSampleSize $ spawn 5 >> sprinkler >> spawn 3
+manySize =
+  sampleIOfixed (weightedSampleSize $ spawn 5 >> sprinkler >> spawn 3)
 
 sprinkler :: MonadInfer m => m Bool
 sprinkler = Sprinkler.soft
 
 sprinklerExact :: [(Bool, Double)]
-sprinklerExact = enumerate Sprinkler.soft
-
---all_check = (mass (Population.all id (spawn 2 >> sprinkler)) True) ~== 0.09
+sprinklerExact = enumerator Sprinkler.soft
 
 transCheck1 :: Bool
 transCheck1 =
-  enumerate (collapse sprinkler)
+  enumerator (collapse sprinkler)
     ~== sprinklerExact
 
 transCheck2 :: Bool
 transCheck2 =
-  enumerate (collapse (spawn 2 >> sprinkler))
+  enumerator (collapse (spawn 2 >> sprinkler))
     ~== sprinklerExact
 
 resampleCheck :: Int -> Bool
 resampleCheck n =
-  (enumerate . collapse . resampleMultinomial) (spawn n >> sprinkler)
+  (enumerator . collapse . resampleMultinomial) (spawn n >> sprinkler)
     ~== sprinklerExact
 
 popAvgCheck :: Bool
diff --git a/test/TestSampler.hs b/test/TestSampler.hs
new file mode 100644
--- /dev/null
+++ b/test/TestSampler.hs
@@ -0,0 +1,15 @@
+module TestSampler where
+
+import qualified Control.Foldl as Fold
+import Control.Monad (replicateM)
+import Control.Monad.Bayes.Class (MonadSample (normal))
+import Control.Monad.Bayes.Sampler.Strict (sampleSTfixed)
+import Control.Monad.ST (ST, runST)
+
+testMeanAndVariance :: Bool
+testMeanAndVariance = isDiff
+  where
+    m = runST (sampleSTfixed (foldWith Fold.mean (normal 2 4)))
+    v = runST (sampleSTfixed (foldWith Fold.variance (normal 2 4)))
+    foldWith f = fmap (Fold.fold f) . replicateM 100000
+    isDiff = abs (2 - m) < 0.01 && abs (16 - v) < 0.1
diff --git a/test/TestSequential.hs b/test/TestSequential.hs
--- a/test/TestSequential.hs
+++ b/test/TestSequential.hs
@@ -1,10 +1,14 @@
-module TestSequential where
+module TestSequential (twoSync, finishedTwoSync, checkTwoSync, checkPreserve, pFinished, isFinished, checkSync) where
 
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Enumerator as Dist
-import Control.Monad.Bayes.Sequential
-import Data.AEq
-import Sprinkler
+  ( MonadInfer,
+    MonadSample (uniformD),
+    factor,
+  )
+import Control.Monad.Bayes.Enumerator as Dist (enumerator, mass)
+import Control.Monad.Bayes.Sequential.Coroutine (advance, finish, finished)
+import Data.AEq (AEq ((~==)))
+import Sprinkler (soft)
 
 twoSync :: MonadInfer m => m Int
 twoSync = do
@@ -18,7 +22,7 @@
 finishedTwoSync n = finished (run n twoSync)
   where
     run 0 d = d
-    run k d = run (k -1) (advance d)
+    run k d = run (k - 1) (advance d)
 
 checkTwoSync :: Int -> Bool
 checkTwoSync 0 = mass (finishedTwoSync 0) False ~== 1
@@ -30,7 +34,7 @@
 sprinkler = Sprinkler.soft
 
 checkPreserve :: Bool
-checkPreserve = enumerate (finish sprinkler) ~== enumerate sprinkler
+checkPreserve = enumerator (finish sprinkler) ~== enumerator sprinkler
 
 pFinished :: Int -> Double
 pFinished 0 = 0.8267716535433071
@@ -42,7 +46,7 @@
 isFinished n = finished (run n sprinkler)
   where
     run 0 d = d
-    run k d = run (k -1) (advance d)
+    run k d = run (k - 1) (advance d)
 
 checkSync :: Int -> Bool
 checkSync n = mass (isFinished n) True ~== pFinished n
diff --git a/test/TestStormerVerlet.hs b/test/TestStormerVerlet.hs
new file mode 100644
--- /dev/null
+++ b/test/TestStormerVerlet.hs
@@ -0,0 +1,98 @@
+module TestStormerVerlet
+  ( passed1,
+  )
+where
+
+import Control.Lens
+import Control.Monad.ST
+import Data.Maybe (fromJust)
+import qualified Data.Vector as V
+import qualified Linear as L
+import Linear.V
+import Math.Integrators.StormerVerlet
+
+gConst :: Double
+gConst = 6.67384e-11
+
+nStepsTwoPlanets :: Int
+nStepsTwoPlanets = 44
+
+stepTwoPlanets :: Double
+stepTwoPlanets = 24 * 60 * 60 * 100
+
+sunMass, jupiterMass :: Double
+sunMass = 1.9889e30
+jupiterMass = 1.8986e27
+
+jupiterPerihelion :: Double
+jupiterPerihelion = 7.405736e11
+
+jupiterV :: [Double]
+jupiterV = [-1.0965244901087316e02, -1.3710001990210707e04, 0.0]
+
+jupiterQ :: [Double]
+jupiterQ = [negate jupiterPerihelion, 0.0, 0.0]
+
+sunV :: [Double]
+sunV = [0.0, 0.0, 0.0]
+
+sunQ :: [Double]
+sunQ = [0.0, 0.0, 0.0]
+
+tm :: V.Vector Double
+tm = V.enumFromStepN 0 stepTwoPlanets nStepsTwoPlanets
+
+keplerP :: L.V2 (L.V3 Double) -> L.V2 (L.V3 Double)
+keplerP (L.V2 p1 p2) = L.V2 dHdP1 dHdP2
+  where
+    dHdP1 = p1 / pure jupiterMass
+    dHdP2 = p2 / pure sunMass
+
+keplerQ :: L.V2 (L.V3 Double) -> L.V2 (L.V3 Double)
+keplerQ (L.V2 q1 q2) = L.V2 dHdQ1 dHdQ2
+  where
+    r = q2 L.^-^ q1
+    ri = r `L.dot` r
+    rr = ri * (sqrt ri)
+    q1' = pure gConst * r / pure rr
+    q2' = negate q1'
+    dHdQ1 = q1' * pure sunMass * pure jupiterMass
+    dHdQ2 = q2' * pure sunMass * pure jupiterMass
+
+listToV3 :: [a] -> L.V3 a
+listToV3 [x, y, z] = fromV . fromJust . fromVector . V.fromList $ [x, y, z]
+listToV3 xs = error $ "Only supply 3 elements not: " ++ show (length xs)
+
+initPQ2s :: L.V2 (L.V2 (L.V3 Double))
+initPQ2s =
+  L.V2
+    (L.V2 (listToV3 jupiterQ) (listToV3 sunQ))
+    (L.V2 (pure jupiterMass * listToV3 jupiterV) (pure sunMass * listToV3 sunV))
+
+result2 :: V.Vector (L.V2 (L.V2 (L.V3 Double)))
+result2 = runST $ integrateV (\h -> stormerVerlet2H (pure h) keplerQ keplerP) initPQ2s tm
+
+energy :: (L.V2 (L.V2 (L.V3 Double))) -> Double
+energy x = keJ + keS + peJ + peS
+  where
+    qs = x ^. _1
+    ps = x ^. _2
+    qJ = qs ^. _1
+    qS = qs ^. _2
+    pJ = ps ^. _1
+    pS = ps ^. _2
+    keJ = (* 0.5) $ (/ jupiterMass) $ sum $ fmap (^ 2) pJ
+    keS = (* 0.5) $ (/ sunMass) $ sum $ fmap (^ 2) pS
+    r = qJ L.^-^ qS
+    ri = r `L.dot` r
+    peJ = 0.5 * gConst * sunMass * jupiterMass / (sqrt ri)
+    peS = 0.5 * gConst * sunMass * jupiterMass / (sqrt ri)
+
+energies :: V.Vector Double
+energies = fmap energy result2
+
+diffs :: V.Vector Double
+diffs = V.zipWith (\x y -> abs (x - y) / x) energies (V.tail energies)
+
+passed1 :: IO Bool
+passed1 = return $ V.all (< 1.0e-3) diffs
diff --git a/test/TestWeighted.hs b/test/TestWeighted.hs
--- a/test/TestWeighted.hs
+++ b/test/TestWeighted.hs
@@ -1,14 +1,19 @@
 {-# LANGUAGE TypeFamilies #-}
 
-module TestWeighted where
+module TestWeighted (check, passed, result, model) where
 
 import Control.Monad.Bayes.Class
-import Control.Monad.Bayes.Sampler
-import Control.Monad.Bayes.Weighted
-import Control.Monad.State
-import Data.AEq
+  ( MonadInfer,
+    MonadSample (normal, uniformD),
+    factor,
+  )
+import Control.Monad.Bayes.Sampler.Strict (sampleIOfixed)
+import Control.Monad.Bayes.Weighted (weighted)
+import Control.Monad.State (unless, when)
+import Data.AEq (AEq ((~==)))
 import Data.Bifunctor (second)
-import Numeric.Log
+import Numeric.Log (Log (Exp, ln))
+import System.Random.Stateful (mkStdGen, newIOGenM)
 
 model :: MonadInfer m => m (Int, Double)
 model = do
@@ -19,7 +24,7 @@
   return (n, x)
 
 result :: MonadSample m => m ((Int, Double), Double)
-result = second (exp . ln) <$> runWeighted model
+result = second (exp . ln) <$> weighted model
 
 passed :: IO Bool
 passed = fmap check (sampleIOfixed result)