Add bus example and necessary Poisson pdf (#232)

* Add bus example and necessary Poisson pdf

* Remove cruft

* Fix lint complaint

* Add more context to the example

src/Control/Monad/Bayes/Class.hs

@@ -58,6 +58,7 @@ module Control.Monad.Bayes.Class
     discrete,
     normalPdf,
     Bayesian (..),
+    poissonPdf,
     posterior,
     priorPredictive,
     posteriorPredictive,
@@ -96,7 +97,7 @@ import Data.Vector.Generic as VG (Vector, map, mapM, null, sum, (!))
 import Numeric.Log (Log (..))
 import Statistics.Distribution
   ( ContDistr (logDensity, quantile),
-    DiscreteDistr (probability),
+    DiscreteDistr (logProbability, probability),
   )
 import Statistics.Distribution.Beta (betaDistr)
 import Statistics.Distribution.Gamma (gammaDistr)
@@ -285,6 +286,9 @@ normalPdf ::
   Log Double
 normalPdf mu sigma x = Exp $ logDensity (normalDistr mu sigma) x
 
+poissonPdf :: Double -> Integer -> Log Double
+poissonPdf rate n = Exp $ logProbability (Poisson.poisson rate) (fromIntegral n)
+
 -- | multivariate normal
 mvNormal :: MonadDistribution m => V.Vector Double -> Matrix Double -> m (V.Vector Double)
 mvNormal mu bigSigma = do

src/Control/Monad/Bayes/Traced/Static.hs

@@ -82,8 +82,35 @@ mhStep (Traced m d) = Traced m d'
   where
     d' = d >>= mhTransFree m
 
+-- $setup
+-- >>> import Control.Monad.Bayes.Class
+-- >>> import Control.Monad.Bayes.Sampler.Strict
+-- >>> import Control.Monad.Bayes.Weighted
+
 -- | Full run of the Trace Metropolis-Hastings algorithm with a specified
 -- number of steps. Newest samples are at the head of the list.
+--
+-- For example:
+--
+-- * I have forgotten what day it is.
+-- * There are ten buses per hour in the week and three buses per hour at the weekend.
+-- * I observe four buses in a given hour.
+-- * What is the probability that it is the weekend?
+--
+-- >>> :{
+--  let
+--    bus = do x <- bernoulli (2/7)
+--             let rate = if x then 3 else 10
+--             factor $ poissonPdf rate 4
+--             return x
+--    mhRunBusSingleObs = do
+--      let nSamples = 2
+--      sampleIOfixed $ unweighted $ mh nSamples bus
+--  in mhRunBusSingleObs
+-- :}
+-- [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)
   where