Add example of distribution driven by an ODE

[adolgert]

This would be a distribution that is a delta function, where the time of the distribution comes from an ordinary differential equation (ODE) that depends on the state of the system.

[slwu89]

The simplest example is something like this "gene regulation with positive feedback" from the R package pdmpsim: https://github.com/CharlotteJana/pdmpsim?tab=readme-ov-file#a-simple-example

There are also examples in the Julia package https://rveltz.github.io/PiecewiseDeterministicMarkovProcesses.jl/latest/

[slwu89]

@adolgert may I suggest that we take this issue off the list of things to do before the initial release? There is actually some delicate coding required to get PDMP simulation working correctly, and I think that any example would be too overwhelming, and actually begin to exceed the complexity for a "tutorial". In short, I think it would start to become another module or even package itself. The most recent review of simulation algorithms for PDMPs is: https://arxiv.org/abs/1504.06873

First I introduce some notation, largely from that paper. Let $X$ be the PDMP, any particular state of the process is given by $x$. The jump instants are an increasing sequence $T_n$ of times. In between jumps, the dynamics of $X$ are governed by a differential equation whose right hand side is $F(x)$, written $\dot{\phi(x)} = F(\phi(x))$. There is also a function $R(x)$ which gives the hazard of the clock(s) that are enabled, which may depend on the continuous states by $R(\phi(x))$. Functions $F,R$ also depend on the discrete part of the state (and on time), but we suppress that notation because it gets too overwhelming. Given the time of the last jump $T_{n-1}$, there are 3 basic algorithms to sample the next jump time $T_n$.

The first type is classic inversion.

1. Draw $S\sim Exp(1)$
2. Solve for $t$ such that (inversion): $S = \int_{T_{n-1}}^{t} R(\phi(x)) ds$, which involves solving the differential equation given by $F$
3. Let $T_n$ be equal to $t$

The second type is known as change-in-variable, and augments the state of the continuous state such that an inversion problem does not have to be solved. It involves changing the continuous flow between jumps to be:

$$ \dot y(s)= F(y(s)) / R(y(s)) $$

$$ \dot\tau(s)=1 / {R(y(s))} $$

Subject to initial conditions $y(0)=X(T_{n-1}) \tau(0)=T_{n-1}$ Then given $S\sim Exp(1)$, the continuous state at the next jump time and the next jump time are given by $(y(S_n),\tau(S_n))$

Finally, there is a rejection method. It needs a rate for the Poisson process which will be used to draw the potential jump times, $\lambda$.

1. Draw $S\sim Exp(\lambda)$
2. To decide the next jump time $T_{n-1}$, accept with probability $\frac{R(X(T_n + S))}{\lambda}$ (note this still requires integration of the differential equations, up to the time of the potential jump).

[adolgert]

Hi Sean. That's a thoughtful suggestion, and I agree. And your explanation makes sense to me, but it will take some work. I'll take it off the list in the planner, if I can figure out how, but it's off the list.