A Poisson processes can be thought of as a stochastic process that increments a counter variable (taking integer values) at random times. The increments are assumed to be exponentially distributed. The master equation, which describes the change in probability of being in a state $n$ at time $t$, is:
\[ \dot{p}_n(t) = q(p_{n-t} - p_n) \]
where $q$ is the transition rate (in jumps per unit time) for changing state.
For this example, we will assume that $n$ ranges between 0 and 10. We are interested in the probability distribution of first-passage times to reach $n = 10$ when starting at $n = 0$, that is, how long does it take this system to reach its maximum counter value? We assume that state $n = 10$ is an absorbing state from which the system cannot escape once it reaches it. The remaining states are transient states, which the system is free to enter or leave.
We assume that $\dot{p}_0 = -qp_0$ and that $\dot{p}_{10} = qp_9$, which creates a system of coupled differential equations. To solve these, we define the initial probability distribution over states to be
p0 = zeros(9) p0[1] = 1.0
1.0
We only define the initial conditions for the transient states, the states that are not absorbing, $p_n, n\in [1, \dots, 9]$.
As shown in, e.g., van Kampen (2007, Stochastic processes in physics and chemistry, Ch. 12), the solution for the probability distribution of the time to reach state $n = 10$ given that the system started at $n = 0$ is
\[ p(t) = \mathbf{A}\exp\left(\mathbf{T}t\right)\mathbf{p}(0) \]
where $\mathbf{p}(0)$ is a column vector of initial probability distribution over states defined above, $\mathbf{T}$ is the matrix of transition rates between transient states, and $\mathbf{A}$ is the matrix of transition rates from the transient states into the absorbing state:
9×9 Matrix{Float64}: -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0
1×9 Matrix{Float64}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
The non-zero entries are set to one for now, but they can simply be rescaled by $q$, which we will do below.
Say we have a data set $y_i$ of $n_{\text{data}}$ data points that we assume come are generated by this absorbing Poisson process with an unknown transition rate parameter $q$. How can we fit $q$ to the data?
We start by loading the packages we need. We'll using Turing.jl
to implement Bayesian parameter fitting of $q$. FirstPassageTools.jl
directly implements the matrix exponential function for the first-passage time distribution in the equation above, and currently, only the Zygote.jl
package implements the necessary automatic differentiation procedures need to use the NUTS() sampler in Turing.
using Distributions using Turing using Zygote Turing.setadbackend(:zygote) using Plots, StatsPlots using Pkg Pkg.activate("../../FirstPassageTools.jl/") using FirstPassageTools
We will assume that $q = 7.5$ transitions per second is the true value that generates the data. We can now set up the first-passage time distribution:
250-element Vector{Float64}: 0.9947187540463421 1.014096157078322 0.5874830550912131 0.9921950439813819 1.1867765065646365 1.6524274459410035 0.9825300343606603 0.8982433198039106 1.6519696611649775 0.6945559768077825 ⋮ 0.8673186612838045 0.6547792730950205 1.5244870123948495 1.1491937149210338 0.6658814142406704 1.3352308340424899 1.5333444239354312 0.8800907884007108 1.578802372336248
histogram(data, xlabel="First-passage times", normalize=:true, label="Random data") plot!(fp, label="True density")
We will assume the following model (priors and likelihood):
@model function mod(y, T=T, A=A, p0=p0) # Set prior on q q ~ truncated(Normal(5, 1), 0.001, Inf) # Likelihood y ~ filldist(fpdistribution(q*T, q*A, p0), length(y)) end
mod (generic function with 5 methods)
We assume a truncated normal distribution (SD = 5) for $q$. The data are distributed according to the first-passage time distribution defined above, but rescaled by $q$.
We can now sample from the posterior distribution of $q$:
posterior = sample(mod(data), NUTS(500, 0.65), MCMCThreads(), 1000, 4)
Chains MCMC chain (1000×13×4 Array{Float64, 3}): Iterations = 501:1:1500 Number of chains = 4 Samples per chain = 1000 Wall duration = 250.13 seconds Compute duration = 854.25 seconds parameters = q internals = lp, n_steps, is_accept, acceptance_rate, log_density, h amiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size Summary Statistics parameters mean std naive_se mcse ess rhat ⋯ Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯ q 7.4030 0.1525 0.0024 0.0033 2107.7598 0.9998 ⋯ 1 column om itted Quantiles parameters 2.5% 25.0% 50.0% 75.0% 97.5% Symbol Float64 Float64 Float64 Float64 Float64 q 7.1044 7.3033 7.4029 7.5037 7.7068
Let's have a look at the trace plot and the posterior density of $q$:
It looks like we recovered the true value of $q$ pretty well.