Skip to content

Commit

Permalink
Lazyppl (#175)
Browse files Browse the repository at this point in the history
* analytic

* docs

* docs

* docs

* docs

* inference

* executable docs

* more comments, added BayesianModel

* Remove MaybeT and use ExceptT

* remove vscode and update gitignore

* delete additional changes: Failure

* delete broken files

* update test for enumerator

* update models

* tips on inference

* density

* analytic

* analytic

* more pipes

* update docs

* linear outliers as in Gen tutorial

* update docs

* cleanup SMC.hs

* more cleanup, and explicit imports

* remove Rejection

* cleanup

* enumerate

* binning

* comment

* explore

* nested inference example

* found bug in Gamma

* debugging

* port lazyppl monad into monad-bayes

* mh

* add conjugacy tests

* working branch

* tried symbolic

* no symbolic

* remove reflection

* midway: need to finish hmm and enumerator:

* update docs

* build docs manually

* insegel theme

* stanford theme

* lint

* lint

* ormulo

* nix ormolu fixes:
git push

* docs

* trying to fix mathjax

* remove empirical

* notebooks

* notebooks: sampling

* remove helpers

* models: lts

* models: lts

* weighted sampling

* formatting

* test enumerator

* first attempt to fix normalization bug

* update tests

* working tests: normalNormal goes awry if observations are huge

* cleanup

* Bayesian in class, and remove LinearWithOutliers

* model updates

* yaml

* prepare for PR

* address warnings

* review

* more warnings fixed

* less warnings

* update ghc version stated in cabal

* remove warnings

* fix build

* version bump

* remove safe

* fix merge

* fix rmsmc bug you introduced

* integrator as MonadInfer

* updates

* updates

* updates

* fix the integrator

* inference

* rope in analyticT

* remove integratorT

* nix notebooks

* notebooks

* rmsmc

* hoistFirst bug again

* hoistFirst bug again

* diagrams

* fp notebook

* remove integratorT

* presentation

* notebooks

* cleanup

* docs

* docs

* presentation

* update talk

* notebook

* works, but readd basic

* first working example of marginalization

* renaming

* continue cleanup

* tests: not yet integrated in

* tests: integrated in

* bayesian pmmh

* switch over to new smc

* switch over to new smc

* improve resampleGeneric sum

* independent in Class

* parser

* nbs

* make clean

* fix tests

* fix tests

* fix tests

* fix tests

* pmmh

* withParticles

* withParticles

* Sequential.Free

* Sequential.Free

* docs

* update nix

* cabal format

* update nix

* fix build

* bump changelog and version. bump ghc version to 9.2.3

* update docs and change mcmc so unweighted is not needed

* remove sequentialCata

* nbs

* merge

* updating notebooks

* mhTrans'

* fix mh

* infinite populations

* simple density

* merge

* lint

* lint

* lazy sampler

* clean up

Co-authored-by: Reuben Cohn-Gordon <reubencohn-gordon@Reubens-MacBook-Air.local>
Co-authored-by: System Administrator <root@Reubens-MacBook-Air.local>
Co-authored-by: Dominic Steinitz <dominic@steinitz.org>
  • Loading branch information
4 people authored Aug 30, 2022
1 parent 099dbea commit d8d3beb
Show file tree
Hide file tree
Showing 81 changed files with 99,193 additions and 59,125 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
venv
_cache
docs/build
*.csv
Expand Down
3 changes: 3 additions & 0 deletions .netlify/state.json
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
{
"siteId": "12c0eb5b-a921-45cf-af85-b903b65b801c"
}
4 changes: 2 additions & 2 deletions benchmark/SSM.hs
Original file line number Diff line number Diff line change
Expand Up @@ -6,8 +6,8 @@ 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, runPopulation)
import Control.Monad.Bayes.Sampler (sampleIO, sampleIOfixed, sampleWith)
import Control.Monad.Bayes.Population (population, resampleMultinomial)
import Control.Monad.Bayes.Sampler.Strict (sampleIO, sampleIOfixed, sampleWith)
import Control.Monad.Bayes.Weighted (unweighted)
import Control.Monad.IO.Class (MonadIO (liftIO))
import NonlinearSSM (generateData, model, param)
Expand Down
28 changes: 21 additions & 7 deletions benchmark/Single.hs
Original file line number Diff line number Diff line change
@@ -1,21 +1,36 @@
{-# LANGUAGE DerivingStrategies #-}
{-# LANGUAGE ImportQualifiedPost #-}

import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Class (MonadInfer)
import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (..), Proposal (SingleSiteMH))
import Control.Monad.Bayes.Inference.RMSMC
import Control.Monad.Bayes.Inference.RMSMC (rmsmcBasic)
import Control.Monad.Bayes.Inference.SMC
( SMCConfig (SMCConfig, numParticles, numSteps, resampler),
smc,
)
import Control.Monad.Bayes.Population
import Control.Monad.Bayes.Population (population)
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Sampler.Strict
import Control.Monad.Bayes.Traced
import Control.Monad.Bayes.Weighted
import Control.Monad.ST (runST)
import Data.Time
import Data.Time (diffUTCTime, getCurrentTime)
import HMM qualified
import LDA qualified
import LogReg qualified
import Options.Applicative
( Applicative (liftA2),
ParserInfo,
auto,
execParser,
fullDesc,
help,
info,
long,
maybeReader,
option,
short,
)
import System.Random.MWC (GenIO, createSystemRandom)

data Model = LR Int | HMM Int | LDA (Int, Int)
Expand All @@ -42,7 +57,7 @@ getModel model = (size model, program model)
data Alg = SMC | MH | RMSMC
deriving stock (Read, Show)

runAlg :: Model -> Alg -> Sampler GenIO IO String
runAlg :: Model -> Alg -> SamplerIO String
runAlg model alg =
case alg of
SMC ->
Expand All @@ -61,8 +76,7 @@ runAlg model alg =

infer :: Model -> Alg -> IO ()
infer model alg = do
g <- createSystemRandom
x <- sampleWith (runAlg model alg) g
x <- sampleIOfixed (runAlg model alg)
print x

opts :: ParserInfo (Model, Alg)
Expand Down
4 changes: 2 additions & 2 deletions benchmark/Speed.hs
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,8 @@ import Control.Monad.Bayes.Class (MonadInfer, MonadSample)
import Control.Monad.Bayes.Inference.MCMC (MCMCConfig (MCMCConfig, numBurnIn, numMCMCSteps, proposal), Proposal (SingleSiteMH))
import Control.Monad.Bayes.Inference.RMSMC (rmsmcDynamic)
import Control.Monad.Bayes.Inference.SMC (SMCConfig (SMCConfig, numParticles, numSteps, resampler), smc)
import Control.Monad.Bayes.Population (population, resampleSystematic, runPopulation)
import Control.Monad.Bayes.Sampler (SamplerIO, sampleIOfixed, sampleWith)
import Control.Monad.Bayes.Population (population, resampleSystematic)
import Control.Monad.Bayes.Sampler.Strict (SamplerIO, sampleIOfixed)
import Control.Monad.Bayes.Traced (mh)
import Control.Monad.Bayes.Weighted (unweighted)
import Criterion.Main
Expand Down
24 changes: 22 additions & 2 deletions docs/source/probprog.md
Original file line number Diff line number Diff line change
Expand Up @@ -261,7 +261,7 @@ which gives
[([1,2,3,4],0.5),([2,3,4,5],0.5)]
```

### Near exact inference for continuous distributions
## Near exact inference for continuous distributions

Monad-Bayes does not currently support exact inference (via symbolic solving) for continuous distributions. However, it *does* support numerical integration. For example, for the distribution defined by

Expand Down Expand Up @@ -289,7 +289,7 @@ model = do

we must first `normalize` the model, as in `probability (0, 0.1) (normalize model)`.

### Independent forward sampling
## Independent forward sampling

For any probabilistic program `p` without any `condition` or `factor` statements, we may do `sampler p` or `sampleIOfixed p` (to run with a fixed seed) to obtain a sample in an ancestral fashion. For example, consider:

Expand Down Expand Up @@ -335,6 +335,26 @@ run = (sampler . weighted) example
is an IO operation which when run, will display either `(False, 0.0)` or `(True, 1.0)`


## Lazy sampling

If you want to forward sample from an infinite program, just as a distribution over infinite lists, you can use monad-bayes's lazy sampler, which is based on LazyPPL. For example,

```haskell
import qualified Control.Monad.Bayes.Sampler.Lazy as Lazy

example :: MonadSample m => m [Double]
example = do
x <- random
fmap (x:) example

infiniteList <- Lazy.sampler example
take 4 infiniteList
```

To perform weighted sampling, use `lwis` from `Control.Monad.Bayes.Inference.Lazy.WIS` as in `lwis 10 example`. This takes 10 weighted samples, and produces an infinite stream of samples, regarding those 10 as an empirical distribution.

LazyPPL's `mh` implementation is also available.

## Markov Chain Monte Carlo

There are several versions of metropolis hastings MCMC defined in monad-bayes. The standard version is found in `Control.Monad.Bayes.Inference.MCMC`. You can use it as follows:
Expand Down
65 changes: 31 additions & 34 deletions docs/source/usage.md
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,7 @@ Summary of key info on `Sequential`:
- `instance MonadSample m => instance MonadSample (Sequential m)`
- `instance MonadCond m => instance MonadCond (Sequential m)`


```haskell
newtype Sequential m a =
Sequential {runSequential :: Coroutine (Await ()) m a}
Expand Down Expand Up @@ -474,56 +475,53 @@ hoistFirst :: (forall x. m x -> m x) -> Sequential m a -> Sequential m a
hoistFirst f = Sequential . Coroutine . f . resume . runSequential
```

<!-- As an example, consider:
TODO: Enumerator example with `trace` -->

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).




### FreeSampler
### Density

Summary of key info on `FreeSampler`:
Summary of key info on `Density`:

- `FreeSampler :: (Type -> Type) -> (Type -> Type)`
- `instance MonadSample (FreeSampler m)`
- `Density :: (Type -> Type) -> (Type -> Type)`
- `instance MonadSample (Density m)`
- **No** instance for `MonadCond`

`FreeSampler m` is not often going to be used on its own, but instead as part of the `Traced` type, defined below. A `FreeSampler m a` represents a reified execution of the program.
A *trace* of a program of type `MonadSample 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.

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.

`FreeSampler m` is best understood if you're familiar with the standard use of a free monad to construct a domain specific language. For probability in particular, see this [blog post](https://jtobin.io/simple-probabilistic-programming). Here's the definition:
The former is relatively straightforward: the `MonadSample` 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`.

The latter is best understood if you're familiar with the standard use of a free monad to construct a domain specific language. For probability in particular, see this [blog post](https://jtobin.io/simple-probabilistic-programming). Here's the definition:

```haskell
newtype SamF a = Random (Double -> a)
newtype FreeSampler m a =
FreeSampler {runFreeSampler :: FT SamF m a}
newtype Density m a =
Density {density :: FT SamF m a}

instance Monad m => MonadSample (FreeSampler m) where
random = FreeSampler $ liftF (Random id)
instance Monad m => MonadSample (Density m) where
random = Density $ 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.

If you unpack the definition, you get:

```haskell
FreeSampler m a ~ m (Either a (Double -> (FreeSampler m a)))
Density m a ~ m (Either a (Double -> (Density m a)))
```

As you can see, this is rather like `Coroutine`, except to "resume", you must provide a new `Double`, corresponding to the value of some particular random choice.

Since `FreeT` is a transformer, we can use `lift` to get a `MonadSample` instance.


A *trace* of a program of type `MonadSample 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`.

Given a probabilistic program interpreted in `FreeSampler m`, we can "run" it to produce a program in the underlying monad `m`. For simplicity, consider the case of a program `bernoulli 0.5 :: FreeSampler SamplerIO Bool`. We can then use the following function:
`density` is then defined using the folding pattern `iterFT`, which interprets `SamF` in the appropriate way:

```haskell
withPartialRandomness :: MonadSample m => [Double] -> FreeSampler m a -> m (a, [Double])
withPartialRandomness randomness (FreeSampler m) =
density :: MonadSample m => [Double] -> Density m a -> m (a, [Double])
density randomness (Density m) =
runWriterT $ evalStateT (iterTM f $ hoistFT lift m) randomness
where
f (Random k) = do
Expand All @@ -538,7 +536,7 @@ withPartialRandomness randomness (FreeSampler m) =
k x
```

This takes a list of `Double`s (a representation of a trace), and a probabilistic program like `example`, and gives back a `SamplerIO (Bool, [Double])`. At each call to `random` in `example`, the next double in the list is used. If the list of doubles runs out, calls are made to `random` using the underlying monad, which in our example is `SamplerIO`. Hence "with*Partial*Randomness".
This takes a list of `Double`s (a representation of a trace), and a probabilistic program like `example`, and gives back a `SamplerIO (Bool, [Double])`. At each call to `random` in `example`, the next double in the list is used. If the list of doubles runs out, calls are made to `random` using the underlying monad.

<!-- The intuition here is that given a list of doubles in $[0,1]$, you can evaluate any probabilistic program. If your list of numbers is shorter than the number of calls to `random` in the program, the remaining calls are made in the underlying `MonadSample` instance `m`. -->

Expand All @@ -554,7 +552,7 @@ Summary of key info on `Traced`:
- `instance MonadSample m => MonadSample (Traced m)`
- `instance MonadCond m => MonadCond (Traced m)`

`Traced m` is actually several related interpretations, each built on top of `FreeSampler`. These range in complexity.
`Traced m` is actually several related interpretations, each built on top of `Density`. These range in complexity.



Expand All @@ -576,12 +574,12 @@ data Trace a = Trace
}
```

We also need a specification of the probabilistic program in question, free of any particular interpretation. That is precisely what `FreeSampler` is for.
We also need a specification of the probabilistic program in question, free of any particular interpretation. That is precisely what `Density` is for.

The simplest version of `Traced` is in `Control.Monad.Bayes.Traced.Basic`

```haskell
Traced m a ~ (FreeSampler Identity a, Log Double), m (Trace a))
Traced m a ~ (Density 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.
Expand Down Expand Up @@ -707,7 +705,7 @@ A single step in this chain (in Metropolis Hasting MCMC) looks like this:

```haskell
mhTrans :: MonadSample m =>
Weighted (FreeSampler m) a -> Trace a -> m (Trace a)
Weighted (Density m) a -> Trace a -> m (Trace a)
mhTrans m t@Trace {variables = us, density = p} = do
let n = length us
us' <- do
Expand All @@ -717,15 +715,14 @@ mhTrans m t@Trace {variables = us, density = p} = do
(xs, _ : ys) -> return $ xs ++ (u' : ys)
_ -> error "impossible"
((b, q), vs) <-
runWriterT $ weighted
$ Weighted.hoist (WriterT . withPartialRandomness us') m
runWriterT $ weighted $ Weighted.hoist (WriterT . density us') m
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 (FreeSampler m) a`, which is an instance of `MonadInfer`. 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 `Weighted (Density m) a`, which is an instance of `MonadInfer`. 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.

Expand All @@ -736,18 +733,18 @@ MH is then easily defined as taking steps with this kernel, in the usual fashion
This is provided by

```haskell
sis ::
sequentially ::
Monad m =>
-- | transformation
(forall x. m x -> m x) ->
-- | number of time steps
Int ->
Sequential m a ->
m a
sis f k = finish . composeCopies k (advance . hoistFirst f)
sequentially f k = finish . composeCopies k (advance . hoistFirst f)
```

in Control.Monad.Bayes.Sequential. You provide a natural transformation in the underlying monad `m`, and `sis` applies that natural transformation at each point of conditioning in your program. The main use case is in defining `smc`, below, but here is a nice alternative use case:
in `Control.Monad.Bayes.Sequential.Coroutine`. You provide a natural transformation in the underlying monad `m`, and `sequentially` applies that natural transformation at each point of conditioning in your program. The main use case is in defining `smc`, below, but here is a nice didactic use case:

Consider the program:

Expand Down
Loading

0 comments on commit d8d3beb

Please sign in to comment.