Place Holder for the Fellowship Work

[idontgetoutmuch]

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.

[idontgetoutmuch]

I've been looking through some of the branches on monad-bayes on Github (including your integration of Random-fu and multivariate normals), and I noticed that Adam had written a Random-fu sampler in this branch: https://github.com/tweag/monad-bayes/blob/kernel/src/Control/Monad/Bayes/Experimental/RandomFu.hs

I wasn't sure if you were aware of this, so I thought it was maybe worth mentioning. I don't understand enough about sampling to follow the lazy sampling part, but it looks like he made RVarT an instance of MonadSample (the typeclass has a different name here because it's an old

[idontgetoutmuch]

It does look somewhat similar and no I wasn’t aware of it. System.Random has changed quite a bit since then so I don’t know how much is now transferable.

https://github.com/tweag/monad-bayes/pull/112/files

One thing that we should do is make this independent of the RNG so that the user can choose: the standard RNG, MWC, Mersenne Twister or some else. Also we really need to do performance measures without the “noise” from my changes to random-fu.

[idontgetoutmuch]

I just discovered that you can pass a resampler as an argument to the sequential Monte Carlo function: https://monad-bayes.netlify.app/usage.html#sequential-monte-carlo - cool.

[idontgetoutmuch]

Here's a stratified resampler that presumably needs modifying to pass to the sequential Monte Carlo function:

resampleStratified :: (UniformRange d, Ord d, Fractional d) => [d] -> [Int]
resampleStratified weights = catMaybes $ unfoldr coalg (0, 0)
  where
    bigN      = length weights
    positions = map (/ (fromIntegral bigN)) $
                -- FIXME: the generator should be configurable
                zipWith (+) (take bigN . unfoldr (Just . uniformR (0.0, 1.0)) $ mkStdGen 23)
                            (map fromIntegral [0 .. bigN - 1])
    cumulativeSum = scanl (+) 0.0 weights
    coalg (i, j) | i < bigN =
                     if (positions!!i) < (cumulativeSum!!j)
                     then Just (Just j, (i + 1, j))
                     else Just (Nothing, (i, j + 1))
                 | otherwise =
                     Nothing

[reubenharry]

Yeah, the modification needed would presumably be to make the type more generic. The resamplers that SMC can take, as the type shows, have to be (forall x. Population m x -> Population m x) i.e. to be agnostic as to the type of x (or actually if you look at the internals, any (forall x. m x -> m x) will do). Unclear to me as yet how to modify the above, since it seems to involve doing something particular with numbers, whereas the other resamplers in the code are suitably type agnostic.

[reubenharry]

I might need a pointer to good materials on random samplers, as I don't know too much about how they work and why some are better than others.

Re. making the RNG independent, perhaps I can add a reader layer to the Sampler monad, so that it gets an RNG that the user can then supply at their choice. Will all RNGs have the same type?

[idontgetoutmuch]

I think you can use https://github.com/tweag/monad-bayes/blob/master/src/Control/Monad/Bayes/Population.hs#L76 but it's not exported. Let me try exporting it and then using my stratified resample and compare the results.

[idontgetoutmuch]

If I export this function then this seems to work (it probably should have more tests before I submit it as a PR):

module SmcExample (run, run') where

import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Population
import Control.Monad.Bayes.Inference.SMC
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Sequential (Sequential)

import Numeric.Log (Log)
import Data.List
import qualified Data.Vector as V
import Data.Vector ((!))
import Data.Maybe (catMaybes)

example :: MonadInfer m => m Bool
example = do
  x <- bernoulli 0.5
  condition x
  return x

run :: IO [(Bool, Log Double)]
run = (sampleIO . runPopulation. smcSystematic 4 4) example

stratified :: MonadSample m => V.Vector Double -> m [Int]
stratified weights = do
  let bigN = V.length weights
  dithers <- V.replicateM bigN (uniform 0.0 1.0)
  let positions = V.map (/ (fromIntegral bigN)) $
                  V.zipWith (+) dithers (V.map fromIntegral $ V.fromList [0 .. bigN - 1])
      cumulativeSum = V.scanl (+) 0.0 weights
      coalg (i, j) | i < bigN =
                       if (positions!i) < (cumulativeSum!j)
                       then Just (Just j, (i + 1, j))
                       else Just (Nothing, (i, j + 1))
                   | otherwise =
                       Nothing
  return $ map (\i -> i - 1) $ catMaybes $ unfoldr coalg (0, 0)

resampleStratified ::
  (MonadSample m) =>
  Population m a ->
  Population m a
resampleStratified = resampleGeneric stratified

smcStratified ::
  MonadSample m =>
  Int ->
  Int ->
  Sequential (Population m) a ->
  Population m a
smcStratified = sir resampleStratified

run' :: IO [(Bool, Log Double)]
run' = (sampleIO . runPopulation. smcStratified 4 4) example

[idontgetoutmuch]

A useful reference for resampling methods https://bisite.usal.es/archivos/resampling_methods_for_particle_filtering_classification_implementation_and_strategies.pdf