turion / rhine

Haskell Functional Reactive Programming framework with type-level clocks
http://hackage.haskell.org/package/rhine
117 stars 21 forks source link

Variable rate observations #243

Open reubenharry opened 11 months ago

reubenharry commented 11 months ago

It would be nice to have a demo of object tracking where the observations arrive on a Poisson process clock. This shouldn't be super hard for me to do, right? Something like:

observationModel @@ PoissonClock   >--resBuf--> particleFilter posterior @@ PeriodicClock
turion commented 5 months ago

There is multiple things one might want to do. Let's use an example, maybe it's easier to reason about which parts we want to do inference.

The physical system consists of a radioactive material that randomly emits energetic particles, and an object that is heated up by them. The temperature of the object is measured by a noisy sensor.

We can model this system as follows:

avgEnergy = 10

-- Assuming some positive distribution for the energy
material = constMCl $ gamma 1 avgEnergy

sensorNoise = 1

sensor = arrM $ \energy -> (energy +) <$> normal 0 sensorNoise

modelRhine poissonClock = material @@ poissonClock >-- collect --> arr sum >-> sumS >-> sensor @@ SensorClock

poissonClock :: PoissonClock ticks for every Poisson event, with some constant density. SensorClock right now is completely formal, it might e.g. be a fixed rate.

Now the real world observation will be more like:

constM getMeasurements @@ sensorClockIO

sensorClockIO corresponds to SensorClock, but is not the same. One might have type SensorClockIO = RescaledClockM IO SensorClock UtcTime and then implementing sensorClockIO :: SensorClockIO by waiting until a given time and then polling on the actual sensor in IO.

What might one want to estimate here? You could estimate e.g. the rate of the Poisson process or the average energy using PMMH or so, but we have not specified where the prior for the clock enters the model. In the WIP PR for Poisson clocks, the rate is part of the clock value, so it has to be drawn before building up the whole rhine, i.e. we have something like

mkPoissonClock :: MonadDistribution m => m PoissonClock
mkPoissonClock = PoissonClock <$> somePrior

runModel = sampleIO $ do
  poissonClock <- mkPoissonClock
  flow $ modelRhine poissonClock

So to estimate the parameter, one would have to do PMMH (#242) or similar. But I think the inference part would run on SensorClockIO, it doesn't need knowledge of the clock structure of the model. Probably, one would do clock erasure on the model, essentially turning it into an MSF, and then do a synchronous particle filter on it.

This does not work anymore if the Poisson clock or a corresponding clock (say, a real world Geiger counter accessed via IO) is part of the Rhine where you want to do inference in. I don't know yet what one would do there.

Now let's say that, instead of a parameter, you want to estimate an internal state of the object, say, the total energy. You'd still need to sample the parameters in some way (or jointly do inference on them), and I guess you'd still do clock erasure on the whole model, and only score with the sensor data.

This is different when you'd like to estimate some state of the running clock, say, the number of radioactive particles emitted.