Idea sketch: FRP and sequential inference

[turion]

Recently @reubenharry asked whether it makes sense to combine Functional Reactive Programming (FRP), stochastic processes, and real-time inference:

By the way, I noticed from your website that you've done a lot of FRP (e.g. dunai). This is a bit of a tangent, so maybe I can switch it to a different issue, but do you have thoughts about implementing real-time inference (as in https://arxiv.org/abs/1908.07563) using Yampa/dunai? The idea of modelling a stochastic process or real-time inference using FRP seems fitting.

And yes, I believe this makes sense. I want to try and explain the correspondence between sequential inference and arrowized FRP with effects as embodied by dunai/bearriver/essence-of-live-coding/rhine. Comments are warmly welcome :) Hopefully, we'll end up with enough material to write a tutorial on FRP & real-time inference.

CC @dataopt

Arrowized FRP with effects

MSFs

The fundamental building block of effectful AFRP is the monadic stream function:

data MSF m a b = MSF { unMSF :: a -> m (b, MSF m a b) }

An MSF runs in steps: At each step, it consumes an a, performs a side effect in m, outputs a b, and continues.

This is implemented in dunai, essence-of-live-coding (as initial algebra though), varying, netwire, and probably other libraries I'm not aware of. It's not yet FRP because we have no notion of time, but let us for the sake of simplicity assume that there are discrete and evenly spaced time steps in which we will execute our program. (If that's not satisfying to you, look at rhine and see how we recover continuity, events, classic FRP, and so on.)

One can create data (by reading data from a file, a hardware sensor, a network socket, a random number generator, ...) like this:

constM :: Functor m => m a -> MSF m () a
constM action = go
  where
    go = MSF $ \() -> ( , go) <$> action

Equally, one can consume data (by writing to a file, creating a plot, ...) with an easy to define function arrM :: Monad m => (a -> m b) -> MSF m a b. MSFs form a category, that is, they can be composed:

(>>>) :: MSF m a b -> MSF m b c -> MSF m a c

So we can add consumers and producers to our pipeline and finally arrive at a "main" MSF with no inputs and outputs that can be run:

reactimate :: MSF m () () -> m ()

(Phew, this is a really terse primer, read more about all this in dunai or essence-of-live-coding.)

Stochastic processes

Definition: A stochastic process is an MSF in a MonadSample.

In other words, a stochastic process is an MSF where we are allowed to introduce some randomness at every step of the computation.

For example, this is a Gaussian random walk (caution, I didn't typecheck or test):

type Mean = Double
type StdDev = Double
gaussianRandomWalkFrom :: Mean -> MSF m StdDev Double
gaussianRandomWalkFrom mean = MSF $ \stddev -> do
  mean' <- normal mean stddev
  return (mean', gaussianRandomWalkFrom mean')

(If you find the syntax MSF $\... daunting, have a look at dunai, which cleans this up a lot.)

Note in particular that MSFs can consume a value drawn from a prior (in this case the standard deviation), on every step. Since MSF form a category (they can be composed), you can easily build hierarchical models with prior processes, hyperprior processes and so on.

Bayes filter

I'm carrying this idea around since a few years already, but since I couldn't find any write-up of it, I'll just dump it here so people can refer to it later:

Probability monads and monadic stream functions naturally give rise to Bayes filters. MSFs here play the role of the hidden Markov model.
It basically works like this:

bayesFilter :: (MonadInfer m, Eq sensor) =>
  -- | model
  MSF m input (sensor, state) ->
  -- | external sensor, data source
  MSF m input sensor ->
  MSF m input state
bayesFilter model sensor = proc input -> do
  (estimatedOutput, estimatedState) <- model -< input
  measuredOutput <- sensor -< input
  arrM condition -< estimatedOutput == measuredOutput
  returnA -< estimatedState

In other words: At every step, draw a sample from the model and the measured data, condition on the observation being equal, and return the current state estimation. This is pretty neat and in principle works for simple MonadInfers! But you have to be careful with naive implementations like Enumerator because you'll easily end up in an exponential explosion of possible paths.

Termination

With an arbitrary monad m, there is no telling whether an MSF will ever terminate. So let's add an effect to the stack:

type MSFExcept a b m r = MSF (ExceptT r m) a b

This is the type of MSFs that can, at every step, terminate, and return a result r. This may sound unmotivated right now, but we'll need it when we discuss Coroutines.

Coroutines

As implemented in https://hackage.haskell.org/package/monad-coroutine-0.9.2/docs/Control-Monad-Coroutine.html and used by monad-bayes, a coroutine is this:

newtype Coroutine s m r = Coroutine {
   resume :: m (Either (s (Coroutine s m r)) r)
   }

A Coroutine decides (based on an action in m) at every step whether it is done (returning an r) or whether it needs to perform a further computation in s. So it is sort of a generalization to an MSFExcept, where s means "consume an a and output a b".

Lemma

newtype InOut a b x = InOut { getInOut :: a -> (b, x) }
Coroutine (InOut a b) m r == ExceptT r m (MSFExcept a b m r)

We need to understand one special case that is used in monad-bayes:

newtype Sequential m a = Sequential {runSequential :: Coroutine (Await ()) m a}

Luckily, this Await () functor is isomorphic to InOut () (), i.e. we have this isomorphism:

Sequential m a == ExceptT r m (MSFExcept () () m r)

So a Sequential is nearly the same as an MSF, except for two things:

1. The monad is specialised to ExceptT r m, so it can finish at any step
2. There is one layer of the monad around the MSF. In Sequential, this has the meaning of the part of the computation that has already been performed, while the "inner" part are the computation steps that are yet to come.

Remember the function reactimate :: Monad m => MSF m () () -> m () to execute an MSF? It specialises to reactimate :: Monad m => MSFExcept () () m r -> ExceptT r m (). The return type is isomorphic to m (Either r ()) which in our case is in principle the same as m r.

The same such function exists for Sequential, it is called finish :: Monad m => Sequential m a -> m a.

Sequential importance (re-)sampling

To summarize so far:

1. Stochastic processes that eventually return can be modelled with MonadSample m => MSFExcept a b m r
2. Composing all stochastic processes with priors etc., leads to something isomorphic to Sequential m r

This brings us to this function:

-- | Sequential importance sampling.
-- Applies a given transformation after each time step.
sis ::
  Monad m =>
  -- | transformation
  (forall x. m x -> m x) ->
  -- | number of time steps
  Int ->
  Sequential m a ->
  m a

I banged my head on this one for a while because here is where the straightforward compatibility ends. sis does two things:

1. Apply the transformation at the first n time steps (the Int parameter). (Question: Why not for all timesteps? After all, the information of how many timesteps there are is encoded in Sequential m a!)
2. Run the sequentially built up model in the m monad.

sis is the basis for all the sequential Monte Carlo methods. There is one simple strategy we can combine FRP with this:

1. Build all your models with stochastic processes as MSFs (or, if you want some more high-level description, use rhine).
2. Make sure your model eventually terminates after some number of steps. (This is awkward if we don't want to batch-process a limited dataset, but instead do live inference)
3. Transform it to a Sequential m a.
4. Use an inference method like sis, sir, smcMultinomial, ..., feeding it a sufficiently high number of time steps (again, why do we have to do this?)
5. Arrive at one (1) result.

This is all nice and dandy, and something similar is described for the pipes library in https://monad-bayes-site.netlify.app/smc, but it's not satisfying! We don't want to run the whole program once and get one result. An FRP program is an ongoing process that constantly consumes and produces data. What we actually want is to use the intermediate inference result all the time while we receive new data!

The conundrum

In terms of monad-bayes, the type signature should be something like:

liveSir :: Monad m =>
  -- | resampler
  (forall x. Population m x -> Population m x) ->
  -- | population size
  Int ->
  -- | model
  Sequential (Population m) a ->
  Sequential (Population m) a

Note that it returns a Sequential again.

A more general, analogous dunai type signature would be:

frpSir :: Monad m =>
  -- | resampler
  (forall x. Population m x -> Population m x) ->
  -- | population size
  Int ->
  -- | model
  MSF (Population m) a b ->
  MSF (Population m) a b

But it's surprisingly hard to implement this! Just following the types one would like to use hoist, but that is wrong, I believe, because it applies the resampler only on the current step, and not the whole model until then.

I think part of the difficulty is that the resampler type is higher order in the effect, but I haven't been able to figure out how to make it first order yet. Population is a composition of StateT and ListT, and while the StateT effects are all first order, I don't know how to make ListT first order.

Where to go now

I think I'll try and write a little rhine backend that:

• implements the bayes filter, for simple applications
• transforms rhine programs into Sequential, so one can at least batch-run them

I'll think a bit more about how particle resampling should work. If you have ideas, please feel free to let me know!

[idontgetoutmuch]

I just wanted to add some of my thinking (which you may already have mentioned but I didn't notice).

I think of a probabilistic program as a stochastic process (if we confine ourselves to programs that sample from a fixed number of times).

Sequential Monte Carlo is a Markov process in which the the paths are modified by a potential function given by the data. This is really close to a program being a stochastic process.

I will re-read your thoughts more carefully.

Maybe this helps

A sequence of random variables $X_{0:T}$ from $(\mathbb{X}, {\mathcal{X}})$ to $(\mathbb{X}, {\mathcal{X}})$ with joint distribution given by

$$ P_T(X_{0:T} \in {\mathrm d}x_{0:T}) = P_0(\mathrm{d}x_0)\prod_{s = 1}^T K_s(x_{s - 1}, \mathrm{d}x_s) $$

where $K_t$ are a sequence of probability kernels is called a (discrete-time) Markov process. The measure so given is a path measure.

And a program is a composition of Markov kernels. Actually I need to generalise the definition of a Markov kernel but you get the idea.

Some more notes (Feynman-Kac measures):

$$ d Q_t=\frac{1}{\mathcal{Z}_t} \exp \left[\int_0^t V_s\left(X_s\right) d s\right] d \mathbb{P} $$

$$ \mathbb{Q}_n\left(d\left(x_0, \ldots, x_n\right)\right)=\frac{1}{\mathcal{Z}_n} \exp \left[\sum_0^{n-1} V_p\left(x_p\right)\right] \mathbb{P}_n\left(d\left(x_0, \ldots, x_n\right)\right) $$

$$ Q = \frac{1}{Z_n} {\prod_{p=0}^{n-1} G_p\left(x_p\right)} \mathbb{P}_n\left(d\left(x_0, \ldots, x_n\right)\right) $$

And some more notes:

Happy to be corrected here but the maps are just translations: either you don't move or you move by $\epsilon$ so are definitely measure preserving. The densities of $\eta$ and $\epsilon$ are given so that you integrate over the (parameterised) paths, so something like

$$ \int \epsilon \left[ \int \mathbb{I} \left( \eta < \frac{f(q + \epsilon)}{f(q)}\right) \mathrm{d}\eta\right] \mathrm{d}\epsilon = \int \epsilon \left[ 1 \wedge \frac{f(q + \epsilon)}{f(q)} \right] \mathrm{d}\epsilon $$

which we can write as

$$ \int \alpha(z, z')p(z,z')\mathrm{d}z' $$

where

$$ \alpha(z, z') = 1 \wedge r(z, z') \quad \textrm{and} \quad r(z, z') = \frac{f(z')p(z',z)}{f(z)p(z, z')} $$

And this is the familiar Metropolis Hastings kernel with a Gaussian random walk.

Apologies for missing out some steps and not considering the case where the chain stays where it is.

[reubenharry]

Thanks for writing this up! I'll read through it more carefully soon and respond in more detail with my ideas/plans (I have been moving this week, so things are a bit chaotic). But basically, I have been thinking on very similar lines! A few things that may or may not be relevant:

• on a local branch (which I can push) I tried implementing "live" inference using pipes and/or streamly. The idea was to work in the IO monad, but unfold a stream where the state is a Population, and resample it at each step. It seemed to work, but there was very clearly a leak of some kind, because it got slow quickly. Moreover, I'd rather do this with an FRP library than a streaming library
• you asked: "Question: Why not for all timesteps? After all, the information of how many timesteps there are is encoded in Sequential m a". I also thought this, and spent a week trying to write a version of sis that could do all timesteps. I didn't succeed, although this certainly doesn't mean it's impossible. I have a suspicion that there's a conceptual problem, because the number of timesteps depends on the stochastic choices in the program, which can differ between different particles (recall that in sir, a population is a population of coroutines).
• are there ideas we can borrow directly from https://arxiv.org/abs/1908.07563 and implement in dunai/rhine? Or even maybe in Yampa?

[reubenharry]

So in my own attempts to write a pipes based version of a similar idea, I ended up with the type

Pipe observation (Population SamplerIO latent) IO ()

For example, I was trying to write a system which would receive a stream of Bools, each of which was known to be generated from some distribution P(b| latent), and live infer latent using SMC. Given some predicate on the expectation of the posterior over latent, I would perform some IO action that would change the nature of the input stream. This would run continuously.

Since the core monad is IO, the pipe is able to describe both reading in data from an input source and writing out actions at a given step, based on beliefs. The probability monad doesn't reside in the pipe per se, but just in the output type.

Is your suggestion to have a more functional and pure version of the same idea by using MSFs? Or to put it another way, can you say more about

frpSir :: Monad m =>
  -- | resampler
  (forall x. Population m x -> Population m x) ->
  -- | population size
  Int ->
  -- | model
  MSF (Population m) a b ->
  MSF (Population m) a b

What does the input signal function represent? And what is the output one? Can it handle the use case described above?

[turion]

But basically, I have been thinking on very similar lines!

Awesome :) I think we'll arrive at something cool in the end.

* on a local branch (which I can push) I tried implementing "live" inference using `pipes` and/or `streamly`. The idea was to work in the `IO` monad, but unfold a stream where the state is a `Population`, and resample it at each step. It _seemed_ to work, but there was very clearly a leak of some kind, because it got slow quickly.

Yes, arrowized FRP was designed partly to deal with some classes of leaks. One reason I don't like streamly a lot is because it is easy to run into leaks when trying to work with the internals. I guess with pipes one should be able to get something to work, but I haven't looked into it.

Moreover, I'd rather do this with an FRP library than a streaming library

The libraries I usually work with are basically streaming libraries (dunai) extended to be time-aware and sample-rate-independent (rhine)

* you asked: "Question: Why not for all timesteps? After all, the information of how many timesteps there are is encoded in Sequential m a". I also thought this, and spent a week trying to write a version of `sis` that could do all timesteps. I didn't succeed, although this certainly doesn't mean it's impossible.

Yes, I think it's possible somehow. But I ended up not being able to use sis, and instead rewriting it in a reactive way. See below.

I have a suspicion that there's a conceptual problem, because the number of timesteps depends on the stochastic choices in the program, which can differ between different particles (recall that in sir, a population is a population of coroutines).

You mean, one particle does two factors, and the other one three?

* are there ideas we can borrow directly from https://arxiv.org/abs/1908.07563

That's a fascinating article, but I still haven't finished reading it. Definitely a lot that can be implemented.

and implement in dunai/rhine?

Yes, certainly!

Or even maybe in Yampa?

Probably not, Yampa knows no monads/effects.

So in my own attempts to write a pipes based version of a similar idea, I ended up with the type

Pipe observation (Population SamplerIO latent) IO ()

So does that mean that Population is the output here? If yes we arrived at something similar.

For example, I was trying to write a system which would receive a stream of Bools, each of which was known to be generated from some distribution P(b| latent), and live infer latent using SMC. Given some predicate on the expectation of the posterior over latent, I would perform some IO action that would change the nature of the input stream. This would run continuously.

Since the core monad is IO, the pipe is able to describe both reading in data from an input source and writing out actions at a given step, based on beliefs. The probability monad doesn't reside in the pipe per se, but just in the output type.

That's the part that I don't quite understand yet. How can you "extract" the randomness to the output?

Is your suggestion to have a more functional and pure version of the same idea by using MSFs? Or to put it another way, can you say more about

frpSir :: Monad m =>
  -- | resampler
  (forall x. Population m x -> Population m x) ->
  -- | population size
  Int ->
  -- | model
  MSF (Population m) a b ->
  MSF (Population m) a b

The type signature I proposed doesn't make much sense it turns out, and instead we can do the following which works:

runPopulationS :: forall m a b . Monad m =>
  -- | Number of particles
  Int ->
  -- | Resampler
  (forall x . Population m x -> Population m x)
  -> MSF (Population m) a b
  -> MSF m a [(b, Log Double)]
runPopulationS nParticles resampler msf = runPopulationCl' $ spawn nParticles $> msf
  where
    runPopulationCl' :: Monad m => Population m (MSF (Population m) a b) -> MSF m a [(b, Log Double)]
    runPopulationCl' msfs = MSF $ \a -> do
      bAndMSFs <- runPopulation $ flip unMSF a =<< msfs
      let (currentPopulation, continuations) = unzip $ (\((b, msf), weight) -> ((b, weight), (msf, weight))) <$> bAndMSFs
      -- FIXME This normalizes, which introduces bias, whatever that means
      return (currentPopulation, runPopulationCl' $ normalize $ resampler $ fromWeightedList $ return continuations)

An MSF m a b where MonadSample m is a stochastic process. If we add particles to the monad stack, we get MonadInfer (Population m) and can condition in the MSF, for example to use a bayes filter.

What I'm doing here is something similar like you were doing: Resample every step of the MSF and output the current particle population of the output values explicitly. It works really well for me and I have a nice interactive graphical example that I'll polish and upload these days. What I don't understand is whether it is ok to normalize each step. I need to do that, otherwise I loose probability mass pretty quickly (10^-100 in a thousand steps or so, reaching the precision of Doubles). But it says that "this introduces bias" in the docs, and I don't understand what that means in this context.

Another work-in-progress question: I could simplify the type of the bayes filter:

-- | Condition on one output of a distribution.
--
--   p(x,y | theta) ~> p(x | y, theta)
bayesFilter :: (MonadInfer m, Eq sensor) =>
  MSF m input (sensor, latent) ->
  -- | external sensor, data source
  MSF m (input, sensor) latent

One gets an additional input, which is used to condition the sensor. I found that in practice, this doesn't work well with particles because hard equality comparison throws away too many particles. One usually uses something called an "importance function" I believe, and you also use a normal distribution in the SMC example, instead of a hard condition $ x == y. I turned this into a type class:

class SoftEq a where
  -- | `similarity a1 a2 == 1` if they are exactly equal, and 0 if they are completely different.
  similarity :: a -> a -> Log Double

  -- | Scores the similarity of the two inputs
  (=~) :: MonadInfer m => a -> a -> m ()
  a1 =~ a2 = score $ similarity a1 a2

Does that make sense to you?

[turion]

See https://github.com/turion/dunai-bayes and https://github.com/turion/rhine/tree/dev_monad_bayes for work in progress

[reubenharry]

This is really great, way cooler than my pipes attempt. I'm reading through your code and in general it makes sense, but here are some things that are currently confusing me:

• Is there a way I can get this up and running locally? I would love to play around with it. Git cloning dunai-bayes and doing cabal repl (is that the right command? I'm used to stack) gave "Could not resolve dependencies". EDIT: nvm, I got it up and running using stack.

• it seems that there are imports going both ways from dunai-bayes to rhine and vice versa. Is that intentional?

• for the type type MySmallMonad = Sequential (GlossConcT SamplerIO), why have Sequential? I thought that the FRP was taking care of that stuff, instead of Sequential. EDIT: I tried replacing Sequential with IdentityT and making the appropriate changes, and things seemed to work.

Btw, I did try using MMorph in monad-bayes, but there was some problem I can't remember. I think one of the hoists for some monad wasn't actually writeable as an MMorph instance. I can try again though.

[reubenharry]

I turned this into a type class:
Does that make sense to you?

Yes, it does.

Re normalization, are you just saying that the absolute values of the weights in the population become too small? That confuses me, because in resampleGeneric, the weights get normalized automatically. So I think maybe I'm misunderstanding.

You mean, one particle does two factors, and the other one three?

Yep, because a particle is in fact a distribution over continuations of the program, whose fate depends on the stochastic choices made at each step.

[reubenharry]

(I'll keep posting questions as I come up with them)

An MSF m a b where MonadSample m is a stochastic process

To make sure I'm understanding, wouldn't that be true in particular when a is (). I would have thought MSF m a b would be a map from one stochastic process to another. Or am I getting that wrong?

[reubenharry]

Can we write runPopulationCl' using morphGS?

[reubenharry]

By the way, the sort of application I have in mind here (and it seems like you're not that far off), is a GUI app where a user can produce observations live (e.g. by clicking points on the screen) and the system will update its posterior in real-time. One could even imagine changing the population size or resampling rate live.

[reubenharry]

Yet another question: is there an easy way using rhine to get the current time as a signal? (I read your paper on Rhine btw, and it seems very cool). It would be cool to define a Gaussian process by just expressing the jump to the next state after time t as a Gaussian parametrized by t (I mainly ask this because I have always found stochastic processes confusing, and writing them this way would be pleasantly clarifying for me)

EDIT: I think this works:

model4 :: forall cl m td. (Float ~ Time cl, MonadSample m) => ClSF m cl () Pos2
model4 = feedback (0,0) model4' 
model4' = MSF \((), (x,y)) -> do
  time <- DunaiReader.asks sinceLast -- @(TimeInfo cl) @ReaderT
  x' <- normal x (float2Double (time*2))
  y' <- normal y (float2Double (time*2))
  let p = (x',y')
  return ((p,p), model4')

[turion]

* Is there a way I can get this up and running locally? I would love to play around with it. Git cloning dunai-bayes and doing `cabal repl`

Yes, that's the right way. But you need to make sure that rhine has a cabal.project.local where you specify packages: path/to/dunai-bayes/dunai-bayes.cabal.

* it seems that there are imports going both ways from dunai-bayes to rhine and vice versa. Is that intentional?

No, dunai-bayes imports from dunai and monad-bayes.

* for the type `type MySmallMonad = Sequential (GlossConcT SamplerIO)`, why have `Sequential`? I thought that the FRP was taking care of that stuff, instead of `Sequential`. EDIT: I tried replacing `Sequential` with `IdentityT` and making the appropriate changes, and things seemed to work.

I had Sequential in there at first because I thought I might be able to reuse the sis etc. functions. Now its only advantage is that it does a suspension every time we use condition. But yes, we might be able to get rid of it now.

Btw, I did try using MMorph in monad-bayes, but there was some problem I can't remember. I think one of the hoists for some monad wasn't actually writeable as an MMorph instance. I can try again though.

Yes, do try, I'll try and help if it doesn't work!

Btw. about naming: Sequential should rather be named SequentialT, and so on for several other transformers.

Re normalization, are you just saying that the absolute values of the weights in the population become too small?

Yes.

That confuses me, because in resampleGeneric, the weights get normalized automatically. So I think maybe I'm misunderstanding.

Makes sense what you're saying, yet still if I remove normalize, I quickly lose the probability mass. So maybe there is a bug in resampleGeneric?

In fact, I don't understand this line here:

monad-bayes/src/Control/Monad/Bayes/Population.hs

Line 117 in 598464d

return $ map (,z / fromIntegral n) offsprings

Why do we multiply with the previous mass? That way, we guarantee that the probability mass is exactly the same as before resampling. Maybe this is intended (to avoid "bias"?).

You mean, one particle does two factors, and the other one three?

Yep, because a particle is in fact a distribution over continuations of the program, whose fate depends on the stochastic choices made at each step.

Ok, then it makes more sense to remove Sequential completely in the FRP approach.

An MSF m a b where MonadSample m is a stochastic process

To make sure I'm understanding, wouldn't that be true in particular when a is (). I would have thought MSF m a b would be a map from one stochastic process to another. Or am I getting that wrong?

Yes, I guess you're right.

Can we write runPopulationCl' using morphGS?

I'm pretty sure that not. The issue is that we want to keep joining the population effects "in the background" over all timesteps, but morphGS only transforms pointwise per timestamp. But that's not to say that there isn't some other general MSF handling principle that runPopulationS is a special case of.

By the way, the sort of application I have in mind here (and it seems like you're not that far off), is a GUI app where a user can produce observations live (e.g. by clicking points on the screen) and the system will update its posterior in real-time. One could even imagine changing the population size or resampling rate live.

Yes, that's sort of what I wanted to do. Spawning or killing particles should be done on user input. The resampling rate is equal to the simulation speed, which is currently tied to the graphics FPS rate, but there is no a priori reason why graphics and simulation has to be synchronous. (In fact rhine is about being able to decouple these two, it's just that I haven't done that so far.)

I'm not so convinced that the user should produce observations. It's super cool to have user input as measurements (or some other input like a physical sensor or so), but for the model to make sense we have to have some kind of rule how the latent value produces an observable output. In the current model, there is a fixed normal error, and we cannot expect the user to play a noisy sensor and produce roughly this noise distribution.

We could put a prior on the noise, and then tell the user (as some kind of game) to try to click onto the latent position, and thus fit the noise that the user themself introduces by inaccurate clicking. Maybe this is a fun minigame. Or maybe it's too convoluted, what do you think?

Yet another question: is there an easy way using rhine to get the current time as a signal? (I read your paper on Rhine btw, and it seems very cool). It would be cool to define a Gaussian process by just expressing the jump to the next state after time t as a Gaussian parametrized by t (I mainly ask this because I have always found stochastic processes confusing, and writing them this way would be pleasantly clarifying for me)

EDIT: I think this works:

model4 :: forall cl m td. (Float ~ Time cl, MonadSample m) => ClSF m cl () Pos2
model4 = feedback (0,0) model4' 
model4' = MSF \((), (x,y)) -> do
  time <- DunaiReader.asks sinceLast -- @(TimeInfo cl) @ReaderT
  x' <- normal x (float2Double (time*2))
  y' <- normal y (float2Double (time*2))
  let p = (x',y')
  return ((p,p), model4')

Yes, that's what rhine (and bearriver) are for 😀 you are probably looking for https://hackage.haskell.org/package/rhine-0.8.0.1/docs/FRP-Rhine.html#v:sinceLastS.

I think I'd try it like this:

model4 :: forall cl m td. (Double ~ Time cl, MonadSample m) => ClSF m cl () Pos2
model4 = gaussianProcess &&& gaussianProcess

gaussianProcess :: forall cl m td. (Double ~ Time cl, MonadSample m) => ClSF m cl () Double
gaussianProcess = feedback (0,0) $ proc ((), x) -> do
  timeSinceLast <- sinceLastS -< ()
  x' <- arrMCl $ uncurry normal -< (x, (time*2))
  return (x', x')

Comments: Floats are stupid and I only use them because of gloss, but I'm trying to push them as far away as possible. &&& combines two independent processes. Prefer not to use the MSF constructor. (I personally like arrow style, matter of taste.)

Workflow

I think we've collected enough material so that we're sure there is something here that can work. How do we proceed now? I'd propose:

• This issue here should be about developing a tutorial on how to use rhine-bayes and dunai-bayes. We can go on developing, but shouldn't merge any material before rhine-bayes and dunai-bayes are merged.
• Questions & improvements about non-time aware stream processing are issues & PR on https://github.com/turion/dunai-bayes
• Everything time aware (e.g. Gaussian process, little standard library of stochastic processes) and in particular the demo app involving real time graphics and user interaction land as comments on Dev monad bayes turion/rhine#186. Once it's merged, issues & PRs on https://github.com/turion/rhine/.

Unfortunately I don't have as much time as I'd like, I have to squeeze this project in between/after my work hours. But this is super exciting, and I'll continue whenever I can!

Since you're already in it and digging, does it make sense if I keep developing rhine-bayes and dunai-bayes and you develop a tutorial (if this is at all interesting for you) while asking all sorts of questions and requesting improvements?

[reubenharry]

Yes, do try, I'll try and help if it doesn't work!

I've made an issue, and will try when free.

Btw. about naming: Sequential should rather be named SequentialT, and so on for several other transformers.

Also made an issue.

Why do we multiply with the previous mass? That way, we guarantee that the probability mass is exactly the same as before resampling.

I think that's intentional. I also found this confusing, but I think the point is to keep the weights from exploding/vanishing. I'm pretty sure it's correct, but we should try to work out why you have to normalize to make things work. Once I understand the FRP code more, I can hopefully work that out.

Prefer not to use the MSF constructor. (I personally like arrow style, matter of taste.)

I've never used the arrow notation before, but this is a pretty good advert for it :) And thanks, sinceLastS was exactly what I wanted, I just didn't find it.

I'm not so convinced that the user should produce observations. It's super cool to have user input as measurements (or some other input like a physical sensor or so), but for the model to make sense we have to have some kind of rule how the latent value produces an observable output

True, but on the other hand, the nice thing about being in Bayesian land is we can specify models that somewhat accommodate this. For example, we could have a model which tries to guess whether the input is coming from a sin wave with noise, or from uncorrelated random samples (an ok description of user input). Or just to do (non-parametric) polynomial regression and show the best fit. Reading sensor data live would be really nice too though. Anyway, that's all down the road...

Workflow

I think we've collected enough material so that we're sure there is something here that can work

Agreed!

This issue here should be about developing a tutorial on how to use rhine-bayes and dunai-bayes. We can go on developing, but shouldn't merge any material before rhine-bayes and dunai-bayes are merged.

Are you planning to have rhine-bayes and dunai-bayes be separate hackage packages from rhine and dunai? Am assuming yes, but just checking.

Everything time aware (e.g. Gaussian process, little standard library of stochastic processes) and in particular the demo app involving real time graphics and user interaction land as comments on turion/rhine#186. Once it's merged, issues & PRs on https://github.com/turion/rhine/.

Yep, agreed.

Unfortunately I don't have as much time as I'd like, I have to squeeze this project in between/after my work hours. But this is super exciting, and I'll continue whenever I can!

That's pretty much where I'm at too. One of the reasons I find this exciting is that it feels like Haskell has a unique advantage at tackling this problem (because of mature FRP libraries, arrow notation, and monad-bayes) in a really declarative way.

Since you're already in it and digging, does it make sense if I keep developing rhine-bayes and dunai-bayes and you develop a tutorial (if this is at all interesting for you) while asking all sorts of questions and requesting improvements?

Yes, that sounds great because a) I like writing tutorials and b) it will help me understand all the moving pieces properly, which I don't quite yet.

Btw, are we envisioning that there are changes needed to monad-bayes to make this work better? Perhaps generalization from double, for example. But if there are other things, let me know. One interesting point to consider is that we could easily extend to more complex inference algorithms, e.g. doing a step of MH every so often (or even on the user's request) to rejuvenate the population. The non-reactive version of that is called RMSMC, and is basically a one-liner in monad-bayes.

[reubenharry]

OK, closing this issue as activity is happening in dunai-bayes and rhine-bayes (see above for links). tldr: combining reactive and probabilistic programming is a great idea, and works.