From a582f755d744e88eed058d931894b1095e72c9bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Manuel=20B=C3=A4renz?= Date: Mon, 17 Jul 2023 16:42:14 +0200 Subject: [PATCH 1/2] Make transformer naming idiomatic. Fixes #189. Fixes #242. --- benchmark/SSM.hs | 2 +- benchmark/Speed.hs | 6 +- docs/docs/probprog.md | 117 ++++--- docs/docs/usage.md | 320 +++++++++--------- models/Helper.hs | 4 +- models/NonlinearSSM/Algorithms.hs | 10 +- notebooks/Implementation.ipynb | 12 +- notebooks/examples/Diagrams.ipynb | 4 +- notebooks/examples/Histogram.ipynb | 2 +- notebooks/examples/Parsing.ipynb | 2 +- notebooks/examples/Streaming.ipynb | 2 +- notebooks/tutorials/AdvancedSampling.ipynb | 6 +- notebooks/tutorials/Bayesian.ipynb | 8 +- notebooks/tutorials/Lazy.ipynb | 2 +- notebooks/tutorials/MCMC.ipynb | 2 +- notebooks/tutorials/SMC.ipynb | 2 +- notebooks/tutorials/Sampling.ipynb | 12 +- src/Control/Monad/Bayes/Density/Free.hs | 40 +-- src/Control/Monad/Bayes/Density/State.hs | 26 +- src/Control/Monad/Bayes/Inference/Lazy/MH.hs | 10 +- src/Control/Monad/Bayes/Inference/Lazy/WIS.hs | 8 +- src/Control/Monad/Bayes/Inference/MCMC.hs | 16 +- src/Control/Monad/Bayes/Inference/PMMH.hs | 18 +- src/Control/Monad/Bayes/Inference/RMSMC.hs | 16 +- src/Control/Monad/Bayes/Inference/SMC.hs | 10 +- src/Control/Monad/Bayes/Inference/SMC2.hs | 16 +- src/Control/Monad/Bayes/Inference/TUI.hs | 4 +- src/Control/Monad/Bayes/Integrator.hs | 36 +- src/Control/Monad/Bayes/Population.hs | 92 +++-- src/Control/Monad/Bayes/Sampler/Lazy.hs | 26 +- src/Control/Monad/Bayes/Sampler/Strict.hs | 34 +- .../Monad/Bayes/Sequential/Coroutine.hs | 40 +-- src/Control/Monad/Bayes/Traced/Basic.hs | 50 +-- src/Control/Monad/Bayes/Traced/Common.hs | 20 +- src/Control/Monad/Bayes/Traced/Dynamic.hs | 58 ++-- src/Control/Monad/Bayes/Traced/Static.hs | 54 +-- src/Control/Monad/Bayes/Weighted.hs | 44 ++- test/Spec.hs | 2 +- test/TestAdvanced.hs | 10 +- test/TestInference.hs | 26 +- test/TestIntegrator.hs | 6 +- test/TestPopulation.hs | 8 +- test/TestWeighted.hs | 4 +- 43 files changed, 590 insertions(+), 597 deletions(-) diff --git a/benchmark/SSM.hs b/benchmark/SSM.hs index c57c8dcb..ec873a0a 100644 --- a/benchmark/SSM.hs +++ b/benchmark/SSM.hs @@ -7,7 +7,7 @@ import Control.Monad.Bayes.Inference.RMSMC (rmsmcDynamic) import Control.Monad.Bayes.Inference.SMC import Control.Monad.Bayes.Inference.SMC2 as SMC2 (smc2) import Control.Monad.Bayes.Population -import Control.Monad.Bayes.Population (population, resampleMultinomial) +import Control.Monad.Bayes.Population (resampleMultinomial, runPopulationT) import Control.Monad.Bayes.Sampler.Strict (sampleIO, sampleIOfixed, sampleWith) import Control.Monad.Bayes.Weighted (unweighted) import Control.Monad.IO.Class (MonadIO (liftIO)) diff --git a/benchmark/Speed.hs b/benchmark/Speed.hs index c9b6f7b2..0d72808c 100644 --- a/benchmark/Speed.hs +++ b/benchmark/Speed.hs @@ -8,7 +8,7 @@ import Control.Monad.Bayes.Class (MonadMeasure) import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (MCMCConfig, numBurnIn, numMCMCSteps, proposal), Proposal (SingleSiteMH)) import Control.Monad.Bayes.Inference.RMSMC (rmsmcDynamic) import Control.Monad.Bayes.Inference.SMC (SMCConfig (SMCConfig, numParticles, numSteps, resampler), smc) -import Control.Monad.Bayes.Population (population, resampleSystematic) +import Control.Monad.Bayes.Population (resampleSystematic, runPopulationT) import Control.Monad.Bayes.Sampler.Strict (SamplerIO, sampleIOfixed) import Control.Monad.Bayes.Traced (mh) import Control.Monad.Bayes.Weighted (unweighted) @@ -59,10 +59,10 @@ instance Show Alg where runAlg :: Model -> Alg -> SamplerIO String runAlg model (MH n) = show <$> unweighted (mh n (buildModel model)) -runAlg model (SMC n) = show <$> population (smc SMCConfig {numSteps = (modelLength model), numParticles = n, resampler = resampleSystematic} (buildModel model)) +runAlg model (SMC n) = show <$> runPopulationT (smc SMCConfig {numSteps = (modelLength model), numParticles = n, resampler = resampleSystematic} (buildModel model)) runAlg model (RMSMC n t) = show - <$> population + <$> runPopulationT ( rmsmcDynamic MCMCConfig {numMCMCSteps = t, numBurnIn = 0, proposal = SingleSiteMH} SMCConfig {numSteps = modelLength model, numParticles = n, resampler = resampleSystematic} diff --git a/docs/docs/probprog.md b/docs/docs/probprog.md index c9661f4f..ea153882 100644 --- a/docs/docs/probprog.md +++ b/docs/docs/probprog.md @@ -20,7 +20,7 @@ $$ $$ $$ -y_{n}=\alpha+\beta x_{n}+\epsilon_{n} +y_{n}=\alpha+\beta x_{n}+\epsilon_{n} $$ but in code as: @@ -48,12 +48,12 @@ This is the *model*. To perform *inference* , suppose we have a data set `xsys` To run the model -```haskell -mhRunsRegression = sampler - $ mcmc MCMCConfig - {numMCMCSteps = 1500, - proposal = SingleSiteMH, - numBurnIn = 500} +```haskell +mhRunsRegression = sampler + $ mcmc MCMCConfig + {numMCMCSteps = 1500, + proposal = SingleSiteMH, + numBurnIn = 500} (regression xsys) ``` @@ -61,7 +61,7 @@ This yields 1000 samples from an MCMC walk using an MH kernel. `mh n` produces a ![](/images/regress.png) -Monad-bayes provides a variety of MCMC and SMC methods, and methods arising from the composition of the two. +Monad-bayes provides a variety of MCMC and SMC methods, and methods arising from the composition of the two. - + As a concrete example, here is a probabilistic program: ```haskell -(sampler . population . smc SMCConfig {numSteps = 2, numParticles = 2, resampler = resampleMultinomial} random) - :: Sequential (Population SamplerIO) a -> IO [(a, Numeric.Log.Log Double)] +(sampler . runPopulationT . smc SMCConfig {numSteps = 2, numParticles = 2, resampler = resampleMultinomial} random) + :: SequentialT (PopulationT SamplerIO) a -> IO [(a, Numeric.Log.Log Double)] ``` --> - + As a concrete example, here is a probabilistic program: @@ -430,13 +430,13 @@ example = do return x ``` -And here is the inference: +And here is the inference: ```haskell run :: IO [(Bool, Log Double)] -run = (sampler . population . smc (SMCConfig { - numSteps = 2, - numParticles = 4, +run = (sampler . runPopulationT . smc (SMCConfig { + numSteps = 2, + numParticles = 4, resampler = resampleMultinomial})) example ``` @@ -462,8 +462,8 @@ rmsmcBasic :: MCMCConfig -> SMCConfig m -> -- | model - Sequential (TrBas.Traced (Population m)) a -> - Population m a + SequentialT (TrBas.TracedT (PopulationT m)) a -> + PopulationT m a ``` For instance: @@ -471,11 +471,11 @@ For instance: ```haskell run :: IO [(Bool, Log Double)] run = ( - sampler . - population . - rmsmcBasic + sampler . + runPopulationT . + rmsmcBasic MCMCConfig {numMCMCSteps = 4, proposal = SingleSiteMH, numBurnIn = 0} - SMCConfig {numParticles = 4, numSteps = 4}) + SMCConfig {numParticles = 4, numSteps = 4}) example ``` @@ -483,14 +483,14 @@ What this returns is a population of samples, just like plain `SMC`. More MCMC s ## Particle Marginal Metropolis Hastings -This inference method takes a prior and a model separately, so it only applies to a (large) subset of probabilistic programs. +This inference method takes a prior and a model separately, so it only applies to a (large) subset of probabilistic programs. Run it like this: ```haskell -sampler $ pmmh - mcmcConfig {numMCMCSteps = 100} - smcConfig {numSteps = 0, numParticles = 100} +sampler $ pmmh + mcmcConfig {numMCMCSteps = 100} + smcConfig {numSteps = 0, numParticles = 100} random (normal 0) ``` @@ -551,7 +551,7 @@ betaBernoulli n = do # Interoperating with other Haskell code -Probabilistic programs in monad-bayes are Haskell programs. This contrasts to many probabilistic programming languages, which are deeply embedded and cannot smoothly interact with their host language. +Probabilistic programs in monad-bayes are Haskell programs. This contrasts to many probabilistic programming languages, which are deeply embedded and cannot smoothly interact with their host language. For example, we can use ordinary monadic combinators, as in: @@ -562,7 +562,7 @@ example = do return x ``` -or +or ```haskell example = whileM (bernoulli 0.99) (normal 0 1) @@ -598,7 +598,7 @@ mixture1 point = do is a piece of code to infer whether an observed point was generated from a Gaussian of mean $1$ or $5$. That is, `mixture1` is a conditional Bernoulli distribution over the mean given an observation. You're not going to be able to do much with `mixture1` though. Exact inference is impossible because of the sample from the normal, and as for sampling, there is zero probability of sampling the normal to exactly match the observed point, which is what the `condition` requires. -However, the same conditional distribution is represented by +However, the same conditional distribution is represented by ```haskell mixture2 point = do @@ -619,7 +619,7 @@ yields [(1.0, 0.98...), (5.0, 0.017...)] ``` -as well as sampling methods. +as well as sampling methods. The local lesson here is that you shouldn't `condition` on samples from a continuous distribution and expect a sampling based inference method to work. But the more general lesson is that you aren't exempted from thinking about inference when specifying your model. Alas. @@ -646,7 +646,7 @@ incremental = do Like in the previous example, `allAtOnce` and `incremental` denote the same distribution, namely `[(True, 1.0), (False, 0.0)]`. However, any inference algorithm for `allAtOnce` will have to explore all 4 possible variable assignments (`[(True,True), (True, False), (False, True), (False, False)]`). Meanwhile, `incremental` opens the possibility for inference algorithms to first determine the value of `x` and then of `y`. -In this example, the performance difference is negligible, but it's easy to extend this to models where it's the difference between something tractable and something intractable. +In this example, the performance difference is negligible, but it's easy to extend this to models where it's the difference between something tractable and something intractable. @@ -674,9 +674,8 @@ There is a variety of universal probabilistic programming libraries and/or langu A lot of engineering work has been put into the above libraries and languages to make them practical for real-world problems. While monad-bayes' core is very nice, it doesn't come with a lot of the batteries you might want. (The author's PhD thesis contains this relevant paragraph: "our library implements basic versions of advanced sampling algorithms. However, their successful application in practice requires incorporating established heuristics, such as: adaptive proposal distributions, controlling resampling with effective sample size, tuning rejuvenation kernels based on population in SMC2, and so on.") -**What Monad-Bayes has that is unique**: +**What Monad-Bayes has that is unique**: Monad-Bayes is just a library, unlike almost all other PPLs, which are separate languages written inside another language. As such, probabilistic programs in monad-bayes are first class programs in Haskell with no new special syntax or keywords. This allows all of Haskell's expressive power to be brought to bear. You can write distributions over any datatype (lists, trees, functions, histograms, JSON files, graphs, diagrams, etc). You can use powerful libraries like `pipes`, `lens` and `Parsec`. Everything is pure. Everything is strongly typed. You can make use of laziness. Models are monadic and inference is modular. Complex inference algorithms like RMSMC or PMMH are built out of simple composable pieces, and so are expressable extraordinarily simply. - diff --git a/docs/docs/usage.md b/docs/docs/usage.md index 52bc96dd..d6cd9e9e 100644 --- a/docs/docs/usage.md +++ b/docs/docs/usage.md @@ -16,7 +16,7 @@ The library relies on two core typeclasses `MonadDistribution` and `MonadFactor` (MonadDistribution m, MonadFactor m) => MonadMeasure m ``` -You can find these in `Control.Monad.Bayes.Class`. +You can find these in `Control.Monad.Bayes.Class`. ### MonadDistribution @@ -44,7 +44,7 @@ normal m s = fmap (quantile (normalDistr m s)) random `normalDistr` comes from a separate library `Statistics.Distribution.Normal` and `quantile (normalDistr m s) :: Double -> Double` is the inverse CDF of the normal, a deterministic function. -Again, to emphasize: **all of our randomness can be reduced to draws from a uniform distribution over the interval $[0,1]$**. +Again, to emphasize: **all of our randomness can be reduced to draws from a uniform distribution over the interval $[0,1]$**. So we now have a way of constructing distributions in a monadic fashion. As a simple example: @@ -56,9 +56,9 @@ example = do return (x + y > 1.5) ``` -Think of this as the procedure of first sampling uniformly from $[0,1]$, then from $[0,x]$, and then returning the Boolean $x + y > 1.5$. More precisely, this is the **marginal** probability of $x + y > 1.5$. +Think of this as the procedure of first sampling uniformly from $[0,1]$, then from $[0,x]$, and then returning the Boolean $x + y > 1.5$. More precisely, this is the **marginal** probability of $x + y > 1.5$. -**Technical note**: `MonadDistribution` actually contains a number of other distributions beyond `random`, which by default are defined in terms of `random`, but allow for different definitions when desired. For example, `Sampler` (an instance of `MonadDistribution` in Control.Monad.Sampler) defines `normal` and other distributions independently of `random`. +**Technical note**: `MonadDistribution` actually contains a number of other distributions beyond `random`, which by default are defined in terms of `random`, but allow for different definitions when desired. For example, `SamplerT` (an instance of `MonadDistribution` in Control.Monad.Sampler) defines `normal` and other distributions independently of `random`. @@ -67,7 +67,7 @@ Think of this as the procedure of first sampling uniformly from $[0,1]$, then fr ### MonadFactor -Here's `MonadFactor`. +Here's `MonadFactor`. ```haskell class Monad m => MonadFactor m where @@ -84,7 +84,7 @@ class Monad m => MonadFactor m where Now core idea of monad-bayes is that various monads will be made to be instances of `MonadDistribution`, `MonadFactor` or both (i.e. an instance of `MonadMeasure`), and different inference algorithms will be written using these instances. This separates the specification of the model (which happens abstractly in `MonadMeasure`) from the inference algorithm, which takes place in on of the concrete instances. The clever part of monad-bayes is that it allows this instances to be constructed in a modular way, using monad transformers. In the paper, these are termed *inference transformers* to emphasize that it doesn't really matter whether they satisfy the monad laws. -For example, to run weighted rejection sampling on a probabilistic program `p`, we can write `(sampler . weighted) p`. Here, `(sampler . weighted) :: Weighted SamplerIO a -> IO a`. So `p` gets understood as being of type `Weighted SamplerIO`, a type we'll encounter soon. +For example, to run weighted rejection sampling on a probabilistic program `p`, we can write `(sampler . runWeightedT) p`. Here, `(sampler . runWeightedT) :: WeightedT SamplerIO a -> IO a`. So `p` gets understood as being of type `WeightedT SamplerIO`, a type we'll encounter soon. Some of these transformers are easy to understand (like `StateT Double`, while others (like the Church transformed Free monad transformer) lie on the more advanced side of things. The following tour of these types goes from easy to hard. @@ -99,7 +99,7 @@ Summary of key info: `Enumerator` is in `Control.Monad.Bayes.Enumerator`, defined as follows: ```haskell -newtype Enumerator a = +newtype Enumerator a = Enumerator (WriterT (Product (Log Double)) [] a) ``` @@ -134,7 +134,7 @@ But for this to be well-typed, we need `Enumerator` to be an instance of `MonadM `Enumerator` is a monad automatically, because `WriterT a m` is a monad for `a` a monoid and `m` a monad. As needed, `[]` is a monad and `Log Double` is a monoid. But what does that monad actually **do**? -For instance, if I have a weighted list `l` like `[(False,0.8914),(True,0.1086)]` and a function `f :: Bool -> Enumerator a`, what is `l >>= f`? It takes each element in `l`, and passes it through `f`, to obtain a weighted list of weighted lists. Then it flattens it, by including each element in the inner lists in the resulting list, with the weight of the list multiplied by the weight of the element. As for `return x`, it gives the singleton list containing the pair of `x` and the unit of the product monoid, i.e. `[(x, 1.0)]`. +For instance, if I have a weighted list `l` like `[(False,0.8914),(True,0.1086)]` and a function `f :: Bool -> Enumerator a`, what is `l >>= f`? It takes each element in `l`, and passes it through `f`, to obtain a weighted list of weighted lists. Then it flattens it, by including each element in the inner lists in the resulting list, with the weight of the list multiplied by the weight of the element. As for `return x`, it gives the singleton list containing the pair of `x` and the unit of the product monoid, i.e. `[(x, 1.0)]`. This is the essence of propagating probability forward. @@ -153,11 +153,11 @@ instance MonadDistribution Enumerator where The first thing to notice is that `random` is actually not defined for `Enumerator`, and consequently, you can't handle any continuous distributions. This makes sense, because you can't represent continuous distributions (whose support is uncountably infinite) as a list of samples. -So really `Enumerator` isn't very general purpose. It's best understood as a didactic tool, rather than a real world inference algorithm. +So really `Enumerator` isn't very general purpose. It's best understood as a didactic tool, rather than a real world inference algorithm. But it can handle discrete distributions. It does this via `bernoulli` (as well as `categorical`, which I've omitted). `bernoulli p` constructs the weighted list corresponding to a `bernoulli` distribution with parameter `p`. -And the `MonadFactor` instance: +And the `MonadFactor` instance: ```haskell instance MonadFactor Enumerator where @@ -201,13 +201,13 @@ Summary of key info on `SamplerIO`: Monad-Bayes actually provides a more general constructor: ```haskell -newtype Sampler g m a = Sampler (ReaderT g m a) +newtype SamplerT g m a = SamplerT (ReaderT g m a) ``` `SamplerIO` just specializes `g`, which is the random number generator, and `m`, the IO-handling monad: ```haskell -type SamplerIO = Sampler (IOGenM StdGen) IO +type SamplerIO = SamplerT (IOGenM StdGen) IO ``` But you can specialize many other ways (see the `random` package), depending on your use case. @@ -215,13 +215,13 @@ But you can specialize many other ways (see the `random` package), depending on As the names suggest, `SamplerIO` instantiates an abstract distribution as a sampling procedure. We have ```haskell -instance MonadDistribution (Sampler g m) where - random = Sampler (ReaderT uniformDouble01M) +instance MonadDistribution (SamplerT g m) where + random = SamplerT (ReaderT uniformDouble01M) - uniform a b = Sampler (ReaderT $ uniformRM (a, b)) - normal m s = Sampler (ReaderT (MWC.normal m s)) - gamma shape scale = Sampler (ReaderT $ MWC.gamma shape scale) - beta a b = Sampler (ReaderT $ MWC.beta a b) + uniform a b = SamplerT (ReaderT $ uniformRM (a, b)) + normal m s = SamplerT (ReaderT (MWC.normal m s)) + gamma shape scale = SamplerT (ReaderT $ MWC.gamma shape scale) + beta a b = SamplerT (ReaderT $ MWC.beta a b) ... ``` @@ -235,34 +235,34 @@ We can unpack values from `SamplerIO a` using `sampler :: SamperIO a -> IO a`. S print =<< sampler random ``` -will print a sample from `random`. Similarly for other distributions. +will print a sample from `random`. Similarly for other distributions. -But note that `SamplerIO` only has a `MonadDistribution` instance. This means that `sampler sprinkler` does not type check and gives the type error: +But note that `SamplerIO` only has a `MonadDistribution` instance. This means that `sampler sprinkler` does not type check and gives the type error: ```haskell No instance for (MonadMeasure SamplerIO) ``` -This is to be expected. `SamplerIO` has no instance for `MonadFactor`. +This is to be expected. `SamplerIO` has no instance for `MonadFactor`. -### Weighted m +### WeightedT m -Summary of key info on `Weighted`: +Summary of key info on `WeightedT`: -- `Weighted :: (Type -> Type) -> (Type -> Type)` -- `instance MonadDistribution m => MonadMeasure (Weighted m)` -- `instance MonadFactor (Weighted m)` +- `WeightedT :: (Type -> Type) -> (Type -> Type)` +- `instance MonadDistribution m => MonadMeasure (WeightedT m)` +- `instance MonadFactor (WeightedT m)` ```haskell -newtype Weighted m a = Weighted (StateT (Log Double) m a) +newtype WeightedT m a = WeightedT (StateT (Log Double) m a) ``` -A **key** difference to `Enumerator` and `SamplerIO` is the **kind** of `Weighted`: it takes a type `m :: Type -> Type` as an argument. +A **key** difference to `Enumerator` and `SamplerIO` is the **kind** of `WeightedT`: it takes a type `m :: Type -> Type` as an argument. -`Weighted m` is isomorphic to: +`WeightedT m` is isomorphic to: ```haskell Log Double -> m (a, Log Double) @@ -271,87 +271,87 @@ Log Double -> m (a, Log Double) Some intuition for what this means comes from the `MonadFactor` instance: ```haskell -instance Monad m => MonadFactor (Weighted m) where - score w = Weighted (modify (* w)) +instance Monad m => MonadFactor (WeightedT m) where + score w = WeightedT (modify (* w)) ``` So if we write: ```haskell -example :: MonadDistribution m => Weighted m Bool +example :: MonadDistribution m => WeightedT m Bool example = do x <- bernoulli 0.5 score (if b then 1 else 0) return x ``` -then the result is that first, we draw a sample from a Bernoulli distribution from the **underlying** distribution `m`, and then multiply the state (which is a `Log Double`) by a number which depends on that sample. For convenience, we write `condition b = score (if b then 1 else 0)`. +then the result is that first, we draw a sample from a Bernoulli distribution from the **underlying** distribution `m`, and then multiply the state (which is a `Log Double`) by a number which depends on that sample. For convenience, we write `condition b = score (if b then 1 else 0)`. -To unpack from `Weighted m a`, we use: +To unpack from `WeightedT m a`, we use: ```haskell -weighted :: Weighted m a -> m (a, Log Double) -weighted (Weighted m) = runStateT m 1 +runWeightedT :: WeightedT m a -> m (a, Log Double) +runWeightedT (WeightedT m) = runStateT m 1 ``` -`Weighted m` is not an instance of `MonadDistribution`, but only as instance of `MonadFactor` (and that, only when `m` is an instance of `Monad`). However, since `StateT` is a monad transformer, there is a function `lift :: m Double -> Weighted m Double`. +`WeightedT m` is not an instance of `MonadDistribution`, but only as instance of `MonadFactor` (and that, only when `m` is an instance of `Monad`). However, since `StateT` is a monad transformer, there is a function `lift :: m Double -> WeightedT m Double`. -So if we take a `MonadDistribution` instance like `SamplerIO`, then `Weighted SamplerIO` is an instance of both `MonadDistribution` and `MonadFactor`. Which means it is an instance of `MonadMeasure`. +So if we take a `MonadDistribution` instance like `SamplerIO`, then `WeightedT SamplerIO` is an instance of both `MonadDistribution` and `MonadFactor`. Which means it is an instance of `MonadMeasure`. -So we can successfully write `(sampler . weighted) sprinkler` and get a program of type `IO (Bool, Log Double)`. When run, this will draw a sample from `sprinkler` along with an **unnormalized** density for that sample. +So we can successfully write `(sampler . runWeightedT) sprinkler` and get a program of type `IO (Bool, Log Double)`. When run, this will draw a sample from `sprinkler` along with an **unnormalized** density for that sample. -It's worth stopping here to remark on what's going on. What has happened is that the `m` in `forall m. MonadMeasure m => m Bool` has been *instantiated* as `Weighted SamplerIO`. This is an example of how the interpreters for inference can be composed in modular ways. +It's worth stopping here to remark on what's going on. What has happened is that the `m` in `forall m. MonadMeasure m => m Bool` has been *instantiated* as `WeightedT SamplerIO`. This is an example of how the interpreters for inference can be composed in modular ways. -Finally, there's a function +Finally, there's a function ```haskell -hoist :: (forall x. m x -> n x) -> Weighted m a -> Weighted n a +hoist :: (forall x. m x -> n x) -> WeightedT m a -> WeightedT n a ``` -This takes a natural transformation `m ~> n` and lifts it into a natural transformation `Weighted m ~> Weighted n`. Most of the inference transformers have an analogous `hoist` function. +This takes a natural transformation `m ~> n` and lifts it into a natural transformation `WeightedT m ~> WeightedT n`. Most of the inference transformers have an analogous `hoist` function. - +In a similar vein, `PopulationT` and `SequentialT` go together to handle Sequential Monte Carlo (SMC), and `TracedT` is associated with Markov Chain Monte Carlo (MCMC). --> -### Population +### PopulationT -Summary of key info on `Population`: +Summary of key info on `PopulationT`: -- `Population :: (Type -> Type) -> (Type -> Type)` -- `instance MonadDistribution m => instance MonadDistribution (Population m)` -- `instance MonadFactor m => instance MonadFactor (Population m)` +- `PopulationT :: (Type -> Type) -> (Type -> Type)` +- `instance MonadDistribution m => instance MonadDistribution (PopulationT m)` +- `instance MonadFactor m => instance MonadFactor (PopulationT m)` ```haskell -newtype Population m a = Population (Weighted (ListT m) a) +newtype PopulationT m a = PopulationT (WeightedT (ListT m) a) ``` So: ```haskell -Population m a ~ m [Log Double -> (a, Log Double)] +PopulationT m a ~ m [Log Double -> (a, Log Double)] ``` Note that while `ListT` isn't in general a valid monad transformer, we're not requiring it to be one here. -`Population` is used to represent a collection of particles (in the statistical sense), along with their weights. +`PopulationT` is used to represent a collection of particles (in the statistical sense), along with their weights. There are several useful functions associated with it: ```haskell -spawn :: Monad m => Int -> Population m () +spawn :: Monad m => Int -> PopulationT m () spawn n = fromWeightedList $ pure $ replicate n ((), 1 / fromIntegral n) ``` `spawn` spawns new particles. As an example: ```haskell -enumerator $ population (spawn 2) +enumerator $ runPopulationT (spawn 2) ``` gives @@ -360,16 +360,16 @@ gives [([((),0.5),((),0.5)],1.0)] ``` -Observe how here we have interpreted `(spawn 2)` as of type `Population Enumerator ()`. +Observe how here we have interpreted `(spawn 2)` as of type `PopulationT Enumerator ()`. -`resampleGeneric` takes a function to probabilistically select a set of indices from a vector, and makes a new population by selecting those indices. +`resampleGeneric` takes a function to probabilistically select a set of indices from a vector, and makes a new population by selecting those indices. ```haskell resampleGeneric :: MonadDistribution m => (V.Vector Double -> m [Int]) -> - Population m a -> - Population m a + PopulationT m a -> + PopulationT m a ``` `pushEvidence`, to quote the API docs, "normalizes the weights in the population, while at the same time incorporating the sum of the weights as a score in m." @@ -377,24 +377,24 @@ resampleGeneric :: ```haskell pushEvidence :: MonadFactor m => - Population m a -> - Population m a + PopulationT m a -> + PopulationT m a ``` -In other words, `pushEvidence` takes a `Population m a` where `m` is a `MonadFactor` instance. It takes the sum of the weights, divides the weights by it, and then factors by the sum in `m`. +In other words, `pushEvidence` takes a `PopulationT m a` where `m` is a `MonadFactor` instance. It takes the sum of the weights, divides the weights by it, and then factors by the sum in `m`. -### Sequential +### SequentialT -Summary of key info on `Sequential`: +Summary of key info on `SequentialT`: -- `Sequential :: (Type -> Type) -> (Type -> Type)` -- `instance MonadDistribution m => instance MonadDistribution (Sequential m)` -- `instance MonadFactor m => instance MonadFactor (Sequential m)` +- `SequentialT :: (Type -> Type) -> (Type -> Type)` +- `instance MonadDistribution m => instance MonadDistribution (SequentialT m)` +- `instance MonadFactor m => instance MonadFactor (SequentialT m)` ```haskell -newtype Sequential m a = - Sequential {runSequential :: Coroutine (Await ()) m a} +newtype SequentialT m a = + SequentialT {runSequentialT :: Coroutine (Await ()) m a} ``` This is a wrapper for the `Coroutine` type applied to the `Await` constructor from `Control.Monad.Coroutine`, which is defined thus: @@ -409,11 +409,11 @@ newtype Await x y = Await (x -> y) Unpacking that: -```haskell -Sequential m a ~ m (Either (() -> Sequential m a) a) +```haskell +SequentialT m a ~ m (Either (() -> SequentialT m a) a) ``` -As usual, `m` is going to be some other probability monad, so understand `Sequential m a` as representing a program which, after making a random choice or doing conditioning, we either obtain an `a` value, or a paused computation, which when resumed gets us back to a new `Sequential m a`. +As usual, `m` is going to be some other probability monad, so understand `SequentialT m a` as representing a program which, after making a random choice or doing conditioning, we either obtain an `a` value, or a paused computation, which when resumed gets us back to a new `SequentialT m a`. (For more on coroutines, see the final article in: https://themonadreader.files.wordpress.com/2011/10/issue19.pdf.) @@ -439,7 +439,7 @@ The `MonadFactor` instance has less trivial content: ```haskell -- | Execution is 'suspend'ed after each 'score'. -instance MonadFactor m => MonadFactor (Sequential m) where +instance MonadFactor m => MonadFactor (SequentialT m) where score w = lift (score w) >> suspend ``` @@ -447,16 +447,16 @@ First you take a `score` in the underlying `MonadFactor` instance, and then you ```haskell -- | A point where the computation is paused. -suspend :: Monad m => Sequential m () -suspend = Sequential (Coroutine (return (Left (Await return)))) +suspend :: Monad m => SequentialT m () +suspend = SequentialT (Coroutine (return (Left (Await return)))) ``` We can move to the next suspension point with: ```haskell -advance :: Monad m => Sequential m a -> Sequential m a -advance = Sequential . bounce extract . runSequential +advance :: Monad m => SequentialT m a -> SequentialT m a +advance = SequentialT . bounce extract . runSequentialT ``` @@ -464,36 +464,36 @@ and move through all with: ```haskell -- | Remove the remaining suspension points. -finish :: Monad m => Sequential m a -> m a -finish = pogoStick extract . runSequential +finish :: Monad m => SequentialT m a -> m a +finish = pogoStick extract . runSequentialT ``` But most importantly, we can apply a natural transformation over the underlying monad `m` to only the current suspension. ```haskell -hoistFirst :: (forall x. m x -> m x) -> Sequential m a -> Sequential m a -hoistFirst f = Sequential . Coroutine . f . resume . runSequential +hoistFirst :: (forall x. m x -> m x) -> SequentialT m a -> SequentialT m a +hoistFirst f = SequentialT . Coroutine . f . resume . runSequentialT ``` -When `m` is `Population n` for some other `n`, then `resampleGeneric` gives us one example of the natural transformation we want. In other words, operating in `Sequential (Population n)` works, and not only works but does something statistically interesting: particle filtering (aka SMC). +When `m` is `PopulationT n` for some other `n`, then `resampleGeneric` gives us one example of the natural transformation we want. In other words, operating in `SequentialT (PopulationT n)` works, and not only works but does something statistically interesting: particle filtering (aka SMC). -**Note**: the running of `Sequential`, i.e. getting from `Sequential m a` to `m a` is surprisingly subtle, and there are many incorrect ways to do it, such as plain folds of the recursive structure. These can result in a semantics in which the transformation gets applied an exponentially large number of times. +**Note**: the running of `SequentialT`, i.e. getting from `SequentialT m a` to `m a` is surprisingly subtle, and there are many incorrect ways to do it, such as plain folds of the recursive structure. These can result in a semantics in which the transformation gets applied an exponentially large number of times. -### Density +### DensityT -Summary of key info on `Density`: +Summary of key info on `DensityT`: -- `Density :: (Type -> Type) -> (Type -> Type)` -- `instance MonadDistribution (Density m)` +- `DensityT :: (Type -> Type) -> (Type -> Type)` +- `instance MonadDistribution (DensityT m)` - **No** instance for `MonadFactor` A *trace* of a program of type `MonadDistribution m => m a` is an execution of the program, so a choice for each of the random values. Recall that `random` underlies all of the random values in a program, so a trace for a program is fully specified by a list of `Double`s, giving the value of each call to `random`. -With this in mind, a `Density m a` is an interpretation of a probabilistic program as a function from a trace to the *density* of that execution of the program. +With this in mind, a `DensityT m a` is an interpretation of a probabilistic program as a function from a trace to the *density* of that execution of the program. -Monad-bayes offers two implementations, in `Control.Monad.Bayes.Density.State` and `Control.Monad.Bayes.Density.Free`. The first is slow but easy to understand, the second is more sophisticated, but faster. +Monad-bayes offers two implementations, in `Control.Monad.Bayes.DensityT.State` and `Control.Monad.Bayes.DensityT.Free`. The first is slow but easy to understand, the second is more sophisticated, but faster. The former is relatively straightforward: the `MonadDistribution` instance implements `random` as `get`ting the trace (using `get` from `MonadState`), using (and removing) the first element (`put` from `MonadState`), and writing that element to the output (using `tell` from `MonadWriter`). If the trace is empty, the `random` from the underlying monad is used, but the result is still written with `tell`. @@ -501,28 +501,28 @@ The latter is best understood if you're familiar with the standard use of a free ```haskell newtype SamF a = Random (Double -> a) -newtype Density m a = - Density {density :: FT SamF m a} - -instance Monad m => MonadDistribution (Density m) where - random = Density $ liftF (Random id) +newtype DensityT m a = + DensityT {getDensityT :: FT SamF m a} + +instance Monad m => MonadDistribution (DensityT m) where + random = DensityT $ liftF (Random id) ``` -The monad-bayes implementation uses a more efficient implementation of `FreeT`, namely `FT` from the `free` package, known as the *Church transformed Free monad*. This is a technique explained in https://begriffs.com/posts/2016-02-04-difference-lists-and-codennsity.html. But that only changes the operational semantics - performance aside, it works just the same as the standard `FreeT` datatype. +The monad-bayes implementation uses a more efficient implementation of `FreeT`, namely `FT` from the `free` package, known as the *Church transformed Free monad*. This is a technique explained in https://begriffs.com/posts/2016-02-04-difference-lists-and-codennsity.html. But that only changes the operational semantics - performance aside, it works just the same as the standard `FreeT` datatype. If you unpack the definition, you get: ```haskell -Density m a ~ m (Either a (Double -> (Density m a))) +DensityT m a ~ m (Either a (Double -> (DensityT m a))) ``` Since `FreeT` is a transformer, we can use `lift` to get a `MonadDistribution` instance. -`density` is then defined using the folding pattern `iterFT`, which interprets `SamF` in the appropriate way: +`runDensityT` is then defined using the folding pattern `iterFT`, which interprets `SamF` in the appropriate way: ```haskell -density :: MonadDistribution m => [Double] -> Density m a -> m (a, [Double]) -density randomness (Density m) = +runDensityT :: MonadDistribution m => [Double] -> DensityT m a -> m (a, [Double]) +runDensityT randomness (DensityT m) = runWriterT $ evalStateT (iterTM f $ hoistFT lift m) randomness where f (Random k) = do @@ -545,15 +545,15 @@ This takes a list of `Double`s (a representation of a trace), and a probabilisti -### Traced +### TracedT -Summary of key info on `Traced`: +Summary of key info on `TracedT`: -- `Traced :: (Type -> Type) -> (Type -> Type)` -- `instance MonadDistribution m => MonadDistribution (Traced m)` -- `instance MonadFactor m => MonadFactor (Traced m)` +- `TracedT :: (Type -> Type) -> (Type -> Type)` +- `instance MonadDistribution m => MonadDistribution (TracedT m)` +- `instance MonadFactor m => MonadFactor (TracedT m)` -`Traced m` is actually several related interpretations, each built on top of `Density`. These range in complexity. +`TracedT m` is actually several related interpretations, each built on top of `DensityT`. These range in complexity. @@ -568,24 +568,24 @@ It's convenient to specify a trace not just as a `[Double]` but also with the re ```haskell data Trace a = Trace - { + { variables :: [Double], output :: a, - density :: Log Double + runDensityT :: Log Double } ``` -We also need a specification of the probabilistic program in question, free of any particular interpretation. That is precisely what `Density` is for. +We also need a specification of the probabilistic program in question, free of any particular interpretation. That is precisely what `DensityT` is for. -The simplest version of `Traced` is in `Control.Monad.Bayes.Traced.Basic` +The simplest version of `TracedT` is in `Control.Monad.Bayes.TracedT.Basic` ```haskell -Traced m a ~ (Density Identity a, Log Double), m (Trace a)) +TracedT m a ~ (DensityT Identity a, Log Double), m (Trace a)) ``` -A `Traced` interpretation of a model is a particular run of the model with its corresponding probability, alongside a distribution over `Trace` info, which records: the value of each call to `random`, the value of the final output, and the density of this program trace. +A `TracedT` interpretation of a model is a particular run of the model with its corresponding probability, alongside a distribution over `Trace` info, which records: the value of each call to `random`, the value of the final output, and the density of this program trace. -This machinery allows us to implement MCMC (see inference methods below for details). +This machinery allows us to implement MCMC (see inference methods below for details). ### Integrator @@ -602,8 +602,8 @@ newtype Integrator a = Integrator {getCont :: Cont Double a} This `MonadDistribution` instance interprets a probabilistic program as a numerical integrator. For a nice explanation, see [this blog post](https://jtobin.io/giry-monad-implementation). -`Integrator a` is isomorphic to `(a -> Double) -> Double`. -A program `model` of type `Integrator a` will take a function `f` and calculate $E_{p}[f] = \int f(x)*p(x)$ where $p$ is the density of `model`. +`Integrator a` is isomorphic to `(a -> Double) -> Double`. +A program `model` of type `Integrator a` will take a function `f` and calculate $E_{p}[f] = \int f(x)*p(x)$ where $p$ is the density of `model`. The integral for the expectation is performed by quadrature, using the tanh-sinh approach. For example, `random :: Integrator Double` is the program which takes a function `f` and integrates `f` over the $(0,1)$ range. @@ -624,9 +624,9 @@ example = replicateM 100 $ do return x ``` -Doing `enumerator example` will create a list of $2^{100}$ entries, all but one of which have $0$ mass. (See below for a way to perform this inference efficiently). +Doing `enumerator example` will create a list of $2^{100}$ entries, all but one of which have $0$ mass. (See below for a way to perform this inference efficiently). -The main purpose of `Enumerator` is didactic, as a way to understand simple discrete distributions in full. In addition, you can use it in concert with transformers like `Weighted`, to get a sense of how they work. For example, consider: +The main purpose of `Enumerator` is didactic, as a way to understand simple discrete distributions in full. In addition, you can use it in concert with transformers like `WeightedT`, to get a sense of how they work. For example, consider: ```haskell example = do @@ -635,13 +635,13 @@ example = do return x ``` -`(enumerator . weighted) example` gives `[((False,0.0),0.5),((True,1.0),0.5)]`. This is quite edifying for understanding `(sampler . weighted) example`. What it says is that there are precisely two ways the program will run, each with equal probability: either you get `False` with weight `0.0` or `True` with weight `1.0`. +`(enumerator . runWeightedT) example` gives `[((False,0.0),0.5),((True,1.0),0.5)]`. This is quite edifying for understanding `(sampler . runWeightedT) example`. What it says is that there are precisely two ways the program will run, each with equal probability: either you get `False` with weight `0.0` or `True` with weight `1.0`. ### Quadrature As described on the section on `Integrator`, we can interpret our probabilistic program of type `MonadDistribution m => m a` as having concrete type `Integrator a`. This views our program as an integrator, allowing us to calculate expectations, probabilities and so on via quadrature (i.e. numerical approximation of an integral). -This can also handle programs of type `MonadMeasure m => m a`, that is, programs with `factor` statements. For these cases, a function `normalize :: Weighted Integrator a -> Integrator a` is employed. For example, +This can also handle programs of type `MonadMeasure m => m a`, that is, programs with `factor` statements. For these cases, a function `normalize :: WeightedT Integrator a -> Integrator a` is employed. For example, ```haskell model :: MonadMeasure m => m Double @@ -652,15 +652,15 @@ model = do return var ``` -is really an unnormalized measure, rather than a probability distribution. `normalize` views it as of type `Weighted Integrator Double`, which is isomorphic to `(Double -> (Double, Log Double) -> Double)`. This can be used to compute the normalization constant, and divide the integrator's output by it, all within `Integrator`. +is really an unnormalized measure, rather than a probability distribution. `normalize` views it as of type `WeightedT Integrator Double`, which is isomorphic to `(Double -> (Double, Log Double) -> Double)`. This can be used to compute the normalization constant, and divide the integrator's output by it, all within `Integrator`. ### Independent forward sampling -For any program of type `p = MonadDistribution m => m a`, we may do `sampler p` or `runST $ sampleSTfixed p`. Note that if there are any calls to `factor` in the program, then it cannot have type `MonadDistribution m`. +For any program of type `p = MonadDistribution m => m a`, we may do `sampler p` or `runST $ sampleSTfixed p`. Note that if there are any calls to `factor` in the program, then it cannot have type `MonadDistribution m`. ### Independent weighted sampling -Consider +Consider ```haskell example :: MonadMeasure m => m Bool @@ -670,15 +670,15 @@ example = do return x ``` -`(weighted . sampler) example :: IO (Bool, Log Double)` returns a tuple of a truth value and a probability mass (or more generally density). How does this work? Types are clarifying: +`(runWeightedT . sampler) example :: IO (Bool, Log Double)` returns a tuple of a truth value and a probability mass (or more generally density). How does this work? Types are clarifying: ```haskell -run = +run = (sampler :: SamplerIO (a, Log Double) -> IO (a, Log Double) ) - . (weighted :: Weighted SamplerIO a -> SamplerIO (a, Log Double) + . (runWeightedT :: WeightedT SamplerIO a -> SamplerIO (a, Log Double) ``` -In other words, the program is being interpreted in the `Weighted SamplerIO` instance of `MonadMeasure`. +In other words, the program is being interpreted in the `WeightedT SamplerIO` instance of `MonadMeasure`. ### Metropolis Hastings MCMC @@ -688,7 +688,7 @@ The version of MCMC in monad-bayes performs its random walk on program traces of A single step in this chain (in Metropolis Hasting MCMC) looks like this: ```haskell -mhTrans :: MonadDistribution m => (Weighted (State.Density m)) a -> Trace a -> m (Trace a) +mhTrans :: MonadDistribution m => (WeightedT (State.DensityT m)) a -> Trace a -> m (Trace a) mhTrans m t@Trace {variables = us, probDensity = p} = do let n = length us us' <- do @@ -697,13 +697,13 @@ mhTrans m t@Trace {variables = us, probDensity = p} = do case splitAt i us of (xs, _ : ys) -> return $ xs ++ (u' : ys) _ -> error "impossible" - ((b, q), vs) <- State.density (weighted m) us' + ((b, q), vs) <- State.runDensityT (runWeightedT m) us' let ratio = (exp . ln) $ min 1 (q * fromIntegral n / (p * fromIntegral (length vs))) accept <- bernoulli ratio return $ if accept then Trace vs b q else t ``` -Our probabilistic program is interpreted in the type `Weighted (Density m) a`, which is an instance of `MonadMeasure`. We use this to define our kernel on traces. We begin by perturbing the list of doubles contained in the trace by selecting a random position in the list and resampling there. We could do this *proposal* in a variety of ways, but here, we do so by choosing a double from the list at random and resampling it (hence, *single site* trace MCMC). We then run the program on this new list of doubles; `((b,q), vs)` is the outcome, probability, and result of all calls to `random`, respectively (recalling that the list of doubles may be shorter than the number of calls to `random`). The value of these is probabilistic in the underlying monad `m`. We then use the MH criterion to decide whether to accept the new list of doubles as our trace. +Our probabilistic program is interpreted in the type `WeightedT (DensityT m) a`, which is an instance of `MonadMeasure`. We use this to define our kernel on traces. We begin by perturbing the list of doubles contained in the trace by selecting a random position in the list and resampling there. We could do this *proposal* in a variety of ways, but here, we do so by choosing a double from the list at random and resampling it (hence, *single site* trace MCMC). We then run the program on this new list of doubles; `((b,q), vs)` is the outcome, probability, and result of all calls to `random`, respectively (recalling that the list of doubles may be shorter than the number of calls to `random`). The value of these is probabilistic in the underlying monad `m`. We then use the MH criterion to decide whether to accept the new list of doubles as our trace. MH is then easily defined as taking steps with this kernel, in the usual fashion. Note that it works for any probabilistic program whatsoever. @@ -711,7 +711,7 @@ MH is then easily defined as taking steps with this kernel, in the usual fashion ### Sequential Importance Sampling -This is provided by +This is provided by ```haskell sequentially :: @@ -720,7 +720,7 @@ sequentially :: (forall x. m x -> m x) -> -- | number of time steps Int -> - Sequential m a -> + SequentialT m a -> m a sequentially f k = finish . composeCopies k (advance . hoistFirst f) ``` @@ -738,7 +738,7 @@ example = replicateM 100 $ do Naive enumeration, as in `enumerator example` is enormously and needlessly inefficient, because it will create a $2^{100}$ size list of possible values. What we'd like to do is to throw away values of `x` that are `False` at each condition statement, rather than carrying them along forever. -Suppose we have a function `removeZeros :: Enumerator a -> Enumerator a`, which removes values of the distribution with $0$ mass from `Enumerator`. We can then write `enumerator $ sequentially removeZeros 100 $ model` to run `removeZeros` at each of the 100 `condition` statements, making the algorithm run quickly. +Suppose we have a function `removeZeros :: Enumerator a -> Enumerator a`, which removes values of the distribution with $0$ mass from `Enumerator`. We can then write `enumerator $ sequentially removeZeros 100 $ model` to run `removeZeros` at each of the 100 `condition` statements, making the algorithm run quickly. ### Sequential Monte Carlo @@ -748,14 +748,14 @@ Sequential importance resampling works by doing `sequentially` with a resampler smc :: Monad m => -- | resampler - (forall x. Population m x -> Population m x) -> + (forall x. PopulationT m x -> PopulationT m x) -> -- | number of timesteps Int -> -- | population size Int -> -- | model - Sequential (Population m) a -> - Population m a + SequentialT (PopulationT m) a -> + PopulationT m a smc resampler k n = sequentially resampler k . Seq.hoistFirst (spawn n >>) ``` @@ -765,7 +765,7 @@ Quoting the monad-bayes paper: "PMMH is only applicable to models with a specific structure, namely the probabilistic program needs to decompose to a prior over the global parameters `m param` and the rest of the model `param -> m a`. Combining these using >>= would yield the complete model of type `m a`." -"The idea is to do MH on the parameters of the model. Recall that for MH we need to compute the likelihood for the particular values of parameters but that involves integrating over the remaining random variables in the model which is intractable. Fortunately to obtain valid MH it is sufficient to have an unbiased estimator for the likelihood which is produced by a single sample from `Weighted`. MH with such an estimator is referred to as pseudo-marginal MH. If instead of taking a single weight from `Weighted` we take the sum of weights from `Population` we obtain an unbiased estimator with lower variance. In particular if such a `Population` is a result of smc the resulting algorithm is known as PMMH." +"The idea is to do MH on the parameters of the model. Recall that for MH we need to compute the likelihood for the particular values of parameters but that involves integrating over the remaining random variables in the model which is intractable. Fortunately to obtain valid MH it is sufficient to have an unbiased estimator for the likelihood which is produced by a single sample from `WeightedT`. MH with such an estimator is referred to as pseudo-marginal MH. If instead of taking a single weight from `WeightedT` we take the sum of weights from `PopulationT` we obtain an unbiased estimator with lower variance. In particular if such a `PopulationT` is a result of smc the resulting algorithm is known as PMMH." Because of the modularity of monad-bayes, the implementation is remarkably simple, and directly linked to the algorithm: @@ -774,7 +774,7 @@ pmmh mcmcConf smcConf param model = mcmc mcmcConf ( param - >>= population + >>= runPopulationT . pushEvidence . Pop.hoist lift . smc smcConf @@ -783,29 +783,29 @@ pmmh mcmcConf smcConf param model = ``` -There's a lot to unpack here. Here's the definition with more types. To shorten the signatures, the synonyms: `T = Traced, S = Sequential, P = Population` are used: +There's a lot to unpack here. Here's the definition with more types. To shorten the signatures, the synonyms: `T = TracedT, S = SequentialT, P = PopulationT` are used: ```haskell pmmh :: MonadDistribution m => MCMCConfig -> - SMCConfig (Weighted m) -> - Traced (Weighted m) a1 -> - (a1 -> Sequential (Population (Weighted m)) a2) -> + SMCConfig (WeightedT m) -> + TracedT (WeightedT m) a1 -> + (a1 -> SequentialT (PopulationT (WeightedT m)) a2) -> m [[(a2, Log Double)]] pmmh mcmcConf smcConf param model = (mcmc mcmcConf :: T m [(a, Log Double)] -> m [[(a, Log Double)]]) - ((param :: T m b) >>= - (population :: P (T m) a -> T m [(a, Log Double)]) - . (pushEvidence :: P (T m) a -> P (T m) a) - . Pop.hoist (lift :: forall x. m x -> T m x) - . (smc smcConf :: S (P m) a -> P m a) + ((param :: T m b) >>= + (runPopulationT :: P (T m) a -> T m [(a, Log Double)]) + . (pushEvidence :: P (T m) a -> P (T m) a) + . Pop.hoist (lift :: forall x. m x -> T m x) + . (smc smcConf :: S (P m) a -> P m a) . (model :: b -> S (P m) a)) ``` -(Note that this uses the version of `mcmc` that uses the `Traced` representation from `Control.Monad.Bayes.Traced.Static`.) +(Note that this uses the version of `mcmc` that uses the `TracedT` representation from `Control.Monad.Bayes.TracedT.Static`.) -To understand this, note that the outer algorithm is just `mcmc`. But the probabilistic program that we pass to `mcmc` does the following: run `param` to get values for the parameters and then pass these to `model`. Then run `n` steps of SMC with `k` particles. Then lift the underlying monad `m` into `Traced m`. Then calculate the sum of the weights of the particles and `factor` on this (this is what `pushEvidence` does). This `factor` takes place in `Traced m`, so it is "visible" to `mcmc`. +To understand this, note that the outer algorithm is just `mcmc`. But the probabilistic program that we pass to `mcmc` does the following: run `param` to get values for the parameters and then pass these to `model`. Then run `n` steps of SMC with `k` particles. Then lift the underlying monad `m` into `TracedT m`. Then calculate the sum of the weights of the particles and `factor` on this (this is what `pushEvidence` does). This `factor` takes place in `TracedT m`, so it is "visible" to `mcmc`. ### Resample-Move Sequential Monte Carlo @@ -816,7 +816,7 @@ The paper introducing monad-bayes has this to say about resample-move SMC: "A common problem with particle filters is that of particle degeneracy, where after resampling many particles are the same, effectively reducing the sample size. One way to ameliorate this problem is to introduce rejuvenation moves, where after each resampling we apply a number of MCMC transitions to each particle independently, thus spreading them around the space. If we use an MCMC kernel that preserves the target distribution at a given step, the resulting algorithm is correct" -monad-bayes provides three versions of RMSMC, each of which uses one of the three `Traced` implementations respectively. Here is the simplest, which I have annotated with types. To shorten the signatures, the synonyms: `T = Traced, S = Sequential, P = Population` are used: +monad-bayes provides three versions of RMSMC, each of which uses one of the three `TracedT` implementations respectively. Here is the simplest, which I have annotated with types. To shorten the signatures, the synonyms: `T = TracedT, S = SequentialT, P = PopulationT` are used: ```haskell @@ -829,19 +829,19 @@ rmsmcBasic :: P m a rmsmc (MCMCConfig {..}) (SMCConfig {..}) = (TrBas.marginal :: T (P m) a -> P m a ) - . sequentially + . sequentially ( (composeCopies numMCMCSteps TrBas.mhStep :: T (P m) a -> T (P m) a ) - . TrBas.hoist (resampleSystematic :: P m a -> P m a ) ) + . TrBas.hoist (resampleSystematic :: P m a -> P m a ) ) numSteps . (S.hoistFirst (TrBas.hoist (spawn numParticles >>)) :: S (T (P m)) a -> S (T (P m)) a )) ``` -What is this doing? Recall that `S m a` represents an `m` of coroutines over `a`. Recall that `Traced m a` represents an `m` of traced computations of `a`. Recall that `P m a` represents an `m` of a list of `a`s. +What is this doing? Recall that `S m a` represents an `m` of coroutines over `a`. Recall that `TracedT m a` represents an `m` of traced computations of `a`. Recall that `P m a` represents an `m` of a list of `a`s. This means that an `S (T (P m)) a` is a program "interpreted as a population of traced coroutines". The paper adds that this "allows us to apply MH transitions to partially executed coroutines, which is exactly what we require for the rejuvenation steps." -So the algorithm works by creating `n` particles, and at each of the first `k` calls to `factor`, first resampling the population and then for each particle in the population, doing an MH-MCMC walk for `t` steps to update it. +So the algorithm works by creating `n` particles, and at each of the first `k` calls to `factor`, first resampling the population and then for each particle in the population, doing an MH-MCMC walk for `t` steps to update it. ### Sequential Monte Carlo Squared @@ -856,21 +856,19 @@ There is one slight complication here, related to the fact that the MTL style ef Basic instances of `MonadDistribution`: - `SamplerIO` -- `Traced m` -- `Sequential m` +- `TracedT m` +- `SequentialT m` Basic instances of `MonadFactor`: -- `Weighted m` -- `Sequential m` +- `WeightedT m` +- `SequentialT m` Basic instances of `MonadMeasure`: - `Enumerator` (discrete distributions only) Composed instances of `MonadMeasure`: -- `Weighted SamplerIO` (used in weighted sampling) -- `Traced (Weighted SamplerIO)` (used in MCMC) -- `Sequential (Population (Weighted SamplerIO)` (used in SMC) +- `WeightedT SamplerIO` (used in weighted sampling) +- `TracedT (WeightedT SamplerIO)` (used in MCMC) +- `SequentialT (PopulationT (WeightedT SamplerIO)` (used in SMC) - more --> - - diff --git a/models/Helper.hs b/models/Helper.hs index 827c0614..3741e4ce 100644 --- a/models/Helper.hs +++ b/models/Helper.hs @@ -55,7 +55,7 @@ runAlg model alg = SMC -> let n = 100 (k, m) = getModel model - in show <$> population (smc SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic} m) + in show <$> runPopulationT (smc SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic} m) MH -> let t = 100 (_, m) = getModel model @@ -64,7 +64,7 @@ runAlg model alg = let n = 10 t = 1 (k, m) = getModel model - in show <$> population (rmsmcBasic MCMCConfig {numMCMCSteps = t, numBurnIn = 0, proposal = SingleSiteMH} (SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic}) m) + in show <$> runPopulationT (rmsmcBasic MCMCConfig {numMCMCSteps = t, numBurnIn = 0, proposal = SingleSiteMH} (SMCConfig {numSteps = k, numParticles = n, resampler = resampleSystematic}) m) runAlgFixed :: Model -> Alg -> IO String runAlgFixed model alg = sampleIOfixed $ runAlg model alg diff --git a/models/NonlinearSSM/Algorithms.hs b/models/NonlinearSSM/Algorithms.hs index 4a893495..5b053781 100644 --- a/models/NonlinearSSM/Algorithms.hs +++ b/models/NonlinearSSM/Algorithms.hs @@ -23,24 +23,24 @@ t = 5 -- FIXME refactor such that it can be reused in ssm benchmark runAlgFixed :: (MonadDistribution m) => SSMData -> Alg -> m String -runAlgFixed ys SMC = fmap show $ population $ smc SMCConfig {numSteps = t, numParticles = 10, resampler = resampleMultinomial} (param >>= model ys) +runAlgFixed ys SMC = fmap show $ runPopulationT $ smc SMCConfig {numSteps = t, numParticles = 10, resampler = resampleMultinomial} (param >>= model ys) runAlgFixed ys RMSMC = fmap show $ - population $ + runPopulationT $ rmsmc MCMCConfig {numMCMCSteps = 10, numBurnIn = 0, proposal = SingleSiteMH} SMCConfig {numSteps = t, numParticles = 10, resampler = resampleSystematic} (param >>= model ys) runAlgFixed ys RMSMCBasic = fmap show $ - population $ + runPopulationT $ rmsmcBasic MCMCConfig {numMCMCSteps = 10, numBurnIn = 0, proposal = SingleSiteMH} SMCConfig {numSteps = t, numParticles = 10, resampler = resampleSystematic} (param >>= model ys) runAlgFixed ys RMSMCDynamic = fmap show $ - population $ + runPopulationT $ rmsmcDynamic MCMCConfig {numMCMCSteps = 10, numBurnIn = 0, proposal = SingleSiteMH} SMCConfig {numSteps = t, numParticles = 10, resampler = resampleSystematic} @@ -53,4 +53,4 @@ runAlgFixed ys PMMH = SMCConfig {numSteps = t, numParticles = 3, resampler = resampleSystematic} param (model ys) -runAlgFixed ys SMC2 = fmap show $ population $ smc2 t 3 2 1 param (model ys) +runAlgFixed ys SMC2 = fmap show $ runPopulationT $ smc2 t 3 2 1 param (model ys) diff --git a/notebooks/Implementation.ipynb b/notebooks/Implementation.ipynb index 9986a2b6..3dc62d9f 100644 --- a/notebooks/Implementation.ipynb +++ b/notebooks/Implementation.ipynb @@ -140,9 +140,9 @@ "\n", "the typeclass\n", "\n", - "`Weighted`\n", - "`Population`\n", - "`Sequential`" + "`WeightedT`\n", + "`PopulationT`\n", + "`SequentialT`" ] }, { @@ -178,7 +178,7 @@ } ], "source": [ - "sampleIO $ runWeighted program " + "sampleIO $ runWeightedT program " ] }, { @@ -197,7 +197,7 @@ " return (x+y) \n", " \n", " \n", - "y <- sampleIO $ runPopulation $ do \n", + "y <- sampleIO $ runPopulationT $ do \n", " x <- resampleMultinomial $ (spawn 10 >> d)\n", " factor (normalPdf 0 1 x)\n", " return x\n" @@ -317,7 +317,7 @@ } ], "source": [ - "sampleIO $ replicateM 10 $ runWeighted $ proper $ (spawn 2 >> random)" + "sampleIO $ replicateM 10 $ runWeightedT $ proper $ (spawn 2 >> random)" ] }, { diff --git a/notebooks/examples/Diagrams.ipynb b/notebooks/examples/Diagrams.ipynb index 365f6663..db392867 100644 --- a/notebooks/examples/Diagrams.ipynb +++ b/notebooks/examples/Diagrams.ipynb @@ -110,8 +110,8 @@ "metadata": {}, "outputs": [], "source": [ - "getTrace :: Weighted (Density SamplerIO) a -> SamplerIO ([Double], Log Double)\n", - "getTrace m = fmap (\\((x,y),z) -> (z,y) ) $ runWriterT $ weighted $ Weighted.hoist (WriterT . density []) m" + "getTrace :: WeightedT (DensityT SamplerIO) a -> SamplerIO ([Double], Log Double)\n", + "getTrace m = fmap (\\((x,y),z) -> (z,y) ) $ runWriterT $ runWeightedT $ Weighted.hoist (WriterT . runDensityT []) m" ] }, { diff --git a/notebooks/examples/Histogram.ipynb b/notebooks/examples/Histogram.ipynb index 7942b2ed..e6c84451 100644 --- a/notebooks/examples/Histogram.ipynb +++ b/notebooks/examples/Histogram.ipynb @@ -43,7 +43,7 @@ " return y\n", " \n", "distributionOverWeightedSamples :: Distribution (Double, Log Double)\n", - "distributionOverWeightedSamples = weighted unnormalizedDistribution\n", + "distributionOverWeightedSamples = runWeightedT unnormalizedDistribution\n", " \n", "iidDistribution :: Distribution [(Double, Log Double)]\n", "iidDistribution = replicateM 100000 distributionOverWeightedSamples\n", diff --git a/notebooks/examples/Parsing.ipynb b/notebooks/examples/Parsing.ipynb index 464fd5a4..1fcac334 100644 --- a/notebooks/examples/Parsing.ipynb +++ b/notebooks/examples/Parsing.ipynb @@ -112,7 +112,7 @@ "\n", "runWordParser w = do\n", " x <- sampler \n", - " . population \n", + " . runPopulationT \n", " . smc SMCConfig {numSteps = 5, numParticles = 3000, resampler = resampleMultinomial} \n", " $ run word w\n", " pPrintCustom $ toEmpiricalWeighted x\n", diff --git a/notebooks/examples/Streaming.ipynb b/notebooks/examples/Streaming.ipynb index 1835f184..fab37d41 100644 --- a/notebooks/examples/Streaming.ipynb +++ b/notebooks/examples/Streaming.ipynb @@ -1781,7 +1781,7 @@ } ], "source": [ - "particles <- sampleIOfixed $ runPopulation $ \n", + "particles <- sampleIOfixed $ runPopulationT $ \n", " smc SMCConfig {\n", " numSteps = 100, \n", " numParticles = 1000, \n", diff --git a/notebooks/tutorials/AdvancedSampling.ipynb b/notebooks/tutorials/AdvancedSampling.ipynb index 1c7730c0..fc2f178a 100644 --- a/notebooks/tutorials/AdvancedSampling.ipynb +++ b/notebooks/tutorials/AdvancedSampling.ipynb @@ -33,7 +33,7 @@ "import Control.Monad.Bayes.Class \n", "import Control.Monad.Bayes.Traced\n", "import Control.Monad.Bayes.Weighted\n", - "import Graphics.Vega.VegaLite hiding (density)\n", + "import Graphics.Vega.VegaLite hiding (runDensityT)\n", "import qualified Graphics.Vega.VegaLite as VL\n", "import IHaskell.Display.Hvega (vlShow)\n", "\n", @@ -853,7 +853,7 @@ } ], "source": [ - "sSMC <- sampleIOfixed $ runPopulation $ smc \n", + "sSMC <- sampleIOfixed $ runPopulationT $ smc \n", " SMCConfig {numSteps = numPoints, numParticles = 10000, resampler = resampleSystematic} \n", " $ P.toListM $ model observations\n" ] @@ -51121,7 +51121,7 @@ } ], "source": [ - "sRMSMC <- sampleIOfixed $ runPopulation $ rmsmcBasic\n", + "sRMSMC <- sampleIOfixed $ runPopulationT $ rmsmcBasic\n", " MCMCConfig {numMCMCSteps = 50, numBurnIn = 0, proposal = SingleSiteMH}\n", " SMCConfig {numSteps = numPoints, numParticles = 10, resampler = resampleSystematic} \n", " $ P.toListM $ model observations" diff --git a/notebooks/tutorials/Bayesian.ipynb b/notebooks/tutorials/Bayesian.ipynb index 57e933a9..311d9612 100644 --- a/notebooks/tutorials/Bayesian.ipynb +++ b/notebooks/tutorials/Bayesian.ipynb @@ -630,7 +630,7 @@ ], "source": [ "observations1 = [-3..3]\n", - "sampler $ toEmpiricalWeighted <$> replicateM 100000 (weighted $ posterior model observations1)" + "sampler $ toEmpiricalWeighted <$> replicateM 100000 (runWeightedT $ posterior model observations1)" ] }, { @@ -659,7 +659,7 @@ ], "source": [ "observations2 = [-2.3, -2, 0, 0.3, 2, 2.3]\n", - "sampler $ toEmpiricalWeighted <$> replicateM 100000 (weighted $ posterior model observations2)" + "sampler $ toEmpiricalWeighted <$> replicateM 100000 (runWeightedT $ posterior model observations2)" ] }, { @@ -707,7 +707,7 @@ "source": [ "sampler $ \n", " plot . fmap (bimap (T.pack . show) (ln . exp)) . toEmpiricalWeighted <$> \n", - " replicateM 1000 (weighted $ posterior model observations2)" + " replicateM 1000 (runWeightedT $ posterior model observations2)" ] }, { @@ -1166,7 +1166,7 @@ } ], "source": [ - "sampler $ plot . histogramToList . histogram 100 <$> replicateM 100000 (weighted $ posteriorPredictive model observations2)" + "sampler $ plot . histogramToList . histogram 100 <$> replicateM 100000 (runWeightedT $ posteriorPredictive model observations2)" ] }, { diff --git a/notebooks/tutorials/Lazy.ipynb b/notebooks/tutorials/Lazy.ipynb index d6961ef5..5ee2fa9a 100644 --- a/notebooks/tutorials/Lazy.ipynb +++ b/notebooks/tutorials/Lazy.ipynb @@ -90,7 +90,7 @@ "id": "daee5898-3b3d-4876-ad21-14ae16f2fe7e", "metadata": {}, "source": [ - "`Control.Monad.Bayes.Sampler.Lazy` is simply a port of LazyPPL, with some refactoring. In particular, what LazyPPL calls the `Prob` monad, we now call `Sampler`, and `Meas` is built modularly as `Weighted Sampler`. \n", + "`Control.Monad.Bayes.Sampler.Lazy` is simply a port of LazyPPL, with some refactoring. In particular, what LazyPPL calls the `Prob` monad, we now call `SamplerT`, and `Meas` is built modularly as `WeightedT SamplerT`. \n", "\n", "LazyPPL also comes with a number of inference algorithms that make special use of the sampler's laziness, including weighted importance sampling and Metropolis-Hastings. Monad-Bayes has these also, as shown below:" ] diff --git a/notebooks/tutorials/MCMC.ipynb b/notebooks/tutorials/MCMC.ipynb index 11202e85..5722d674 100644 --- a/notebooks/tutorials/MCMC.ipynb +++ b/notebooks/tutorials/MCMC.ipynb @@ -17,7 +17,7 @@ "import qualified Data.Text as T\n", "import Control.Arrow (first,second)\n", "import Control.Monad\n", - "import Graphics.Vega.VegaLite hiding (density)\n", + "import Graphics.Vega.VegaLite hiding (runDensityT)\n", "import qualified Graphics.Vega.VegaLite as VL\n", "import IHaskell.Display.Hvega (vlShow)\n", "\n", diff --git a/notebooks/tutorials/SMC.ipynb b/notebooks/tutorials/SMC.ipynb index 5afd2651..be4973b9 100644 --- a/notebooks/tutorials/SMC.ipynb +++ b/notebooks/tutorials/SMC.ipynb @@ -1885,7 +1885,7 @@ } ], "source": [ - "particles <- sampleIOfixed $ runPopulation $ \n", + "particles <- sampleIOfixed $ runPopulationT $ \n", " smc SMCConfig {\n", " numSteps = 100, \n", " numParticles = 1000, \n", diff --git a/notebooks/tutorials/Sampling.ipynb b/notebooks/tutorials/Sampling.ipynb index c5432bcd..5ed42a37 100644 --- a/notebooks/tutorials/Sampling.ipynb +++ b/notebooks/tutorials/Sampling.ipynb @@ -2242,7 +2242,7 @@ "source": [ "sampler $ \n", " plot . histogramToList . histogram 200\n", - " <$> replicateM 100000 (weighted model3)\n" + " <$> replicateM 100000 (runWeightedT model3)\n" ] }, { @@ -2278,7 +2278,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Here, we use `weighted` to convert `model4` into a distribution over pairs of numbers and their weights. We then sample from that. To view " + "Here, we use `runWeightedT` to convert `model4` into a distribution over pairs of numbers and their weights. We then sample from that. To view " ] }, { @@ -3117,7 +3117,7 @@ "source": [ "sampler $ \n", " plot . histogramToList . histogram 200\n", - " <$> replicateM 100000 (weighted model4)" + " <$> replicateM 100000 (runWeightedT model4)" ] }, { @@ -3961,7 +3961,7 @@ } ], "source": [ - "sampler . unweighted $ plot . histogramToList . histogram 200 <$> replicateM 100000 (weighted $ unweighted model4)\n" + "sampler . unweighted $ plot . histogramToList . histogram 200 <$> replicateM 100000 (runWeightedT $ unweighted model4)\n" ] }, { @@ -3987,7 +3987,7 @@ } ], "source": [ - "draws4 <- sampler $ replicateM 10 $ weighted model4\n", + "draws4 <- sampler $ replicateM 10 $ runWeightedT model4\n", "draws4" ] }, @@ -4833,7 +4833,7 @@ "\n", "sampler $ \n", " plot . histogramToList . histogram 200\n", - " <$> replicateM 100000 (weighted model5)" + " <$> replicateM 100000 (runWeightedT model5)" ] } ], diff --git a/src/Control/Monad/Bayes/Density/Free.hs b/src/Control/Monad/Bayes/Density/Free.hs index e67b2a78..d554056f 100644 --- a/src/Control/Monad/Bayes/Density/Free.hs +++ b/src/Control/Monad/Bayes/Density/Free.hs @@ -12,13 +12,13 @@ -- Stability : experimental -- Portability : GHC -- --- 'Density' is a free monad transformer over random sampling. +-- 'DensityT' is a free monad transformer over random sampling. module Control.Monad.Bayes.Density.Free - ( Density, + ( DensityT (..), hoist, interpret, withRandomness, - density, + runDensityT, traced, ) where @@ -36,40 +36,40 @@ 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} +newtype DensityT m a = DensityT {getDensityT :: FT SamF m a} deriving newtype (Functor, Applicative, Monad, MonadTrans) -instance MonadFree SamF (Density m) where - wrap = Density . wrap . fmap runDensity +instance MonadFree SamF (DensityT m) where + wrap = DensityT . wrap . fmap getDensityT -instance (Monad m) => MonadDistribution (Density m) where - random = Density $ liftF (Random id) +instance (Monad m) => MonadDistribution (DensityT m) where + random = DensityT $ liftF (Random id) --- | Hoist 'Density' through a monad transform. -hoist :: (Monad m, Monad n) => (forall x. m x -> n x) -> Density m a -> Density n a -hoist f (Density m) = Density (hoistFT f m) +-- | Hoist 'DensityT' through a monad transform. +hoist :: (Monad m, Monad n) => (forall x. m x -> n x) -> DensityT m a -> DensityT n a +hoist f (DensityT m) = DensityT (hoistFT f m) -- | Execute random sampling in the transformed monad. -interpret :: (MonadDistribution m) => Density m a -> m a -interpret (Density m) = iterT f m +interpret :: (MonadDistribution m) => DensityT m a -> m a +interpret (DensityT m) = iterT f m where f (Random k) = random >>= k -- | Execute computation with supplied values for random choices. -withRandomness :: (Monad m) => [Double] -> Density m a -> m a -withRandomness randomness (Density m) = evalStateT (iterTM f m) randomness +withRandomness :: (Monad m) => [Double] -> DensityT m a -> m a +withRandomness randomness (DensityT m) = evalStateT (iterTM f m) randomness where f (Random k) = do xs <- get case xs of - [] -> error "Density: the list of randomness was too short" + [] -> error "DensityT: the list of randomness was too short" y : ys -> put ys >> k y -- | Execute computation with supplied values for a subset of random choices. -- Return the output value and a record of all random choices used, whether -- taken as input or drawn using the transformed monad. -density :: (MonadDistribution m) => [Double] -> Density m a -> m (a, [Double]) -density randomness (Density m) = +runDensityT :: (MonadDistribution m) => [Double] -> DensityT m a -> m (a, [Double]) +runDensityT randomness (DensityT m) = runWriterT $ evalStateT (iterTM f $ hoistFT lift m) randomness where f (Random k) = do @@ -84,5 +84,5 @@ density randomness (Density m) = k x -- | Like 'density', but use an arbitrary sampling monad. -traced :: (MonadDistribution m) => [Double] -> Density Identity a -> m (a, [Double]) -traced randomness m = density randomness $ hoist (return . runIdentity) m +traced :: (MonadDistribution m) => [Double] -> DensityT Identity a -> m (a, [Double]) +traced randomness m = runDensityT randomness $ hoist (return . runIdentity) m diff --git a/src/Control/Monad/Bayes/Density/State.hs b/src/Control/Monad/Bayes/Density/State.hs index 8bc514ab..800d7e38 100644 --- a/src/Control/Monad/Bayes/Density/State.hs +++ b/src/Control/Monad/Bayes/Density/State.hs @@ -13,21 +13,21 @@ import Control.Monad.Bayes.Class (MonadDistribution (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) +newtype DensityT m a = DensityT {getDensityT :: WriterT [Double] (StateT [Double] m) a} deriving newtype (Functor, Applicative, Monad) -instance MonadTrans Density where - lift = Density . lift . lift +instance MonadTrans DensityT where + lift = DensityT . lift . lift -instance (Monad m) => MonadState [Double] (Density m) where - get = Density $ lift $ get - put = Density . lift . put +instance (Monad m) => MonadState [Double] (DensityT m) where + get = DensityT $ lift $ get + put = DensityT . lift . put -instance (Monad m) => MonadWriter [Double] (Density m) where - tell = Density . tell - listen = Density . listen . runDensity - pass = Density . pass . runDensity +instance (Monad m) => MonadWriter [Double] (DensityT m) where + tell = DensityT . tell + listen = DensityT . listen . getDensityT + pass = DensityT . pass . getDensityT -instance (MonadDistribution m) => MonadDistribution (Density m) where +instance (MonadDistribution m) => MonadDistribution (DensityT m) where random = do trace <- get x <- case trace of @@ -36,5 +36,5 @@ instance (MonadDistribution m) => MonadDistribution (Density m) where tell [x] pure x -density :: (Monad m) => Density m b -> [Double] -> m (b, [Double]) -density (Density m) = evalStateT (runWriterT m) +runDensityT :: (Monad m) => DensityT m b -> [Double] -> m (b, [Double]) +runDensityT (DensityT m) = evalStateT (runWriterT m) diff --git a/src/Control/Monad/Bayes/Inference/Lazy/MH.hs b/src/Control/Monad/Bayes/Inference/Lazy/MH.hs index fe9d2c52..f25bb7d4 100644 --- a/src/Control/Monad/Bayes/Inference/Lazy/MH.hs +++ b/src/Control/Monad/Bayes/Inference/Lazy/MH.hs @@ -7,18 +7,18 @@ module Control.Monad.Bayes.Inference.Lazy.MH where import Control.Monad.Bayes.Class (Log (ln)) import Control.Monad.Bayes.Sampler.Lazy - ( Sampler (runSampler), + ( SamplerT (runSamplerT), Tree (..), Trees (..), randomTree, ) -import Control.Monad.Bayes.Weighted (Weighted, weighted) +import Control.Monad.Bayes.Weighted (WeightedT, runWeightedT) import Control.Monad.Extra (iterateM) import Control.Monad.State.Lazy (MonadState (get, put), runState) import System.Random (RandomGen (split), getStdGen, newStdGen) import System.Random qualified as R -mh :: forall a. Double -> Weighted Sampler a -> IO [(a, Log Double)] +mh :: forall a. Double -> WeightedT SamplerT a -> IO [(a, Log Double)] mh p m = do -- Top level: produce a stream of samples. -- Split the random number generator in two @@ -27,7 +27,7 @@ mh p m = do g <- newStdGen >> getStdGen let (g1, g2) = split g let t = randomTree g1 - let (x, w) = runSampler (weighted m) t + let (x, w) = runSamplerT (runWeightedT m) t -- Now run step over and over to get a stream of (tree,result,weight)s. let (samples, _) = runState (iterateM step (t, x, w)) g2 -- The stream of seeds is used to produce a stream of result/weight pairs. @@ -50,7 +50,7 @@ mh p m = do let t' = mutateTree p g1 t -- Rerun the model with the new tree, to get a new -- weight w'. - let (x', w') = runSampler (weighted m) t' + let (x', w') = runSamplerT (runWeightedT m) t' -- MH acceptance ratio. This is the probability of either -- returning the new seed or the old one. let ratio = w' / w diff --git a/src/Control/Monad/Bayes/Inference/Lazy/WIS.hs b/src/Control/Monad/Bayes/Inference/Lazy/WIS.hs index 76a10651..ba76a41b 100644 --- a/src/Control/Monad/Bayes/Inference/Lazy/WIS.hs +++ b/src/Control/Monad/Bayes/Inference/Lazy/WIS.hs @@ -1,15 +1,15 @@ module Control.Monad.Bayes.Inference.Lazy.WIS where -import Control.Monad.Bayes.Sampler.Lazy (Sampler, weightedsamples) -import Control.Monad.Bayes.Weighted (Weighted) +import Control.Monad.Bayes.Sampler.Lazy (SamplerT, weightedsamples) +import Control.Monad.Bayes.Weighted (WeightedT) import Numeric.Log (Log (Exp)) import System.Random (Random (randoms), getStdGen, newStdGen) --- | Weighted Importance Sampling +-- | WeightedT Importance Sampling -- | Likelihood weighted importance sampling first draws n weighted samples, -- and then samples a stream of results from that regarded as an empirical distribution -lwis :: Int -> Weighted Sampler a -> IO [a] +lwis :: Int -> WeightedT SamplerT a -> IO [a] lwis n m = do xws <- weightedsamples m let xws' = take n $ accumulate xws 0 diff --git a/src/Control/Monad/Bayes/Inference/MCMC.hs b/src/Control/Monad/Bayes/Inference/MCMC.hs index 36f1f674..56c2a9b4 100644 --- a/src/Control/Monad/Bayes/Inference/MCMC.hs +++ b/src/Control/Monad/Bayes/Inference/MCMC.hs @@ -21,7 +21,7 @@ import Control.Monad.Bayes.Traced.Common ) import Control.Monad.Bayes.Traced.Dynamic qualified as Dynamic import Control.Monad.Bayes.Traced.Static qualified as Static -import Control.Monad.Bayes.Weighted (Weighted, unweighted) +import Control.Monad.Bayes.Weighted (WeightedT, unweighted) import Pipes ((>->)) import Pipes qualified as P import Pipes.Prelude qualified as P @@ -33,25 +33,25 @@ data MCMCConfig = MCMCConfig {proposal :: Proposal, numMCMCSteps :: Int, numBurn defaultMCMCConfig :: MCMCConfig defaultMCMCConfig = MCMCConfig {proposal = SingleSiteMH, numMCMCSteps = 1, numBurnIn = 0} -mcmc :: (MonadDistribution m) => MCMCConfig -> Static.Traced (Weighted m) a -> m [a] +mcmc :: (MonadDistribution m) => MCMCConfig -> Static.TracedT (WeightedT m) a -> m [a] mcmc (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Static.mh numMCMCSteps m -mcmcBasic :: (MonadDistribution m) => MCMCConfig -> Basic.Traced (Weighted m) a -> m [a] +mcmcBasic :: (MonadDistribution m) => MCMCConfig -> Basic.TracedT (WeightedT m) a -> m [a] mcmcBasic (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Basic.mh numMCMCSteps m -mcmcDynamic :: (MonadDistribution m) => MCMCConfig -> Dynamic.Traced (Weighted m) a -> m [a] +mcmcDynamic :: (MonadDistribution m) => MCMCConfig -> Dynamic.TracedT (WeightedT m) a -> m [a] mcmcDynamic (MCMCConfig {..}) m = burnIn numBurnIn $ unweighted $ Dynamic.mh numMCMCSteps m -- -- | draw iid samples until you get one that has non-zero likelihood -independentSamples :: (Monad m) => Static.Traced m a -> P.Producer (MHResult a) m (Trace a) -independentSamples (Static.Traced _w d) = +independentSamples :: (Monad m) => Static.TracedT m a -> P.Producer (MHResult a) m (Trace a) +independentSamples (Static.TracedT _w d) = P.repeatM d >-> P.takeWhile' ((== 0) . probDensity) >-> P.map (MHResult False) -- | convert a probabilistic program into a producer of samples -mcmcP :: (MonadDistribution m) => MCMCConfig -> Static.Traced m a -> P.Producer (MHResult a) m () -mcmcP MCMCConfig {..} m@(Static.Traced w _) = do +mcmcP :: (MonadDistribution m) => MCMCConfig -> Static.TracedT m a -> P.Producer (MHResult a) m () +mcmcP MCMCConfig {..} m@(Static.TracedT w _) = do initialValue <- independentSamples m >-> P.drain ( P.unfoldr (fmap (Right . (\k -> (k, trace k))) . mhTransWithBool w) initialValue >-> P.drop numBurnIn diff --git a/src/Control/Monad/Bayes/Inference/PMMH.hs b/src/Control/Monad/Bayes/Inference/PMMH.hs index f63f9848..a58b1946 100644 --- a/src/Control/Monad/Bayes/Inference/PMMH.hs +++ b/src/Control/Monad/Bayes/Inference/PMMH.hs @@ -22,13 +22,13 @@ import Control.Monad.Bayes.Class (Bayesian (generative), MonadDistribution, Mona import Control.Monad.Bayes.Inference.MCMC (MCMCConfig, mcmc) import Control.Monad.Bayes.Inference.SMC (SMCConfig (), smc) import Control.Monad.Bayes.Population as Pop - ( Population, + ( PopulationT, hoist, - population, pushEvidence, + runPopulationT, ) -import Control.Monad.Bayes.Sequential.Coroutine (Sequential) -import Control.Monad.Bayes.Traced.Static (Traced) +import Control.Monad.Bayes.Sequential.Coroutine (SequentialT) +import Control.Monad.Bayes.Traced.Static (TracedT) import Control.Monad.Bayes.Weighted import Control.Monad.Trans (lift) import Numeric.Log (Log) @@ -37,15 +37,15 @@ import Numeric.Log (Log) pmmh :: (MonadDistribution m) => MCMCConfig -> - SMCConfig (Weighted m) -> - Traced (Weighted m) a1 -> - (a1 -> Sequential (Population (Weighted m)) a2) -> + SMCConfig (WeightedT m) -> + TracedT (WeightedT m) a1 -> + (a1 -> SequentialT (PopulationT (WeightedT m)) a2) -> m [[(a2, Log Double)]] pmmh mcmcConf smcConf param model = mcmc mcmcConf ( param - >>= population + >>= runPopulationT . pushEvidence . Pop.hoist lift . smc smcConf @@ -56,7 +56,7 @@ pmmh mcmcConf smcConf param model = pmmhBayesianModel :: (MonadMeasure m) => MCMCConfig -> - SMCConfig (Weighted m) -> + SMCConfig (WeightedT m) -> (forall m'. (MonadMeasure m') => Bayesian m' a1 a2) -> m [[(a2, Log Double)]] pmmhBayesianModel mcmcConf smcConf bm = pmmh mcmcConf smcConf (prior bm) (generative bm) diff --git a/src/Control/Monad/Bayes/Inference/RMSMC.hs b/src/Control/Monad/Bayes/Inference/RMSMC.hs index 7e5d77d7..d9bc8a9b 100644 --- a/src/Control/Monad/Bayes/Inference/RMSMC.hs +++ b/src/Control/Monad/Bayes/Inference/RMSMC.hs @@ -24,7 +24,7 @@ import Control.Monad.Bayes.Class (MonadDistribution) import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (..)) import Control.Monad.Bayes.Inference.SMC import Control.Monad.Bayes.Population - ( Population, + ( PopulationT, spawn, withParticles, ) @@ -33,7 +33,7 @@ import Control.Monad.Bayes.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, + ( TracedT, marginal, mhStep, ) @@ -46,8 +46,8 @@ rmsmc :: MCMCConfig -> SMCConfig m -> -- | model - Sequential (Traced (Population m)) a -> - Population m a + SequentialT (TracedT (PopulationT m)) a -> + PopulationT m a rmsmc (MCMCConfig {..}) (SMCConfig {..}) = marginal . S.sequentially (composeCopies numMCMCSteps mhStep . TrStat.hoist resampler) numSteps @@ -60,8 +60,8 @@ rmsmcBasic :: MCMCConfig -> SMCConfig m -> -- | model - Sequential (TrBas.Traced (Population m)) a -> - Population m a + SequentialT (TrBas.TracedT (PopulationT m)) a -> + PopulationT m a rmsmcBasic (MCMCConfig {..}) (SMCConfig {..}) = TrBas.marginal . S.sequentially (composeCopies numMCMCSteps TrBas.mhStep . TrBas.hoist resampler) numSteps @@ -75,8 +75,8 @@ rmsmcDynamic :: MCMCConfig -> SMCConfig m -> -- | model - Sequential (TrDyn.Traced (Population m)) a -> - Population m a + SequentialT (TrDyn.TracedT (PopulationT m)) a -> + PopulationT m a rmsmcDynamic (MCMCConfig {..}) (SMCConfig {..}) = TrDyn.marginal . S.sequentially (TrDyn.freeze . composeCopies numMCMCSteps TrDyn.mhStep . TrDyn.hoist resampler) numSteps diff --git a/src/Control/Monad/Bayes/Inference/SMC.hs b/src/Control/Monad/Bayes/Inference/SMC.hs index 22d77d71..3f3a30b2 100644 --- a/src/Control/Monad/Bayes/Inference/SMC.hs +++ b/src/Control/Monad/Bayes/Inference/SMC.hs @@ -22,14 +22,14 @@ where import Control.Monad.Bayes.Class (MonadDistribution, MonadMeasure) import Control.Monad.Bayes.Population - ( Population, + ( PopulationT, pushEvidence, withParticles, ) import Control.Monad.Bayes.Sequential.Coroutine as Coroutine data SMCConfig m = SMCConfig - { resampler :: forall x. Population m x -> Population m x, + { resampler :: forall x. PopulationT m x -> PopulationT m x, numSteps :: Int, numParticles :: Int } @@ -39,8 +39,8 @@ data SMCConfig m = SMCConfig smc :: (MonadDistribution m) => SMCConfig m -> - Coroutine.Sequential (Population m) a -> - Population m a + Coroutine.SequentialT (PopulationT m) a -> + PopulationT m a smc SMCConfig {..} = Coroutine.sequentially resampler numSteps . Coroutine.hoistFirst (withParticles numParticles) @@ -49,5 +49,5 @@ smc SMCConfig {..} = -- Weights are normalized at each timestep and the total weight is pushed -- as a score into the transformed monad. smcPush :: - (MonadMeasure m) => SMCConfig m -> Coroutine.Sequential (Population m) a -> Population m a + (MonadMeasure m) => SMCConfig m -> Coroutine.SequentialT (PopulationT m) a -> PopulationT m a smcPush config = smc config {resampler = (pushEvidence . resampler config)} diff --git a/src/Control/Monad/Bayes/Inference/SMC2.hs b/src/Control/Monad/Bayes/Inference/SMC2.hs index f4b877e9..5570a2ba 100644 --- a/src/Control/Monad/Bayes/Inference/SMC2.hs +++ b/src/Control/Monad/Bayes/Inference/SMC2.hs @@ -27,17 +27,17 @@ import Control.Monad.Bayes.Class import Control.Monad.Bayes.Inference.MCMC import Control.Monad.Bayes.Inference.RMSMC (rmsmc) import Control.Monad.Bayes.Inference.SMC (SMCConfig (SMCConfig, numParticles, numSteps, resampler), smcPush) -import Control.Monad.Bayes.Population as Pop (Population, population, resampleMultinomial) -import Control.Monad.Bayes.Sequential.Coroutine (Sequential) +import Control.Monad.Bayes.Population as Pop (PopulationT, resampleMultinomial, runPopulationT) +import Control.Monad.Bayes.Sequential.Coroutine (SequentialT) import Control.Monad.Bayes.Traced import Control.Monad.Trans (MonadTrans (..)) import Numeric.Log (Log) -- | Helper monad transformer for preprocessing the model for 'smc2'. -newtype SMC2 m a = SMC2 (Sequential (Traced (Population m)) a) +newtype SMC2 m a = SMC2 (SequentialT (TracedT (PopulationT m)) a) deriving newtype (Functor, Applicative, Monad) -setup :: SMC2 m a -> Sequential (Traced (Population m)) a +setup :: SMC2 m a -> SequentialT (TracedT (PopulationT m)) a setup (SMC2 m) = m instance MonadTrans SMC2 where @@ -63,12 +63,12 @@ smc2 :: -- | number of MH transitions Int -> -- | model parameters - Sequential (Traced (Population m)) b -> + SequentialT (TracedT (PopulationT m)) b -> -- | model - (b -> Sequential (Population (SMC2 m)) a) -> - Population m [(a, Log Double)] + (b -> SequentialT (PopulationT (SMC2 m)) a) -> + PopulationT m [(a, Log Double)] smc2 k n p t param m = rmsmc MCMCConfig {numMCMCSteps = t, proposal = SingleSiteMH, numBurnIn = 0} SMCConfig {numParticles = p, numSteps = k, resampler = resampleMultinomial} - (param >>= setup . population . smcPush (SMCConfig {numSteps = k, numParticles = n, resampler = resampleMultinomial}) . m) + (param >>= setup . runPopulationT . smcPush (SMCConfig {numSteps = k, numParticles = n, resampler = resampleMultinomial}) . m) diff --git a/src/Control/Monad/Bayes/Inference/TUI.hs b/src/Control/Monad/Bayes/Inference/TUI.hs index 4d155337..56addd34 100644 --- a/src/Control/Monad/Bayes/Inference/TUI.hs +++ b/src/Control/Monad/Bayes/Inference/TUI.hs @@ -18,7 +18,7 @@ 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 (TracedT) import Control.Monad.Bayes.Traced.Common hiding (burnIn) import Control.Monad.Bayes.Weighted import Data.Scientific (FPFormat (Exponent), formatScientific, fromFloatDigits) @@ -130,7 +130,7 @@ theMap = (attrName "highlight", fg yellow) ] -tui :: (Show a) => Int -> Traced (Weighted SamplerIO) a -> ([a] -> Widget ()) -> IO () +tui :: (Show a) => Int -> TracedT (WeightedT SamplerIO) a -> ([a] -> Widget ()) -> IO () tui burnIn distribution visualizer = void do eventChan <- B.newBChan 10 initialVty <- buildVty diff --git a/src/Control/Monad/Bayes/Integrator.hs b/src/Control/Monad/Bayes/Integrator.hs index 5376f7b5..a38e4d8c 100644 --- a/src/Control/Monad/Bayes/Integrator.hs +++ b/src/Control/Monad/Bayes/Integrator.hs @@ -35,7 +35,7 @@ import Control.Applicative (Applicative (..)) import Control.Foldl (Fold) import Control.Foldl qualified as Foldl import Control.Monad.Bayes.Class (MonadDistribution (bernoulli, random, uniformD)) -import Control.Monad.Bayes.Weighted (Weighted, weighted) +import Control.Monad.Bayes.Weighted (WeightedT, runWeightedT) import Control.Monad.Cont ( Cont, ContT (ContT), @@ -49,12 +49,14 @@ 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} +newtype Integrator a = Integrator {getIntegrator :: Cont Double a} deriving newtype (Functor, Applicative, Monad) -integrator, runIntegrator :: (a -> Double) -> Integrator a -> Double -integrator f (Integrator a) = runCont a f -runIntegrator = integrator +runIntegrator :: (a -> Double) -> Integrator a -> Double +runIntegrator f (Integrator a) = runCont a f + +integrator :: ((a -> Double) -> Double) -> Integrator a +integrator = Integrator . cont instance MonadDistribution Integrator where random = fromDensityFunction $ Statistics.density $ Statistics.uniformDistr 0 1 @@ -85,28 +87,28 @@ empirical = Integrator . cont . flip weightedAverage averageFold = (/) <$> Foldl.sum <*> Foldl.genericLength expectation :: Integrator Double -> Double -expectation = integrator id +expectation = runIntegrator id variance :: Integrator Double -> Double -variance nu = integrator (^ 2) nu - expectation nu ^ 2 +variance nu = runIntegrator (^ 2) nu - expectation nu ^ 2 momentGeneratingFunction :: Integrator Double -> Double -> Double -momentGeneratingFunction nu t = integrator (\x -> exp (t * x)) nu +momentGeneratingFunction nu t = runIntegrator (\x -> exp (t * x)) nu cumulantGeneratingFunction :: Integrator Double -> Double -> Double cumulantGeneratingFunction nu = log . momentGeneratingFunction nu -normalize :: Weighted Integrator a -> Integrator a +normalize :: WeightedT Integrator a -> Integrator a normalize m = - let m' = weighted m - z = integrator (ln . exp . snd) m' + let m' = runWeightedT m + z = runIntegrator (ln . exp . snd) m' in do - (x, d) <- weighted m + (x, d) <- runWeightedT m Integrator $ cont $ \f -> (f () * (ln $ exp d)) / z return x cdf :: Integrator Double -> Double -> Double -cdf nu x = integrator (negativeInfinity `to` x) nu +cdf nu x = runIntegrator (negativeInfinity `to` x) nu where negativeInfinity :: Double negativeInfinity = negate (1 / 0) @@ -117,7 +119,7 @@ cdf nu x = integrator (negativeInfinity `to` x) nu | otherwise = 0 volume :: Integrator Double -> Double -volume = integrator (const 1) +volume = runIntegrator (const 1) containing :: (Num a, Eq b) => [b] -> b -> a containing xs x @@ -133,12 +135,12 @@ instance (Num a) => Num (Integrator a) where fromInteger = pure . fromInteger probability :: (Ord a) => (a, a) -> Integrator a -> Double -probability (lower, upper) = integrator (\x -> if x < upper && x >= lower then 1 else 0) +probability (lower, upper) = runIntegrator (\x -> if x < upper && x >= lower then 1 else 0) enumeratorWith :: (Ord a) => Set a -> Integrator a -> [(a, Double)] enumeratorWith ls meas = [ ( val, - integrator + runIntegrator (\x -> if x == val then 1 else 0) meas ) @@ -149,7 +151,7 @@ histogram :: (Enum a, Ord a, Fractional a) => Int -> a -> - Weighted Integrator a -> + WeightedT Integrator a -> [(a, Double)] histogram nBins binSize model = do x <- take nBins [1 ..] diff --git a/src/Control/Monad/Bayes/Population.hs b/src/Control/Monad/Bayes/Population.hs index 4a351cde..2a384f77 100644 --- a/src/Control/Monad/Bayes/Population.hs +++ b/src/Control/Monad/Bayes/Population.hs @@ -13,11 +13,10 @@ -- Stability : experimental -- Portability : GHC -- --- 'Population' turns a single sample into a collection of weighted samples. +-- 'PopulationT' turns a single sample into a collection of weighted samples. module Control.Monad.Bayes.Population - ( Population, - population, - runPopulation, + ( PopulationT (..), + runPopulationT, explicitPopulation, fromWeightedList, spawn, @@ -47,11 +46,11 @@ import Control.Monad.Bayes.Class factor, ) import Control.Monad.Bayes.Weighted - ( Weighted, + ( WeightedT, applyWeight, extractWeight, - weighted, - withWeight, + runWeightedT, + weightedT, ) import Control.Monad.List (ListT (..), MonadIO, MonadTrans (..)) import Data.List (unfoldr) @@ -64,46 +63,43 @@ import 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) +newtype PopulationT m a = PopulationT {getPopulationT :: WeightedT (ListT m) a} deriving newtype (Functor, Applicative, Monad, MonadIO, MonadDistribution, MonadFactor, MonadMeasure) -instance MonadTrans Population where - lift = Population . lift . lift +instance MonadTrans PopulationT where + lift = PopulationT . lift . lift -- | Explicit representation of the weighted sample with weights in the log -- domain. -population, runPopulation :: Population m a -> m [(a, Log Double)] -population (Population m) = runListT $ weighted m - --- | deprecated synonym -runPopulation = population +runPopulationT :: PopulationT m a -> m [(a, Log Double)] +runPopulationT = runListT . runWeightedT . getPopulationT -- | Explicit representation of the weighted sample. -explicitPopulation :: (Functor m) => Population m a -> m [(a, Double)] -explicitPopulation = fmap (map (second (exp . ln))) . population +explicitPopulation :: (Functor m) => PopulationT m a -> m [(a, Double)] +explicitPopulation = fmap (map (second (exp . ln))) . runPopulationT --- | Initialize 'Population' with a concrete weighted sample. -fromWeightedList :: (Monad m) => m [(a, Log Double)] -> Population m a -fromWeightedList = Population . withWeight . ListT +-- | Initialize 'PopulationT' with a concrete weighted sample. +fromWeightedList :: (Monad m) => m [(a, Log Double)] -> PopulationT m a +fromWeightedList = PopulationT . weightedT . ListT -- | Increase the sample size by a given factor. -- The weights are adjusted such that their sum is preserved. -- It is therefore safe to use 'spawn' in arbitrary places in the program -- without introducing bias. -spawn :: (Monad m) => Int -> Population m () +spawn :: (Monad m) => Int -> PopulationT m () spawn n = fromWeightedList $ pure $ replicate n ((), 1 / fromIntegral n) -withParticles :: (Monad m) => Int -> Population m a -> Population m a +withParticles :: (Monad m) => Int -> PopulationT m a -> PopulationT m a withParticles n = (spawn n >>) resampleGeneric :: (MonadDistribution m) => -- | resampler (V.Vector Double -> m [Int]) -> - Population m a -> - Population m a + PopulationT m a -> + PopulationT m a resampleGeneric resampler m = fromWeightedList $ do - pop <- population m + pop <- runPopulationT m let (xs, ps) = unzip pop let n = length xs let z = Log.sum ps @@ -150,8 +146,8 @@ systematic u ps = f 0 (u / fromIntegral n) 0 0 [] -- The total weight is preserved. resampleSystematic :: (MonadDistribution m) => - Population m a -> - Population m a + PopulationT m a -> + PopulationT m a resampleSystematic = resampleGeneric (\ps -> (`systematic` ps) <$> random) -- | Stratified sampler. @@ -192,8 +188,8 @@ stratified weights = do -- The total weight is preserved. resampleStratified :: (MonadDistribution m) => - Population m a -> - Population m a + PopulationT m a -> + PopulationT m a resampleStratified = resampleGeneric stratified -- | Multinomial sampler. Sample from \(0, \ldots, n - 1\) \(n\) @@ -206,18 +202,18 @@ multinomial ps = replicateM (V.length ps) (categorical ps) -- The total weight is preserved. resampleMultinomial :: (MonadDistribution m) => - Population m a -> - Population m a + PopulationT m a -> + PopulationT m a resampleMultinomial = resampleGeneric multinomial --- | Separate the sum of weights into the 'Weighted' transformer. +-- | Separate the sum of weights into the 'WeightedT' transformer. -- Weights are normalized after this operation. extractEvidence :: (Monad m) => - Population m a -> - Population (Weighted m) a + PopulationT m a -> + PopulationT (WeightedT m) a extractEvidence m = fromWeightedList $ do - pop <- lift $ population m + pop <- lift $ runPopulationT m let (xs, ps) = unzip pop let z = sum ps let ws = map (if z > 0 then (/ z) else const (1 / fromIntegral (length ps))) ps @@ -228,26 +224,26 @@ extractEvidence m = fromWeightedList $ do -- Weights are normalized after this operation. pushEvidence :: (MonadFactor m) => - Population m a -> - Population m a + PopulationT m a -> + PopulationT m a pushEvidence = hoist applyWeight . extractEvidence -- | A properly weighted single sample, that is one picked at random according -- to the weights, with the sum of all weights. proper :: (MonadDistribution m) => - Population m a -> - Weighted m a + PopulationT m a -> + WeightedT m a proper m = do - pop <- population $ extractEvidence m + pop <- runPopulationT $ extractEvidence m let (xs, ps) = unzip pop index <- logCategorical $ V.fromList ps let x = xs !! index return x -- | Model evidence estimator, also known as pseudo-marginal likelihood. -evidence :: (Monad m) => Population m a -> m (Log Double) -evidence = extractWeight . population . extractEvidence +evidence :: (Monad m) => PopulationT m a -> m (Log Double) +evidence = extractWeight . runPopulationT . extractEvidence -- | Picks one point from the population and uses model evidence as a 'score' -- in the transformed monad. @@ -255,12 +251,12 @@ evidence = extractWeight . population . extractEvidence -- introducing bias. collapse :: (MonadMeasure m) => - Population m a -> + PopulationT m a -> m a collapse = applyWeight . proper --- | Population average of a function, computed using unnormalized weights. -popAvg :: (Monad m) => (a -> Double) -> Population m a -> m Double +-- | PopulationT average of a function, computed using unnormalized weights. +popAvg :: (Monad m) => (a -> Double) -> PopulationT m a -> m Double popAvg f p = do xs <- explicitPopulation p let ys = map (\(x, w) -> f x * w) xs @@ -271,6 +267,6 @@ popAvg f p = do hoist :: (Monad n) => (forall x. m x -> n x) -> - Population m a -> - Population n a -hoist f = fromWeightedList . f . population + PopulationT m a -> + PopulationT n a +hoist f = fromWeightedList . f . runPopulationT diff --git a/src/Control/Monad/Bayes/Sampler/Lazy.hs b/src/Control/Monad/Bayes/Sampler/Lazy.hs index aa878ab4..c12e874b 100644 --- a/src/Control/Monad/Bayes/Sampler/Lazy.hs +++ b/src/Control/Monad/Bayes/Sampler/Lazy.hs @@ -8,7 +8,7 @@ module Control.Monad.Bayes.Sampler.Lazy where import Control.Monad (ap) import Control.Monad.Bayes.Class (MonadDistribution (random)) -import Control.Monad.Bayes.Weighted (Weighted, weighted) +import Control.Monad.Bayes.Weighted (WeightedT, runWeightedT) import Numeric.Log (Log (..)) import System.Random ( RandomGen (split), @@ -35,7 +35,7 @@ data Trees = 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} +newtype SamplerT a = SamplerT {runSamplerT :: Tree -> a} deriving (Functor) -- | Two key things to do with trees: @@ -51,29 +51,29 @@ 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 +instance Applicative SamplerT where + pure = SamplerT . const (<*>) = ap -- | probabilities for a monad. -- | Sequencing is done by splitting the tree -- | and using different bits for different computations. -instance Monad Sampler where +instance Monad SamplerT where return = pure - (Sampler m) >>= f = Sampler \g -> + (SamplerT m) >>= f = SamplerT \g -> let (g1, g2) = splitTree g - (Sampler m') = f (m g1) + (SamplerT m') = f (m g1) in m' g2 -instance MonadDistribution Sampler where - random = Sampler \(Tree r _) -> r +instance MonadDistribution SamplerT where + random = SamplerT \(Tree r _) -> r -sampler :: Sampler a -> IO a -sampler m = newStdGen *> (runSampler m . randomTree <$> getStdGen) +sampler :: SamplerT a -> IO a +sampler m = newStdGen *> (runSamplerT m . randomTree <$> getStdGen) independent :: (Monad m) => m a -> m [a] independent = sequence . repeat -- | 'weightedsamples' runs a probability measure and gets out a stream of (result,weight) pairs -weightedsamples :: Weighted Sampler a -> IO [(a, Log Double)] -weightedsamples = sampler . independent . weighted +weightedsamples :: WeightedT SamplerT a -> IO [(a, Log Double)] +weightedsamples = sampler . independent . runWeightedT diff --git a/src/Control/Monad/Bayes/Sampler/Strict.hs b/src/Control/Monad/Bayes/Sampler/Strict.hs index 000aee54..f5f9b645 100644 --- a/src/Control/Monad/Bayes/Sampler/Strict.hs +++ b/src/Control/Monad/Bayes/Sampler/Strict.hs @@ -15,7 +15,7 @@ -- 'SamplerIO' and 'SamplerST' are instances of 'MonadDistribution'. Apply a 'MonadFactor' -- transformer to obtain a 'MonadMeasure' that can execute probabilistic models. module Control.Monad.Bayes.Sampler.Strict - ( Sampler, + ( SamplerT (..), SamplerIO, SamplerST, sampleIO, @@ -48,27 +48,27 @@ import System.Random.Stateful (IOGenM (..), STGenM, StatefulGen, StdGen, initStd -- | The sampling interpretation of a probabilistic program -- Here m is typically IO or ST -newtype Sampler g m a = Sampler (ReaderT g m a) deriving (Functor, Applicative, Monad, MonadIO) +newtype SamplerT g m a = SamplerT {runSamplerT :: ReaderT g m a} deriving (Functor, Applicative, Monad, MonadIO) --- | convenient type synonym to show specializations of Sampler +-- | convenient type synonym to show specializations of SamplerT -- to particular pairs of monad and RNG -type SamplerIO = Sampler (IOGenM StdGen) IO +type SamplerIO = SamplerT (IOGenM StdGen) IO --- | convenient type synonym to show specializations of Sampler +-- | convenient type synonym to show specializations of SamplerT -- to particular pairs of monad and RNG -type SamplerST s = Sampler (STGenM StdGen s) (ST s) +type SamplerST s = SamplerT (STGenM StdGen s) (ST s) -instance (StatefulGen g m) => MonadDistribution (Sampler g m) where - random = Sampler (ReaderT uniformDouble01M) +instance (StatefulGen g m) => MonadDistribution (SamplerT g m) where + random = SamplerT (ReaderT uniformDouble01M) - uniform a b = Sampler (ReaderT $ uniformRM (a, b)) - normal m s = Sampler (ReaderT (MWC.normal m s)) - gamma shape scale = Sampler (ReaderT $ MWC.gamma shape scale) - beta a b = Sampler (ReaderT $ MWC.beta a b) + uniform a b = SamplerT (ReaderT $ uniformRM (a, b)) + normal m s = SamplerT (ReaderT (MWC.normal m s)) + gamma shape scale = SamplerT (ReaderT $ MWC.gamma shape scale) + beta a b = SamplerT (ReaderT $ MWC.beta a b) - bernoulli p = Sampler (ReaderT $ MWC.bernoulli p) - categorical ps = Sampler (ReaderT $ MWC.categorical ps) - geometric p = Sampler (ReaderT $ MWC.geometric0 p) + bernoulli p = SamplerT (ReaderT $ MWC.bernoulli p) + categorical ps = SamplerT (ReaderT $ MWC.categorical ps) + geometric p = SamplerT (ReaderT $ MWC.geometric0 p) -- | Sample with a random number generator of your choice e.g. the one -- from `System.Random`. @@ -77,8 +77,8 @@ instance (StatefulGen g m) => MonadDistribution (Sampler g m) where -- >>> import System.Random.Stateful hiding (random) -- >>> newIOGenM (mkStdGen 1729) >>= sampleWith random -- 4.690861245089605e-2 -sampleWith :: Sampler g m a -> g -> m a -sampleWith (Sampler m) = runReaderT m +sampleWith :: SamplerT g m a -> g -> m a +sampleWith (SamplerT m) = runReaderT m -- | initialize random seed using system entropy, and sample sampleIO, sampler :: SamplerIO a -> IO a diff --git a/src/Control/Monad/Bayes/Sequential/Coroutine.hs b/src/Control/Monad/Bayes/Sequential/Coroutine.hs index 43614ffe..926c3db1 100644 --- a/src/Control/Monad/Bayes/Sequential/Coroutine.hs +++ b/src/Control/Monad/Bayes/Sequential/Coroutine.hs @@ -11,9 +11,9 @@ -- Stability : experimental -- Portability : GHC -- --- 'Sequential' represents a computation that can be suspended. +-- 'SequentialT' represents a computation that can be suspended. module Control.Monad.Bayes.Sequential.Coroutine - ( Sequential, + ( SequentialT, suspend, finish, advance, @@ -48,55 +48,55 @@ import Data.Either (isRight) -- useful for implementation of Sequential Monte Carlo related methods. -- All the probabilistic effects are lifted from the transformed monad, but -- also `suspend` is inserted after each `factor`. -newtype Sequential m a = Sequential {runSequential :: Coroutine (Await ()) m a} +newtype SequentialT m a = SequentialT {runSequentialT :: Coroutine (Await ()) m a} deriving newtype (Functor, Applicative, Monad, MonadTrans, MonadIO) extract :: Await () a -> a extract (Await f) = f () -instance (MonadDistribution m) => MonadDistribution (Sequential m) where +instance (MonadDistribution m) => MonadDistribution (SequentialT m) where random = lift random bernoulli = lift . bernoulli categorical = lift . categorical -- | Execution is 'suspend'ed after each 'score'. -instance (MonadFactor m) => MonadFactor (Sequential m) where +instance (MonadFactor m) => MonadFactor (SequentialT m) where score w = lift (score w) >> suspend -instance (MonadMeasure m) => MonadMeasure (Sequential m) +instance (MonadMeasure m) => MonadMeasure (SequentialT m) -- | A point where the computation is paused. -suspend :: (Monad m) => Sequential m () -suspend = Sequential await +suspend :: (Monad m) => SequentialT m () +suspend = SequentialT await -- | Remove the remaining suspension points. -finish :: (Monad m) => Sequential m a -> m a -finish = pogoStick extract . runSequential +finish :: (Monad m) => SequentialT m a -> m a +finish = pogoStick extract . runSequentialT -- | Execute to the next suspension point. -- If the computation is finished, do nothing. -- -- > finish = finish . advance -advance :: (Monad m) => Sequential m a -> Sequential m a -advance = Sequential . bounce extract . runSequential +advance :: (Monad m) => SequentialT m a -> SequentialT m a +advance = SequentialT . bounce extract . runSequentialT -- | Return True if no more suspension points remain. -finished :: (Monad m) => Sequential m a -> m Bool -finished = fmap isRight . resume . runSequential +finished :: (Monad m) => SequentialT m a -> m Bool +finished = fmap isRight . resume . runSequentialT -- | Transform the inner monad. -- This operation only applies to computation up to the first suspension. -hoistFirst :: (forall x. m x -> m x) -> Sequential m a -> Sequential m a -hoistFirst f = Sequential . Coroutine . f . resume . runSequential +hoistFirst :: (forall x. m x -> m x) -> SequentialT m a -> SequentialT m a +hoistFirst f = SequentialT . Coroutine . f . resume . runSequentialT -- | Transform the inner monad. -- The transformation is applied recursively through all the suspension points. hoist :: (Monad m, Monad n) => (forall x. m x -> n x) -> - Sequential m a -> - Sequential n a -hoist f = Sequential . mapMonad f . runSequential + SequentialT m a -> + SequentialT n a +hoist f = SequentialT . mapMonad f . runSequentialT -- | Apply a function a given number of times. composeCopies :: Int -> (a -> a) -> (a -> a) @@ -111,7 +111,7 @@ sequentially, (forall x. m x -> m x) -> -- | number of time steps Int -> - Sequential m a -> + SequentialT m a -> m a sequentially f k = finish . composeCopies k (advance . hoistFirst f) diff --git a/src/Control/Monad/Bayes/Traced/Basic.hs b/src/Control/Monad/Bayes/Traced/Basic.hs index 0db46590..d2b727d1 100644 --- a/src/Control/Monad/Bayes/Traced/Basic.hs +++ b/src/Control/Monad/Bayes/Traced/Basic.hs @@ -9,7 +9,7 @@ -- Stability : experimental -- Portability : GHC module Control.Monad.Bayes.Traced.Basic - ( Traced, + ( TracedT, hoist, marginal, mhStep, @@ -23,7 +23,7 @@ import Control.Monad.Bayes.Class MonadFactor (..), MonadMeasure, ) -import Control.Monad.Bayes.Density.Free (Density) +import Control.Monad.Bayes.Density.Free (DensityT) import Control.Monad.Bayes.Traced.Common ( Trace (..), bind, @@ -31,56 +31,56 @@ import Control.Monad.Bayes.Traced.Common scored, singleton, ) -import Control.Monad.Bayes.Weighted (Weighted) +import Control.Monad.Bayes.Weighted (WeightedT) import Data.Functor.Identity (Identity) import Data.List.NonEmpty as NE (NonEmpty ((:|)), toList) -- | Tracing monad that records random choices made in the program. -data Traced m a = Traced +data TracedT m a = TracedT { -- | Run the program with a modified trace. - model :: Weighted (Density Identity) a, + model :: WeightedT (DensityT Identity) a, -- | Record trace and output. traceDist :: m (Trace a) } -instance (Monad m) => Functor (Traced m) where - fmap f (Traced m d) = Traced (fmap f m) (fmap (fmap f) d) +instance (Monad m) => Functor (TracedT m) where + fmap f (TracedT m d) = TracedT (fmap f m) (fmap (fmap f) d) -instance (Monad m) => Applicative (Traced m) where - pure x = Traced (pure x) (pure (pure x)) - (Traced mf df) <*> (Traced mx dx) = Traced (mf <*> mx) (liftA2 (<*>) df dx) +instance (Monad m) => Applicative (TracedT m) where + pure x = TracedT (pure x) (pure (pure x)) + (TracedT mf df) <*> (TracedT mx dx) = TracedT (mf <*> mx) (liftA2 (<*>) df dx) -instance (Monad m) => Monad (Traced m) where - (Traced mx dx) >>= f = Traced my dy +instance (Monad m) => Monad (TracedT m) where + (TracedT mx dx) >>= f = TracedT my dy where my = mx >>= model . f dy = dx `bind` (traceDist . f) -instance (MonadDistribution m) => MonadDistribution (Traced m) where - random = Traced random (fmap singleton random) +instance (MonadDistribution m) => MonadDistribution (TracedT m) where + random = TracedT random (fmap singleton random) -instance (MonadFactor m) => MonadFactor (Traced m) where - score w = Traced (score w) (score w >> pure (scored w)) +instance (MonadFactor m) => MonadFactor (TracedT m) where + score w = TracedT (score w) (score w >> pure (scored w)) -instance (MonadMeasure m) => MonadMeasure (Traced m) +instance (MonadMeasure m) => MonadMeasure (TracedT m) -hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a -hoist f (Traced m d) = Traced m (f d) +hoist :: (forall x. m x -> m x) -> TracedT m a -> TracedT m a +hoist f (TracedT m d) = TracedT m (f d) -- | Discard the trace and supporting infrastructure. -marginal :: (Monad m) => Traced m a -> m a -marginal (Traced _ d) = fmap output d +marginal :: (Monad m) => TracedT m a -> m a +marginal (TracedT _ d) = fmap output d -- | A single step of the Trace Metropolis-Hastings algorithm. -mhStep :: (MonadDistribution m) => Traced m a -> Traced m a -mhStep (Traced m d) = Traced m d' +mhStep :: (MonadDistribution m) => TracedT m a -> TracedT m a +mhStep (TracedT m d) = TracedT m d' where d' = d >>= mhTrans' m -- | Full run of the Trace Metropolis-Hastings algorithm with a specified -- number of steps. -mh :: (MonadDistribution m) => Int -> Traced m a -> m [a] -mh n (Traced m d) = fmap (map output . NE.toList) (f n) +mh :: (MonadDistribution m) => Int -> TracedT m a -> m [a] +mh n (TracedT m d) = fmap (map output . NE.toList) (f n) where f k | k <= 0 = fmap (:| []) d diff --git a/src/Control/Monad/Bayes/Traced/Common.hs b/src/Control/Monad/Bayes/Traced/Common.hs index 0daff0da..a112c5f1 100644 --- a/src/Control/Monad/Bayes/Traced/Common.hs +++ b/src/Control/Monad/Bayes/Traced/Common.hs @@ -26,10 +26,10 @@ import Control.Monad.Bayes.Class ) import Control.Monad.Bayes.Density.Free qualified as Free import Control.Monad.Bayes.Density.State qualified as State -import Control.Monad.Bayes.Weighted as Weighted - ( Weighted, +import Control.Monad.Bayes.Weighted as WeightedT + ( WeightedT, hoist, - weighted, + runWeightedT, ) import Control.Monad.Writer (WriterT (WriterT, runWriterT)) import Data.Functor.Identity (Identity (runIdentity)) @@ -81,7 +81,7 @@ bind dx f = do return $ t2 {variables = variables t1 ++ variables t2, probDensity = probDensity t1 * probDensity t2} -- | A single Metropolis-corrected transition of single-site Trace MCMC. -mhTrans :: (MonadDistribution m) => (Weighted (State.Density m)) a -> Trace a -> m (Trace a) +mhTrans :: (MonadDistribution m) => (WeightedT (State.DensityT m)) a -> Trace a -> m (Trace a) mhTrans m t@Trace {variables = us, probDensity = p} = do let n = length us us' <- do @@ -90,16 +90,16 @@ mhTrans m t@Trace {variables = us, probDensity = p} = do case splitAt i us of (xs, _ : ys) -> return $ xs ++ (u' : ys) _ -> error "impossible" - ((b, q), vs) <- State.density (weighted m) us' + ((b, q), vs) <- State.runDensityT (runWeightedT m) us' let ratio = (exp . ln) $ min 1 (q * fromIntegral n / (p * fromIntegral (length vs))) accept <- bernoulli ratio return $ if accept then Trace vs b q else t -mhTransFree :: (MonadDistribution m) => Weighted (Free.Density m) a -> Trace a -> m (Trace a) +mhTransFree :: (MonadDistribution m) => WeightedT (Free.DensityT m) a -> Trace a -> m (Trace a) mhTransFree m t = trace <$> mhTransWithBool m t -- | A single Metropolis-corrected transition of single-site Trace MCMC. -mhTransWithBool :: (MonadDistribution m) => Weighted (Free.Density m) a -> Trace a -> m (MHResult a) +mhTransWithBool :: (MonadDistribution m) => WeightedT (Free.DensityT m) a -> Trace a -> m (MHResult a) mhTransWithBool m t@Trace {variables = us, probDensity = p} = do let n = length us us' <- do @@ -108,14 +108,14 @@ mhTransWithBool m t@Trace {variables = us, probDensity = p} = do case splitAt i us of (xs, _ : ys) -> return $ xs ++ (u' : ys) _ -> error "impossible" - ((b, q), vs) <- runWriterT $ weighted $ Weighted.hoist (WriterT . Free.density us') m + ((b, q), vs) <- runWriterT $ runWeightedT $ WeightedT.hoist (WriterT . Free.runDensityT us') m let ratio = (exp . ln) $ min 1 (q * fromIntegral n / (p * fromIntegral (length vs))) accept <- bernoulli ratio return if accept then MHResult True (Trace vs b q) else MHResult False t -- | A variant of 'mhTrans' with an external sampling monad. -mhTrans' :: (MonadDistribution m) => Weighted (Free.Density Identity) a -> Trace a -> m (Trace a) -mhTrans' m = mhTransFree (Weighted.hoist (Free.hoist (return . runIdentity)) m) +mhTrans' :: (MonadDistribution m) => WeightedT (Free.DensityT Identity) a -> Trace a -> m (Trace a) +mhTrans' m = mhTransFree (WeightedT.hoist (Free.hoist (return . runIdentity)) m) -- | burn in an MCMC chain for n steps (which amounts to dropping samples of the end of the list) burnIn :: (Functor m) => Int -> m [a] -> m [a] diff --git a/src/Control/Monad/Bayes/Traced/Dynamic.hs b/src/Control/Monad/Bayes/Traced/Dynamic.hs index 4609394c..58ef231b 100644 --- a/src/Control/Monad/Bayes/Traced/Dynamic.hs +++ b/src/Control/Monad/Bayes/Traced/Dynamic.hs @@ -9,7 +9,7 @@ -- Stability : experimental -- Portability : GHC module Control.Monad.Bayes.Traced.Dynamic - ( Traced, + ( TracedT, hoist, marginal, freeze, @@ -24,7 +24,7 @@ import Control.Monad.Bayes.Class MonadFactor (..), MonadMeasure, ) -import Control.Monad.Bayes.Density.Free (Density) +import Control.Monad.Bayes.Density.Free (DensityT) import Control.Monad.Bayes.Traced.Common ( Trace (..), bind, @@ -32,75 +32,75 @@ import Control.Monad.Bayes.Traced.Common scored, singleton, ) -import Control.Monad.Bayes.Weighted (Weighted) +import Control.Monad.Bayes.Weighted (WeightedT) import Control.Monad.Trans (MonadTrans (..)) import Data.List.NonEmpty as NE (NonEmpty ((:|)), toList) -- | A tracing monad where only a subset of random choices are traced and this -- subset can be adjusted dynamically. -newtype Traced m a = Traced {runTraced :: m (Weighted (Density m) a, Trace a)} +newtype TracedT m a = TracedT {runTraced :: m (WeightedT (DensityT m) a, Trace a)} -pushM :: (Monad m) => m (Weighted (Density m) a) -> Weighted (Density m) a +pushM :: (Monad m) => m (WeightedT (DensityT m) a) -> WeightedT (DensityT m) a pushM = join . lift . lift -instance (Monad m) => Functor (Traced m) where - fmap f (Traced c) = Traced $ do +instance (Monad m) => Functor (TracedT m) where + fmap f (TracedT c) = TracedT $ do (m, t) <- c let m' = fmap f m let t' = fmap f t return (m', t') -instance (Monad m) => Applicative (Traced m) where - pure x = Traced $ pure (pure x, pure x) - (Traced cf) <*> (Traced cx) = Traced $ do +instance (Monad m) => Applicative (TracedT m) where + pure x = TracedT $ pure (pure x, pure x) + (TracedT cf) <*> (TracedT cx) = TracedT $ do (mf, tf) <- cf (mx, tx) <- cx return (mf <*> mx, tf <*> tx) -instance (Monad m) => Monad (Traced m) where - (Traced cx) >>= f = Traced $ do +instance (Monad m) => Monad (TracedT m) where + (TracedT cx) >>= f = TracedT $ do (mx, tx) <- cx let m = mx >>= pushM . fmap fst . runTraced . f t <- return tx `bind` (fmap snd . runTraced . f) return (m, t) -instance MonadTrans Traced where - lift m = Traced $ fmap ((,) (lift $ lift m) . pure) m +instance MonadTrans TracedT where + lift m = TracedT $ fmap ((,) (lift $ lift m) . pure) m -instance (MonadDistribution m) => MonadDistribution (Traced m) where - random = Traced $ fmap ((,) random . singleton) random +instance (MonadDistribution m) => MonadDistribution (TracedT m) where + random = TracedT $ fmap ((,) random . singleton) random -instance (MonadFactor m) => MonadFactor (Traced m) where - score w = Traced $ fmap (score w,) (score w >> pure (scored w)) +instance (MonadFactor m) => MonadFactor (TracedT m) where + score w = TracedT $ fmap (score w,) (score w >> pure (scored w)) -instance (MonadMeasure m) => MonadMeasure (Traced m) +instance (MonadMeasure m) => MonadMeasure (TracedT m) -hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a -hoist f (Traced c) = Traced (f c) +hoist :: (forall x. m x -> m x) -> TracedT m a -> TracedT m a +hoist f (TracedT c) = TracedT (f c) -- | Discard the trace and supporting infrastructure. -marginal :: (Monad m) => Traced m a -> m a -marginal (Traced c) = fmap (output . snd) c +marginal :: (Monad m) => TracedT m a -> m a +marginal (TracedT c) = fmap (output . snd) c -- | Freeze all traced random choices to their current values and stop tracing -- them. -freeze :: (Monad m) => Traced m a -> Traced m a -freeze (Traced c) = Traced $ do +freeze :: (Monad m) => TracedT m a -> TracedT m a +freeze (TracedT c) = TracedT $ do (_, t) <- c let x = output t return (return x, pure x) -- | A single step of the Trace Metropolis-Hastings algorithm. -mhStep :: (MonadDistribution m) => Traced m a -> Traced m a -mhStep (Traced c) = Traced $ do +mhStep :: (MonadDistribution m) => TracedT m a -> TracedT m a +mhStep (TracedT c) = TracedT $ do (m, t) <- c t' <- mhTransFree m t return (m, t') -- | Full run of the Trace Metropolis-Hastings algorithm with a specified -- number of steps. -mh :: (MonadDistribution m) => Int -> Traced m a -> m [a] -mh n (Traced c) = do +mh :: (MonadDistribution m) => Int -> TracedT m a -> m [a] +mh n (TracedT c) = do (m, t) <- c let f k | k <= 0 = return (t :| []) diff --git a/src/Control/Monad/Bayes/Traced/Static.hs b/src/Control/Monad/Bayes/Traced/Static.hs index a635f327..fc99b327 100644 --- a/src/Control/Monad/Bayes/Traced/Static.hs +++ b/src/Control/Monad/Bayes/Traced/Static.hs @@ -10,7 +10,7 @@ -- Stability : experimental -- Portability : GHC module Control.Monad.Bayes.Traced.Static - ( Traced (..), + ( TracedT (..), hoist, marginal, mhStep, @@ -24,7 +24,7 @@ import Control.Monad.Bayes.Class MonadFactor (..), MonadMeasure, ) -import Control.Monad.Bayes.Density.Free (Density) +import Control.Monad.Bayes.Density.Free (DensityT) import Control.Monad.Bayes.Traced.Common ( Trace (..), bind, @@ -32,7 +32,7 @@ import Control.Monad.Bayes.Traced.Common scored, singleton, ) -import Control.Monad.Bayes.Weighted (Weighted) +import Control.Monad.Bayes.Weighted (WeightedT) import Control.Monad.Trans (MonadTrans (..)) import Data.List.NonEmpty as NE (NonEmpty ((:|)), toList) @@ -40,45 +40,45 @@ import Data.List.NonEmpty as NE (NonEmpty ((:|)), toList) -- -- The random choices that are not to be traced should be lifted from the -- transformed monad. -data Traced m a = Traced - { model :: Weighted (Density m) a, +data TracedT m a = TracedT + { model :: WeightedT (DensityT m) a, traceDist :: m (Trace a) } -instance (Monad m) => Functor (Traced m) where - fmap f (Traced m d) = Traced (fmap f m) (fmap (fmap f) d) +instance (Monad m) => Functor (TracedT m) where + fmap f (TracedT m d) = TracedT (fmap f m) (fmap (fmap f) d) -instance (Monad m) => Applicative (Traced m) where - pure x = Traced (pure x) (pure (pure x)) - (Traced mf df) <*> (Traced mx dx) = Traced (mf <*> mx) (liftA2 (<*>) df dx) +instance (Monad m) => Applicative (TracedT m) where + pure x = TracedT (pure x) (pure (pure x)) + (TracedT mf df) <*> (TracedT mx dx) = TracedT (mf <*> mx) (liftA2 (<*>) df dx) -instance (Monad m) => Monad (Traced m) where - (Traced mx dx) >>= f = Traced my dy +instance (Monad m) => Monad (TracedT m) where + (TracedT mx dx) >>= f = TracedT my dy where my = mx >>= model . f dy = dx `bind` (traceDist . f) -instance MonadTrans Traced where - lift m = Traced (lift $ lift m) (fmap pure m) +instance MonadTrans TracedT where + lift m = TracedT (lift $ lift m) (fmap pure m) -instance (MonadDistribution m) => MonadDistribution (Traced m) where - random = Traced random (fmap singleton random) +instance (MonadDistribution m) => MonadDistribution (TracedT m) where + random = TracedT random (fmap singleton random) -instance (MonadFactor m) => MonadFactor (Traced m) where - score w = Traced (score w) (score w >> pure (scored w)) +instance (MonadFactor m) => MonadFactor (TracedT m) where + score w = TracedT (score w) (score w >> pure (scored w)) -instance (MonadMeasure m) => MonadMeasure (Traced m) +instance (MonadMeasure m) => MonadMeasure (TracedT m) -hoist :: (forall x. m x -> m x) -> Traced m a -> Traced m a -hoist f (Traced m d) = Traced m (f d) +hoist :: (forall x. m x -> m x) -> TracedT m a -> TracedT m a +hoist f (TracedT m d) = TracedT m (f d) -- | Discard the trace and supporting infrastructure. -marginal :: (Monad m) => Traced m a -> m a -marginal (Traced _ d) = fmap output d +marginal :: (Monad m) => TracedT m a -> m a +marginal (TracedT _ d) = fmap output d -- | A single step of the Trace Metropolis-Hastings algorithm. -mhStep :: (MonadDistribution m) => Traced m a -> Traced m a -mhStep (Traced m d) = Traced m d' +mhStep :: (MonadDistribution m) => TracedT m a -> TracedT m a +mhStep (TracedT m d) = TracedT m d' where d' = d >>= mhTransFree m @@ -111,8 +111,8 @@ mhStep (Traced m d) = Traced m d' -- [True,True,True] -- -- Of course, it will need to be run more than twice to get a reasonable estimate. -mh :: (MonadDistribution m) => Int -> Traced m a -> m [a] -mh n (Traced m d) = fmap (map output . NE.toList) (f n) +mh :: (MonadDistribution m) => Int -> TracedT m a -> m [a] +mh n (TracedT m d) = fmap (map output . NE.toList) (f n) where f k | k <= 0 = fmap (:| []) d diff --git a/src/Control/Monad/Bayes/Weighted.hs b/src/Control/Monad/Bayes/Weighted.hs index 8c7e9ae5..a1bbe44a 100644 --- a/src/Control/Monad/Bayes/Weighted.hs +++ b/src/Control/Monad/Bayes/Weighted.hs @@ -11,17 +11,16 @@ -- Stability : experimental -- Portability : GHC -- --- 'Weighted' is an instance of 'MonadFactor'. Apply a 'MonadDistribution' transformer to +-- 'WeightedT' is an instance of 'MonadFactor'. Apply a 'MonadDistribution' transformer to -- obtain a 'MonadMeasure' that can execute probabilistic models. module Control.Monad.Bayes.Weighted - ( Weighted, - withWeight, - weighted, + ( WeightedT, + weightedT, extractWeight, unweighted, applyWeight, hoist, - runWeighted, + runWeightedT, ) where @@ -35,46 +34,45 @@ import Control.Monad.State (MonadIO, MonadTrans, StateT (..), lift, mapStateT, m import Numeric.Log (Log) -- | Execute the program using the prior distribution, while accumulating likelihood. -newtype Weighted m a = Weighted (StateT (Log Double) m a) +newtype WeightedT m a = WeightedT (StateT (Log Double) m a) -- StateT is more efficient than WriterT deriving newtype (Functor, Applicative, Monad, MonadIO, MonadTrans, MonadDistribution) -instance (Monad m) => MonadFactor (Weighted m) where - score w = Weighted (modify (* w)) +instance (Monad m) => MonadFactor (WeightedT m) where + score w = WeightedT (modify (* w)) -instance (MonadDistribution m) => MonadMeasure (Weighted m) +instance (MonadDistribution m) => MonadMeasure (WeightedT m) -- | Obtain an explicit value of the likelihood for a given value. -weighted, runWeighted :: Weighted m a -> m (a, Log Double) -weighted (Weighted m) = runStateT m 1 -runWeighted = weighted +runWeightedT :: WeightedT m a -> m (a, Log Double) +runWeightedT (WeightedT m) = runStateT m 1 -- | Compute the sample and discard the weight. -- -- This operation introduces bias. -unweighted :: (Functor m) => Weighted m a -> m a -unweighted = fmap fst . weighted +unweighted :: (Functor m) => WeightedT m a -> m a +unweighted = fmap fst . runWeightedT -- | Compute the weight and discard the sample. -extractWeight :: (Functor m) => Weighted m a -> m (Log Double) -extractWeight = fmap snd . weighted +extractWeight :: (Functor m) => WeightedT m a -> m (Log Double) +extractWeight = fmap snd . runWeightedT -- | Embed a random variable with explicitly given likelihood. -- --- > weighted . withWeight = id -withWeight :: (Monad m) => m (a, Log Double) -> Weighted m a -withWeight m = Weighted $ do +-- > runWeightedT . weightedT = id +weightedT :: (Monad m) => m (a, Log Double) -> WeightedT m a +weightedT m = WeightedT $ do (x, w) <- lift m modify (* w) return x -- | Use the weight as a factor in the transformed monad. -applyWeight :: (MonadFactor m) => Weighted m a -> m a +applyWeight :: (MonadFactor m) => WeightedT m a -> m a applyWeight m = do - (x, w) <- weighted m + (x, w) <- runWeightedT m factor w return x -- | Apply a transformation to the transformed monad. -hoist :: (forall x. m x -> n x) -> Weighted m a -> Weighted n a -hoist t (Weighted m) = Weighted $ mapStateT t m +hoist :: (forall x. m x -> n x) -> WeightedT m a -> WeightedT n a +hoist t (WeightedT m) = WeightedT $ mapStateT t m diff --git a/test/Spec.hs b/test/Spec.hs index 6557feba..dcdb756c 100644 --- a/test/Spec.hs +++ b/test/Spec.hs @@ -57,7 +57,7 @@ main = hspec do -- Because of rounding issues, require the variance to be a bit bigger than 0 -- See https://github.com/tweag/monad-bayes/issues/275 var > 0.1 ==> property $ TestIntegrator.normalVariance mean (sqrt var) ~== var - describe "Sampler mean and variance" do + describe "SamplerT mean and variance" do it "gets right mean and variance" $ TestSampler.testMeanAndVariance `shouldBe` True describe "Integrator Volume" do diff --git a/test/TestAdvanced.hs b/test/TestAdvanced.hs index bb7569dc..41222efb 100644 --- a/test/TestAdvanced.hs +++ b/test/TestAdvanced.hs @@ -23,16 +23,16 @@ 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 + sample <- sampleIOfixed $ runPopulationT $ smc (SMCConfig {numSteps = 0, numParticles = 10000, resampler = resampleMultinomial}) random return $ close 0.5 sample passed3 = do - sample <- sampleIOfixed $ population $ rmsmcDynamic mcmcConfig smcConfig random + sample <- sampleIOfixed $ runPopulationT $ rmsmcDynamic mcmcConfig smcConfig random return $ close 0.5 sample passed4 = do - sample <- sampleIOfixed $ population $ rmsmcBasic mcmcConfig smcConfig random + sample <- sampleIOfixed $ runPopulationT $ rmsmcBasic mcmcConfig smcConfig random return $ close 0.5 sample passed5 = do - sample <- sampleIOfixed $ population $ rmsmc mcmcConfig smcConfig random + sample <- sampleIOfixed $ runPopulationT $ rmsmc mcmcConfig smcConfig random return $ close 0.5 sample passed6 = do sample <- @@ -48,7 +48,7 @@ passed6 = do close :: Double -> [(Double, Log Double)] -> Bool passed7 = do - sample <- fmap join $ sampleIOfixed $ fmap (fmap (\(x, y) -> fmap (second (* y)) x)) $ population $ smc2 0 100 100 100 random (normal 0) + sample <- fmap join $ sampleIOfixed $ fmap (fmap (\(x, y) -> fmap (second (* y)) x)) $ runPopulationT $ smc2 0 100 100 100 random (normal 0) return $ close 0.0 sample close n sample = abs (n - (expectation id $ fromList $ toEmpiricalWeighted sample)) < 0.01 diff --git a/test/TestInference.hs b/test/TestInference.hs index d09eecba..64340b75 100644 --- a/test/TestInference.hs +++ b/test/TestInference.hs @@ -20,10 +20,10 @@ 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.Strict (Sampler, sampleIOfixed) +import Control.Monad.Bayes.Sampler.Strict (SamplerT, sampleIOfixed) import Control.Monad.Bayes.Sampler.Strict qualified as Sampler -import Control.Monad.Bayes.Weighted (Weighted) -import Control.Monad.Bayes.Weighted qualified as Weighted +import Control.Monad.Bayes.Weighted (WeightedT) +import Control.Monad.Bayes.Weighted qualified as WeightedT import Data.AEq (AEq ((~==))) import Numeric.Log (Log) import Sprinkler (soft) @@ -35,18 +35,18 @@ sprinkler = Sprinkler.soft -- | Count the number of particles produced by SMC checkParticles :: Int -> Int -> IO Int checkParticles observations particles = - sampleIOfixed (fmap length (population $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleMultinomial} Sprinkler.soft)) + sampleIOfixed (fmap length (runPopulationT $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleMultinomial} Sprinkler.soft)) checkParticlesSystematic :: Int -> Int -> IO Int checkParticlesSystematic observations particles = - sampleIOfixed (fmap length (population $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleSystematic} Sprinkler.soft)) + sampleIOfixed (fmap length (runPopulationT $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleSystematic} Sprinkler.soft)) checkParticlesStratified :: Int -> Int -> IO Int checkParticlesStratified observations particles = - sampleIOfixed (fmap length (population $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleStratified} Sprinkler.soft)) + sampleIOfixed (fmap length (runPopulationT $ smc SMCConfig {numSteps = observations, numParticles = particles, resampler = resampleStratified} Sprinkler.soft)) checkTerminateSMC :: IO [(Bool, Log Double)] -checkTerminateSMC = sampleIOfixed (population $ smc SMCConfig {numSteps = 2, numParticles = 5, resampler = resampleMultinomial} sprinkler) +checkTerminateSMC = sampleIOfixed (runPopulationT $ smc SMCConfig {numSteps = 2, numParticles = 5, resampler = resampleMultinomial} sprinkler) checkPreserveSMC :: Bool checkPreserveSMC = @@ -54,8 +54,8 @@ checkPreserveSMC = ~== enumerator sprinkler expectationNearNumeric :: - Weighted Integrator.Integrator Double -> - Weighted Integrator.Integrator Double -> + WeightedT Integrator.Integrator Double -> + WeightedT Integrator.Integrator Double -> Double expectationNearNumeric x y = let e1 = Integrator.expectation $ normalize x @@ -63,12 +63,12 @@ expectationNearNumeric x y = in (abs (e1 - e2)) expectationNearSampling :: - Weighted (Sampler (IOGenM StdGen) IO) Double -> - Weighted (Sampler (IOGenM StdGen) IO) Double -> + WeightedT (SamplerT (IOGenM StdGen) IO) Double -> + WeightedT (SamplerT (IOGenM StdGen) IO) Double -> IO Double expectationNearSampling x y = do - e1 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ Weighted.weighted x - e2 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ Weighted.weighted y + e1 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ WeightedT.runWeightedT x + e2 <- sampleIOfixed $ fmap Sampler.sampleMean $ replicateM 10 $ WeightedT.runWeightedT y return (abs (e1 - e2)) testNormalNormal :: [Double] -> IO Bool diff --git a/test/TestIntegrator.hs b/test/TestIntegrator.hs index 481b8e9a..efe54bb6 100644 --- a/test/TestIntegrator.hs +++ b/test/TestIntegrator.hs @@ -13,7 +13,7 @@ import Control.Monad.Bayes.Class ) import Control.Monad.Bayes.Integrator import Control.Monad.Bayes.Sampler.Strict -import Control.Monad.Bayes.Weighted (weighted) +import Control.Monad.Bayes.Weighted (runWeightedT) import Control.Monad.ST (runST) import Data.AEq (AEq ((~==))) import Data.List (sortOn) @@ -98,7 +98,7 @@ passed7 = passed8 = 1 == ( volume $ - fmap (ln . exp . snd) $ weighted do + fmap (ln . exp . snd) $ runWeightedT do x <- bernoulli 0.5 factor $ if x then 0.2 else 0.1 return x @@ -137,7 +137,7 @@ passed13 = ~== 1 -- sampler and integrator agree on a non-trivial model passed14 = - let sample = runST $ sampleSTfixed $ fmap sampleMean $ replicateM 10000 $ weighted $ model1 + let sample = runST $ sampleSTfixed $ fmap sampleMean $ replicateM 10000 $ runWeightedT $ model1 quadrature = expectation $ normalize $ model1 in abs (sample - quadrature) < 0.01 diff --git a/test/TestPopulation.hs b/test/TestPopulation.hs index 33ee289d..ccf86ffc 100644 --- a/test/TestPopulation.hs +++ b/test/TestPopulation.hs @@ -3,20 +3,20 @@ module TestPopulation (weightedSampleSize, popSize, manySize, sprinkler, sprinkl import Control.Monad.Bayes.Class (MonadDistribution, MonadMeasure) import Control.Monad.Bayes.Enumerator (enumerator, expectation) import Control.Monad.Bayes.Population as Population - ( Population, + ( PopulationT, collapse, popAvg, - population, pushEvidence, resampleMultinomial, + runPopulationT, spawn, ) import Control.Monad.Bayes.Sampler.Strict (sampleIOfixed) import Data.AEq (AEq ((~==))) import Sprinkler (soft) -weightedSampleSize :: (MonadDistribution m) => Population m a -> m Int -weightedSampleSize = fmap length . population +weightedSampleSize :: (MonadDistribution m) => PopulationT m a -> m Int +weightedSampleSize = fmap length . runPopulationT popSize :: IO Int popSize = diff --git a/test/TestWeighted.hs b/test/TestWeighted.hs index 31cb16f0..3e068e0f 100644 --- a/test/TestWeighted.hs +++ b/test/TestWeighted.hs @@ -8,7 +8,7 @@ import Control.Monad.Bayes.Class factor, ) import Control.Monad.Bayes.Sampler.Strict (sampleIOfixed) -import Control.Monad.Bayes.Weighted (weighted) +import Control.Monad.Bayes.Weighted (runWeightedT) import Control.Monad.State (unless, when) import Data.AEq (AEq ((~==))) import Data.Bifunctor (second) @@ -23,7 +23,7 @@ model = do return (n, x) result :: (MonadDistribution m) => m ((Int, Double), Double) -result = second (exp . ln) <$> weighted model +result = second (exp . ln) <$> runWeightedT model passed :: IO Bool passed = fmap check (sampleIOfixed result) From b73c22d7086daff82ca459d5a534df65863bb962 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Manuel=20B=C3=A4renz?= Date: Wed, 4 Oct 2023 09:49:37 +0200 Subject: [PATCH 2/2] Update notebook htmls --- docs/docs/notebooks/Diagrams.html | 4 ++-- docs/docs/notebooks/Histogram.html | 2 +- docs/docs/notebooks/Parsing.html | 2 +- docs/docs/notebooks/Streaming.html | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/docs/notebooks/Diagrams.html b/docs/docs/notebooks/Diagrams.html index bbc54d16..49779584 100644 --- a/docs/docs/notebooks/Diagrams.html +++ b/docs/docs/notebooks/Diagrams.html @@ -14729,8 +14729,8 @@

Viewing tracesIn [3]: