diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,3 @@
+# 0.1.0.0 (2020-02-17)
+
+Initial release.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,22 @@
+The MIT License (MIT)
+
+Copyright (c) 2015-2020 Adam Scibior
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/benchmark/SSM.hs b/benchmark/SSM.hs
new file mode 100644
--- /dev/null
+++ b/benchmark/SSM.hs
@@ -0,0 +1,35 @@
+module Main where
+
+import Control.Monad.IO.Class
+
+import Control.Monad.Bayes.Sampler
+import Control.Monad.Bayes.Weighted
+import Control.Monad.Bayes.Population
+import Control.Monad.Bayes.Inference.SMC
+import Control.Monad.Bayes.Inference.RMSMC
+import Control.Monad.Bayes.Inference.PMMH as PMMH
+import Control.Monad.Bayes.Inference.SMC2 as SMC2
+
+import NonlinearSSM
+
+main :: IO ()
+main = sampleIO $ do
+  let t = 5
+  dat <- generateData t
+  let ys = map snd dat
+
+  liftIO $ print "SMC"
+  smcRes <- runPopulation $ smcMultinomial t 10 (param >>= model ys)
+  liftIO $ print $ show smcRes
+
+  liftIO $ print "RM-SMC"
+  smcrmRes <- runPopulation $ rmsmcLocal t 10 10 (param >>= model ys)
+  liftIO $ print $ show smcrmRes
+
+  liftIO $ print "PMMH"
+  pmmhRes <- prior $ pmmh 2 t 3 param (model ys)
+  liftIO $ print $ show pmmhRes
+
+  liftIO $ print "SMC2"
+  smc2Res <- runPopulation $ smc2 t 3 2 1 param (model ys)
+  liftIO $ print $ show smc2Res
diff --git a/benchmark/Single.hs b/benchmark/Single.hs
new file mode 100644
--- /dev/null
+++ b/benchmark/Single.hs
@@ -0,0 +1,87 @@
+import System.Random.MWC (createSystemRandom, GenIO)
+import Data.Time
+import Options.Applicative
+import Data.Semigroup ((<>))
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Sampler
+import Control.Monad.Bayes.Weighted
+import Control.Monad.Bayes.Inference.SMC
+import Control.Monad.Bayes.Inference.RMSMC
+import Control.Monad.Bayes.Population
+import Control.Monad.Bayes.Sequential
+import Control.Monad.Bayes.Traced
+
+import qualified HMM
+import qualified LogReg
+import qualified LDA
+
+data Model = LR Int | HMM Int | LDA (Int, Int)
+  deriving(Show, Read)
+parseModel :: String -> Maybe Model
+parseModel s =
+  case s of
+    'L':'R':n -> Just $ LR (read n)
+    'H':'M':'M':n -> Just $ HMM (read n)
+    'L':'D':'A':n -> Just $ LDA (5, read n)
+    _ -> Nothing
+
+getModel :: MonadInfer m => Model ->  (Int, m String)
+getModel model = (size model, program model) where
+  size (LR n) = n
+  size (HMM n) = n
+  size (LDA (d,w)) = d*w
+  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
+
+data Alg = SMC | MH | RMSMC
+  deriving(Read,Show)
+
+runAlg :: Model -> Alg -> SamplerIO String
+runAlg model alg =
+  case alg of
+    SMC ->
+      let n = 100
+          (k, m) = getModel model
+      in show <$> (runPopulation $ smcSystematic k n m)
+    MH  ->
+      let t = 100
+          (_, m) = getModel model
+      in show <$> (prior $ mh t m)
+    RMSMC ->
+      let n = 10
+          t = 1
+          (k, m) = getModel model
+      in show <$> (runPopulation $ rmsmcBasic k n t m)
+
+
+infer :: Model -> Alg -> IO ()
+infer model alg = do
+  g <- createSystemRandom
+  x <- sampleIOwith (runAlg model alg) g
+  print x
+
+opts :: ParserInfo (Model, Alg)
+opts = flip info fullDesc $ liftA2 (,) model alg where
+  model = option (maybeReader parseModel)
+          ( long "model"
+          <> short 'm'
+          <> help "Model")
+  alg = option auto
+          ( long "alg"
+         <> short 'a'
+         <> help "Inference algorithm")
+
+main :: IO ()
+main = do
+  (model, alg) <- execParser opts
+
+  startTime <- getCurrentTime
+
+  infer model alg
+
+  endTime <- getCurrentTime
+  print (diffUTCTime endTime startTime)
diff --git a/benchmark/Speed.hs b/benchmark/Speed.hs
new file mode 100644
--- /dev/null
+++ b/benchmark/Speed.hs
@@ -0,0 +1,253 @@
+import Criterion.Main
+import Criterion.Types
+import System.Random.MWC (createSystemRandom, GenIO)
+import Control.Monad (replicateM)
+
+import Debug.Trace
+
+import GHC.IO.Handle
+import System.Exit
+import System.Process hiding (env)
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Sampler
+import Control.Monad.Bayes.Weighted
+import Control.Monad.Bayes.Inference.SMC
+import Control.Monad.Bayes.Inference.RMSMC
+import Control.Monad.Bayes.Population
+import Control.Monad.Bayes.Traced
+
+-- import NonlinearSSM
+import qualified HMM
+import qualified LogReg
+import qualified LDA
+
+-- | 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 lein model = do
+  anglican lein ["(use 'nstools.ns)\n"]
+  anglican lein ["(ns+ " ++ anglicanModelName model ++ ")\n"]
+  case model of
+    LR dataset -> do
+      let (xs, labels) = unzip dataset
+      anglican lein ["(ns-unmap *ns* 'xs)\n"]
+      anglican lein ["(def xs " ++ clojureShowVector xs ++ ")\n", "nil\n"]
+      anglican lein ["(ns-unmap *ns* 'labels)\n"]
+      anglican lein ["(def labels " ++ clojureBoolVector labels ++ ")\n", "nil\n"]
+    HMM observations -> do
+      anglican lein ["(ns-unmap *ns* 'observations)\n"]
+      anglican lein ["(def observations " ++ clojureShowVector observations ++ ")\n", "nil\n"]
+    LDA docs -> do
+      anglican lein ["(ns-unmap *ns* 'docs)\n"]
+      anglican lein ["(def docs " ++ clojureVector (map clojureShowVector docs) ++ ")\n", "nil\n"]
+  anglican lein ["(use 'nstools.ns)\n"]
+  anglican lein ["(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"
+    if l == "nil" then
+      return ()
+    else
+      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 lein = LeinProc input output process
+  anglican lein ["(use 'nstools.ns)\n"]
+  return lein
+
+
+-- | 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)
+
+webpplModelName :: Model -> String
+webpplModelName (LR _) = "logisticRegression"
+webpplModelName (HMM _) = "hmm"
+webpplModelName (LDA _) = "lda"
+
+-- | Environment to execute benchmarks in.
+data Env = Env {rng :: GenIO, lein :: LeinProc}
+
+data ProbProgSys = MonadBayes | Anglican | WebPPL
+  deriving(Show)
+
+data Model = LR [(Double,Bool)] | HMM [Double] | LDA [[String]]
+instance Show Model where
+  show (LR xs) = "LR" ++ show (length xs)
+  show (HMM xs) = "HMM" ++ show (length xs)
+  show (LDA xs) = "LDA" ++ show (length $ head xs)
+buildModel :: MonadInfer m => Model -> m String
+buildModel (LR dataset) = show <$> LogReg.logisticRegression dataset
+buildModel (HMM dataset) = show <$> HMM.hmm dataset
+buildModel (LDA dataset) = show <$> LDA.lda dataset
+modelLength :: Model -> Int
+modelLength (LR xs) = length xs
+modelLength (HMM xs) = length xs
+modelLength (LDA xs) = sum (map length xs)
+
+data Alg = MH Int | SMC Int | RMSMC Int Int
+instance Show Alg where
+  show (MH n) = "MH" ++ show n
+  show (SMC n) = "SMC" ++ show n
+  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))
+
+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"
+
+prepareBenchmark :: Env -> ProbProgSys -> Model -> Alg -> Benchmark
+prepareBenchmark e 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 = do
+    anglican (lein e) $ ["(time (m! " ++ anglicanModelName model ++ " " ++ algString alg ++ "))\n"]
+prepareBenchmark _ WebPPL model alg = bench name $ whnfIO run 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 else 0) labels) ++ "'"
+  dataString (HMM obs) = "--obs='" ++ javascriptList obs ++ "'"
+  dataString (LDA docs) = unwords $ map (\(i,doc) -> "--doc" ++ show i ++ "='" ++ unwords doc ++ "'") (zip [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
+
+
+-- | Checks if the requested benchmark is implemented.
+supported :: (ProbProgSys, Model, Alg) -> Bool
+supported (Anglican, _, RMSMC _ _) = False
+supported _ = True
+
+systems = [
+            MonadBayes
+            -- Anglican,
+            -- WebPPL
+          ]
+
+lengthBenchmarks e lrData hmmData ldaData = benchmarks where
+  lrLengths = [10] ++ map (*100) [1..10]
+  hmmLengths = [10] ++ map (*100) [1..10]
+  ldaLengths = [5] ++ map (*50) [1..10]
+  models =
+    map (LR . (`take` lrData)) lrLengths ++
+    map (HMM . (`take` hmmData)) hmmLengths ++
+    map (\n -> LDA $ map (take n) ldaData) ldaLengths
+  algs = [
+    MH 100,
+    SMC 100,
+    RMSMC 10 1
+    ]
+  benchmarks = map (uncurry3 (prepareBenchmark e)) $ filter supported xs where
+        uncurry3 f (x,y,z) = f x y z
+        xs = do
+          m <- models
+          s <- systems
+          a <- algs
+          return (s,m,a)
+
+samplesBenchmarks e lrData hmmData ldaData = benchmarks where
+  lrLengths = [50]
+  hmmLengths = [20]
+  ldaLengths = [10]
+  models = map (LR . (`take` lrData)) lrLengths ++
+           map (HMM . (`take` hmmData)) hmmLengths ++
+           map (\n -> LDA $ map (take n) ldaData) ldaLengths
+  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 where
+        uncurry3 f (x,y,z) = f x y z
+        xs = do
+          a <- algs
+          s <- systems
+          m <- models
+          return (s,m,a)
+
+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
+
+  let configLength = defaultConfig{csvFile = Just "speed-length.csv", rawDataFile = Just "raw.dat"}
+  defaultMainWith configLength (lengthBenchmarks e lrData hmmData ldaData)
+
+  let configSamples = defaultConfig{csvFile = Just "speed-samples.csv", rawDataFile = Just "raw.dat"}
+  defaultMainWith configSamples (samplesBenchmarks e lrData hmmData ldaData)
diff --git a/models/Dice.hs b/models/Dice.hs
new file mode 100644
--- /dev/null
+++ b/models/Dice.hs
@@ -0,0 +1,32 @@
+
+
+module Dice 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.Monad.Bayes.Class
+
+-- | A toss of a six-sided die.
+die :: MonadSample m => m Int
+die = uniformD [1..6]
+
+-- | 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))
+
+-- | Toss of two dice where the output is greater than 4.
+dice_hard :: MonadInfer m => m Int
+dice_hard = do
+  result <- dice 2
+  condition (result > 4)
+  return result
+
+-- | Toss of two dice with an artificial soft constraint.
+dice_soft :: MonadInfer m => m Int
+dice_soft = do
+  result <- dice 2
+  score (1 / fromIntegral result)
+  return result
diff --git a/models/HMM.hs b/models/HMM.hs
new file mode 100644
--- /dev/null
+++ b/models/HMM.hs
@@ -0,0 +1,55 @@
+-- HMM from Anglican (https://bitbucket.org/probprog/anglican-white-paper)
+
+{-# LANGUAGE
+ FlexibleContexts,
+ TypeFamilies
+ #-}
+
+module HMM (
+  values,
+  hmm,
+  syntheticData
+  ) where
+
+--Hidden Markov Models
+
+import Control.Monad (replicateM)
+import Data.Vector (fromList)
+
+import Control.Monad.Bayes.Class
+
+-- | Observed values
+values :: [Double]
+values = [0.9,0.8,0.7,0,-0.025,-5,-2,-0.1,0,
+          0.13,0.45,6,0.2,0.3,-1,-1]
+
+-- | The transition model.
+trans :: MonadSample m => Int -> m Int
+trans 0 = categorical $ fromList [0.1, 0.4, 0.5]
+trans 1 = categorical $ fromList [0.2, 0.6, 0.2]
+trans 2 = categorical $ fromList [0.15,0.7,0.15]
+
+-- | The emission model.
+emissionMean :: Int -> Double
+emissionMean x = mean x where
+  mean 0 = -1
+  mean 1 = 1
+  mean 2 = 0
+
+-- | Initial state distribution
+start :: MonadSample m => m Int
+start = uniformD [0,1,2]
+
+-- | Example HMM from http://dl.acm.org/citation.cfm?id=2804317
+hmm :: (MonadInfer m) => [Double] -> m [Int]
+hmm dataset = f dataset (const . return) where
+  expand x y = do
+    x' <- trans x
+    factor $ normalPdf (emissionMean x') 1 y
+    return x'
+  f [] k = start >>= k []
+  f (y:ys) k = f ys (\xs x -> expand x y >>= k (x:xs))
+
+syntheticData :: MonadSample m => Int -> m [Double]
+syntheticData n = replicateM n syntheticPoint where
+    syntheticPoint = uniformD [0,1,2]
diff --git a/models/LDA.hs b/models/LDA.hs
new file mode 100644
--- /dev/null
+++ b/models/LDA.hs
@@ -0,0 +1,56 @@
+-- LDA model from Anglican
+-- (https://bitbucket.org/probprog/anglican-white-paper)
+
+module LDA where
+
+import Numeric.Log
+import qualified Control.Monad as List (replicateM)
+import Data.Vector as V hiding (length, zip, mapM, mapM_)
+import qualified Data.Map as Map
+
+import Control.Monad.Bayes.Class
+
+vocabluary :: [String]
+vocabluary = ["bear", "wolf", "python", "prolog"]
+topics :: [String]
+topics = ["topic1", "topic2"]
+
+docs :: [[String]]
+docs = [
+  words "bear wolf bear wolf bear wolf python wolf bear wolf",
+  words "python prolog python prolog python prolog python prolog python prolog",
+  words "bear wolf bear wolf bear wolf bear wolf bear wolf",
+  words "python prolog python prolog python prolog python prolog python prolog",
+  words "bear wolf bear python bear wolf bear wolf bear wolf"
+  ]
+
+word_dist_prior :: MonadSample m => m (Vector Double)
+word_dist_prior = dirichlet $ V.replicate (length vocabluary) 1
+
+topic_dist_prior :: MonadSample m => m (Vector Double)
+topic_dist_prior = dirichlet $ V.replicate (length topics) 1
+
+word_index :: Map.Map String Int
+word_index = Map.fromList $ zip vocabluary [0..]
+
+lda :: MonadInfer m => [[String]] -> m [Int]
+lda docs = do
+  word_dist_for_topic <- do
+    ts <- mapM (const word_dist_prior) [0 .. length topics]
+    return $ Map.fromList $ zip [0 .. length topics] ts
+
+  let obs doc = do
+        topic_dist <- fmap categorical topic_dist_prior
+        let f word = do
+              topic <- topic_dist
+              factor $ (Exp . log) $ (word_dist_for_topic Map.! topic) V.! (word_index Map.! word)
+        mapM_ f doc
+
+  mapM_ obs docs
+
+  -- return samples since Discrete is not NFData
+  mapM (categorical . snd) $ Map.toList word_dist_for_topic
+
+syntheticData :: MonadSample m => Int -> Int -> m [[String]]
+syntheticData d w = List.replicateM d (List.replicateM w syntheticWord) where
+  syntheticWord = uniformD vocabluary
diff --git a/models/LogReg.hs b/models/LogReg.hs
new file mode 100644
--- /dev/null
+++ b/models/LogReg.hs
@@ -0,0 +1,35 @@
+-- Logistic regression model from Anglican
+-- (https://bitbucket.org/probprog/anglican-white-paper)
+
+module LogReg where
+
+import Control.Monad (replicateM)
+
+import Numeric.Log
+import Control.Monad.Bayes.Class
+
+xs :: [Double]
+xs = [-10, -5, 2, 6, 10]
+
+labels :: [Bool]
+labels = [False, False, True, True, True]
+
+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))
+      obs x label = do
+        p <- sigmoid x
+        factor $ (Exp . log) $ if label then p else 1 - p
+  mapM_ (uncurry obs) dat
+  sigmoid 8
+
+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)
diff --git a/models/NonlinearSSM.hs b/models/NonlinearSSM.hs
new file mode 100644
--- /dev/null
+++ b/models/NonlinearSSM.hs
@@ -0,0 +1,53 @@
+module NonlinearSSM where
+
+import Control.Monad.Bayes.Class
+
+param :: MonadSample m => m (Double, Double)
+param = do
+  let a = 0.01
+  let b = 0.01
+  precX <- gamma a b
+  let sigmaX = 1 / sqrt precX
+  precY <- gamma a b
+  let sigmaY = 1 / sqrt precY
+  return (sigmaX, sigmaY)
+
+-- | A nonlinear series model from Doucet et al. (2000)
+-- "On sequential Monte Carlo sampling methods" section VI.B
+model :: (MonadInfer m)
+      => [Double]  -- ^ observed data
+      -> (Double, Double) -- ^ prior on the parameters
+      -> m [Double] -- ^ list of latent states from t=1
+model obs (sigmaX, sigmaY) = do
+  let sq x = x * x
+      simulate [] _ acc = return acc
+      simulate (y:ys) x acc = do
+        let n = length acc
+        let mean = 0.5 * x + 25 * x / (1 + sq x) +
+                   8 * cos (1.2 * fromIntegral n)
+        x' <- normal mean sigmaX
+        factor $ normalPdf (sq x' / 20) sigmaY y
+        simulate ys x' (x':acc)
+
+  x0 <- normal 0 (sqrt 5)
+  xs <- simulate obs x0 []
+  return $ reverse xs
+
+generateData :: MonadSample m
+      => Int  -- ^ T
+      -> m [(Double,Double)] -- ^ list of latent and observable states from t=1
+generateData t = do
+  (sigmaX, sigmaY) <- param
+  let sq x = x * x
+      simulate 0 _ acc = return acc
+      simulate k x acc = do
+        let n = length acc
+        let mean = 0.5 * x + 25 * x / (1 + sq x) +
+                   8 * cos (1.2 * fromIntegral n)
+        x' <- normal mean sigmaX
+        y' <- normal (sq x' / 20) sigmaY
+        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
new file mode 100644
--- /dev/null
+++ b/models/Sprinkler.hs
@@ -0,0 +1,24 @@
+module Sprinkler where
+
+import Control.Monad (when)
+
+import Control.Monad.Bayes.Class
+
+hard :: MonadInfer m => m Bool
+hard = do
+  rain <- bernoulli 0.3
+  sprinkler <- bernoulli $ if rain then 0.1 else 0.4
+  wet <- bernoulli $ case (rain,sprinkler) of (True,True) -> 0.98
+                                              (True,False) -> 0.8
+                                              (False,True) -> 0.9
+                                              (False,False) -> 0.0
+  condition (not wet)
+  return rain
+
+soft :: MonadInfer m => m Bool
+soft = do
+  rain <- bernoulli 0.3
+  when rain (factor 0.2)
+  sprinkler <- bernoulli $ if rain then 0.1 else 0.4
+  when sprinkler (factor 0.1)
+  return rain
diff --git a/monad-bayes.cabal b/monad-bayes.cabal
new file mode 100644
--- /dev/null
+++ b/monad-bayes.cabal
@@ -0,0 +1,141 @@
+cabal-version:       2.0
+name:                monad-bayes
+version:             0.1.0.0
+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.
+homepage:            http://github.com/tweag/monad-bayes#readme
+license:             MIT
+license-file:        LICENSE
+author:              Adam Scibior <adscib@gmail.com>
+maintainer:          leonhard.markert@tweag.io
+copyright:           2015-2020 Adam Scibior
+bug-reports:         https://github.com/tweag/monad-bayes/issues
+stability:           experimental
+category:            Statistics
+build-type:          Simple
+
+extra-source-files:
+  CHANGELOG.md
+
+executable example
+  hs-source-dirs:      benchmark, models
+  main-is:             Single.hs
+  build-depends:       base
+                     , monad-bayes
+                     , log-domain
+                     , vector
+                     , containers
+                     , mwc-random
+                     , time
+                     , optparse-applicative
+  default-language:    Haskell2010
+  other-modules:       LogReg
+                     , HMM
+                     , LDA
+                     , Dice
+
+library
+  hs-source-dirs:      src
+  exposed-modules:     Control.Monad.Bayes.Class
+                     , Control.Monad.Bayes.Free
+                     , Control.Monad.Bayes.Sampler
+                     , Control.Monad.Bayes.Weighted
+                     , Control.Monad.Bayes.Sequential
+                     , Control.Monad.Bayes.Population
+                     , Control.Monad.Bayes.Traced.Static
+                     , Control.Monad.Bayes.Traced.Dynamic
+                     , Control.Monad.Bayes.Traced.Basic
+                     , Control.Monad.Bayes.Traced
+                     , Control.Monad.Bayes.Enumerator
+                     , Control.Monad.Bayes.Helpers
+                     , Control.Monad.Bayes.Inference.SMC
+                     , Control.Monad.Bayes.Inference.RMSMC
+                     , Control.Monad.Bayes.Inference.PMMH
+                     , Control.Monad.Bayes.Inference.SMC2
+  other-modules:       Control.Monad.Bayes.Traced.Common
+  -- See MAINTAINERS.md for when and how to update the version bounds.
+  build-depends:       base >= 4.10.1 && < 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
+  ghc-options:         -Wall -fno-warn-redundant-constraints
+  default-language:    Haskell2010
+  default-extensions:  RankNTypes
+                     , GeneralizedNewtypeDeriving
+                     , StandaloneDeriving
+                     , TypeFamilies
+                     , FlexibleContexts
+                     , FlexibleInstances
+                     , TupleSections
+                     , MultiParamTypeClasses
+                     , GADTs
+  other-extensions:    ScopedTypeVariables
+                     , DeriveFunctor
+
+test-suite monad-bayes-test
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test, models
+  main-is:             Spec.hs
+  build-depends:       base
+                     , monad-bayes
+                     , hspec
+                     , QuickCheck
+                     , ieee754
+                     , mtl
+                     , math-functions
+                     , transformers
+                     , log-domain
+                     , vector
+  ghc-options:         -Wall -threaded -rtsopts "-with-rtsopts=-N -M1g -K1g"
+  default-language:    Haskell2010
+  other-modules:       Sprinkler,
+                       TestEnumerator,
+                       TestPopulation,
+                       TestInference,
+                       TestSequential,
+                       TestWeighted
+
+benchmark ssm-bench
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      models
+                     , benchmark
+  build-depends:       base
+                     , monad-bayes
+  default-language:    Haskell2010
+  default-extensions:  RankNTypes
+  main-is:             SSM.hs
+  other-modules:       NonlinearSSM
+
+benchmark speed-bench
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      models
+                     , benchmark
+  build-depends:       base
+                     , monad-bayes
+                     , criterion
+                     , abstract-par
+                     , log-domain
+                     , vector
+                     , mwc-random
+                     , containers
+                     , process
+  ghc-options:         -rtsopts "-with-rtsopts=-M1g -K1g"
+  default-language:    Haskell2010
+  default-extensions:  RankNTypes
+  main-is:             Speed.hs
+  other-modules:       LogReg
+                     , HMM
+                     , LDA
+
+source-repository head
+  type:     git
+  location: https://github.com/tweag/monad-bayes.git
diff --git a/src/Control/Monad/Bayes/Class.hs b/src/Control/Monad/Bayes/Class.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Class.hs
@@ -0,0 +1,296 @@
+{-|
+Module      : Control.Monad.Bayes.Class
+Description : Types for probabilistic modelling
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+This module defines 'MonadInfer', which can be used to represent a simple model
+like the following:
+
+@
+import Control.Monad (when)
+import Control.Monad.Bayes.Class
+
+model :: MonadInfer m => m Bool
+model = do
+  rain <- bernoulli 0.3
+  sprinkler <-
+    bernoulli $
+    if rain
+      then 0.1
+      else 0.4
+  let wetProb =
+    case (rain, sprinkler) of
+      (True,  True)  -> 0.98
+      (True,  False) -> 0.80
+      (False, True)  -> 0.90
+      (False, False) -> 0.00
+  score wetProb
+  return rain
+@
+-}
+
+module Control.Monad.Bayes.Class (
+  MonadSample,
+  random,
+  uniform,
+  normal,
+  gamma,
+  beta,
+  bernoulli,
+  categorical,
+  logCategorical,
+  uniformD,
+  geometric,
+  poisson,
+  dirichlet,
+  MonadCond,
+  score,
+  factor,
+  condition,
+  MonadInfer,
+  normalPdf
+) where
+
+import Control.Monad.Trans.Class
+import Control.Monad.Trans.Identity
+import Control.Monad.Trans.Maybe
+import Control.Monad.Trans.State
+import Control.Monad.Trans.Writer
+import Control.Monad.Trans.Reader
+import Control.Monad.Trans.RWS hiding (tell)
+import Control.Monad.Trans.List
+import Control.Monad.Trans.Cont
+
+import Statistics.Distribution
+import Statistics.Distribution.Uniform (uniformDistr)
+import Statistics.Distribution.Normal (normalDistr)
+import Statistics.Distribution.Gamma (gammaDistr)
+import Statistics.Distribution.Beta (betaDistr)
+import Statistics.Distribution.Geometric (geometric0)
+import qualified Statistics.Distribution.Poisson as Poisson
+
+import Numeric.Log
+
+import Data.Vector.Generic as VG
+import qualified Data.Vector as V
+import Control.Monad (when)
+
+-- | Monads that can draw random variables.
+class Monad m => MonadSample m where
+  -- | Draw from a uniform distribution.
+  random :: m Double -- ^ \(\sim \mathcal{U}(0, 1)\)
+
+  -- | Draw from a uniform distribution.
+  uniform ::
+       Double -- ^ lower bound a
+    -> Double -- ^ upper bound b
+    -> m Double -- ^ \(\sim \mathcal{U}(a, b)\).
+  uniform a b = draw (uniformDistr a b)
+
+  -- | Draw from a normal distribution.
+  normal ::
+       Double -- ^ mean μ
+    -> Double -- ^ standard deviation σ
+    -> m Double -- ^ \(\sim \mathcal{N}(\mu, \sigma^2)\)
+  normal m s = draw (normalDistr m s)
+
+  -- | Draw from a gamma distribution.
+  gamma ::
+       Double -- ^ shape k
+    -> Double -- ^ scale θ
+    -> m Double -- ^ \(\sim \Gamma(k, \theta)\)
+  gamma shape scale = draw (gammaDistr shape scale)
+
+  -- | Draw from a beta distribution.
+  beta ::
+       Double -- ^ shape α
+    -> Double -- ^ shape β
+    -> m Double -- ^ \(\sim \mathrm{Beta}(\alpha, \beta)\)
+  beta a b = draw (betaDistr a b)
+
+  -- | Draw from a Bernoulli distribution.
+  bernoulli ::
+       Double -- ^ probability p
+    -> m Bool -- ^ \(\sim \mathrm{B}(1, p)\)
+  bernoulli p = fmap (< p) random
+
+  -- | Draw from a categorical distribution.
+  categorical ::
+       Vector v Double
+    => v Double -- ^ event probabilities
+    -> m Int -- ^ outcome category
+  categorical ps = fromPMF (ps !)
+
+  -- | Draw from a categorical distribution in the log domain.
+  logCategorical ::
+       (Vector v (Log Double), Vector v Double)
+    => v (Log Double) -- ^ event probabilities
+    -> m Int -- ^ outcome category
+  logCategorical = categorical . VG.map (exp . ln)
+
+  -- | Draw from a discrete uniform distribution.
+  uniformD ::
+       [a] -- ^ observable outcomes @xs@
+    -> m a -- ^ \(\sim \mathcal{U}\{\mathrm{xs}\}\)
+  uniformD xs = do
+    let n = Prelude.length xs
+    i <- categorical $ V.replicate n (1 / fromIntegral n)
+    return (xs !! i)
+
+  -- | Draw from a geometric distribution.
+  geometric ::
+       Double -- ^ success rate p
+    -> m Int -- ^ \(\sim\) number of failed Bernoulli trials with success probability p before first success
+  geometric = discrete . geometric0
+
+  -- | Draw from a Poisson distribution.
+  poisson ::
+       Double -- ^ parameter λ
+    -> m Int -- ^ \(\sim \mathrm{Pois}(\lambda)\)
+  poisson = discrete . Poisson.poisson
+
+  -- | Draw from a Dirichlet distribution.
+  dirichlet ::
+       Vector v Double
+    => v Double -- ^ concentration parameters @as@
+    -> m (v Double) -- ^ \(\sim \mathrm{Dir}(\mathrm{as})\)
+  dirichlet as = do
+    xs <- VG.mapM (`gamma` 1) as
+    let s = VG.sum xs
+    let ys = VG.map (/ s) xs
+    return ys
+
+-- | Draw from a continuous distribution using the inverse cumulative density
+-- function.
+draw :: (ContDistr d, MonadSample m) => d -> m Double
+draw d = fmap (quantile d) random
+
+-- | Draw from a discrete distribution using a sequence of draws from
+-- Bernoulli.
+fromPMF :: MonadSample m => (Int -> Double) -> m Int
+fromPMF p = f 0 1 where
+  f i r = do
+    when (r < 0) $ error "fromPMF: total PMF above 1"
+    let q = p i
+    when (q < 0 || q > 1) $ error "fromPMF: invalid probability value"
+    b <- bernoulli (q / r)
+    if b then pure i else f (i+1) (r-q)
+
+-- | Draw from a discrete distributions using the probability mass function.
+discrete :: (DiscreteDistr d, MonadSample m) => d -> m Int
+discrete = fromPMF . probability
+
+-- | Monads that can score different execution paths.
+class Monad m => MonadCond m where
+  -- | Record a likelihood.
+  score ::
+       Log Double -- ^ likelihood of the execution path
+    -> m ()
+
+-- | Synonym for 'score'.
+factor ::
+     MonadCond m
+  => Log Double -- ^ likelihood of the execution path
+  -> m ()
+factor = score
+
+-- | Hard conditioning.
+condition :: MonadCond m => Bool -> m ()
+condition b = score $ if b then 1 else 0
+
+-- | Monads that support both sampling and scoring.
+class (MonadSample m, MonadCond m) => MonadInfer m
+
+-- | Probability density function of the normal distribution.
+normalPdf ::
+     Double -- ^ mean μ
+  -> Double -- ^ standard deviation σ
+  -> Double -- ^ sample x
+  -> Log Double -- ^ relative likelihood of observing sample x in \(\mathcal{N}(\mu, \sigma^2)\)
+normalPdf mu sigma x = Exp $ logDensity (normalDistr mu sigma) x
+
+----------------------------------------------------------------------------
+-- Instances that lift probabilistic effects to standard tranformers.
+
+instance MonadSample m => MonadSample (IdentityT m) where
+  random = lift random
+  bernoulli = lift . bernoulli
+
+instance MonadCond m => MonadCond (IdentityT m) where
+  score = lift . score
+
+instance MonadInfer m => MonadInfer (IdentityT m)
+
+
+instance MonadSample m => MonadSample (MaybeT m) where
+  random = lift random
+
+instance MonadCond m => MonadCond (MaybeT m) where
+  score = lift . score
+
+instance MonadInfer m => MonadInfer (MaybeT m)
+
+
+instance MonadSample m => MonadSample (ReaderT r m) where
+  random = lift random
+  bernoulli = lift . bernoulli
+
+instance MonadCond m => MonadCond (ReaderT r m) where
+  score = lift . score
+
+instance MonadInfer m => MonadInfer (ReaderT r m)
+
+
+instance (Monoid w, MonadSample m) => MonadSample (WriterT w m) where
+  random = lift random
+  bernoulli = lift . bernoulli
+  categorical = lift . categorical
+
+instance (Monoid w, MonadCond m) => MonadCond (WriterT w m) where
+  score = lift . score
+
+instance (Monoid w, MonadInfer m) => MonadInfer (WriterT w m)
+
+
+instance MonadSample m => MonadSample (StateT s m) where
+  random = lift random
+  bernoulli = lift . bernoulli
+  categorical = lift . categorical
+
+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
+  bernoulli = lift . bernoulli
+  categorical = lift . categorical
+
+instance MonadCond m => MonadCond (ListT m) where
+  score = lift . score
+
+instance MonadInfer m => MonadInfer (ListT m)
+
+
+instance MonadSample m => MonadSample (ContT r m) where
+  random = lift random
+
+instance MonadCond m => MonadCond (ContT r m) where
+  score = lift . score
+
+instance MonadInfer m => MonadInfer (ContT r m)
diff --git a/src/Control/Monad/Bayes/Enumerator.hs b/src/Control/Monad/Bayes/Enumerator.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Enumerator.hs
@@ -0,0 +1,116 @@
+{-|
+Module      : Control.Monad.Bayes.Enumerator
+Description : Exhaustive enumeration of discrete random variables
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+-}
+
+module Control.Monad.Bayes.Enumerator (
+    Enumerator,
+    logExplicit,
+    explicit,
+    evidence,
+    mass,
+    compact,
+    enumerate,
+    expectation,
+    normalForm
+            ) where
+
+import Data.AEq (AEq, (===), (~==))
+import Control.Applicative (Applicative, Alternative)
+import Control.Monad (MonadPlus)
+import Control.Arrow (second)
+import qualified Data.Map as Map
+import qualified Data.Vector.Generic as V
+import Numeric.Log as Log
+import Control.Monad.Trans.Writer
+import Data.Monoid
+import Data.Maybe
+
+import Control.Monad.Bayes.Class
+
+
+-- | 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)
+
+instance MonadSample Enumerator where
+  random = error "Infinitely supported random variables not supported in Enumerator"
+  bernoulli p = fromList [(True, (Exp . log) p), (False, (Exp . log) (1-p))]
+  categorical v = fromList $ zip [0..] $ map (Exp . log) (V.toList v)
+
+instance MonadCond Enumerator where
+  score w = fromList [((), w)]
+
+instance MonadInfer Enumerator
+
+-- | Construct Enumerator from a list of values and associated weights.
+fromList :: [(a, Log Double)] -> Enumerator a
+fromList = Enumerator . WriterT . map (second Product)
+
+-- | Returns the posterior as a list of weight-value pairs without any post-processing,
+-- such as normalization or aggregation
+logExplicit :: Enumerator a -> [(a, Log Double)]
+logExplicit (Enumerator m) = map (second getProduct) $ runWriterT m
+
+-- | Same as `toList`, only weights are converted from log-domain.
+explicit :: Enumerator a -> [(a,Double)]
+explicit = map (second (exp . ln)) . logExplicit
+
+-- | Returns the model evidence, that is sum of all weights.
+evidence :: Enumerator a -> Log Double
+evidence = Log.sum . map snd . logExplicit
+
+-- | Normalized probability mass of a specific value.
+mass :: Ord a => Enumerator a -> a -> Double
+mass d = f where
+  f a = fromMaybe 0 $ lookup a m
+  m = enumerate 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 (+)
+
+-- | 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) where
+  (xs, ws) = second (map (exp . ln) . normalize) $ unzip (logExplicit d)
+
+-- | 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 xs = map (/ z) xs where
+  z = Log.sum xs
+
+-- | Divide all weights by their sum.
+normalizeWeights :: [(a, Log Double)] -> [(a, Log Double)]
+normalizeWeights ls = zip xs ps where
+  (xs, ws) = unzip ls
+  ps = normalize ws
+
+-- | 'compact' followed by removing values with zero weight.
+normalForm :: Ord a => Enumerator a -> [(a, Double)]
+normalForm = filter ((/= 0) . snd) . compact . explicit
+
+instance Ord a => Eq (Enumerator a) where
+  p == q = normalForm p == normalForm q
+
+instance Ord a => AEq (Enumerator a) where
+  p === q = xs == ys && ps === qs where
+    (xs,ps) = unzip (normalForm p)
+    (ys,qs) = unzip (normalForm q)
+  p ~== q = xs == ys && ps ~== qs where
+    (xs,ps) = unzip $ filter (not . (~== 0) . snd) $ normalForm p
+    (ys,qs) = unzip $ filter (not . (~== 0) . snd) $ normalForm q
diff --git a/src/Control/Monad/Bayes/Free.hs b/src/Control/Monad/Bayes/Free.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Free.hs
@@ -0,0 +1,90 @@
+{-|
+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 Data.Functor.Identity
+
+import Control.Monad.Trans
+import Control.Monad.Writer
+import Control.Monad.State
+import Control.Monad.Trans.Free.Church
+
+import Control.Monad.Bayes.Class
+
+-- | 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 (FT SamF m a)
+  deriving(Functor,Applicative,Monad,MonadTrans)
+
+runFreeSampler :: FreeSampler m a -> FT SamF m a
+runFreeSampler (FreeSampler m) = m
+
+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
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Helpers.hs
@@ -0,0 +1,67 @@
+{-|
+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.Weighted as Weighted
+import Control.Monad.Bayes.Population as Pop
+import Control.Monad.Bayes.Sequential as Seq
+import Control.Monad.Bayes.Free as Free
+import Control.Monad.Bayes.Traced as Tr
+
+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/PMMH.hs b/src/Control/Monad/Bayes/Inference/PMMH.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/PMMH.hs
@@ -0,0 +1,38 @@
+{-|
+Module      : Control.Monad.Bayes.Inference.PMMH
+Description : Particle Marginal Metropolis-Hastings (PMMH)
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+Particle Marginal Metropolis-Hastings (PMMH) sampling.
+
+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
+)  where
+
+import Numeric.Log
+
+import Control.Monad.Trans (lift)
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Sequential
+import Control.Monad.Bayes.Population as Pop
+import Control.Monad.Bayes.Traced
+import Control.Monad.Bayes.Inference.SMC
+
+-- | Particle Marginal Metropolis-Hastings sampling.
+pmmh :: MonadInfer m
+     => Int -- ^ number of Metropolis-Hastings steps
+     -> Int -- ^ number of time steps
+     -> Int -- ^ number of particles
+     -> Traced m b -- ^ model parameters prior
+     -> (b -> Sequential (Population m) a) -- ^ model
+     -> m [[(a, Log Double)]]
+pmmh t k n param model =
+  mh t (param >>= runPopulation . pushEvidence . Pop.hoist lift . smcSystematic k n . model)
diff --git a/src/Control/Monad/Bayes/Inference/RMSMC.hs b/src/Control/Monad/Bayes/Inference/RMSMC.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/RMSMC.hs
@@ -0,0 +1,70 @@
+{-|
+Module      : Control.Monad.Bayes.Inference.RMSMC
+Description : Resample-Move Sequential Monte Carlo (RM-SMC)
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+Resample-move Sequential Monte Carlo (RM-SMC) sampling.
+
+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,
+  rmsmcBasic
+) where
+
+import Control.Monad.Bayes.Class
+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.Dynamic as TrDyn
+import qualified Control.Monad.Bayes.Traced.Basic as TrBas
+import Control.Monad.Bayes.Helpers
+
+-- | Resample-move Sequential Monte Carlo.
+rmsmc :: MonadSample m
+      => Int -- ^ number of timesteps
+      -> Int -- ^ number of particles
+      -> Int -- ^ number of Metropolis-Hastings transitions after each resampling
+      -> Sequential (Traced (Population m)) a -- ^ model
+      -> Population m a
+rmsmc k n t =
+  marginal .
+  sis (composeCopies t mhStep . hoistT resampleSystematic) k .
+  hoistS (hoistT (spawn n >>))
+
+-- | Resample-move Sequential Monte Carlo with a more efficient
+-- tracing representation.
+rmsmcBasic :: MonadSample m
+      => Int -- ^ number of timesteps
+      -> Int -- ^ number of particles
+      -> Int -- ^ number of Metropolis-Hastings transitions after each resampling
+      -> Sequential (TrBas.Traced (Population m)) a -- ^ model
+      -> Population m a
+rmsmcBasic k n t =
+  TrBas.marginal .
+  sis (composeCopies t TrBas.mhStep . TrBas.hoistT resampleSystematic) k .
+  hoistS (TrBas.hoistT (spawn n >>))
+
+-- | A variant of resample-move Sequential Monte Carlo
+-- where only random variables since last resampling are considered
+-- for rejuvenation.
+rmsmcLocal :: MonadSample m
+           => Int -- ^ number of timesteps
+           -> Int -- ^ number of particles
+           -> Int -- ^ number of Metropolis-Hastings transitions after each resampling
+           -> Sequential (TrDyn.Traced (Population m)) a -- ^ model
+           -> Population m a
+rmsmcLocal k n t =
+  TrDyn.marginal .
+  sis (TrDyn.freeze . composeCopies t TrDyn.mhStep . TrDyn.hoistT resampleSystematic) k .
+  hoistS (TrDyn.hoistT (spawn n >>))
+
+-- | Apply a function a given number of times.
+composeCopies :: Int -> (a -> a) -> (a -> a)
+composeCopies k f = foldr (.) id (replicate k f)
diff --git a/src/Control/Monad/Bayes/Inference/SMC.hs b/src/Control/Monad/Bayes/Inference/SMC.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/SMC.hs
@@ -0,0 +1,73 @@
+{-|
+Module      : Control.Monad.Bayes.Inference.SMC
+Description : Sequential Monte Carlo (SMC)
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+Sequential Monte Carlo (SMC) sampling.
+
+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
+) where
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Population
+import Control.Monad.Bayes.Sequential as Seq
+
+-- | Sequential importance resampling.
+-- Basically an SMC template that takes a custom resampler.
+sir :: Monad m
+    => (forall x. Population m x -> Population m x) -- ^ resampler
+    -> Int -- ^ number of timesteps
+    -> Int -- ^ population size
+    -> Sequential (Population m) a -- ^ model
+    -> 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
+               => Int -- ^ number of timesteps
+               -> Int -- ^ number of particles
+               -> Sequential (Population m) a -- ^ model
+               -> Population m a
+smcMultinomial = sir resampleMultinomial
+
+-- | Sequential Monte Carlo with systematic resampling at each timestep.
+-- Weights are not normalized.
+smcSystematic  :: MonadSample m
+               => Int -- ^ number of timesteps
+               -> Int -- ^ number of particles
+               -> Sequential (Population m) a -- ^ model
+               -> Population m a
+smcSystematic = sir resampleSystematic
+
+-- | 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
+                   => Int -- ^ number of timesteps
+                   -> Int -- ^ number of particles
+                   -> Sequential (Population m) a -- ^ model
+                   -> 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
+                   => Int -- ^ number of timesteps
+                   -> Int -- ^ number of particles
+                   -> Sequential (Population m) a -- ^ model
+                   -> Population m a
+smcSystematicPush = sir (pushEvidence . resampleSystematic)
diff --git a/src/Control/Monad/Bayes/Inference/SMC2.hs b/src/Control/Monad/Bayes/Inference/SMC2.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Inference/SMC2.hs
@@ -0,0 +1,55 @@
+{-|
+Module      : Control.Monad.Bayes.Inference.SMC2
+Description : Sequential Monte Carlo squared (SMC²)
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+Sequential Monte Carlo squared (SMC²) sampling.
+
+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
+) where
+
+import Numeric.Log
+import Control.Monad.Trans
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Population as Pop
+import Control.Monad.Bayes.Inference.SMC
+import Control.Monad.Bayes.Inference.RMSMC
+import Control.Monad.Bayes.Helpers
+
+-- | Helper monad transformer for preprocessing the model for 'smc2'.
+newtype SMC2 m a = SMC2 (S (T (P m)) a)
+  deriving(Functor,Applicative,Monad)
+setup :: SMC2 m a -> S (T (P m)) a
+setup (SMC2 m) = m
+
+instance MonadTrans SMC2 where
+  lift = SMC2 . lift . lift . lift
+
+instance MonadSample m => MonadSample (SMC2 m) where
+  random = lift random
+
+instance Monad m => MonadCond (SMC2 m) where
+  score = SMC2 . score
+
+instance MonadSample m => MonadInfer (SMC2 m)
+
+-- | Sequential Monte Carlo squared.
+smc2 :: MonadSample m
+     => Int -- ^ number of time steps
+     -> Int -- ^ number of inner particles
+     -> Int -- ^ number of outer particles
+     -> Int -- ^ number of MH transitions
+     -> S (T (P m)) b -- ^ model parameters
+     -> ( b -> S (P (SMC2 m)) a) -- ^ model
+     -> P m [(a, Log Double)]
+smc2 k n p t param model =
+  rmsmc k p t (param >>= setup . runPopulation . smcSystematicPush k n . model)
diff --git a/src/Control/Monad/Bayes/Population.hs b/src/Control/Monad/Bayes/Population.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Population.hs
@@ -0,0 +1,192 @@
+{-|
+Module      : Control.Monad.Bayes.Population
+Description : Representation of distributions using multiple samples
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+'Population' turns a single sample into a collection of weighted samples.
+-}
+
+module Control.Monad.Bayes.Population (
+    Population,
+    runPopulation,
+    explicitPopulation,
+    fromWeightedList,
+    spawn,
+    resampleMultinomial,
+    resampleSystematic,
+    extractEvidence,
+    pushEvidence,
+    proper,
+    evidence,
+    collapse,
+    mapPopulation,
+    normalize,
+    popAvg,
+    flatten,
+    hoist
+                 ) where
+
+import Prelude hiding (sum, all)
+
+import Control.Arrow (second)
+import Control.Monad.Trans
+import Control.Monad.Trans.List
+import Control.Monad (replicateM)
+
+import qualified Data.List
+import qualified Data.Vector as V
+
+import Numeric.Log
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Weighted hiding (flatten, hoist)
+
+-- | A collection of weighted samples, or particles.
+newtype Population m a = Population (Weighted (ListT m) a)
+  deriving(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
+
+-- | Explicit representation of the weighted sample.
+explicitPopulation :: Functor m => Population m a -> m [(a, Double)]
+explicitPopulation = fmap (map (second (exp . ln))) . runPopulation
+
+-- | Initialize 'Population' with a concrete weighted sample.
+fromWeightedList :: Monad m => m [(a,Log Double)] -> Population m a
+fromWeightedList = Population . withWeight . ListT
+
+-- | Increase the sample size by a given factor.
+-- The weights are adjusted such that their sum is preserved.
+-- It is therefore safe to use 'spawn' in arbitrary places in the program
+-- without introducing bias.
+spawn :: Monad m => Int -> Population m ()
+spawn n = fromWeightedList $ pure $ replicate n ((), 1 / fromIntegral n)
+
+resampleGeneric :: MonadSample m
+                => (V.Vector Double -> m [Int]) -- ^ resampler
+                -> Population m a -> Population m a
+resampleGeneric resampler m = fromWeightedList $ do
+  pop <- runPopulation m
+  let (xs, ps) = unzip pop
+  let n = length xs
+  let z = sum ps
+  if z > 0 then do
+    let weights = V.fromList (map (exp . ln . (/z)) ps)
+    ancestors <- resampler weights
+    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
+    return pop
+
+-- | Systematic resampling helper.
+systematic :: Double -> V.Vector Double -> [Int]
+systematic u ps = f 0 (u / fromIntegral n) 0 0 [] where
+  prob i = ps V.! i
+  n = length ps
+  inc = 1 / fromIntegral n
+  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)
+    else
+      f i v (j + 1) (q + prob j) acc
+
+-- | Resample the population using the underlying monad and a systematic resampling scheme.
+-- The total weight is preserved.
+resampleSystematic :: (MonadSample m)
+         => Population m a -> Population m a
+resampleSystematic = resampleGeneric (\ps -> (`systematic` ps) <$> random)
+
+-- | Multinomial resampler.
+multinomial :: MonadSample m => V.Vector Double -> m [Int]
+multinomial ps = replicateM (V.length ps) (categorical ps)
+
+-- | Resample the population using the underlying monad and a multinomial resampling scheme.
+-- The total weight is preserved.
+resampleMultinomial :: (MonadSample m)
+         => Population m a -> Population m a
+resampleMultinomial = resampleGeneric multinomial
+
+-- | Separate the sum of weights into the 'Weighted' transformer.
+-- Weights are normalized after this operation.
+extractEvidence :: Monad m
+                => Population m a -> Population (Weighted m) a
+extractEvidence m = fromWeightedList $ do
+  pop <- lift $ runPopulation 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
+  factor z
+  return $ zip xs ws
+
+-- | Push the evidence estimator as a score to the transformed monad.
+-- Weights are normalized after this operation.
+pushEvidence :: MonadCond m
+           => Population m a -> Population m a
+pushEvidence = hoist applyWeight . extractEvidence
+
+-- | A properly weighted single sample, that is one picked at random according
+-- to the weights, with the sum of all weights.
+proper :: (MonadSample m)
+       => Population m a -> Weighted m a
+proper m = do
+  pop <- runPopulation $ extractEvidence m
+  let (xs, ps) = unzip pop
+  index <- logCategorical $ V.fromList ps
+  let x = xs !! index
+  return x
+
+-- | Model evidence estimator, also known as pseudo-marginal likelihood.
+evidence :: (Monad m) => Population m a -> m (Log Double)
+evidence = extractWeight . runPopulation . extractEvidence
+
+-- | Picks one point from the population and uses model evidence as a 'score'
+-- in the transformed monad.
+-- This way a single sample can be selected from a population without
+-- introducing bias.
+collapse :: (MonadInfer m)
+         => Population m a -> 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
+  xs <- explicitPopulation p
+  let ys = map (\(x,w) -> f x * w) xs
+  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)
+      => (forall x. m x -> n x) -> Population m a -> Population n a
+hoist f = fromWeightedList . f . runPopulation
diff --git a/src/Control/Monad/Bayes/Sampler.hs b/src/Control/Monad/Bayes/Sampler.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Sampler.hs
@@ -0,0 +1,110 @@
+{-|
+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.ST (ST, runST, stToIO)
+import System.Random.MWC
+import qualified System.Random.MWC.Distributions as MWC
+import Control.Monad.State (State, state)
+import Control.Monad.Trans (lift, MonadIO)
+import Control.Monad.Trans.Reader (ReaderT, runReaderT, ask, mapReaderT)
+
+import Control.Monad.Bayes.Class
+
+-- | 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/Sequential.hs b/src/Control/Monad/Bayes/Sequential.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Sequential.hs
@@ -0,0 +1,94 @@
+{-|
+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.Trans
+import Control.Monad.Coroutine hiding (suspend)
+import Control.Monad.Coroutine.SuspensionFunctors
+import Data.Either
+
+import Control.Monad.Bayes.Class
+
+-- | 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
+    => (forall x. m x -> m x) -- ^ transformation
+    -> Int -- ^ number of time steps
+    -> Sequential m a
+    -> m a
+sis f k = finish . composeCopies k (advance . hoistFirst f)
diff --git a/src/Control/Monad/Bayes/Traced.hs b/src/Control/Monad/Bayes/Traced.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Traced.hs
@@ -0,0 +1,16 @@
+{-|
+Module      : Control.Monad.Bayes.Traced
+Description : Distributions on execution traces
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+-}
+
+module Control.Monad.Bayes.Traced (
+  module Control.Monad.Bayes.Traced.Static
+) where
+
+import Control.Monad.Bayes.Traced.Static
diff --git a/src/Control/Monad/Bayes/Traced/Basic.hs b/src/Control/Monad/Bayes/Traced/Basic.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Traced/Basic.hs
@@ -0,0 +1,81 @@
+{-|
+Module      : Control.Monad.Bayes.Traced.Basic
+Description : Distributions on full execution traces of full programs
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+-}
+
+module Control.Monad.Bayes.Traced.Basic (
+  Traced,
+  hoistT,
+  marginal,
+  mhStep,
+  mh
+) where
+
+import Data.Functor.Identity
+import Control.Applicative (liftA2)
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Weighted as Weighted
+import Control.Monad.Bayes.Free as FreeSampler
+
+import Control.Monad.Bayes.Traced.Common
+
+-- | Tracing monad that records random choices made in the program.
+-- The first component is used to run the program with a modified trace,
+-- while the second records a trace and an output value from a run.
+data Traced m a = Traced (Weighted (FreeSampler Identity) a) (m (Trace a))
+
+traceDist :: Traced m a -> m (Trace a)
+traceDist (Traced _ d) = d
+
+model :: Traced m a -> Weighted (FreeSampler Identity) a
+model (Traced m _) = m
+
+instance Monad m => Functor (Traced m) where
+  fmap f (Traced m d) = Traced (fmap f m) (fmap (fmap f) d)
+
+instance Monad m => Applicative (Traced m) where
+  pure x = Traced (pure x) (pure (pure x))
+  (Traced mf df) <*> (Traced mx dx) = Traced (mf <*> mx) (liftA2 (<*>) df dx)
+
+instance Monad m => Monad (Traced m) where
+  (Traced mx dx) >>= f = Traced my dy where
+    my = mx >>= model . f
+    dy = dx `bind` (traceDist . f)
+
+instance MonadSample m => MonadSample (Traced m) where
+  random = Traced random (fmap singleton random)
+
+instance MonadCond m => MonadCond (Traced m) where
+  score w = Traced (score w) (score w >> pure (scored w))
+
+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)
+
+-- | Discard the trace and supporting infrastructure.
+marginal :: Monad m => Traced m a -> m a
+marginal (Traced _ d) = fmap output d
+
+-- | A single step of the Trace MH algorithm.
+mhStep :: MonadSample m => Traced m a -> Traced m a
+mhStep (Traced m d) = Traced m d' where
+  d' = d >>= mhTrans' m
+
+-- | Full run of the Trace MH algorithm with a specified
+-- number of steps.
+mh :: MonadSample m => Int -> Traced m a -> m [a]
+mh n (Traced m d) = fmap (map output) t where
+  t = f n
+  f 0 = fmap (:[]) d
+  f k = 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
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Traced/Common.hs
@@ -0,0 +1,81 @@
+{-|
+Module      : Control.Monad.Bayes.Traced.Common
+Description : Numeric code for Trace MCMC
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+-}
+
+module Control.Monad.Bayes.Traced.Common (
+  Trace,
+  singleton,
+  output,
+  scored,
+  bind,
+  mhTrans,
+  mhTrans'
+) where
+
+import Control.Monad.Trans.Writer
+import qualified Data.Vector as V
+import Data.Functor.Identity
+
+import Numeric.Log (Log, ln)
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Weighted as Weighted
+import Control.Monad.Bayes.Free as FreeSampler
+
+data Trace a =
+  Trace {
+    variables :: [Double],
+    output :: a,
+    density :: 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}
+  tf <*> tx = Trace {variables = variables tf ++ variables tx, output = output tf (output tx), density = density tf * density tx}
+
+instance Monad Trace where
+  t >>= f =
+    let t' = f (output t) in
+    t' {variables = variables t ++ variables t', density = density t * density t'}
+
+singleton :: Double -> Trace Double
+singleton u = Trace {variables = [u], output = u, density = 1}
+
+scored :: Log Double -> Trace ()
+scored w = Trace {variables = [], output = (), density = 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}
+
+-- | 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 = do
+  let us = variables t
+      p = density t
+  us' <- do
+    let n = length us
+    i <- categorical $ V.replicate n (1 / fromIntegral n)
+    u' <- random
+    let (xs, _:ys) = splitAt i us
+    return $ xs ++ (u':ys)
+  ((b, q), vs) <- runWriterT $ runWeighted $ Weighted.hoist (WriterT . withPartialRandomness us') m
+  let ratio = (exp . ln) $ min 1 (q * fromIntegral (length us) / (p * fromIntegral (length vs)))
+  accept <- bernoulli ratio
+  return $ if accept then Trace vs b q else 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)
diff --git a/src/Control/Monad/Bayes/Traced/Dynamic.hs b/src/Control/Monad/Bayes/Traced/Dynamic.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Traced/Dynamic.hs
@@ -0,0 +1,101 @@
+{-|
+Module      : Control.Monad.Bayes.Traced.Dynamic
+Description : Distributions on execution traces that can be dynamically frozen
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+-}
+
+module Control.Monad.Bayes.Traced.Dynamic (
+  Traced,
+  hoistT,
+  marginal,
+  freeze,
+  mhStep,
+  mh
+) where
+
+import Control.Monad (join)
+import Control.Monad.Trans
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Weighted as Weighted
+import Control.Monad.Bayes.Free as FreeSampler
+
+import Control.Monad.Bayes.Traced.Common
+
+-- | A tracing monad where only a subset of random choices are traced
+-- and this subset can be adjusted dynamically.
+newtype Traced m a = Traced (m (Weighted (FreeSampler m) a, Trace a))
+runTraced :: Traced m a -> m (Weighted (FreeSampler m) a, Trace a)
+runTraced (Traced c) = c
+
+pushM :: Monad m => m (Weighted (FreeSampler m) a) -> Weighted (FreeSampler m) a
+pushM = join . lift . lift
+
+instance Monad m => Functor (Traced m) where
+  fmap f (Traced c) = Traced $ do
+    (m, t) <- c
+    let m' = fmap f m
+    let t' = fmap f t
+    return (m', t')
+
+instance Monad m => Applicative (Traced m) where
+  pure x = Traced $ pure (pure x, pure x)
+  (Traced cf) <*> (Traced cx) = Traced $ do
+    (mf, tf) <- cf
+    (mx, tx) <- cx
+    return (mf <*> mx, tf <*> tx)
+
+instance Monad m => Monad (Traced m) where
+  (Traced cx) >>= f = Traced $ do
+    (mx, tx) <- cx
+    let m = mx >>= pushM . fmap fst . runTraced . f
+    t <- return tx `bind` (fmap snd . runTraced . f)
+    return (m, t)
+
+instance MonadTrans Traced where
+  lift m = Traced $ fmap ((,) (lift $ lift m) . pure) m
+
+instance MonadSample m => MonadSample (Traced m) where
+  random = Traced $ fmap ((,) random . singleton) random
+
+instance MonadCond m => MonadCond (Traced m) where
+  score w = Traced $ fmap (score w,) (score w >> pure (scored w))
+
+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)
+
+marginal :: Monad m => Traced m a -> m a
+marginal (Traced c) = fmap (output . snd) c
+
+-- | Freeze all traced random choices to their current
+-- values and stop tracing them.
+freeze :: Monad m => Traced m a -> Traced m a
+freeze (Traced c) = Traced $ do
+  (_, t) <- c
+  let x = output t
+  return (return x, pure x)
+
+mhStep :: MonadSample m => Traced m a -> Traced m a
+mhStep (Traced c) = Traced $ do
+  (m, t) <- c
+  t' <- mhTrans m t
+  return (m, t')
+
+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
diff --git a/src/Control/Monad/Bayes/Traced/Static.hs b/src/Control/Monad/Bayes/Traced/Static.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Traced/Static.hs
@@ -0,0 +1,80 @@
+{-|
+Module      : Control.Monad.Bayes.Traced.Static
+Description : Distributions on execution traces of full programs
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+-}
+
+module Control.Monad.Bayes.Traced.Static (
+  Traced,
+  hoistT,
+  marginal,
+  mhStep,
+  mh
+) where
+
+import Control.Monad.Trans
+import Control.Applicative (liftA2)
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Weighted as Weighted
+import Control.Monad.Bayes.Free as FreeSampler
+
+import Control.Monad.Bayes.Traced.Common
+
+-- | 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 (Weighted (FreeSampler m) a) (m (Trace a))
+
+traceDist :: Traced m a -> m (Trace a)
+traceDist (Traced _ d) = d
+
+model :: Traced m a -> Weighted (FreeSampler m) a
+model (Traced m _) = m
+
+instance Monad m => Functor (Traced m) where
+  fmap f (Traced m d) = Traced (fmap f m) (fmap (fmap f) d)
+
+instance Monad m => Applicative (Traced m) where
+  pure x = Traced (pure x) (pure (pure x))
+  (Traced mf df) <*> (Traced mx dx) = Traced (mf <*> mx) (liftA2 (<*>) df dx)
+
+instance Monad m => Monad (Traced m) where
+  (Traced mx dx) >>= f = Traced my dy where
+    my = mx >>= model . f
+    dy = dx `bind` (traceDist . f)
+
+instance MonadTrans Traced where
+  lift m = Traced (lift $ lift m) (fmap pure m)
+
+instance MonadSample m => MonadSample (Traced m) where
+  random = Traced random (fmap singleton random)
+
+instance MonadCond m => MonadCond (Traced m) where
+  score w = Traced (score w) (score w >> pure (scored w))
+
+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)
+
+marginal :: Monad m => Traced m a -> m a
+marginal (Traced _ d) = fmap output d
+
+mhStep :: MonadSample m => Traced m a -> Traced m a
+mhStep (Traced m d) = Traced m d' where
+  d' = d >>= mhTrans m
+
+mh :: MonadSample m => Int -> Traced m a -> m [a]
+mh n (Traced m d) = fmap (map output) t where
+  t = f n
+  f 0 = fmap (:[]) d
+  f k = do
+    ~(x:xs) <- f (k-1)
+    y <- mhTrans m x
+    return (y:x:xs)
diff --git a/src/Control/Monad/Bayes/Weighted.hs b/src/Control/Monad/Bayes/Weighted.hs
new file mode 100644
--- /dev/null
+++ b/src/Control/Monad/Bayes/Weighted.hs
@@ -0,0 +1,76 @@
+{-|
+Module      : Control.Monad.Bayes.Weighted
+Description : Probability monad accumulating the likelihood
+Copyright   : (c) Adam Scibior, 2015-2020
+License     : MIT
+Maintainer  : leonhard.markert@tweag.io
+Stability   : experimental
+Portability : GHC
+
+'Weighted' is an instance of 'MonadCond'. Apply a 'MonadSample' transformer to
+obtain a 'MonadInfer' that can execute probabilistic models.
+-}
+
+module Control.Monad.Bayes.Weighted (
+    Weighted,
+    withWeight,
+    runWeighted,
+    extractWeight,
+    prior,
+    flatten,
+    applyWeight,
+    hoist,
+                  ) where
+
+import Control.Monad.Trans
+import Control.Monad.Trans.State
+
+import Numeric.Log
+import Control.Monad.Bayes.Class
+
+-- | 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)
+
+instance Monad m => MonadCond (Weighted m) where
+  score w = Weighted (modify (* w))
+
+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
+
+-- | Compute the weight and discard the sample.
+extractWeight :: Functor m => Weighted m a -> m (Log Double)
+extractWeight m = snd <$> runWeighted m
+
+-- | Embed a random variable with explicitly given likelihood.
+--
+-- > runWeighted . withWeight = id
+withWeight :: (Monad m) => m (a, Log Double) -> Weighted m a
+withWeight m = Weighted $ do
+  (x,w) <- lift m
+  modify (* w)
+  return x
+
+-- | Discard the weight.
+-- This operation introduces bias.
+prior :: (Functor m) => Weighted m a -> m a
+prior = fmap fst . runWeighted
+
+-- | 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
+  factor w
+  return x
+
+-- | Apply a transformation to the transformed monad.
+hoist :: (forall x. m x -> n x) -> Weighted m a -> Weighted n a
+hoist t (Weighted m) = Weighted $ mapStateT t m
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,68 @@
+import Test.Hspec
+import Test.Hspec.QuickCheck
+import Test.QuickCheck
+
+import qualified TestWeighted
+import qualified TestEnumerator
+import qualified TestPopulation
+import qualified TestSequential
+import qualified TestInference
+
+
+main :: IO ()
+main = hspec $ do
+  describe "Weighted" $
+    it "accumulates likelihood correctly" $ do
+      passed <- TestWeighted.passed
+      passed `shouldBe` True
+  describe "Dist" $ 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 expicitly altered" $ do
+        popSize <- TestPopulation.popSize
+        popSize `shouldBe` 5
+      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
+        TestPopulation.transCheck1 `shouldBe` True
+        TestPopulation.transCheck2 `shouldBe` True
+      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
+      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
+      TestSequential.checkSync 0 `shouldBe` True
+      TestSequential.checkSync 1 `shouldBe` True
+      TestSequential.checkSync 2 `shouldBe` True
+  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
+          checkParticles <- TestInference.checkParticles observations particles
+          return $ checkParticles == particles
+  describe "SMC with systematic resampling" $ do
+    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
diff --git a/test/TestEnumerator.hs b/test/TestEnumerator.hs
new file mode 100644
--- /dev/null
+++ b/test/TestEnumerator.hs
@@ -0,0 +1,31 @@
+module TestEnumerator where
+
+import Data.AEq
+import qualified Data.Vector as V
+
+import Numeric.Log
+import Control.Monad.Bayes.Enumerator
+import Control.Monad.Bayes.Class
+import Sprinkler
+
+unnorm :: MonadSample m => m Int
+unnorm = categorical $ V.fromList [0.5,0.8]
+
+passed1 :: Bool
+passed1 = (exp . ln) (evidence unnorm) ~== 1
+
+agg :: MonadSample m => m Int
+agg = do
+  x <- uniformD [0,1]
+  y <- uniformD [2,1]
+  return (x+y)
+
+passed2 :: Bool
+passed2 = enumerate agg ~== [(1,0.25), (2,0.5), (3,0.25)]
+
+passed3 :: Bool
+passed3 = enumerate Sprinkler.hard ~== enumerate Sprinkler.soft
+
+passed4 :: Bool
+passed4 =
+ expectation (^ (2 :: Int)) (fmap (fromIntegral . (+1)) $ categorical $ V.fromList [0.5, 0.5]) ~== 2.5
diff --git a/test/TestInference.hs b/test/TestInference.hs
new file mode 100644
--- /dev/null
+++ b/test/TestInference.hs
@@ -0,0 +1,35 @@
+{-# LANGUAGE
+  Rank2Types,
+  TypeFamilies
+ #-}
+
+module TestInference where
+
+import Data.AEq
+import Numeric.Log
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Enumerator
+import Control.Monad.Bayes.Sampler
+import Control.Monad.Bayes.Population
+import Control.Monad.Bayes.Inference.SMC
+import Sprinkler
+
+sprinkler :: MonadInfer m => m Bool
+sprinkler = Sprinkler.soft
+
+-- | Count the number of particles produced by SMC
+checkParticles :: Int -> Int -> IO Int
+checkParticles observations particles =
+  sampleIOfixed (fmap length (runPopulation $ smcMultinomial observations particles Sprinkler.soft))
+
+checkParticlesSystematic :: Int -> Int -> IO Int
+checkParticlesSystematic observations particles =
+  sampleIOfixed (fmap length (runPopulation $ smcSystematic observations particles Sprinkler.soft))
+
+checkTerminateSMC :: IO [(Bool, Log Double)]
+checkTerminateSMC = sampleIOfixed (runPopulation $ smcMultinomial 2 5 sprinkler)
+
+checkPreserveSMC :: Bool
+checkPreserveSMC = (enumerate . collapse . smcMultinomial 2 2) sprinkler ~==
+                      enumerate sprinkler
diff --git a/test/TestPopulation.hs b/test/TestPopulation.hs
new file mode 100644
--- /dev/null
+++ b/test/TestPopulation.hs
@@ -0,0 +1,43 @@
+module TestPopulation where
+
+import Data.AEq
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Enumerator
+import Control.Monad.Bayes.Sampler
+import Control.Monad.Bayes.Population as Population
+import Sprinkler
+
+weightedSampleSize :: MonadSample m => Population m a -> m Int
+weightedSampleSize = fmap length . runPopulation
+
+popSize :: IO Int
+popSize = sampleIOfixed $ weightedSampleSize $ spawn 5 >> sprinkler
+
+manySize :: IO Int
+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
+
+transCheck1 :: Bool
+transCheck1 = enumerate (collapse sprinkler) ~==
+               sprinklerExact
+transCheck2 :: Bool
+transCheck2 = enumerate (collapse (spawn 2 >> sprinkler)) ~==
+               sprinklerExact
+
+resampleCheck :: Int -> Bool
+resampleCheck n =
+  (enumerate . collapse . resampleMultinomial) (spawn n >> sprinkler) ~==
+  sprinklerExact
+
+popAvgCheck :: Bool
+popAvgCheck = expectation f Sprinkler.soft ~== expectation id (popAvg f $ pushEvidence Sprinkler.soft) where
+  f True = 10
+  f False = 4
diff --git a/test/TestSequential.hs b/test/TestSequential.hs
new file mode 100644
--- /dev/null
+++ b/test/TestSequential.hs
@@ -0,0 +1,47 @@
+module TestSequential where
+
+import Data.AEq
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Enumerator as Dist
+import Control.Monad.Bayes.Sequential
+import Sprinkler
+
+twoSync :: MonadInfer m => m Int
+twoSync = do
+  x <- uniformD[0,1]
+  factor (fromIntegral x)
+  y <- uniformD[0,1]
+  factor (fromIntegral y)
+  return (x+y)
+
+finishedTwoSync :: MonadInfer m => Int -> m Bool
+finishedTwoSync n = finished (run n twoSync) where
+  run 0 d = d
+  run k d = run (k-1) (advance d)
+
+checkTwoSync :: Int -> Bool
+checkTwoSync 0 = mass (finishedTwoSync 0) False ~== 1
+checkTwoSync 1 = mass (finishedTwoSync 1) False ~== 1
+checkTwoSync 2 = mass (finishedTwoSync 2) True  ~== 1
+checkTwoSync _ = error "Unexpected argument"
+
+sprinkler :: MonadInfer m => m Bool
+sprinkler = Sprinkler.soft
+
+checkPreserve :: Bool
+checkPreserve = enumerate (finish sprinkler) ~== enumerate sprinkler
+
+pFinished :: Int -> Double
+pFinished 0 = 0.8267716535433071
+pFinished 1 = 0.9988062077198566
+pFinished 2 = 1
+pFinished _ = error "Unexpected argument"
+
+isFinished :: MonadInfer m => Int -> m Bool
+isFinished n = finished (run n sprinkler) where
+  run 0 d = d
+  run k d = run (k-1) (advance d)
+
+checkSync :: Int -> Bool
+checkSync n = mass (isFinished n) True ~== pFinished n
diff --git a/test/TestWeighted.hs b/test/TestWeighted.hs
new file mode 100644
--- /dev/null
+++ b/test/TestWeighted.hs
@@ -0,0 +1,34 @@
+{-# LANGUAGE
+  TypeFamilies
+ #-}
+
+module TestWeighted where
+
+import Data.AEq
+import Control.Monad.State
+import Data.Bifunctor (second)
+import Numeric.Log
+
+import Control.Monad.Bayes.Class
+import Control.Monad.Bayes.Sampler
+import Control.Monad.Bayes.Weighted
+
+model :: MonadInfer m => m (Int,Double)
+model = do
+  n <- uniformD [0,1,2]
+  unless (n == 0) (factor 0.5)
+  x <- if n == 0 then return 1 else normal 0 1
+  when (n == 2) (factor $ (Exp . log) (x*x))
+  return (n,x)
+
+result :: MonadSample m => m ((Int,Double), Double)
+result = second (exp . ln) <$> runWeighted model
+
+passed :: IO Bool
+passed = fmap check (sampleIOfixed result)
+
+check :: ((Int,Double), Double) -> Bool
+check ((0,1),1) = True
+check ((1,_),y) =  y ~== 0.5
+check ((2,x),y) =  y ~== 0.5 * x * x
+check _ = False