tweag/monad-bayes

A library for probabilistic programming in Haskell.

Monad-Bayes

A library for probabilistic programming in Haskell.

See the docs for a user guide, notebook-style tutorials, an example gallery, and a detailed account of the implementation.

Created by Adam Scibior (@adscib), documentation, website and newer features by Reuben, maintained by Tweag.

Project status

Now that monad-bayes has been released on Hackage, and the documentation and the API has been updated, we will focus on adding new features. See the Github issues to get a sense of what is being prepared, and please feel free to make requests.

Background

The basis for the code in this repository is the ICFP 2018 paper [2]. For the code associated with the Haskell2015 paper [1], see the haskell2015 tag.

[1] Adam M. Ścibior, Zoubin Ghahramani, and Andrew D. Gordon. 2015. Practical probabilistic programming with monads. In Proceedings of the 2015 ACM SIGPLAN Symposium on Haskell (Haskell ’15), Association for Computing Machinery, Vancouver, BC, Canada, 165–176.

[2] Adam M. Ścibior, Ohad Kammar, and Zoubin Ghahramani. 2018. Functional programming for modular Bayesian inference. In Proceedings of the ACM on Programming Languages Volume 2, ICFP (July 2018), 83:1–83:29.

[3] Adam M. Ścibior. 2019. Formally justified and modular Bayesian inference for probabilistic programs. Thesis. University of Cambridge.

Hacking

1. Install stack by following these instructions.

2. Clone the repository using one of these URLs:

   git clone [email protected]:tweag/monad-bayes.git
   git clone https://github.com/tweag/monad-bayes.git
   

Now you can use stack build, stack test and stack ghci.

To view the notebooks, go to the website. To use the notebooks interactively:

1. Compile the source: stack build
2. If you do not have nix install it.
3. Run nix develop --system x86_64-darwin --extra-experimental-features nix-command --extra-experimental-features flakes - this should open a nix shell. For Linux use x86_64-linux for --system option instead.
4. Run jupyter-lab from the nix shell to load the notebooks.

Your mileage may vary.

monad-bayes's Issues

*** Exception: System.Random.MWC.Distributions.categorical: bad weights!

{-# OPTIONS_GHC -Wall     #-}

{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE BangPatterns #-}

module Main where

import Control.Monad.Bayes.LogDomain
import Control.Monad.Bayes.Primitive
import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Population
import Control.Monad.Bayes.Conditional
import Control.Monad.Bayes.Inference

import Numeric.GSL.ODE
import Numeric.LinearAlgebra hiding ( step, Vector )

import qualified Data.Vector as V
import Data.Vector ( Vector )
import Control.Monad

import Control.Monad.Bayes.Sampler

import Debug.Trace

-- number of particles used in PMMH
n_particles :: Int
n_particles = 5

-- model data and constants
h :: Double
h = 0.1

k1', b', d', k2', c' :: Double
k1' = 2.0e2  -- Hare carrying capacity
b'  = 2.0e-2 -- Hare death rate per lynx
d'  = 4.0e-1 -- Lynx death rate
k2' = 2.0e1  -- Lynx carrying capacity
c'  = 4.0e-3 -- Lynx birth rate per hare

-- a :: Double
-- a = 0.5

-- type of simulation state
data S = S {p :: Double, z :: Double, log_alpha :: Double}

parameters :: (MonadDist m, CustomReal m ~ Double) => m (Double, Double)
parameters = do
  mu <- uniform 0 1
  sigma <- uniform 0 5
  return (mu,sigma)

-- initial state of the simulation
initial_state :: (MonadBayes m, CustomReal m ~ Double) => (Double, Double) -> m S
initial_state (mu, sigma) = do
  log_p_init <- normal (log 100) 0.2
  log_z_init <- normal (log  50) 0.1
  log_alpha_init <- normal (log mu) sigma
  return (S (exp log_p_init) (exp log_z_init) log_alpha_init)

-- transition model
transition :: (MonadBayes m, CustomReal m ~ Double) => (Double, Double) -> S -> m S
transition params state = do
  w <- normal 0 (sqrt h)
  -- TODO: ODE solver updates state here
  let a = exp (log_alpha state)
      m = solPp a (p state) (z state)
      newP = m ! 1 ! 0
      newZ = m ! 1 ! 1
      sigma = snd params
      newLog_alpha = -sigma * sigma * h / 2 - sigma * w
  return $ S newP newZ newLog_alpha
    where
      ppOde a k1 b d k2 c _t [pp, zz] =
        [
          a * pp * (1 - pp / k1) - b * pp * zz
        , -d * zz * (1 + zz / k2) + c * pp * zz
        ]
      ppOde _a _k1 _b _d _k2 _c _t vars =
        error $ "ppOde called with: " ++ show (length vars) ++ " variables"

      solPp a x y = odeSolve (ppOde a k1' b' d' k2' c')
                             [x, y]
                             (fromList [0.0, h])

-- full simulation given parameters, returns the parameters
model :: (MonadBayes m, CustomReal m ~ Double) => (Double, Double) -> m (Double, Double)
model params = foldl step (initial_state params) obs >> return params where
  step old_state p_obs = do
    state <- old_state
    new_state <- transition params state
    observe (Continuous (Normal (log (p state)) 0.1)) (log p_obs)
    return new_state

-- full model with particle filter, returns posterior over model parameters
full_model :: (MonadBayes m, CustomReal m ~ Double) => m (Double, Double)
full_model = parameters >>= (collapse . smc (length obs) n_particles . model)

-- PMMH transition kernel
-- monad-bayes does not currently have truncated normals
pmmh_kernel :: (MonadDist m, CustomReal m ~ Double) => [Double] -> m [Double]
pmmh_kernel [mu, sigma] = do
  mu' <- normal mu 1
  sigma' <- normal sigma 1
  return [mu', sigma']
pmmh_kernel xs = error $ "pmmh_kernel called with: " ++ show (length xs) ++ " variables"

-- full PMMH transition step
pmmh_step :: (MonadDist m, CustomReal m ~ Double) => [Double] -> m [Double]
pmmh_step params = do
  params' <- trace ("Params: " ++ show params) $ pmmh_kernel params
  let kernel_density  = unsafeContJointDensity (pmmh_kernel params) params'
  let kernel_density' = unsafeContJointDensity (pmmh_kernel params') params
  pm_density  <- pseudoDensity full_model (map Just params , [])
  pm_density' <- pseudoDensity full_model (map Just params', [])
  let mh_ratio = pm_density' * kernel_density' / (pm_density * kernel_density)
  accept <- bernoulli (min 1 (fromLogDomain mh_ratio))
  trace ("Accept: " ++ show accept) $ return (if accept then params' else params)

iterateNM :: Monad m => Int -> (a -> m a) -> a -> m (Vector a)
iterateNM n f x
  | n == 0 = return $ V.singleton x
  | otherwise = do
      y <- f x
      liftM (y `V.cons`) (iterateNM (n - 1) f y)

main :: IO ()
main = do
  ps1 <- sampleIOfixed $ do ps <- parameters
                            foo <- iterateNM 2 pmmh_step [fst ps, snd ps]
                            return foo
  let mus = V.map (!!0) $ V.drop 1 ps1
  putStrLn $ show (sum mus / 1)
  let sigmas = V.map (!!1) $ V.drop 1 ps1
  putStrLn $ show (sum sigmas / 1)

obs :: [Double]
obs = take 3 $ [
  76.2403679754159,
  65.0098784873532,
  72.7307834964011,
  66.9032651694089,
  69.4511465390588,
  64.2803362618095,
  56.2848299444718,
  55.38783776159,
  58.4038970077208,
  65.1111562800144,
  61.8487358181798,
  54.6268190304816,
  61.7832255452353,
  57.2580992086863,
  51.6506692415061,
  53.513024643873,
  53.2516281763045,
  52.8808601383571,
  50.8172783401389,
  56.4015933848598,
  51.9694880143296,
  53.6405052882993,
  54.575606130029,
  59.3492395433371,
  60.0888983492202,
  44.4110818026195,
  54.4982828776438,
  73.4084712052465,
  54.7644338856397,
  65.9215790688941,
  64.4989883914755,
  80.9772033868803,
  74.4092779374009,
  69.1141635364459,
  67.4281988428452,
  73.2048837536519,
  68.120224957282,
  62.0069646761111,
  80.1892530691497,
  61.1742005270423,
  79.8803349040011,
  83.2200370887438,
  77.5918040799712,
  72.5731170839407,
  77.4547803718811,
  76.6217772509249,
  73.676388095968,
  85.6784703161388,
  104.170224855567,
  89.2220461329041,
  97.8557959956232,
  108.537541058624,
  80.7233691710187,
  101.184198386732,
  106.119103176147,
  100.475879494068,
  125.13032727884,
  91.0701758431548,
  88.8812050725806,
  100.679238971492,
  110.589530407881,
  89.6049957679831,
  103.926588273376,
  100.068112237358,
  115.512346913051,
  113.109941769312,
  114.117083726776,
  126.223426393656,
  116.655177064036,
  119.47938528113,
  119.673062865914,
  122.17342562363,
  115.260394834536,
  118.68195366763,
  124.943355707336,
  151.936975691988,
  114.130007656096,
  148.946645952698,
  144.565611912741,
  118.759397178281,
  126.673416199079,
  142.600114603459,
  124.392431912102,
  123.845615416597,
  131.162727200371,
  127.530085947053,
  149.07114126,
  123.535118463953,
  166.415187417087,
  131.719924167561,
  139.991005535832,
  130.070166508874,
  141.353801701055,
  136.634266000686,
  120.874938386323,
  135.843010192177,
  158.44951183627,
  150.728722863145,
  176.881378505555,
  133.774741754552,
  137.445456445669,
  124.34236836753,
  152.292445573619,
  145.969440861177,
  131.636047513673,
  146.963132155529,
  151.695165885035,
  165.977302905119,
  142.031450539543,
  164.846194851848,
  168.013420195697,
  152.767971695764,
  143.851647250068,
  160.436734291349,
  152.302294176293,
  155.067612585515,
  175.247982500084,
  136.897735270247,
  151.97211044291,
  160.887803973557,
  136.486833382118,
  158.883061549789,
  152.458944513575,
  149.023829176969,
  155.340944696967,
  157.829634728021,
  158.633421934879,
  132.036887061947,
  135.101844477751,
  153.29922033573,
  158.052850728581,
  143.369824868863,
  157.714843025493,
  144.325777788735,
  131.0549741648,
  144.337275322618,
  145.026723249829,
  160.030142868999,
  159.542825278056,
  147.87991947625,
  126.381489671728,
  145.226447360054,
  136.894318063363,
  159.268060107793,
  119.239998434002,
  132.788094593534,
  150.663078585921,
  144.430826746593,
  165.173155950077,
  147.929951578988,
  157.760394555343,
  131.053973697495,
  165.580489628044,
  146.718090811465,
  142.773179403759,
  146.392167581097,
  146.025983883888,
  141.47610191321,
  134.214216068753,
  119.248789990169,
  135.031489547166,
  120.561743419339,
  167.557212051444,
  148.168532426848,
  136.547249833344,
  152.416629324237,
  173.5981255812,
  132.895016695342,
  138.669418114316,
  167.112761599289,
  119.007943156896,
  144.496358594722,
  138.515189246643,
  158.554459930596,
  136.975247626486,
  129.513055665373,
  169.020473191163,
  156.116991826441,
  129.524981775373,
  148.85363969582,
  157.801922552141,
  157.577368682676,
  140.430975061626,
  136.205715759791,
  158.040217033728,
  125.493638253065,
  130.114292183725,
  142.750827896278,
  174.071675002185,
  155.168991141177,
  165.5803346601,
  133.967097475562,
  133.930047497004,
  127.601345216493,
  135.602745924639,
  163.280148579137,
  150.634881342163,
  144.898580660782,
  123.221738509521,
  170.394076047434,
  152.596479493751,
  156.695148025635,
  125.317613954204,
  152.957730717693,
  168.012683439549,
  120.716643641628,
  150.531587838906,
  145.104736014324,
  124.182149449396,
  133.269771329542,
  114.976356844067,
  139.860395302311,
  137.801555771354,
  139.678968034677,
  134.434880549483,
  138.011761267802,
  125.439482936179,
  151.74813237231,
  131.847695399538,
  123.319231908423,
  152.547130543284,
  125.491719298048,
  137.096981763213,
  138.768540718737,
  181.5283483391,
  126.098802734704,
  163.482620601302,
  128.342617368436,
  141.586517467223,
  119.474015267941,
  120.9364132138,
  112.984981940618,
  141.341437153727,
  135.160684269048,
  134.509200686011,
  163.319201985027,
  135.415854723738,
  159.063490681458,
  153.148663733748,
  149.026302196215,
  164.826019093328,
  123.071508271513,
  176.387627783159,
  149.483954577518,
  123.02095153944,
  156.650270662224,
  181.195533023262,
  141.238386891745,
  134.701979169655,
  131.340667101923,
  144.548953780858,
  131.588662675565,
  144.541170469088,
  125.578498701588,
  164.74439785325,
  127.766537075735,
  133.273837701266,
  145.93004226061,
  123.374900614912,
  151.596803755421,
  144.894219796883,
  126.998541138479,
  147.035896225008,
  136.476333130728,
  135.608414577445,
  131.807584253206,
  127.366264864742,
  134.964091100405,
  143.421168142046,
  142.054495699452,
  151.438892001945,
  143.05917211411,
  136.712518442789,
  120.003089311654,
  140.869561010692,
  135.058679734824,
  134.081730533486,
  139.43057733862,
  162.26933796043,
  136.471725650913,
  125.873826665898,
  151.097217418369,
  121.162250703819,
  150.767682408018,
  130.606266701801,
  105.20437547654,
  132.359328972528,
  116.989153980971,
  148.954925443656,
  142.105725642533,
  143.553119435633,
  127.751657026012,
  147.148424004414,
  114.85477950561,
  139.694548684757,
  146.171259764561,
  127.789841217549,
  163.168916396005,
  147.231424950466,
  125.058830151811,
  119.074143576494,
  133.536702795857,
  137.262453768502,
  179.672247639064,
  118.207552485961,
  136.618926147412,
  142.352972731303,
  136.479602357591,
  118.750751452668,
  147.986276828413,
  147.411481380629,
  144.052405584174,
  146.153708105283,
  131.557481669184,
  134.453090884824,
  153.443307349816,
  146.984446479185,
  138.421394583867,
  132.707447826652,
  136.516554357634,
  122.191521371155,
  134.80783216139,
  155.447901778778,
  136.926153057977,
  125.184160305696,
  122.367097582219,
  134.726963096563,
  138.72230909606,
  133.531816345683,
  138.817134160604,
  151.680755115041,
  143.081423812855,
  144.507716078662,
  150.54905673819,
  164.395431529991,
  136.035307291772,
  135.963680724784,
  159.088071380666,
  115.544651209282,
  126.449866575417,
  145.861556307966,
  121.299809587301,
  135.041076468273,
  115.420027815205,
  120.649592931961,
  149.222665704902,
  136.787170312266,
  117.577955349667,
  123.09131466115,
  149.613666057669,
  127.542230989168,
  148.826516120672,
  172.024769211347,
  119.163511852413,
  168.846806388603,
  122.577905606408,
  163.691820942384,
  121.152394606119,
  125.364015734805,
  134.053429227603,
  116.102707413675,
  115.932034184824,
  157.856002757677,
  120.417539311501,
  159.448081983552,
  146.523081800403,
  142.244157886188,
  147.49875516021,
  149.794865803213,
  142.424753997512,
  155.377391656754,
  133.457569330548,
  157.979204500375,
  132.867976837153,
  150.766104208022,
  151.574164752913,
  142.79494077867,
  135.300935952747,
  155.087644429195,
  134.339799297657,
  117.780905975319,
  141.325010750032,
  137.141299531074,
  150.973508888167,
  135.037762357837,
  145.031582826828,
  138.551290507779,
  144.969403511833,
  162.819743874685,
  132.441267184392,
  131.380740568251,
  141.527642238539,
  124.238855213983,
  155.588321253241,
  145.068004565977,
  120.849885116319,
  138.600747823775,
  138.767637532808,
  149.658261434865,
  153.218173715505,
  123.173078726806,
  144.980628125895,
  133.601607881188,
  143.600386923789,
  134.99375013592,
  178.108618865062,
  162.589164120014,
  149.801157966121,
  124.574389244282,
  159.749711094236,
  154.462697945675,
  124.745887246537,
  134.619145288027,
  132.829029331846,
  133.515016847945,
  130.370138175619,
  134.058600191345,
  138.967206000683,
  144.844845114948,
  140.150588785047,
  168.155026956341,
  105.663614240897,
  188.849879696026,
  134.066594679321,
  118.875307519085,
  152.36326566797,
  126.208445705652,
  168.297022523569,
  152.27258877451,
  140.382228072818,
  119.424396673087,
  161.465352774379,
  151.655127807313,
  135.515116635163,
  151.532357393466,
  151.810511351233,
  136.361410128331,
  153.741012514592,
  147.431121676277,
  157.91915964734,
  120.327187813985,
  120.171500909668,
  142.428896431138,
  170.204819149259,
  158.214506571837,
  138.733887337899,
  123.808097192822,
  152.903612727097,
  140.769759838882,
  119.667303148295,
  129.45932949513,
  149.16660912873,
  150.254944400129,
  128.273570440445,
  153.278564595965,
  121.864809541511,
  161.542701435884,
  156.946807512148,
  158.881028441797,
  125.569197718704,
  137.940620707105,
  129.90603859364,
  152.757290961669,
  127.74685675702,
  137.183947553826,
  150.462794723172,
  140.929843406962,
  164.975435125998,
  182.344447097492,
  156.919848441299,
  150.750704137777,
  141.693725748231,
  110.676915031473,
  173.230786634625,
  148.248841760484,
  119.716751121111,
  143.947798672895,
  130.440507876832,
  131.964077776465,
  134.646898631033,
  100.778877180885,
  124.348744604234,
  127.180330688462,
  158.146060814462,
  171.358807156152,
  131.864294544909,
  143.085939246005,
  149.586050733912,
  132.2767374166,
  129.632802049237,
  137.02205312678,
  146.243848392623,
  125.713745961482,
  132.545872818923,
  129.355428299295,
  150.515004650034,
  149.481788573074,
  154.478850748833,
  137.329325663629,
  127.9504469047,
  133.350229974433]

Add instructions to use BuildKite's cache

... to speed up the initial Nix build.

Follow-up to #63
Related: https://nixos.org/nix/manual/#conf-substituters

User's guide

We need a user's guide to show people how to use the library. I like the style used by Hspec.

Tentative table of contents:

1. MonadDist and MonadBayes type classes
2. How to write models
3. Existing inference algorithms
4. Inference by composition of monad transformers
5. How to write a new inference algorithm

Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) is a widely used, gradient based, MCMC algorithm, that is the backbone of Stan's inference. I plan to implement it for monad-bayes. Todos (checkboxes indicate things done on branch):

• a means of converting a probabilistic program into a function f :: Floating n => [n] -> n (Or: can I use Numeric.AD.Double?)
• refactor of monad-bayes to move from the concrete type Double to any number, which is needed for automatic differentiation of f
• implementation of Hamiltonian dynamics with a numerical integrator
• Thorough testing

Optional:

• Non-parametric HMC: this allows HMC to work for arbitrary probabilistic programs, not just ones with a fixed number of variables
• NUTS dynamics
• Riemannian HMC

I have the basic pieces working on branch, but the checkboxes indicate further steps. Current roadblocks:

• many distributions are defined using inverse cdfs from the statistics package, which use concrete type Double. I need versions that are in terms of arbitrary Num n => n, for autodiff to work.
• Non-parametric HMC requires some further understanding on my part, in particular if I want to do something beautifully lazy

Classical MCMC with Metropolis-Hasting proposals

Hello,

I fail to implement a classical MCMC method with Metropolis-Hasting (MH) proposals. How can I define the proposal distribution conditioned on the current state? I looked at the code of the MH proposal, and it is hidden in the Trace data type. The acceptance ratio is computed using the likelihood ratio of the proposed vs. current state (and some correction factor which I do not understand; it involves the number of sampled variables?).

To me, this is a standard Metropolis update, and not a Metropolis-Hastings update which includes a proposal distribution dependent on the current state. Can you comment on this?

Thank you for your help!
Dominik

As a side question: Is it possible to run an MCMC with MH proposals without using Trace. I.e., by just saving or outputting the current state of the chain?

Fixing the badges

The badges on the front page say things like the build is failing, which makes the repo look broken (which it isn't).

Remove T from transformer names

We use the T postfix to indicate monad transformers used for inference, for example WeightedT. This is unnecessary, since we never apply those to Identity, so the names without T are available. We should rename WeightedT to Weighted and similarly for other transformers.

Remove weights from Empirical

Currently Empirical is a synonym for Weighted ListT, which means it's a collection of weighted samples. This is a bit misleading, since normally empirical distribution is an unweighted collection of samples. We should change Empirical to be a synonym for ListT, then the current implementation can be recovered as Weighted Empirical, which is both more modular and more intuitive.

Abstract LogFloat

We are currently using LogFloat to represent probabilities. This is a sensible default, but may not always be the best choice. It would be better to parametrise the type classes MonadDist and MonadBayes by the type representing probability.

sampling from distributions over lazy lists?

For educational purposes, I modified your HMM example to work as an even simpler Markov chain. I decided to really live the Haskell lifestyle and try to define a distribution over lazy infinite lists.

states :: [Int]
states = [-1,0,1]

trans :: MonadDist m => Int -> m Int
trans (-1) = categorical $ zip states [0.1, 0.4, 0.5]
trans 0 = categorical $ zip states [0.2, 0.6, 0.2]
trans 1 = categorical $ zip states [0.15, 0.7, 0.15]

step startDist = startDist >>= trans

start :: MonadDist m => m Int
start = uniformD states

mmlazy :: MonadDist m => m Int -> m [Int]
mmlazy startDist = liftM2 (:) startDist (mmlazy (step startDist) )

Ideally, I'd be able to, say sampleIO $ fmap head (mmlazy start) and get a sample of just the first state; likewise use take, etc.. In practice, although this typechecks, the computation just hangs.

Perhaps there are good reasons why I shouldn't even be trying to do this -- in which case, is there any piece of the code, paper, etc. that I might want to be looking at to gain a better understanding of why it won't work?

Applicative Do for improving efficiency

GHC 8 introduces the Applicative Do extension, which can be used to desugar do notation to Applicative operations where possible. The paper introducing it actually mentions specifically that monad-bayes could benefit from this extension. We should therefore modify the existing code to use Applicative instead of Monad interface wherever possible.

Incorporating failure into enumerate leads to exception

I wrote this code:

b :: MaybeT Enumerator Bool
b = do
  x <- uniformD [1,2]
  return (x > 1)

out = enumerate $ runMaybeT b

Running out causes:

*** Exception: Infinitely supported random variables not supported in Enumerator

But the RV would appear to be finitely supported, so I'm not sure what is going on

Abstract numerical types for distributions

The distributions we use are currently fixed to use Int and Double. We should consider extending it to more abstract types, provided that it wouldn't cause too many problems.

Build benchmarks in CI

Benchmarks are not currently built in CI, but they should be.

ssm-bench broke due to recent refactorings. This breakage was not caught in CI:

$ stack --stack-yaml stack-ghc881.yaml build --bench --no-run-benchmarks
[...]
Building benchmark 'ssm-bench' for monad-bayes-0.1.0.0..
[1 of 2] Compiling NonlinearSSM
            
/home/leo/Code/monad-bayes/models/NonlinearSSM.hs:16:36: error:
    • Variable not in scope: sq :: Double -> Double
    • Perhaps you meant ‘seq’ (imported from Prelude)
   |        
16 | mean x n = 0.5 * x + 25 * x / (1 + sq x) + 8 * cos (1.2 * fromIntegral n)
   |                                    ^^
            
            
--  While building package monad-bayes-0.1.0.0 using:
      /home/leo/.stack/setup-exe-cache/x86_64-linux-nix/Cabal-simple_mPHDZzAJ_3.0.0.0_ghc-8.8.1 --builddir=.stack-work/dist/x86_64-linux-nix/Cabal-3.0.0.0 build lib:monad-bayes exe:example bench:speed-bench bench:ssm-bench --ghc-options " -fdiagnostics-color=always"
    Process exited with code: ExitFailure 1

Replace categorical with discrete in Primitive

Currently one constructor for Primitive is Categorical :: [(a,LogFloat)] -> Primitive a with arbitrary type a. This is a little inconvenient, and forces us to make certain restrictions on a which do not make sense in all contexts. It would make more sense to replace Categorical with Discrete :: [LogFloat] -> Primitive Int and then provide default implementation for categorical in terms of discrete.

Publish to Hackage/Stackage

Would it be possible to release on Hackage/Stackage ?
It would be way more convenient.

Thanks

Trace MCMC with proposals

To make trace MCMC actually useful for real-world problems, you really want to be able to design your own traces in a natural way. This is currently hard, because a trace is a list of Doubles, which is not a very interpretable thing. Gen (another probabilistic programming language) has a nice way of doing things, where you name variables in your trace and then pass a "trace map" to your MCMC algorithm, which tells it how to update those named variables.

Using MonadState (and a slightly different base functor from the free monad trace representation that records the random choice), I've experimented with getting a similar capability here. The user experience looks a lot like Gen's, namely:

example :: (MonadInfer m, MonadState String m) => m Bool
example = do
  x <- traced "x" (bernoulli 0.4)
  y <- if x then bernoulli 0.2 else bernoulli 0.9
  condition y
  return x

Because of the nice typeclassy nature of monad-bayes, you just add a constraint to your type MonadState String m and that gives you the ability to name variables.

It seems to work, more or less, but before I add it to the code, I would like to make sure it's fast (it probably isn't). There's also the issue that because Haskell is strongly typed, it's not easy to have a single data structure which records the values of arbitrary random variables of arbitrary types. I'm not going to worry about that for now: getting it to work for bernoulli and uniform is good enough for me.

This paper is also almost certainly relevant to doing things better: https://dl.acm.org/doi/10.1145/3371087

Investigate message passing algorithms

monad-bayes doesn't currently implement anything like belief propagation. I saw https://arxiv.org/pdf/2112.13251.pdf recently, and since Haskell has some nice avenues for functional reactive programming, I would at least like to spend some time thinking about what a message passing implementation would look like.

Alternative and MonadPlus instance for Dist

If you have an Alternative instance for Dist your condition function would be guard.

instance Alternative (Dist r) where
  empty = Dist []
  Dist d1 <|> Dist d2 = Dist (d1 ++ d2)

Other distribution types might also allow an Alternative instance.

Build with "-Werror" in CI

#63 introduced the dev flag for use in development and CI. Currently there are deprecation warnings that would make the build fail if we added -Werror. So let's fix those warnings and then set -Werror to make sure we don't regress.

• Get rid of all warnings
• Add the -Werror GHC flag if dev is set

Automatic Differentiation for density gradient

At the moment we can compute PDF for Primitive distributions, but not its gradient. Backward mode AD would allow us to efficiently compute gradients for advanced inference algorithms.

This is related to #14.

Testing

Our current approach to testing is rather ad-hoc. A more systematic one would be helpful, ideally including property-based testing. Then there's also the issue of testing stochastic programs. Any suggestions welcome.

Probabilities larger than 1.0?

Hello!

I am trying to use your library but fail at a very simple example test. I would like to estimate the mean of a Gaussian from some sample data.

I create the data using this function:

generateDataGaussian :: MonadSample m
      => Int  -- ^ Number of points
      -> m [Double] -- ^ List of latent and observable states from t=1
generateDataGaussian n = replicateM n $ normal 0 0.1

My (very simple) model with a uniform prior for the mean from -0.5 to 0.5 is the
following:

gaussian :: MonadInfer m => [Double] -> m (Double, Double)
gaussian os = do
  c <- uniform (-0.5) 0.5
  mapM_ (score . normalPdf c 0.1) os
  return (c, s)

The main function is:

main :: IO ()
main = sampleIOfixed $ do
  os <- generateDataGaussian 50
  res <- runPopulation $ smcMultinomial 0 10 (gaussian ys)
  liftIO $ putStr $ unlines $ map show res

So, basically I only sample from the prior, without any particle filter step. Then, the result I get is:

((-0.48477485643126317,0.1),4.4371710323320446e-225)
((-0.1477194307293065,0.1),2.08421997614896e-3)
((-0.14943329883414702,0.1),6.543157297094635e-4)
((-0.29494328640792944,0.1),6.883965851944845e-70)
((-9.60734486685807e-2,0.1),3.083465431826563e9)
((-0.3536764046071378,0.1),1.5220519619418593e-109)
((0.4094283707386025,0.1),6.635827036999188e-178)
((-0.16277680281860063,0.1),4.7892419999331464e-8)
((3.550408159678675e-2,0.1),2.0870669318446438e14)
((-0.28782974706078146,0.1),1.3557955847978798e-65)

This means, that two draws have very high values as probabilities (e14 and e9). Why is that?

Furthermore, if I do one resample step in the particle filter, the value with the highest probability vanishes?? Why would that be the case? I.e.,

main :: IO ()
main = sampleIOfixed $ do
  os <- generateDataGaussian 50
  res <- runPopulation $ smcMultinomial 1 10 (gaussian ys)
  liftIO $ putStr $ unlines $ map show res

Has a result of:

((-0.16277680281860063,0.1),2.5237137500055417e-8)
((-9.60734486685807e-2,0.1),2.5386626855022683e9)
((-0.16277680281860063,0.1),2.5237137500055417e-8)
((-0.16277680281860063,0.1),2.5237137500055417e-8)
((-0.14943329883414702,0.1),3.638030875815378e-4)
((-0.1477194307293065,0.1),1.1683592101578812e-3)
((-0.28782974706078146,0.1),1.0265137549975822e-65)
((-0.16277680281860063,0.1),2.5237137500055417e-8)
((-0.28782974706078146,0.1),1.0265137549975822e-65)
((-0.16277680281860063,0.1),2.5237137500055417e-8)

Maybe I am doing something fundamentally wrong, please let me know if this is
the case. Also, if I increase the number of points, particles or particle filter
steps, I do not get conclusive results. In general, a manual with examples would
be of great help!

Updated docs

• documentation (of which a draft is https://monad-bayes.netlify.app/)
• a tutorial, in notebook form, covering independent sampling, MCMC, particle filters, rmsmc, pmmh, smc2 (https://probmods.org/ is an excellent example of how to do this well (in that case the focus is on using probabilistic programming to model cognition) )
• update the haddock docs with more comments (e.g. difference between different resamplers)
• update hackage to link to the docs and tutorials

The broad goal is to make monad-bayes accessible to anyone with a beginner understanding of Haskell and a beginner understanding of Bayesian inference. That will require clear tutorials with simple working examples.

Confusion over the use of pmmh

I am trying to implement a small discrete time SIR model as an HMM and fit this using pmmh.
I am not sure if I am using monad-bayes correctly and, if I am, I am unsure as to why the performance of pmmh on my model is so poor.

The type signature for pmmh returns as output m [[(a, Log Double)]] where a is the output of the model. My normal understanding of PMMH is that it returns samples from the posterior but, when running the implemented version, we receive lists of samples for each MH iteration, one for each particle. Upon inspection, each of these values appear to be the same which follows my understanding of PMMH. Is the information from each particle returned in order to have access to the weights of each particle or do I have a more fundamental misunderstanding?

Furthermore, my interpretation of how I am meant to use the PMMH function to fit parameter is that I am to pass a model that scores against the given data whilst evolving the underlying model, returning the parameters at the end that are reweighted by the scoring procedure. The following model does that, folding over the data, scoring each individually.

scoreEpidemicToData 
    :: MonadInfer m 
    => FixedParams 
    -> Epidemic 
    ->  LatentState 
    -> Params 
    -> m Params 
scoreEpidemicToData fixedParams ys initialState params  = do
    foldM_ (scoreEpidemicToDatum fixedParams params) initialState (unwrapEpidemic ys)
    return params

where scoreEpidemicToDatum takes the following form

scoreEpidemicToDatum 
    :: (MonadSample m, MonadCond m)  
    => FixedParams 
    -> Params 
    -> LatentState 
    -> Int 
    ->  m LatentState 
scoreEpidemicToDatum fixedParams params x datum = do
    let obs lambda y = score (poissonPdf lambda y)
    x' <- transitionModel fixedParams params x 
    obs ((fromIntegral $ inf x') * (rho params)) datum
    return x'

I am performing inference by using pmmh as follows

sampleIO $ do
    pmmhRes <- prior $ pmmh nsteps 14 nparticles paramsPrior (scoreEpidemicToData fixedParams dat initialState)

where my data is [Double] is length 14. Is this the correct way to use pmmh to perform inference over parameters? If so, is there any way to increase the performance of pmmh as it is not feasible to do inference on what is quite a simple model (sampling from binomial and poisson distributions) with enough particles or steps to have a good sense of if the samples converge?

The full code can be found here.

Example Gallery

Monad-bayes currently has examples, in the models folder, but a large subset don't compile with the current code version. And for the ones that do, it would be nice if it was much clearer how to run them. Ideally we'd have runnable examples that produced pretty visualizations for a core set of models.

There's also benchmarking code, and it would be nice if it came with some discussion of the various problems being modelled and the benefits of the different inference algorithms. The goal is that for each inference method that monad-bayes provides, we could see 1) how to use it and 2) when to use it.

New examples to add, which might go in the repo or in a separate monad-bayes-app repo to avoid incurring large dependencies:

• an implementation of odeModelExpanded.pdf
• An Ising model (I've done this on branch)
• A simulation of a simple physical system using the Hamilton library (fun with Monte Carlo)
• Online inference of a mixture model: i.e. as observations come in, you can see the posterior update in real time
• Generation of sentences from a PCFG with some scoring (also done on branch)
• A probability distribution over JSONs, with use of lenses (on branch)
• A distribution over databases
• The Rational Speech Acts model (https://www.problang.org/chapters/01-introduction.html)
• Some sort of example of a decision making agent (maybe using FRP?)
• A combinator parser that incorporates probability in a non-trivial way (some progress here via megaparsec's ParsecT, on branch)

Towards monad-bayes 1.0

Notes from monad-bayes discussion

Date: [2021-05-17 Mon]

Those present: @idontgetoutmuch (Dominic Steinitz), @MMesch (Matthias Meschede)

monad-bayes contains two distinct elements:

1. A method of generating random values according to given probability distributions and then being able to combine these together via a monad. It does this by using a free monad (in fact the Church encoding of such for efficiency).
2. A way of defining statistical models using this framework and then defining inference methods by composition so e.g. Particle Marginal Metropolis Hastings (PMMH) is an almost magically simple composition of Metropolis Hastings using Sequential Monte Carlo as an approximation to the posterior (the composition turns out to actually sample from the posterior not just sample from an approximation).

It uses mwc-random to supply a) a source of randomness and b) distributions over that source.

random-fu also contains a method of generating random values according to given probability distributions. It also uses a free monad via the MonadPrompt (https://hackage.haskell.org/package/MonadPrompt) package - a precursor of free monad packages such as free (https://hackage.haskell.org/package/free).

So the long term plan is to create a dev branch for monad-bayes (and possibly random-fu) that have the following structure:

1. A source of randomness (I have already demonstrated we can replace mwc-random in monad-bayes by random - see #112 for the full details).
2. A way of generating random values from distributions. Maybe we can use the monad-bayes free monad to replace the MonadPrompt in random-fu? And then use random-fu as the way of generating random values from distributions. This would also allow the removal of monad-bayes dependency on the statistics package.
3. A way of composing inference methods and specifying models. I feel that monad-bayes should concentrate on this rather than being a "vertical integration" of related but disparate concepts.
4. In addition, we'd like to be able to change the MH proposal step (make it a parameter?) for inference. This definitely fits in the 3rd layer.
   In addition, we ought to have a model and tests that demonstrate more clearly how to use the package and show that we get convergence the more particles we use and the more MH steps we use.

Stack is not working on stack yaml

Describe the bug
Stack build is getting 404ed(like Test and GHCI stack)
To Reproduce
run Stack build in folder
Expected behavior
Stack works
Environment

• Ubuntu 18.04
• just git cloned

Additional context
james.hennessy@ASR-1335A-U:~/blog-resources/monad-bayes$ stack test Downloading lts-15.6 build plan ...RedownloadFailed Request { host = "raw.githubusercontent.com" port = 443 secure = True requestHeaders = [] path = "/fpco/lts-haskell/master//lts-15.6.yaml" queryString = "" method = "GET" proxy = Nothing rawBody = False redirectCount = 10 responseTimeout = ResponseTimeoutDefault requestVersion = HTTP/1.1 } "/home/AURIS.LOCAL/james.hennessy/.stack/build-plan/lts-15.6.yaml" (Response {responseStatus = Status {statusCode = 404, statusMessage = "Not Found"}, responseVersion = HTTP/1.1, responseHeaders = [("Connection","keep-alive"),("Content-Length","14"),("Content-Security-Policy","default-src 'none'; style-src 'unsafe-inline'; sandbox"),("Strict-Transport-Security","max-age=31536000"),("X-Content-Type-Options","nosniff"),("X-Frame-Options","deny"),("X-XSS-Protection","1; mode=block"),("Content-Type","text/plain; charset=utf-8"),("Via","1.1 varnish (Varnish/6.0), 1.1 varnish"),("X-GitHub-Request-Id","BA62:66B7:1B8749:207C20:5FD6B980"),("Accept-Ranges","bytes"),("Date","Mon, 14 Dec 2020 01:02:19 GMT"),("X-Served-By","cache-pao17445-PAO"),("X-Cache","HIT, HIT"),("X-Cache-Hits","3, 1"),("X-Timer","S1607907740.781343,VS0,VE1"),("Vary","Authorization,Accept-Encoding"),("Access-Control-Allow-Origin","*"),("X-Fastly-Request-ID","2f1eb56d61c037f83dcb38ab41ad9f96438eb481"),("Expires","Mon, 14 Dec 2020 01:07:19 GMT"),("Source-Age","22")], responseBody = (), responseCookieJar = CJ {expose = []}, responseClose' = ResponseClose})

Link to "Practical Probabilistic Programming with Monads"

The article on this library is "Practical Probabilistic Programming with Monads" by Ścibior, Ghahramani, and Gordon. Would a link in the readme be appropriate?

Failed to load interface for ‘Control.Monad.Bayes.Simple’

When I load models/Dice in ghci, I got the following error.

:l models/Dice
[1 of 2] Compiling Control.Monad.Bayes.Class ( /home/vincent/worksapce/monad-bayes/src/Control/Monad/Bayes/Class.hs, interpreted )
[2 of 2] Compiling Dice             ( models/Dice.hs, interpreted )

models/Dice.hs:11:1: error:
    Failed to load interface for ‘Control.Monad.Bayes.Simple’
    Perhaps you meant
      Control.Monad.Bayes.Sampler (needs flag -package-key monad-bayes-0.1.0.0)
      Control.Monad.Bayes.Free (needs flag -package-key monad-bayes-0.1.0.0)
      Control.Monad.Bayes.Class (needs flag -package-key monad-bayes-0.1.0.0)
    Use -v to see a list of the files searched for.
Failed, modules loaded: Control.Monad.Bayes.Class.

how to solve it?

Thank a lot.

Examples for variational/hmc inference

This library is really awesome, and I've been using it to implement some linguistic models recently (along the lines of http://gscontras.github.io/ESSLLI-2016/chapters/1-introduction.html). I'm having trouble getting the fancier inference methods to work (currently I'm just using enumerate and working with discrete variables), so was wondering if there were any simple examples of variational or hmc inference that I could try out, to get the hang of it.

Warnings and error messages

A big part of making monad-bayes usable is to have warnings and error messages that guide users. Issue #121 is a good example: the MCMC chain was never accepting proposals, and it would have been much easier to debug if there had been a warning message like "MCMC chain fails n% of the time. This suggests that the proposal needs improvement."

Other warnings and errors:

• if initial sample in MCMC chain has zero probability, that should probably be an error, but at least a warning
• if MCMC chain doesn't seem to be converging, that should be a warning
• if SMC is used in a model with only a single factor statement, that should raise a warning like: "Are you sure you want to use SMC here? This model only has a single factor, so won't benefit from it."

A somewhat related issue that I'll lump into this is to add some reporting to an MCMC run, like percentage acceptance, time taken, and maybe even a progress bar. (Note: the natural transformation in sis could update a progress bar at each step of SMC)

Automatic MonadDist derivation

We often provide MonadDist instances by lifting distributions from the transformed monad. Currently this is only done for primitive, which means that all the specialised implementations of distributions used by the transformed monad are lost. We should instead lift all of them from the transformed monad, bearing in mind that more distributions can be added to MonadDist in the future.

I'm not sure what's the best way to do it. Metaprogramming and a compiler extension for automatic derivation come to mind.

Weird sampling outputs on a toy example

To try out the library, i implemented the Monty Hall problem. My intention is to infer the possible states of the unobserved rooms, given some observations, e.g., given that the guest chooses room A and Monty chooses room B, what are the probabilities for the prize location. Using monad-bayes (maybe naively) for this somehow gives wrong results.

Here's the code:

module Monty where

import Control.Monad.Bayes.Class (MonadInfer, uniformD, condition)
import Control.Monad (forM_)
import Control.Monad.Bayes.Sampler (sampleIOfixed, sampleIO)
import Control.Monad.Bayes.Weighted (prior)
import Control.Monad.Bayes.Traced (mh)

data Room = A | B | C
  deriving (Eq, Show)
data Observation = Guest Room | Prize Room | Monty Room

condMonty :: MonadInfer m => Room -> Room -> m Room
condMonty guest prize = let options = [A, B, C]
                            validOptions = filter (\x -> x /= guest && x /= prize) options
                        in uniformD validOptions

data GameState = GameState
  { s_guest :: Room
  , s_prize :: Room
  , s_monty :: Room
  }
stateToTuple :: GameState -> (Room, Room, Room)
stateToTuple s = (s_guest s, s_prize s, s_monty s)

montyModel :: MonadInfer m => [Observation] -> m GameState
montyModel observations = do
  guest <- uniformD [A, B, C]
  prize <- uniformD [A, B, C]
  monty <- condMonty guest prize

  forM_ observations (scoreObs guest prize monty)
  return $ GameState guest prize monty
  where
    scoreObs guest _ _ (Guest room) = condition $ guest == room
    scoreObs _ prize _ (Prize room) = condition $ prize == room
    scoreObs _ _ monty (Monty room) = condition $ monty == room

data RoomFreq = RoomFreq
  { freqA :: Double
  , freqB :: Double
  , freqC :: Double
  }
  deriving Show

getRates :: [Room] -> RoomFreq
getRates rs =
  let
    (a, b, c) = foldl addOne (0, 0, 0) rs
    len = fromIntegral $ length rs
  in RoomFreq (a / len) (b / len) (c / len)
  where
    addOne (a, b, c) A = (a + 1, b, c)
    addOne (a, b, c) B = (a, b + 1, c)
    addOne (a, b, c) C = (a, b, c + 1)


runInference :: IO ()
runInference = do
  let nsamples = 5000
      burnin = 500

  -- Given the guest's choice and Monty's choice, i want to infer where the prize is located
  s <- take nsamples <$> (sampleIO $ prior $ mh (nsamples + burnin) $ montyModel [Guest A, Monty B])
  putStrLn "Last 20 samples:"
  print $ map stateToTuple $ take 20 s
  putStrLn ""

  let guestRates = getRates $ s_guest <$> s
      prizeRates = getRates $ s_prize <$> s
      montyRates = getRates $ s_monty <$> s
  putStrLn "Guest room choice frequency:" >> putStr "  " >> print guestRates
  putStrLn "Prize location frequency:" >> putStr "  " >> print prizeRates
  putStrLn "Monty room choice frequency:" >> putStr "  " >> print montyRates

And here are some typical outputs from running it three times:
Run 1

Last 20 samples:
[(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B),(A,C,B)]

Guest room choice frequency:
  RoomFreq {freqA = 1.0, freqB = 0.0, freqC = 0.0}
Prize location frequency:
  RoomFreq {freqA = 0.0, freqB = 0.0, freqC = 1.0}
Monty room choice frequency:
  RoomFreq {freqA = 0.0, freqB = 1.0, freqC = 0.0}

Run 2

Last 20 samples:
[(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C),(B,A,C)]

Guest room choice frequency:
  RoomFreq {freqA = 0.0, freqB = 1.0, freqC = 0.0}
Prize location frequency:
  RoomFreq {freqA = 1.0, freqB = 0.0, freqC = 0.0}
Monty room choice frequency:
  RoomFreq {freqA = 0.0, freqB = 0.0, freqC = 1.0}

Run 3

Last 20 samples:
[(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C),(B,B,C)]

Guest room choice frequency:
  RoomFreq {freqA = 0.0, freqB = 1.0, freqC = 0.0}
Prize location frequency:
  RoomFreq {freqA = 0.0, freqB = 1.0, freqC = 0.0}
Monty room choice frequency:
  RoomFreq {freqA = 0.0, freqB = 0.0, freqC = 1.0}

All of these are wrong and 2 and 3 should have been impossible, given the model specification, i think.

To Reproduce
Run the code above and observe the outputs.

Expected behavior
Given the inputs that are hardcoded in the code – guest chooses A and Monty chooses B –, i'd expect the following output:

Guest room choice frequency:
  RoomFreq {freqA = 1.0, freqB = 0.0, freqC = 0.0}
Prize location frequency:
  RoomFreq {freqA = 0.333, freqB = 0.0, freqC = 0.666}
Monty room choice frequency:
  RoomFreq {freqA = 0.0, freqB = 1.0, freqC = 0.0}

Environment

• OS name + version: Arch Linux + umm... 5.15.11 kernel
• Version of the code: 3f49a8b

Additional context
I don't know whether this is a valid approach to implementing the problem. Originally i wanted to see how easy it would be to convert this Bayesian network example to monad-bayes. My guess is that the MCMC sampler might be diverging or that condMonty shouldn't be used like that to define the conditional probability distribution, or both, but i've no idea how to debug this.

Error in benchmark

On running stack bench:

ssm-bench: Oops!  Entered absent arg $dMonad Monad (Population m)
"RM-SMC"
Benchmark ssm-bench: ERROR
Completed 2 action(s).

CI does not run tests

CI should run tests, but currently does not.

#89 commits a broken test, but the "build and test" step passes: https://buildkite.com/tweag-1/monad-bayes/builds/88.

doesn't install with ghc 9.0.1

Describe the bug
doesn't install with ghc 9.0.1

To Reproduce

 ghc --version
The Glorious Glasgow Haskell Compilation System, version 9.0.1
cabal install monad-bayes --lib
Resolving dependencies...
cabal: Could not resolve dependencies:
[__0] trying: monad-bayes-0.1.1.0 (user goal)
[__1] next goal: base (dependency of monad-bayes)
[__1] rejecting: base-4.15.0.0/installed-4.15.0.0 (conflict: monad-bayes =>
base>=4.11 && <4.14)
[__1] skipping: base-4.16.0.0, base-4.15.0.0, base-4.14.3.0, base-4.14.2.0,
base-4.14.1.0, base-4.14.0.0 (has the same characteristics that caused the
previous version to fail: excluded by constraint '>=4.11 && <4.14' from
'monad-bayes')
[__1] rejecting: base-4.13.0.0, base-4.12.0.0, base-4.11.1.0, base-4.11.0.0,
base-4.10.1.0, base-4.10.0.0, base-4.9.1.0, base-4.9.0.0, base-4.8.2.0,
base-4.8.1.0, base-4.8.0.0, base-4.7.0.2, base-4.7.0.1, base-4.7.0.0,
base-4.6.0.1, base-4.6.0.0, base-4.5.1.0, base-4.5.0.0, base-4.4.1.0,
base-4.4.0.0, base-4.3.1.0, base-4.3.0.0, base-4.2.0.2, base-4.2.0.1,
base-4.2.0.0, base-4.1.0.0, base-4.0.0.0, base-3.0.3.2, base-3.0.3.1
(constraint from non-upgradeable package requires installed instance)
[__1] fail (backjumping, conflict set: base, monad-bayes)
After searching the rest of the dependency tree exhaustively, these were the
goals I've had most trouble fulfilling: base, monad-bayes

Expected behavior
Expect it to install

Environment

• Version of the code: monad-bayes-0.1.1.0

Additional context
Installs with --allow-new

Monad-Coroutine break

Describe the bug
Installation of the package fails from hackage, so does building from source. I received an error during building that said building monad-coroutine-0.9.1 failed with an error about MonadFail not being defined.

To Reproduce
Either running cabal install monad-bayes or building the source code should reproduce.

Expected behavior
Package to install without a problem.

Environment

• OS name + version: Pop OS 20.04
• Version of the code: 0.1.1

Additional context
Fixing the monad-coroutine package to 0.9.0.4 solved the issue, I suspect the problem comes from version 0.9.1. I can open a pull request for this to temporarily fix the version, until monad-coroutine is fixed.

[feature request] Add progress bar capabilities?

Could we please add a monad transformer that adds a progress bar capability to an inference algorithm?

E.g., every time smc takes a step, we can decrease a counter.

Place Holder for the Fellowship Work

Is your feature request related to a problem? Please describe.
A clear and concise description of what the problem is. Ex. I'm always frustrated when [...]

Describe the solution you'd like
A clear and concise description of what you want to happen.

Describe alternatives you've considered
A clear and concise description of any alternative solutions or features you've considered.

Additional context
Add any other context or screenshots about the feature request here.

Variational Inference

Most of the more popular probabilistic programming languages have implementations of variational inference (VI). As such, it's absence in monad-bayes is something of an obstacle to real-world use.

A fast variational inference implementation needs (to my current understanding) a faster autodiff library than ad, which is the most currently accessible Haskell one. There seems to be some work in this direction (vis a vis Accelerate and https://github.com/VMatthijs/CHAD) but it's early days.

In the meantime, I'd like to implement a reference implementation using ad with no concern for speed at all. Previous versions of monad-bayes had something like this, so that will be my starting point. One of the conceptual challenges is that VI only makes sense (at least simple versions) for probabilistic programs with differentiable scoring functions, so it's still unclear to me how/if to express that in the types.

As I understand more about how the implementation should work, I'll update this issue.

Is it possible to implement Gaussian Processes based on monad-bayes?

Is it possible to implement Gaussian Processes based on monad-bayes?
If yes, would you mind giving me some hints how to implement Gaussian Processes based on monad-bayes?

Support for GHC 8

We depend on random-fu which currently doesn't build with GHC 8.

graphical rendering of Bayesian networks ?

Would it be possible to render the "reified" network in graphviz or similar? I suspect that the current, Monad-based implementation is a fundamental limitation for this, but I cannot articulate why.

Build with stack requires manual intervention

Building with stack requires one of the folllowing commands:

ln -s stack-ghc844.yaml stack.yaml
ln -s stack-ghc865.yaml stack.yaml
ln -s stack-ghc881.yaml stack.yaml

I suggest mentioning this in the Readme.

Continuous Integration

We should deploy a CI system and I'm leaning towards Travis. We should at least test GHC versions 7.10 and 8.0, building both with Stack and Cabal.

hident dependency is broken

Describe the bug

error: Package ‘hindent-5.3.1’ in /nix/store/16zk3bb2m3q2xlh1pnc5rhj1m5xvpq6q-nixpkgs-2009-src/pkgs/development/haskell-modules/hackage-packages.nix:125608 is marked as broken, refusing to evaluate.

       a) For `nixos-rebuild` you can set
         { nixpkgs.config.allowBroken = true; }
       in configuration.nix to override this.

       b) For `nix-env`, `nix-build`, `nix-shell` or any other Nix command you can add
         { allowBroken = true; }
       to ~/.config/nixpkgs/config.nix.

To Reproduce
nix-shell

Expected behavior
I should get a nix shell.

Environment

• OS name + version: macOS 11.5.1
• Version of the code:

Additional context
Add any other context about the problem here.

Any way for outsiders to help?

I've been following this library with interest for a while. Recently I've noticed that lots of changes are happening. Is there a big release coming up? Any help needed, e.g. around documentation, examples, or testing?

Benchmarks fail

Describe the bug

[nix-shell:~/monad-bayes]$ cabal bench speed-bench
cabal bench speed-bench
Build profile: -w ghc-8.10.4 -O1
In order, the following will be built (use -v for more details):
 - monad-bayes-0.1.1.0 (bench:speed-bench) (file benchmark/Speed.hs changed)
./monad-bayes.cabal has been changed. Re-configuring with most recently used
options. If this fails, please run configure manually.
Configuring benchmark 'speed-bench' for monad-bayes-0.1.1.0..
Warning: The package has an extraneous version range for a dependency on an
internal library: monad-bayes >=0 && ==0.1.1.0, monad-bayes >=0 && ==0.1.1.0,
monad-bayes >=0 && ==0.1.1.0, monad-bayes >=0 && ==0.1.1.0. This version range
includes the current package but isn't needed as the current package's library
will always be used.
Preprocessing benchmark 'speed-bench' for monad-bayes-0.1.1.0..
Building benchmark 'speed-bench' for monad-bayes-0.1.1.0..
[4 of 4] Compiling Main             ( benchmark/Speed.hs, /Users/dom/monad-bayes/dist-newstyle/build/x86_64-osx/ghc-8.10.4/monad-bayes-0.1.1.0/b/speed-bench/build/speed-bench/speed-bench-tmp/Main.o )
Linking /Users/dom/monad-bayes/dist-newstyle/build/x86_64-osx/ghc-8.10.4/monad-bayes-0.1.1.0/b/speed-bench/build/speed-bench/speed-bench ...
Running 1 benchmarks...
Benchmark speed-bench: RUNNING...
speed-bench: /bin/sh: createProcess: runInteractiveProcess: chdir: does not exist (No such file or directory)
Benchmark speed-bench: ERROR
cabal: Benchmarks failed for bench:speed-bench from monad-bayes-0.1.1.0.

To Reproduce
cabal bench speed-bench

Expected behavior
I should get some performance statistics

Environment

• OS name + version: macOS 11.5.1
• Version of the code: 642ab71 (HEAD -> master, origin/master, origin/HEAD, idontgetoutmuch/master) (build) update haskell.nix, use NixOS 20.09 channel (#110)

Unable to get analytic result to match sampler

models/Gamma.hs contains an old example of a Gamma prior on the precision of a normal, which I am updating. The code contains both a probabilistic program defining the conjugate prior, and an analytic version. I can't get them to agree, which is concerning, and suggests (hopefully) a bug in my code or (worse) a bug in one of the samplers. Here's the code to replicate the problem:

import Control.Monad.Bayes.Class
import Prelude
import Control.Monad.Bayes.Sampler (sampleIO)
import Control.Foldl (mean, fold, variance)
import Control.Applicative (Applicative(liftA2))
import Control.Monad.Bayes.Weighted (prior)
import Control.Monad.Bayes.Traced (mh)
import Numeric.Log (Log(Exp, ln))
import qualified Control.Foldl as F

type GammaParams = (Double, Double)

model, exact :: (MonadInfer m, Foldable t, Functor t) => GammaParams -> t Double -> m Double

model (a,b) points = do
  prec <- gamma a b
  let stddev = sqrt (1 / prec)
  let observe x = factor $ Exp $ log $ nrml 0 stddev x --  $ log $ normalPdf 0 stddev x
  mapM_ observe points
  return prec

exact (a,b) points = gamma a' b'
  where
    a' = a + fromIntegral (length points) / 2
    b' = b + sum (fmap (** 2) points) / 2


run s = fold (liftA2 (,) mean variance) <$> (fmap (take 5000) . sampleIO  . prior . mh 10000) s

-- e.g.
-- >>> run $ model (1,1) [1]
-- >>> run $ exact (1,1) [1]