Jumps
Continuous-time Markov chains by uniformisation, birth–death and telegraph processes, compound Poisson moments and Gillespie's stochastic simulation step.
This module covers processes in continuous time with a discrete state. The central technique is uniformisation: instead of a matrix exponential, the transition probabilities of a generator \(Q\) are written as a Poisson-distributed number of jumps of an ordinary discrete chain. Choose \(\lambda\) at least as large as the fastest exit rate and set \(J = I + Q/\lambda\), a genuine stochastic matrix; then
\[\pi(t) = \sum_{k \ge 0} e^{-\lambda t} \frac{(\lambda t)^k}{k!}\; \pi(0) J^k\]Every continuous-time question becomes a discrete-time one, and the categorical draw a simulated jump needs is the one from Chains. The module also gives the other reading of a generator (exponential holding times and an embedded jump chain), the stationary distribution of a birth–death chain, closed forms for the two-state telegraph process and the compound Poisson process, and one step of Gillespie’s algorithm.
Conventions. A generator q is n×n with non-negative off-diagonal rates \(q_{ij}\) and rows summing to zero, so \(q_{ii} = -\sum_{j \ne i} q_{ij}\). Distributions are vectors n. Uniform draws in \([0, 1)\) are passed in by the caller. The Generator, BirthDeath, Telegraph, Compound and Gillespie sample models declare processes over these functions.
Limitations. transientOf carries a distribution rather than a matrix, so it answers \(\pi(t)\) from one starting distribution, not the whole matrix \(P(t)\). Its series is truncated after the number of turns you give; choose it well above \(\lambda t\) (for example \(\lambda t + 6\sqrt{\lambda t} + 10\)). Gillespie propensities are linear in the species counts, which is exact for first-order reactions such as decay, migration, birth and death, but cannot express mass action between two distinct reactants (a term like \(\beta S I\)).
In the statements below, \(\operatorname{nonzero}(v) = v + [\,\lvert v \rvert < 10^{-300}\,] \cdot 10^{-300}\) (applied entrywise) is the divisor guard from Chains; it moves only entries of magnitude below \(10^{-300}\), in practice exact zeros.
open Prelude.Jumps
Uniformisation
The Poisson weights are computed in logarithms, \(\exp(-\lambda t + k \ln(\lambda t) - \ln\Gamma(k+1))\), because for large \(\lambda t\) the factors \(e^{-\lambda t}\) and \((\lambda t)^k\) each overflow or underflow on their own. The same generator can also be read as a chain that waits an exponential time with mean \(1/(-q_{ii})\) in state \(i\) and then jumps to \(j\) with probability \(q_{ij}/(-q_{ii})\).
exitRate
exitRate(q) · q generator n×n → λ scalar
Let \(Q \in \mathbb{R}^{n \times n}\) (the argument q) be a generator. Then
The fastest exit rate, the smallest valid uniformisation rate. Any larger value also works; it only makes the series slower.
Intuition In state \(i\) the chain leaves at total rate \(-Q_{ii}\), and the largest of these is the busiest clock in the system. Uniformisation runs one Poisson clock at that rate for every state. A slower state lets some ticks pass without moving, and those are the self-loops of the jump matrix.
Use it when you need the rate for jumper, or want to size the series: about \(\lambda t\) ticks happen by time \(t\), so transientOf needs turns well above that. One very fast state makes \(\lambda\) large and the series long for every state. transientOf already calls it for you.
transient-analysissimulation
def exitRate(q) ≔ maximum(0 − diag(q))
jumper
jumper(q, l) · q generator n×n · l uniformisation rate λ scalar → jump matrix n×n
Let \(Q \in \mathbb{R}^{n \times n}\) (the argument q) be a generator and \(\lambda \ge \max_i(-Q_{ii})\) (the argument l) a uniformisation rate. Then
The uniformised jump matrix, which is row-stochastic when \(\lambda\) is at least the fastest exit rate.
Intuition Slice time by the ticks of one rate-\(\lambda\) Poisson clock. At each tick the chain moves from \(i\) to \(j\) with probability \(q_{ij}/\lambda\) and otherwise stays put, which is the diagonal \(1 + q_{ii}/\lambda\). Watched at those Poisson ticks, this ordinary discrete chain behaves exactly like the continuous one.
Use it when you want discrete-chain tools (matrix powers, a categorical draw per tick) on a continuous-time model. If \(\lambda\) is below the fastest exit rate a diagonal entry goes negative and \(J\) is no longer stochastic, so pass exitRate(q) or anything larger. If you want only real moves between states, with no self-loops, use jumpChain.
transient-analysissimulation
def jumper(q, l) ≔ identity(q) + q / nonzero(l)
transientOf
transientOf(start, q, t, turns) · start initial distribution vector n · q generator n×n · t time scalar · turns number of series terms scalar → distribution at time t vector n
Let \(\pi_0 \in \Delta^{n-1}\) (the argument start) be the initial distribution, \(Q \in \mathbb{R}^{n \times n}\) (the argument q) a generator, \(t \ge 0\) a time and \(K \in \mathbb{N}\) (the argument turns) the number of turns. Let \(\lambda = \operatorname{exitRate}(Q)\), \(J = \operatorname{jumper}(Q, \lambda)\) and \(w_k = \operatorname{massAt}(\lambda, t, k)\). Let \(v_0 = \pi_0\), \(a_0 = w_0\, \pi_0\), and for \(k = 0, 1, \dots, K-1\)
Then
\[\operatorname{transientOf}(\pi_0, Q, t, K) = a_K = \sum_{k=0}^{K} w_k \, \pi_0^{\top} J^{k} \;\approx\; \sum_{k=0}^{K} e^{-\lambda t} \frac{(\lambda t)^k}{k!} \, \pi_0^{\top} J^{k} .\]The distribution of the chain at time t, \(\pi(0)P(t)\), from the uniformisation series truncated after turns jumps with \(\lambda\) equal to the fastest exit rate.
Intuition \(\pi(t)\) is an average over “how many ticks happened by \(t\)”: the start distribution pushed \(k\) steps through \(J\), weighted by the Poisson chance of exactly \(k\) ticks. Those weights peak near \(\lambda t\) and vanish a few \(\sqrt{\lambda t}\) beyond it, so a finite number of turns holds almost all the mass.
Use it when you need state probabilities at a given time from a known start: the chance a machine is down at \(t\), a customer’s tier after 30 days, a chain’s behaviour before it settles. Too few turns leaves mass missing, so the result sums to less than 1 (a cheap check). Cost grows with \(\lambda t\). It gives one starting distribution per call. For the long-run answer of a birth–death chain use birthDeathStationary.
transient-analysisreliability-analysis
def transientOf(start, q, t, turns) ≔ seriesLoop(turns, start, massAt(exitRate(q), t, 0) · start, jumper(q, exitRate(q)), exitRate(q), t, 0)
holding
holding(q) · q generator n×n → mean holding times vector n
Let \(Q \in \mathbb{R}^{n \times n}\) (the argument q) be a generator. Then
The mean time spent in each state per visit. An absorbing state gives an effectively infinite value (\(10^{300}\)).
Intuition A stay in state \(i\) ends when the first of its exits fires. Competing exponential exits add their rates, so the stay is exponential with rate \(-Q_{ii}\) and its mean is the reciprocal: many ways out make a short visit.
Use it when you want a generator in human units, such as mean repair time, mean time in a customer tier or how long a molecule stays bound, or you are building an event-by-event simulation together with jumpChain. An absorbing state (a zero row) returns \(10^{300}\), so test for it rather than reading it as a real duration. This is the length of one visit, not the long-run share of time spent in the state.
reliability-analysiswaiting-timessimulation
def holding(q) ≔ 1 / nonzero(0 − diag(q))
jumpChain
jumpChain(q) · q generator n×n → embedded jump chain n×n
Let \(Q \in \mathbb{R}^{n \times n}\) (the argument q) be a generator. Then
The embedded discrete chain: each off-diagonal rate divided by its row’s exit rate, with zeros on the diagonal. Each row is scaled by its own exit rate, so the rows sum to one; an absorbing state gives a row of zeros.
Intuition The other reading of a generator: forget how long the chain waits and record only where it goes next. From \(i\) the exits race at rates \(q_{ij}\), and exit \(j\) wins with probability \(q_{ij}\) over the total exit rate. Holding times from holding plus this chain are exactly how a continuous-time chain runs.
Use it when you simulate event by event (an exponential stay, then the next state from this row), or ask questions about order rather than timing, such as which failure mode happens first or where the chain is finally absorbed. Unlike jumper it has no self-loops, so its powers count real transitions, not clock ticks. Absorbing states come out as zero rows, so put a 1 on their diagonal if a tool expects a stochastic matrix.
simulationreliability-analysisabsorption-analysis
def jumpChain(q) ≔ einsum("ij,i->ij", q − einsum("i,ij->ij", diag(q), identity(q)), holding(q))
Building blocks
massAt(l, t, k)— Let \(\lambda \ge 0\) (the argumentl), \(t \ge 0\) and \(k \in \mathbb{N}\). Then \(\operatorname{massAt}(\lambda, t, k) = \exp\big(-\lambda t + k \ln \operatorname{nonzero}(\lambda t) - \operatorname{lgamma}(k + 1)\big) \approx e^{-\lambda t} (\lambda t)^k / k!\). The Poisson probability of k jumps in time t at rate l, computed in logs.pushed(v, j)— Let \(v \in \mathbb{R}^n\) and \(J \in \mathbb{R}^{n \times n}\). Then \(\operatorname{pushed}(v, J)_t = \sum_s v_s J_{st} = (v^{\top} J)_t\). A distribution moved one jump through the jump matrix.seriesLoop(k, v, acc, j, l, t, i)— Let \(v, a \in \mathbb{R}^n\) (the argumentsvandacc), \(J \in \mathbb{R}^{n \times n}\), \(\lambda, t \ge 0\) and \(k, i \in \mathbb{N}\). Then \(\operatorname{seriesLoop}(0, v, a, J, \lambda, t, i) = a\) and \(\operatorname{seriesLoop}(k, v, a, J, \lambda, t, i) = \operatorname{seriesLoop}\big(k - 1,\, v^{\top} J,\, a + \operatorname{massAt}(\lambda, t, i + 1)\, v^{\top} J,\, J, \lambda, t,\, i + 1\big)\). Adds the next k weighted terms of the uniformisation series; used bytransientOf.
Birth and death
A birth–death chain on states \(0, 1, \dots, n-1\) moves up from \(i-1\) to \(i\) at rate \(\lambda_{i-1}\) and down from \(i\) to \(i-1\) at rate \(\mu_i\). Its stationary distribution is a product of ratios,
\[\pi_i \propto \prod_{j \le i} \frac{\lambda_{j-1}}{\mu_j}\]evaluated as a cumulative sum of logarithms so that long chains do not overflow.
The birth and death vectors are aligned by the state being entered: entry \(i\) of birth is the rate into state \(i\) from below (\(\lambda_{i-1}\)) and entry \(i\) of death is the rate out of state \(i\) downward (\(\mu_i\)). Entry 0 of each is a placeholder: its ratio multiplies every term and cancels in the normalisation, so 1 and 1 is the natural choice.
ratios
ratios(birth, death) · birth up-rates vector n · death down-rates vector n → ratios vector n
Let \(b = (b_0, \dots, b_{n-1})\) (the argument birth, with \(b_i = \lambda_{i-1}\)) and \(d = (d_0, \dots, d_{n-1})\) (the argument death, with \(d_i = \mu_i\)) be non-negative rate vectors, positions counted from zero. Then
The elementwise ratios \(\lambda_{i-1}/\mu_i\); for a queue these are the traffic intensities at each level.
Intuition A birth–death chain only moves between neighbours, so in equilibrium the flow up from level \(i-1\) equals the flow down from level \(i\): \(\pi_{i-1}\lambda_{i-1} = \pi_i \mu_i\). Each ratio is the factor by which occupancy rises (above 1) or falls (below 1) from one level to the next.
Use it when you are inspecting a queue, population or stock-level model before solving it. Ratios that stay below 1 at high levels mean the chain is pulled back down. Ratios of 1 or more mean mass piles up at the top of the states you gave. For a single server with constant rates every ratio is \(\rho\), which offered in Queues gives directly.
congestion-analysispopulation-dynamics
def ratios(birth, death) ≔ birth / nonzero(death)
birthDeathStationary
birthDeathStationary(birth, death) · birth up-rates vector n · death down-rates vector n → stationary distribution vector n
Let \(b, d \in \mathbb{R}^n\) (the arguments birth and death) be positive rate vectors aligned as above, \(\rho = \operatorname{ratios}(b, d)\), and
Then
\[\operatorname{birthDeathStationary}(b, d)_i = \frac{w_i}{\operatorname{nonzero}\big(\sum_{k=0}^{n-1} w_k\big)} .\]The stationary distribution of the finite birth–death chain, normalised to sum to one.
Intuition Balancing each pair of neighbours chains the ratios together: level \(i\) weighs the product of the ratios up to \(i\), and normalising turns the weights into the long-run share of time at each level. Sums of logarithms keep long products from overflowing.
Use it when you need the equilibrium of a finite level process whose rates change with the level: a waiting room where customers balk as it fills, busy servers, stock on hand, a population with a carrying capacity. The closed-form queue formulas cannot handle rates like that. Set entry 0 of both vectors to 1 and keep the down-rates positive. A zero up-rate simply makes the levels above it unreachable. For a truncated infinite chain check settles, and for probabilities at a finite time use transientOf.
capacity-planningcongestion-analysispopulation-dynamics
def birthDeathStationary(birth, death) ≔ unnormalised(birth, death) / nonzero(total(unnormalised(birth, death)))
settles
settles(birth, death) · birth up-rates vector n · death down-rates vector n → flag 1 or 0
Let \(b, d \in \mathbb{R}^n\) (the arguments birth and death) be positive rate vectors and \(w = \operatorname{unnormalised}(b, d)\). Then
1 when the unnormalised weights stay below \(10^{12}\) in total, a practical sign that the ratios fall below one and the chain settles; 0 when they grow without bound over the given states.
Intuition An unbounded birth–death chain has an equilibrium only if the running products of the ratios add up to something finite, which needs the ratios to fall below 1 eventually. Over a finite list of states the sum is always finite, so this checks a practical stand-in: do the weights stay under \(10^{12}\), or blow up?
Use it when you cut an unbounded queue or population off at \(n\) levels and want to know whether the cut is harmless. A 0 means mass piles up at the top: the answer depends on where you cut, and the real system has no steady state. A 1 is only a heuristic, since ratios above 1 over a short list can still stay under the threshold. Add levels and check that birthDeathStationary does not change.
stability-checkcapacity-planning
def settles(birth, death) ≔ total(unnormalised(birth, death)) < 1e12
Building blocks
prefixLogs(v)— Let \(v = (v_0, \dots, v_{n-1}) \in \mathbb{R}^n\) with positive entries, positions counted from zero. Then \(\operatorname{prefixLogs}(v)_i = \sum_{j=0}^{n-1} [\, i \ge j \,] \ln \operatorname{nonzero}(v_j) = \sum_{j=0}^{i} \ln \operatorname{nonzero}(v_j)\). The cumulative sum of ln v, the logarithm of a running product.unnormalised(birth, death)— Let \(b, d \in \mathbb{R}^n\) and \(\rho = \operatorname{ratios}(b, d)\). Then \(\operatorname{unnormalised}(b, d)_i = \exp\big(\operatorname{prefixLogs}(\rho)_i\big) = \prod_{j=0}^{i} \operatorname{nonzero}(\rho_j)\). The running products of the ratios, before normalisation.
The telegraph process
A two-state process switches from off to on at rate \(\alpha\) (up) and from on to off at rate \(\beta\) (down). In the long run it is on with probability \(p = \alpha/(\alpha + \beta)\), and its autocovariance decays at the sum of the two rates, however unequal they are:
telegraphUp
telegraphUp(up, down) · up off-to-on rate α · down on-to-off rate β → probability scalar
Let \(\alpha \ge 0\) (the argument up) and \(\beta \ge 0\) (the argument down) be the switching rates, with \(\alpha + \beta > 0\). Then
The long-run probability of being on.
Intuition A switch that turns on at rate \(\alpha\) and off at rate \(\beta\) spends on average \(1/\alpha\) off and \(1/\beta\) on per cycle. The share of time on is \((1/\beta)/(1/\alpha + 1/\beta)\), which simplifies to \(\alpha/(\alpha+\beta)\).
Use it when something alternates between two states, such as a machine up or down, a gene on or off, a channel busy or idle, or an on-off traffic source, and you need its availability or duty cycle. Pass rates per unit time, not probabilities per step. For more than two states write a generator and use transientOf, or birthDeathStationary when the states form a line.
reliability-analysisintermittent-signals
def telegraphUp(up, down) ≔ up / nonzero(up + down)
telegraphSpread
telegraphSpread(up, down) · up α · down β → variance scalar
Let \(\alpha, \beta \ge 0\) (the arguments up and down) and \(p = \operatorname{telegraphUp}(\alpha, \beta)\). Then
The stationary variance of the 0/1 state.
Intuition At a random moment the switch is a 0/1 coin that lands on 1 with probability \(p\), so its variance is a Bernoulli variance. It peaks at \(1/4\) when the switch is on half the time and nears zero when it is almost always on or almost always off.
Use it when you need the size of the fluctuations of an on/off signal around its average, such as the variability of load from an intermittent source or the noise from a gene switching. If “on” means a level \(h\) rather than 1, multiply by \(h^2\). For how fast those fluctuations are forgotten, use telegraphDecay.
intermittent-signalsuncertainty-quantification
def telegraphSpread(up, down) ≔ telegraphUp(up, down) · (1 − telegraphUp(up, down))
telegraphDecay
telegraphDecay(up, down, t) · up α · down β · t lag scalar → autocovariance scalar
Let \(\alpha, \beta \ge 0\) (the arguments up and down), \(p = \operatorname{telegraphUp}(\alpha, \beta)\) and \(t \ge 0\) a lag. Then
The stationary autocovariance at lag t.
Intuition How much the switch’s state now tells you about its state a time \(t\) later. Either transition erases that information, so the memory fades at the combined rate \(\alpha + \beta\): a switch that rarely turns on but turns off quickly is still forgotten quickly.
Use it when you need the correlation time \(1/(\alpha+\beta)\) of an intermittent signal: to decide how far apart samples must be to count as roughly independent, or to fit the two rates from a measured autocovariance (the decay rate gives \(\alpha + \beta\), the value at lag 0 gives \(p(1-p)\)). At \(t = 0\) it equals telegraphSpread.
intermittent-signalsparameter-estimation
def telegraphDecay(up, down, t) ≔ telegraphSpread(up, down) · exp(0 − (up + down) · t)
The compound Poisson process
The total of a Poisson number of independent jumps, with jump rate \(\lambda\), jump mean \(\mu\) and jump standard deviation \(\sigma\). By Wald’s identities
\[E[S(t)] = \lambda t \mu, \qquad \operatorname{Var} S(t) = \lambda t (\sigma^2 + \mu^2)\]The \(\mu^2\) term in the variance, from the randomness in the number of jumps, is the one a first guess tends to leave out.
compoundMean
compoundMean(rate, t, jumpMean) · rate jump rate λ · t time · jumpMean mean jump size μ → expected total scalar
Let \(\lambda \ge 0\) (the argument rate) be the jump rate, \(t \ge 0\) a time and \(\mu \in \mathbb{R}\) (the argument jumpMean) the mean jump size. Then
The expected total.
Intuition On average \(\lambda t\) events arrive by time \(t\) and each adds \(\mu\), so the expected total is the product. Neither the randomness of the count nor the spread of the sizes moves the mean.
Use it when you need expected aggregate amounts: insurance claims over a policy year, total demand from randomly arriving orders, rainfall from storms. It is the starting point for a premium or a budget and needs only the mean size, not its distribution. It says nothing about risk, so pair it with compoundSpread.
insurance-claimsrisk-measurement
def compoundMean(rate, t, jumpMean) ≔ rate · t · jumpMean
compoundSpread
compoundSpread(rate, t, jumpMean, jumpSpread) · rate jump rate λ · t time · jumpMean mean jump size μ · jumpSpread jump standard deviation σ → variance of the total scalar
Let \(\lambda \ge 0\) (the argument rate), \(t \ge 0\), \(\mu \in \mathbb{R}\) (the argument jumpMean) and \(\sigma \ge 0\) (the argument jumpSpread) the standard deviation of a jump. Then
The variance of the total.
Intuition The total varies for two reasons: each jump’s size is random (\(\sigma^2\)), and so is the number of jumps. A Poisson count has variance equal to its mean, which adds \(\lambda t\,\mu^2\). Together that is \(\lambda t\) times the mean squared jump.
Use it when you set reserves, capital or safety stock for aggregate losses or demand. With many expected events, the mean plus \(z\) standard deviations (the square root of this) is a fair approximation. With few events or heavy-tailed sizes that normal band understates the tail, so simulate instead. jumpSpread is a standard deviation, not a variance.
insurance-claimsrisk-measurementuncertainty-quantification
def compoundSpread(rate, t, jumpMean, jumpSpread) ≔ rate · t · (jumpSpread · jumpSpread + jumpMean · jumpMean)
Gillespie’s algorithm
Each step of the stochastic simulation algorithm needs two draws: when the next reaction happens, and which one. With propensities \(a_r\) and total \(a_0 = \sum_r a_r\), the waiting time is exponential with rate \(a_0\) and reaction \(r\) is chosen with probability \(a_r/a_0\).
Shapes: with \(R\) reactions and \(S\) species, rates is a vector R, reactants an R×S matrix giving each reaction’s reactant abundance as a linear combination of the counts, stoichiometry an R×S matrix of the change each reaction makes to each species, and counts a vector S.
propensities
propensities(rates, reactants, counts) · rates rate constants vector R · reactants reactant map R×S · counts species counts vector S → propensities vector R
Let \(\kappa \in \mathbb{R}^R\) (the argument rates) be the rate constants, \(A \in \mathbb{R}^{R \times S}\) (the argument reactants) the reactant map, and \(c \in \mathbb{R}^S\) (the argument counts) the species counts. Then
Each reaction’s propensity, its rate constant times its reactant abundance.
Intuition A propensity is a reaction’s instantaneous firing rate: its rate constant times how much reactant is around. Ten molecules that can each decay make decay ten times as likely as one molecule does.
Use it when you set up a stochastic simulation of first-order systems such as radioactive or chemical decay, migration between compartments, cell birth and death, or mRNA degradation. On its own it also shows which reaction dominates in the current state. Reactions between two different species (an \(S \cdot I\) infection term, \(A + B \to C\)) are not linear in the counts and cannot be written with this reactant map.
reaction-networkspopulation-dynamicssimulation
def propensities(rates, reactants, counts) ≔ rates · einsum("rs,s->r", reactants, counts)
totalPropensity
totalPropensity(rates, reactants, counts) · rates vector R · reactants R×S · counts vector S → total propensity scalar
Let \(\kappa \in \mathbb{R}^R\), \(A \in \mathbb{R}^{R \times S}\) and \(c \in \mathbb{R}^S\) be as in propensities, and \(a = \operatorname{propensities}(\kappa, A, c)\). Then
The sum of all propensities, the rate at which some reaction occurs.
Intuition With several independent exponential clocks running, the first one rings at the sum of their rates. \(a_0\) measures how busy the whole system is: a large total means short gaps between events.
Use it when you need the rate for waitFor in a Gillespie step, or want to estimate simulation cost: over a stretch where \(a_0\) stays roughly constant, expect about \(a_0\) events per unit of simulated time. A total of zero means no reaction can fire, so the system is frozen and the simulation loop should stop.
reaction-networkssimulation
def totalPropensity(rates, reactants, counts) ≔ total(propensities(rates, reactants, counts))
waitFor
waitFor(a, u) · a rate scalar · u uniform draw scalar → waiting time scalar
Let \(a > 0\) be a rate and \(u \in [0, 1)\) a uniform draw. Then
\[\operatorname{waitFor}(a, u) = -\frac{\ln \operatorname{nonzero}(1 - u)}{\operatorname{nonzero}(a)} .\]An exponential waiting time with rate a.
Intuition Inverse-transform sampling. The chance of no event for a time \(w\) is \(e^{-aw}\). Set that equal to \(1 - u\) and solve for \(w\). Draws near 1 map to the rare long quiet stretches, and draws near 0 to events that follow almost at once.
Use it when you need the time to the next event in any Poisson-driven simulation: a Gillespie step (rate = totalPropensity), customer arrivals, machine failures, claim arrivals. Use a different uniform from the one that picks the reaction, or the timing and the choice become linked. A rate of zero returns an astronomically long wait rather than an error, so check for it.
simulationwaiting-times
def waitFor(a, u) ≔ 0 − ln(nonzero(1 − u)) / nonzero(a)
reactionAt
reactionAt(rates, reactants, counts, u) · rates vector R · reactants R×S · counts vector S · u uniform draw scalar → chosen reaction one-hot vector R
Let \(\kappa \in \mathbb{R}^R\), \(A \in \mathbb{R}^{R \times S}\) and \(c \in \mathbb{R}^S\) be as in propensities, \(a = \operatorname{propensities}(\kappa, A, c)\), \(a_0 = \operatorname{totalPropensity}(\kappa, A, c)\), and \(u \in [0, 1)\) a uniform draw. Then
the one-hot vector with entry \(r\) (positions \(0, \dots, R-1\)) equal to \([\, C_r > u \,]\,[\, C_r - a_r/\operatorname{nonzero}(a_0) \le u \,]\), where \(C_r = \sum_{r' \le r} a_{r'} / \operatorname{nonzero}(a_0)\).
The reaction that fires, drawn in proportion to the propensities.
Intuition Lay the reactions end to end along \([0, 1)\), each taking a slice as wide as its share of the total propensity. The uniform draw lands in one slice and that reaction fires, so faster reactions own wider slices and fire more often.
Use it when you pick which event happens next in a Gillespie step or any competing-risks simulation, such as which of several failure modes strikes or which kind of claim arrives. The result is one-hot, so it can be multiplied straight into a stoichiometry matrix, as afterReaction does. If every propensity is zero the vector is all zeros and nothing fires.
simulationreaction-networks
def reactionAt(rates, reactants, counts, u) ≔ landedOn(propensities(rates, reactants, counts) / nonzero(totalPropensity(rates, reactants, counts)), u)
afterReaction
afterReaction(counts, stoichiometry, rates, reactants, u) · counts vector S · stoichiometry changes R×S · rates vector R · reactants R×S · u uniform draw scalar → new counts vector S
Let \(c \in \mathbb{R}^S\) (the argument counts), \(\nu \in \mathbb{R}^{R \times S}\) (the argument stoichiometry), \(\kappa \in \mathbb{R}^R\) (the argument rates), \(A \in \mathbb{R}^{R \times S}\) (the argument reactants), \(u \in [0, 1)\) a uniform draw, and \(e = \operatorname{reactionAt}(\kappa, A, c, u)\). Then
The species counts after the reaction chosen by u fires. Draw the time of the reaction separately with waitFor(totalPropensity(…), v).
Intuition Firing a reaction changes each species by a fixed amount, which is the reaction’s row of the stoichiometry matrix (\(-1\) A and \(+1\) B for \(A \to B\)). The one-hot choice picks out that row, and the row is added to the counts.
Use it when you advance the state of a stochastic simulation of decay chains, compartment models or linear birth–death systems. A full Gillespie step is two independent draws: waitFor(totalPropensity(…), v) moves the clock, and this function with u updates the counts. Repeat until the time horizon, or until the total propensity is zero. Counts are not clamped at zero, so make sure a species at zero gives its consuming reactions zero propensity.
simulationreaction-networkspopulation-dynamics
def afterReaction(counts, stoichiometry, rates, reactants, u) ≔ counts + einsum("r,rs->s", reactionAt(rates, reactants, counts, u), stoichiometry)
Examples
The distribution of a two-state chain at time 1
The chain leaves state 0 at rate 1 and state 1 at rate 2, and starts in state 0.
open Prelude.Jumps
fn distributionAt(start, q, t, turns) ↦ transientOf(start, q, t, turns)
fn meanHolding(q) ↦ holding(q)
fn embedded(q) ↦ jumpChain(q)
using Examples;
using MLambda.Hilbert.Runtime;
var q = (Tensor)new double[,] { { -1, 1 }, { 2, -2 } };
var atOne = Generators.DistributionAt((Tensor)new[] { 1d, 0d }, q, 1.0, 64.0);
var holds = Generators.MeanHolding(q);
var jumps = Generators.Embedded(q);
Output (verified, rounded to 6 decimal places):
transient = [0.683262, 0.316738]
holding = [1, 0.5]
jumpChain = shape [2, 2] values [0, 1, 1, 0]
For this chain \(p_{00}(t) = \tfrac23 + \tfrac13 e^{-3t}\), which at \(t = 1\) is \(0.683262\).
A queue with a finite waiting room
Arrivals at rate 1, service at rate 2, and room for at most three customers.
open Prelude.Jumps
fn queueLengths(birth, death) ↦ birthDeathStationary(birth, death)
using Examples;
using MLambda.Hilbert.Runtime;
var birth = (Tensor)new[] { 1d, 1, 1, 1 }; // entry 0 is a placeholder
var death = (Tensor)new[] { 1d, 2, 2, 2 };
var lengths = Generators.QueueLengths(birth, death);
Output (verified, rounded to 6 decimal places):
bd = [0.533333, 0.266667, 0.133333, 0.066667]
The weights are \(1, \tfrac12, \tfrac14, \tfrac18\), which sum to \(1.875\).
Moments of the telegraph and compound Poisson processes
open Prelude.Jumps
fn totalClaims(rate, t, jumpMean, jumpSpread) ↦ compoundSpread(rate, t, jumpMean, jumpSpread)
fn flipMemory(up, down, t) ↦ telegraphDecay(up, down, t)
using Examples;
var claimsVariance = Generators.TotalClaims(2.0, 3.0, 1.0, 0.5);
var covariance = Generators.FlipMemory(1.0, 3.0, 0.5);
Output (verified, rounded to 6 decimal places):
compoundSpread = 7.5
telegraph = 0.025375
Claims arriving at rate 2 over 3 years with mean 1 and spread 0.5 have variance \(6 \times (0.25 + 1)\). The telegraph process is on a quarter of the time, so its variance is \(0.1875\), and after half a unit it has decayed by \(e^{-2}\).
One Gillespie step
Species A decays to B at rate 0.5 per molecule, and B is removed at rate 0.1 per molecule, starting from 10 A and 5 B.
open Prelude.Jumps
fn nextWait(a, u) ↦ waitFor(a, u)
fn fired(rates, reactants, counts, u) ↦ reactionAt(rates, reactants, counts, u)
fn afterFiring(counts, stoichiometry, rates, reactants, u) ↦ afterReaction(counts, stoichiometry, rates, reactants, u)
using Examples;
using MLambda.Hilbert.Runtime;
var rates = (Tensor)new[] { 0.5, 0.1 };
var reactants = (Tensor)new double[,] { { 1, 0 }, { 0, 1 } };
var stoich = (Tensor)new double[,] { { -1, 1 }, { 0, -1 } };
var counts = (Tensor)new[] { 10d, 5d };
var wait = Generators.NextWait(5.5, 0.5);
var which = Generators.Fired(rates, reactants, counts, 0.95);
var after = Generators.AfterFiring(counts, stoich, rates, reactants, 0.95);
var other = Generators.AfterFiring(counts, stoich, rates, reactants, 0.5);
Output (verified, rounded to 6 decimal places):
wait = 0.126027
fired = [0, 1]
after = [10, 4]
after2 = [9, 6]
The propensities are 5 and 0.5, total 5.5, so a draw of 0.5 waits \(\ln 2 / 5.5\). The first reaction owns \([0, 0.909)\) of the unit interval: a draw of 0.95 removes a B, and a draw of 0.5 turns an A into a B.