packages feed

monad-bayes 0.1.1.0 → 1.0.0

raw patch · 52 files changed

+2938/−1179 lines, 52 filesdep +brickdep +foldldep +histogram-filldep ~basedep ~containersdep ~ieee754new-uploader

Dependencies added: brick, foldl, histogram-fill, integration, lens, linear, matrix, monad-extras, pipes, pretty-simple, primitive, profunctors, random, scientific, text, typed-process, vty

Dependency ranges changed: base, containers, ieee754, math-functions, monad-coroutine, mtl, mwc-random, safe, statistics, transformers, vector

Files

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