packages feed

monad-bayes (empty) → 0.1.0.0

raw patch · 36 files changed

+2692/−0 lines, 36 filesdep +QuickCheckdep +abstract-pardep +basesetup-changed

Dependencies added: QuickCheck, abstract-par, base, containers, criterion, free, hspec, ieee754, log-domain, math-functions, monad-bayes, monad-coroutine, mtl, mwc-random, optparse-applicative, process, safe, statistics, time, transformers, vector

Files

+ CHANGELOG.md view
@@ -0,0 +1,3 @@+# 0.1.0.0 (2020-02-17)++Initial release.
+ LICENSE view
@@ -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.+
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ benchmark/SSM.hs view
@@ -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
+ benchmark/Single.hs view
@@ -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)
+ benchmark/Speed.hs view
@@ -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)
+ models/Dice.hs view
@@ -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
+ models/HMM.hs view
@@ -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]
+ models/LDA.hs view
@@ -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
+ models/LogReg.hs view
@@ -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)
+ models/NonlinearSSM.hs view
@@ -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
+ models/Sprinkler.hs view
@@ -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
+ monad-bayes.cabal view
@@ -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
+ src/Control/Monad/Bayes/Class.hs view
@@ -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)
+ src/Control/Monad/Bayes/Enumerator.hs view
@@ -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
+ src/Control/Monad/Bayes/Free.hs view
@@ -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
+ src/Control/Monad/Bayes/Helpers.hs view
@@ -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
+ src/Control/Monad/Bayes/Inference/PMMH.hs view
@@ -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)
+ src/Control/Monad/Bayes/Inference/RMSMC.hs view
@@ -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)
+ src/Control/Monad/Bayes/Inference/SMC.hs view
@@ -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)
+ src/Control/Monad/Bayes/Inference/SMC2.hs view
@@ -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)
+ src/Control/Monad/Bayes/Population.hs view
@@ -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
+ src/Control/Monad/Bayes/Sampler.hs view
@@ -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
+ src/Control/Monad/Bayes/Sequential.hs view
@@ -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)
+ src/Control/Monad/Bayes/Traced.hs view
@@ -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
+ src/Control/Monad/Bayes/Traced/Basic.hs view
@@ -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)
+ src/Control/Monad/Bayes/Traced/Common.hs view
@@ -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)
+ src/Control/Monad/Bayes/Traced/Dynamic.hs view
@@ -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
+ src/Control/Monad/Bayes/Traced/Static.hs view
@@ -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)
+ src/Control/Monad/Bayes/Weighted.hs view
@@ -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
+ test/Spec.hs view
@@ -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
+ test/TestEnumerator.hs view
@@ -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
+ test/TestInference.hs view
@@ -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
+ test/TestPopulation.hs view
@@ -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
+ test/TestSequential.hs view
@@ -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
+ test/TestWeighted.hs view
@@ -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