Points
Point processes — Poisson arrivals, thinning, renewal gaps, Hawkes self-excitation, Cox overdispersion, Ripley's K and simple fits.
This module covers point processes, where the randomness is in when events happen rather than in how large they are. A Poisson process is built from exponential gaps added up; thinning keeps each arrival with some probability; an inhomogeneous process is thinning against the peak intensity; a renewal process replaces the exponential gaps with Weibull ones. Hawkes self-exciting processes, Cox (doubly stochastic) processes and Ripley’s K function for points in the plane complete the family, together with moment-based fits for an observed set of arrival times.
Conventions. Draws are passed in: u is a vector of uniforms in \([0, 1)\), and an exponential gap is produced as \(-\ln(1 - u)/\lambda\) (using \(1 - u\) so a draw of exactly zero cannot give an infinite gap). times is a vector of arrival times in increasing order, horizon the length of the observation window, and rate an intensity per unit time. The Hawkes kernel is exponential: each past arrival adds \(\alpha e^{-\beta(t - t_i)}\) to the intensity, where excite is \(\alpha\) and decay is \(\beta\). The Poisson, Inhomogeneous, Renewal, Hawkes, Doubly and Spatial sample models declare processes over these functions.
Limitations. Cumulative sums, the Hawkes intensities at every arrival and Ripley’s K all build n×n matrices, so they suit hundreds to a few thousand events. A Hawkes process cannot be generated here, because that needs thinning against a bound that changes as arrivals occur; its intensity and likelihood can be evaluated on given times, and its parameters are fitted by moments rather than by maximum likelihood.
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.Points
Interarrivals, and the process they make
A Poisson process of rate \(\lambda\) has independent exponential gaps; the arrival times are their running totals, and the count by time \(t\) has mean and variance both equal to \(\lambda t\).
gapsOf
gapsOf(u, rate) · u uniform draws vector n · rate intensity λ scalar → gaps vector n
Let \(u \in [0, 1)^n\) be independent uniform draws and let \(\lambda > 0\) (the argument rate) be the intensity. Then
Exponential interarrival gaps \(-\ln(1 - u_i)/\lambda\), each exponential with rate \(\lambda\).
Intuition Waiting for an event that has no memory: however long you have waited, the rest of the wait looks the same. Such gaps are exponential with mean \(1/\lambda\), mostly short with the occasional long silence. Inverting the exponential distribution turns each uniform draw into one gap.
Use it when events happen independently at a steady rate: customers arriving, calls to a support line, requests to a server, radioactive decays, failures with a constant hazard. For gaps that wear out (more regular) or come in bursts, use weibullGap. For absolute timestamps, use timesOf.
event-arrivalssimulationreliability-analysis
def gapsOf(u, rate) ≔ 0 − ln(nonzero(1 − u)) / nonzero(rate)
timesOf
timesOf(u, rate) · u uniform draws vector n · rate intensity λ scalar → arrival times vector n
Let \(u = (u_0, \dots, u_{n-1}) \in [0, 1)^n\) be independent uniform draws (positions counted from zero), \(\lambda > 0\) the intensity, and \(g = \operatorname{gapsOf}(u, \lambda)\). Then
\[\operatorname{timesOf}(u, \lambda)_k = \sum_{j=0}^{k} g_j , \qquad k = 0, \dots, n-1 .\]The arrival times of a Poisson process: the running totals of gapsOf(u, rate).
Intuition Lay the exponential gaps end to end on a timeline; the tick marks are the arrival times of a Poisson process.
Use it when you generate synthetic event logs: arrivals for a queue or discrete-event simulation, load-test traffic, test data for an anomaly detector. You get exactly one arrival per draw, so you choose the number of events, not the time span; to fill a window \([0, T]\), supply comfortably more than \(\lambda T\) draws and drop times beyond \(T\). The running total builds an n×n mask, which limits it to a few thousand arrivals.
event-arrivalssimulation
def timesOf(u, rate) ≔ runningTotal(gapsOf(u, rate))
countedBy
countedBy(u, rate, t) · u uniform draws vector n · rate intensity λ scalar · t time scalar → count scalar
Let \(u \in [0, 1)^n\) be independent uniform draws, \(\lambda > 0\) the intensity, \(t \ge 0\) a time, and \(\tau = \operatorname{timesOf}(u, \lambda)\). Then
\[\operatorname{countedBy}(u, \lambda, t) = N(t) = \sum_{k=0}^{n-1} [\, \tau_k \le t \,] .\]The number of arrivals at or before time t. The count cannot exceed the number of draws, so supply comfortably more than \(\lambda t\) of them.
Intuition Count the tick marks on the timeline up to time \(t\). Across many streams of draws that count is Poisson distributed with mean \(\lambda t\).
Use it when you simulate the number of events in a window: orders per hour, defects per shift, insurance claims per year, and compare the spread over many runs with poissonMean and poissonSpread. If the draws run out before time \(t\), the count stops at \(n\) and is biased low, so supply well over \(\lambda t\) draws, for example \(\lambda t + 5\sqrt{\lambda t}\).
event-countingsimulation
def countedBy(u, rate, t) ≔ einsum("i->", timesOf(u, rate) ≤ t)
poissonMean
poissonMean(rate, t) · rate intensity λ · t time → expected count scalar
Let \(\lambda > 0\) (the argument rate) be the intensity and \(t \ge 0\) a time. Then
The expected number of arrivals by time \(t\).
Intuition A steady rate adds up linearly: twice the time, twice the expected events.
Use it when you plan capacity (expected arrivals in a shift, failures over a warranty period, claims in a year) or need the baseline an anomaly detector compares a count against. The number is the same as coxMean, because a random rate does not change the average; what differs is the spread, so pair this with poissonSpread or coxSpread.
event-countingforecastingcapacity-planning
def poissonMean(rate, t) ≔ rate · t
poissonSpread
poissonSpread(rate, t) · rate intensity λ · t time → count variance scalar
Let \(\lambda > 0\) (the argument rate) be the intensity and \(t \ge 0\) a time. Then
The variance of the count by time \(t\), also \(\lambda t\). Equal mean and variance is the signature of a Poisson process.
Intuition For Poisson counts the variance equals the mean, so a typical fluctuation is \(\sqrt{\lambda t}\) and the relative noise shrinks like \(1/\sqrt{\lambda t}\): 100 expected events vary by about 10, 10,000 by about 100.
Use it when you set control limits for counts (\(\lambda t \pm 3\sqrt{\lambda t}\)), size a staffing buffer, or judge whether an observed count is unusual. It also tests the model: if the sample variance of real counts is well above their mean, the data are overdispersed and coxSpread or overdispersion fit better.
uncertainty-quantificationevent-countinganomaly-detection
def poissonSpread(rate, t) ≔ rate · t
Thinning, superposition, and a changing rate
Keeping each arrival of a rate-\(\lambda\) process independently with probability \(p\) gives a Poisson process of rate \(p\lambda\), and the thinned process is a subset of the original. The Lewis–Shedler method uses this to draw an inhomogeneous process with intensity \(\lambda(t) \le \lambda^*\): draw a homogeneous process at \(\lambda^*\) and keep an arrival at time \(t\) with probability \(\lambda(t)/\lambda^*\). The module’s example intensity rises and falls once over the horizon \(H\):
\[\lambda(t) = \lambda^* \, \frac{1 + \sin(2\pi t / H)}{2}, \qquad \Lambda(t) = \int_0^t \lambda = \frac{\lambda^* t}{2} + \frac{\lambda^* H \,(1 - \cos(2\pi t/H))}{4\pi}\]keptAt
keptAt(v, p) · v uniform draws vector n · p keep probability scalar → kept flags vector n of 1 or 0
Let \(v \in [0, 1)^n\) be independent uniform draws, one per arrival, and let \(0 \le p \le 1\) be the keep probability. Then
\[\operatorname{keptAt}(v, p)_i = [\, v_i < p \,] , \qquad i = 1, \dots, n .\]1 for each arrival that survives thinning, 0 for each that is removed.
Intuition Each arrival flips its own biased coin to survive. The survivors still form a Poisson process, at rate \(p\lambda\), and they are a subset of the original arrivals rather than a fresh, independent stream.
Use it when you split a stream: visitors who buy, calls that escalate, packets that are dropped, events a sensor detects. Use the flags to select times from timesOf. The v draws must be independent of the draws that made the times. For a keep probability that changes with time, use waveKept as the pattern.
event-arrivalssimulation
def keptAt(v, p) ≔ v < p
thinnedCount
thinnedCount(v, p) · v uniform draws vector n · p keep probability scalar → count scalar
Let \(v \in [0, 1)^n\) be independent uniform draws and \(0 \le p \le 1\) the keep probability. Then
\[\operatorname{thinnedCount}(v, p) = \sum_{i=1}^{n} [\, v_i < p \,] .\]The number of arrivals that survive thinning with probability p.
Intuition The number of coins that came up “keep”. Given \(n\) arrivals it is Binomial\((n, p)\); if \(n\) was itself a Poisson count, the survivors are Poisson with mean \(p\lambda t\).
Use it when you simulate how many events pass a filter: buyers among \(n\) visitors, incidents that get detected, leads that convert. It counts only; to know which arrivals survived, use keptAt.
event-countingsimulation
def thinnedCount(v, p) ≔ einsum("i->", keptAt(v, p))
waveIntensity
waveIntensity(peak, horizon, t) · peak maximum intensity λ* · horizon period H · t time scalar or vector → intensity
Let \(\lambda^* > 0\) (the argument peak) be the maximum intensity, \(H > 0\) (the argument horizon) the period, and \(t\) a time (or a vector of times, taken entrywise). Then
The wave-shaped intensity, which runs between 0 and peak.
Intuition A rate that swells and ebbs once per period: it peaks at peak a quarter of the way through, drops to zero three quarters of the way, and averages half the peak.
Use it when you need a ready-made time-varying demand to test with: a daily traffic cycle, seasonal sales, rush-hour requests. It is one fixed sine shape that reaches zero; real cycles usually keep a floor, so write your own \(\lambda(t) \le\) peak and thin against it the way waveKept does.
event-arrivalsseasonality
def waveIntensity(peak, horizon, t) ≔ peak · (1 + sin(6.283185307179586 · t / nonzero(horizon))) / 2
waveKept
waveKept(peak, horizon, t, v) · peak maximum intensity λ* · horizon period H · t arrival times · v uniform draws → kept flags 1 or 0
Let \(\lambda^* > 0\) (the argument peak), \(H > 0\) (the argument horizon), \(t\) the arrival times of a rate-\(\lambda^*\) process and \(v \in [0, 1)\) independent uniform draws of the same shape. Then, entrywise,
The Lewis–Shedler thinning step: keeps an arrival of the rate-peak process at time t when \(v < \lambda(t)/\lambda^*\).
Intuition Overproduce, then trim. Generate arrivals everywhere at the peak rate, and at each time keep one with probability \(\lambda(t)/\lambda^*\); quiet periods lose most of their arrivals and busy ones keep nearly all.
Use it when you simulate arrivals with a time-of-day or seasonal pattern: website hits, bus passengers, call-centre peaks. Pass times from timesOf(u, peak) and a separate vector of draws v. Over a whole period only half the candidate arrivals survive on average, so generate about twice as many as you need.
event-arrivalssimulationseasonality
def waveKept(peak, horizon, t, v) ≔ v < waveIntensity(peak, horizon, t) / nonzero(peak)
waveExpected
waveExpected(peak, horizon, t) · peak maximum intensity λ* · horizon period H · t time → expected count scalar
Let \(\lambda^* > 0\) (the argument peak), \(H > 0\) (the argument horizon) and \(t \ge 0\). Then
The expected number of arrivals of the wave-intensity process by time \(t\), the integral \(\Lambda(t)\) above.
Intuition The area under the rate curve is the expected number of events: a straight trend \(\lambda^* t/2\) from the average rate, plus a wobble from where in the cycle you stop.
Use it when you check a waveKept simulation against the exact expected count, or plan capacity for part of a cycle. The expected count in a window \([t_1, t_2]\) is the difference of two calls. Because the thinned process is still Poisson, this number is also the variance of the count.
event-countingforecastingmodel-validation
def waveExpected(peak, horizon, t) ≔ peak · t / 2 + peak · horizon · (1 − cos(6.283185307179586 · t / nonzero(horizon))) / 12.566370614359172
Renewal processes
A renewal process has independent, identically distributed gaps. Weibull gaps with shape \(k\) and scale \(s\) cover all three regimes with one parameter: \(k < 1\) is burstier than Poisson, \(k > 1\) more regular, and \(k = 1\) is exactly Poisson. The coefficient of variation of the gaps says which regime applies, the long-run rate is one over the mean gap whatever the distribution (the elementary renewal theorem), and the inspection paradox says that an observer arriving at a random moment lands, on average, in a longer gap than a typical one, and waits until the next arrival for the mean residual life
\[E[\text{residual}] = \mu \, \frac{1 + c^2}{2}, \qquad \mu = s\,\Gamma(1 + 1/k), \quad c = \text{coefficient of variation}\]In the statements below \(\tilde\Gamma = \operatorname{gammaOf}\) is the gamma function as the module computes it.
weibullGap
weibullGap(u, shape, scale) · u uniform draws vector n · shape k · scale s → gaps vector n
Let \(u \in [0, 1)^n\) be independent uniform draws, \(k > 0\) (the argument shape) and \(s > 0\) (the argument scale). Then
Weibull gaps by inversion.
Intuition An exponential gap with its time axis bent by the shape \(k\). With \(k > 1\) the hazard rises with age, so gaps cluster around a typical length (wear-out, regular service); with \(k < 1\) the hazard falls, giving many tiny gaps and rare long silences (bursts); \(k = 1\) is exactly exponential.
Use it when gaps are not memoryless: time between failures of ageing equipment (\(k > 1\)), bursty user activity or infant-mortality failures (\(k < 1\)). Turn the gaps into event times with runningTotal. For steady, memoryless arrivals gapsOf is the simpler choice.
reliability-analysisevent-arrivalssimulation
def weibullGap(u, shape, scale) ≔ scale · (0 − ln(nonzero(1 − u))) ^ (1 / nonzero(shape))
gammaOf
gammaOf(x) · x positive argument → Γ(x)
Let \(x > 0\). Then
\[\operatorname{gammaOf}(x) = \exp\big(\operatorname{lgamma}(x)\big) \approx \Gamma(x) ,\]where lgamma is the Lanczos approximation to \(\ln \Gamma\) in Distributions.
The gamma function.
Intuition The factorial stretched to real numbers: \(\Gamma(n) = (n-1)!\) for whole \(n\), and a smooth curve in between. Here it is computed as the exponential of the Lanczos log-gamma.
Use it when you need moments of distributions built on powers, as the Weibull functions on this page do. Positive arguments only. \(\Gamma(x)\) exceeds the largest double a little above \(x = 171\), so for likelihoods or large arguments stay on the log scale with lgamma from Distributions.
special-functionsdistribution-moments
def gammaOf(x) ≔ exp(lgamma(x))
weibullMean
weibullMean(shape, scale) · shape k · scale s → mean gap scalar
Let \(k > 0\) (the argument shape) and \(s > 0\) (the argument scale). Then
The mean Weibull gap \(s\,\Gamma(1 + 1/k)\).
Intuition The scale sets the unit of time, and the shape adds the factor \(\Gamma(1 + 1/k)\). For \(k \ge 1\) that factor stays between about 0.89 and 1, so the mean gap is close to the scale; as \(k\) drops below 1 it grows fast (2 at \(k = \tfrac12\)), because rare long silences pull the mean up.
Use it when you turn fitted Weibull parameters into a mean time between failures or mean lifetime, schedule maintenance, or feed renewalRate and residualLife. For the spread around it use weibullSpread.
reliability-analysisdistribution-moments
def weibullMean(shape, scale) ≔ scale · gammaOf(1 + 1 / nonzero(shape))
weibullSpread
weibullSpread(shape, scale) · shape k · scale s → gap variance scalar
Let \(k > 0\) (the argument shape), \(s > 0\) (the argument scale), and write \(\kappa = \operatorname{nonzero}(k)\). Then
The variance of a Weibull gap.
Intuition How much individual gaps scatter around the mean, in squared time units. It scales with \(s^2\) and grows quickly as \(k\) falls below 1, when a few very long gaps dominate.
Use it when you need the uncertainty around a lifetime or a gap: safety stock for replacement parts, a band around a maintenance schedule, a check on simulated gaps. To compare regularity across systems on different time scales, use the scale-free weibullVariation instead.
reliability-analysisuncertainty-quantification
def weibullSpread(shape, scale) ≔ scale · scale · (gammaOf(1 + 2 / nonzero(shape)) − gammaOf(1 + 1 / nonzero(shape)) · gammaOf(1 + 1 / nonzero(shape)))
weibullVariation
weibullVariation(shape) · shape k → coefficient of variation scalar
Let \(k > 0\) (the argument shape) and \(\kappa = \operatorname{nonzero}(k)\). Then
The coefficient of variation of a Weibull gap. It does not depend on the scale, and is 1 at \(k = 1\).
Intuition Standard deviation over mean, a pure shape number with the time unit cancelled out. The exponential’s value of 1 is the Poisson yardstick: above 1 the gaps are burstier than random, below 1 more regular than random.
Use it when you classify the regime implied by a shape parameter or compare burstiness across systems that run on different time scales. It is the theoretical twin of gapVariation, which measures the same thing on observed gaps; matching the two is a quick way to choose \(k\).
clustering-of-eventsdistribution-moments
def weibullVariation(shape) ≔ sqrt(max(gammaOf(1 + 2 / nonzero(shape)) / nonzero(gammaOf(1 + 1 / nonzero(shape)) · gammaOf(1 + 1 / nonzero(shape))) − 1, 0))
burstierThan
burstierThan(shape) · shape k → flag 1 or 0
Let \(k > 0\) (the argument shape) and \(c = \operatorname{weibullVariation}(k)\). Then
1 when the gaps are burstier than Poisson (coefficient of variation above one). A tolerance of \(10^{-9}\) keeps shape exactly 1, which is Poisson, from being reported as bursty because of rounding.
Intuition A yes/no on whether the gaps vary more than Poisson gaps would, which for Weibull gaps means a shape below 1. The tiny tolerance keeps rounding from calling the Poisson case bursty.
Use it when you flag a workload as bursty in a report or a rule, since bursty arrivals build longer queues than Poisson-based sizing predicts and need extra headroom. It only tests the bursty side: regular and Poisson shapes both return 0. For a three-way answer from observed data use gapRegime.
clustering-of-eventscapacity-planning
def burstierThan(shape) ≔ weibullVariation(shape) > 1 + 1e-9
renewalRate
renewalRate(shape, scale) · shape k · scale s → long-run rate scalar
Let \(k > 0\) (the argument shape), \(s > 0\) (the argument scale) and \(\mu = \operatorname{weibullMean}(k, s)\). Then
The long-run arrival rate, one over the mean gap.
Intuition Over a long stretch of time, events per unit time are one over the mean gap, however regular or bursty the gaps are. Spread changes the timing of events, not how many arrive in the long run.
Use it when you budget long-run replacements or spare parts, estimate the throughput of a repeating cycle, or approximate the count over a long horizon \(T\) as rate × \(T\). Over short windows the count can differ a lot from this, especially for bursty gaps; simulate with weibullGap if the short run matters.
long-run-behaviourreliability-analysisforecasting
def renewalRate(shape, scale) ≔ 1 / nonzero(weibullMean(shape, scale))
residualLife
residualLife(shape, scale) · shape k · scale s → expected residual life scalar
Let \(k > 0\) (the argument shape), \(s > 0\) (the argument scale), \(\mu = \operatorname{weibullMean}(k, s)\) and \(c = \operatorname{weibullVariation}(k)\). Then
The expected wait until the next arrival for an observer who arrives at a random moment (the mean residual life). The gap that observer lands in is on average twice as long, \(\mu(1 + c^2)\), which is longer than a typical gap. For a Poisson process the wait is a whole mean gap, not half of one.
Intuition The inspection paradox, or why the bus always seems late: arriving at a random moment, you are more likely to land in a long gap than a short one, simply because long gaps take up more of the timeline. The more the gaps vary, the worse it gets, through the \(c^2\) term.
Use it when you estimate the wait of someone who shows up at a random time: a passenger at a stop, an audit that checks at random for the next failure, a user opening an app between notifications. It shows why irregular service hurts waiting times even at the same average rate. It assumes the process has been running long enough that the arrival moment is truly random.
waiting-timesreliability-analysis
def residualLife(shape, scale) ≔ weibullMean(shape, scale) · (1 + weibullVariation(shape) · weibullVariation(shape)) / 2
Hawkes self-excitation
A Hawkes process has intensity
\[\lambda(t) = \mu + \sum_{t_i < t} \alpha\, e^{-\beta (t - t_i)}\]so every arrival raises the chance of more. The branching ratio \(\alpha/\beta\) is the expected number of direct offspring of one event; the process is stable when it is below one, and then fires on average at \(\mu / (1 - \alpha/\beta)\). The log-likelihood over \([0, T]\) is
\[\ell = \sum_i \ln \lambda(t_i) - \mu T - \frac{\alpha}{\beta} \sum_i \big(1 - e^{-\beta (T - t_i)}\big)\]excitedAt
excitedAt(times, background, excite, decay, t) · times past arrivals vector n · background μ · excite α · decay β · t time scalar → intensity scalar
Let \((t_1, \dots, t_n) \in \mathbb{R}^n\) (the argument times) be the arrival times, \(\mu \ge 0\) (the argument background), \(\alpha \ge 0\) (the argument excite), \(\beta > 0\) (the argument decay), and \(t \in \mathbb{R}\). Then
The intensity at time t, counting only arrivals strictly before t.
Intuition Events that breed events. There is a steady background rate \(\mu\), and every past event lights a small fire that adds \(\alpha\) to the rate and dies down at rate \(\beta\). Right after a burst the rate is high, and it drifts back to \(\mu\) as things go quiet, like aftershocks after an earthquake.
Use it when you want the current risk of the next event given the history: aftershocks, clustered trades in a market, retweet cascades, repeat burglaries nearby, bursts of server errors, for example to raise an alert while the rate is elevated. It evaluates one time and ignores arrivals exactly at t; for the rate at every arrival at once use hawkesRates.
clustering-of-eventsrisk-measurementevent-arrivals
def excitedAt(times, background, excite, decay, t) ≔ background + einsum("i->", (times < t) · excite · exp(0 − decay · (t − times)))
branchingRatio
branchingRatio(excite, decay) · excite α · decay β → branching ratio scalar
Let \(\alpha \ge 0\) (the argument excite) and \(\beta > 0\) (the argument decay). Then
The expected number of events directly triggered by one event, \(\alpha/\beta\).
Intuition One event’s fire adds \(\alpha e^{-\beta s}\) to the rate, and the total area under it is \(\alpha/\beta\): the expected number of direct “children”. This is the offspring mean of a branching process (offspringMean in Chains), and in a stable process it is also the fraction of events that were triggered rather than spontaneous.
Use it when you measure how self-reinforcing a stream is: how much market activity is reaction rather than news, how viral posts are, how contagious defaults are. Values near 1 mean a near-critical system in which a small background produces a large amount of activity. Check stability with hawkesStable.
clustering-of-eventsepidemic-modelling
def branchingRatio(excite, decay) ≔ excite / nonzero(decay)
hawkesStable
hawkesStable(excite, decay) · excite α · decay β → flag 1 or 0
Let \(\alpha \ge 0\) (the argument excite) and \(\beta > 0\) (the argument decay). Then
1 when the branching ratio is below one, so the process does not explode.
Intuition If each event has fewer than one child on average, every cascade eventually dies out and the long-run rate stays finite. At one or above, cascades can go on forever and the process explodes.
Use it when you validate fitted or hand-picked parameters before using hawkesRate (which is meaningless otherwise), constrain an optimiser, or flag runaway contagion in a monitored system.
model-validationclustering-of-events
def hawkesStable(excite, decay) ≔ branchingRatio(excite, decay) < 1
hawkesRate
hawkesRate(background, excite, decay) · background μ · excite α · decay β → long-run rate scalar
Let \(\mu \ge 0\) (the argument background), \(\alpha \ge 0\) (the argument excite) and \(\beta > 0\) (the argument decay), with \(n^\ast = \operatorname{branchingRatio}(\alpha, \beta) < 1\). Then
The long-run event rate of a stable process, \(\mu/(1 - \alpha/\beta)\). Meaningful only when hawkesStable is 1.
Intuition Each background event brings its children, grandchildren and so on: \(1 + m + m^2 + \cdots = 1/(1 - m)\) events in total. The long-run rate is the background rate multiplied by that cascade size, so at \(m = 0.9\) activity is ten times the background.
Use it when you plan capacity for a self-exciting stream or split activity into its sources: \(1 - \mu/\text{rate}\) equals \(m\), the share generated from within. Only use it when hawkesStable is 1; otherwise the value is negative or huge and means nothing.
long-run-behaviourevent-countingforecasting
def hawkesRate(background, excite, decay) ≔ background / nonzero(1 − branchingRatio(excite, decay))
hawkesCompensator
hawkesCompensator(times, background, excite, decay, horizon) · times arrivals vector n · background μ · excite α · decay β · horizon T → integrated intensity scalar
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times in \([0, T]\), \(T > 0\) (the argument horizon), \(\mu \ge 0\), \(\alpha \ge 0\) and \(\beta > 0\). Then
The integral of the intensity over \([0, T]\), in closed form.
Intuition The area under the intensity curve: the number of events the model expected over \([0, T]\) given the history it saw. Each arrival contributes up to \(\alpha/\beta\), less if it came so near the end that its fire had no time to burn out.
Use it when you compute the penalty term of hawkesLikelihood, or compare the expected count with the observed one as a quick fit check. Pass only arrivals at or before horizon: later times are not masked out and would add wrong, negative contributions.
likelihood-evaluationmodel-validation
def hawkesCompensator(times, background, excite, decay, horizon) ≔ background · horizon + einsum("i->", excite / nonzero(decay) · (1 − exp(0 − decay · (horizon − times))))
hawkesRates
hawkesRates(times, background, excite, decay) · times arrivals in increasing order vector n · background μ · excite α · decay β → intensities vector n
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times, \(\mu \ge 0\), \(\alpha \ge 0\), \(\beta > 0\), and \(K = \operatorname{hawkesKernel}(t_{1:n}, \alpha, \beta)\). Then
The intensity seen by each arrival, from the arrivals strictly before it, all at once.
Intuition The rate each event “saw” at the moment it happened. Events that arrive during a flurry see a high rate, which marks them as likely triggered; events in quiet times see nearly \(\mu\) and look spontaneous.
Use it when you evaluate the likelihood, or attribute events to causes: \(\mu/\lambda(t_i)\) is the probability that event \(i\) came from the background rather than from an earlier event (stochastic declustering, used for aftershocks and for separating news from reaction in markets). It builds an n×n matrix, so keep to a few thousand events.
clustering-of-eventslikelihood-evaluation
def hawkesRates(times, background, excite, decay) ≔ background + einsum("ij->i", hawkesKernel(times, excite, decay))
hawkesLikelihood
hawkesLikelihood(times, background, excite, decay, horizon) · times arrivals vector n · background μ · excite α · decay β · horizon T → log-likelihood scalar
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times in \([0, T]\), \(T > 0\) (the argument horizon), \(\mu \ge 0\), \(\alpha \ge 0\), \(\beta > 0\), \(\lambda = \operatorname{hawkesRates}(t_{1:n}, \mu, \alpha, \beta)\) and \(\Lambda = \operatorname{hawkesCompensator}(t_{1:n}, \mu, \alpha, \beta, T)\). Then
The log-likelihood of the observed arrivals over \([0, T]\). This is the quantity a full maximum-likelihood fit would maximise.
Intuition A good model puts a high rate where events actually happened (the \(\sum \ln \lambda\) term) without predicting many events overall (the compensator). Parameters that explain the bursts without inflating the expected count score best.
Use it when you compare parameter sets, refine the moment estimates of branchingFrom and exciteFrom with an optimiser, or test for self-excitation: with excite = 0 it is the plain Poisson log-likelihood \(n \ln \mu - \mu T\), so the difference is a likelihood-ratio statistic. Pass only arrivals inside \([0, T]\).
likelihood-evaluationparameter-estimation
def hawkesLikelihood(times, background, excite, decay, horizon) ≔ einsum("i->", ln(nonzero(hawkesRates(times, background, excite, decay)))) − hawkesCompensator(times, background, excite, decay, horizon)
Building blocks
hawkesKernel(times, excite, decay)— Let \((t_1, \dots, t_n) \in \mathbb{R}^n\), \(\alpha \ge 0\) and \(\beta > 0\). Then \(\operatorname{hawkesKernel}(t_{1:n}, \alpha, \beta)_{ij} = [\, t_j < t_i \,]\, \alpha\, e^{-\beta (t_i - t_j)}\) for \(i, j = 1, \dots, n\). The n×n matrix of excitation that arrival j contributes to arrival i, zero unless j came first.
Cox, the doubly stochastic process
A Cox process is a Poisson process whose rate is itself random, with mean \(\bar\lambda\) and standard deviation \(\sigma_\lambda\). Its count keeps the Poisson mean but not the Poisson variance:
\[E[N(t)] = \bar\lambda t, \qquad \operatorname{Var} N(t) = \bar\lambda t + \sigma_\lambda^2 t^2\]The variance-to-mean ratio above one, overdispersion, is the usual sign that count data are not Poisson.
coxMean
coxMean(meanRate, t) · meanRate mean rate · t time → expected count scalar
Let \(\bar\lambda > 0\) (the argument meanRate) be the mean of the random rate and \(t \ge 0\) a time. Then
The expected count.
Intuition When the rate is itself uncertain (busy days and quiet days, risky and safe customers), the expected count simply uses the average rate. Uncertainty about the rate leaves the mean untouched.
Use it when you forecast expected counts for events whose rate varies from unit to unit: insurance claims across policyholders of differing risk, demand under uncertain market conditions, defects across batches. The value equals poissonMean; the Cox model’s distinctive part is its larger variance, so use this together with coxSpread.
event-countingforecasting
def coxMean(meanRate, t) ≔ meanRate · t
coxSpread
coxSpread(meanRate, rateSpread, t) · meanRate mean rate · rateSpread standard deviation of the rate · t time → count variance scalar
Let \(\bar\lambda > 0\) (the argument meanRate) and \(\sigma_\lambda \ge 0\) (the argument rateSpread) be the mean and standard deviation of the random rate, and \(t \ge 0\) a time. Then
The variance of the count.
Intuition Two sources of noise add up: Poisson luck given the rate (\(\bar\lambda t\)) and not knowing the rate itself (\(\sigma_\lambda^2 t^2\)). The second grows with the square of the window, so over long windows the uncertainty about the rate dominates and does not average away.
Use it when you need honest uncertainty on counts: forecast intervals for demand, claims reserves, staffing for a volume whose rate is itself uncertain. The formula treats the rate as one random value held for the whole window. If the observed variance of your counts is close to the mean, poissonSpread is enough.
uncertainty-quantificationevent-countingforecasting
def coxSpread(meanRate, rateSpread, t) ≔ meanRate · t + rateSpread · rateSpread · t · t
overdispersion
overdispersion(meanRate, rateSpread, t) · meanRate mean rate · rateSpread standard deviation of the rate · t time → variance-to-mean ratio scalar
Let \(\bar\lambda > 0\) (the argument meanRate), \(\sigma_\lambda \ge 0\) (the argument rateSpread) and \(t > 0\). Then
The ratio of the count variance to its mean; 1 for a Poisson process, larger for a Cox process.
Intuition The variance-to-mean ratio of the count: exactly 1 for Poisson, and \(1 + \sigma_\lambda^2 t/\bar\lambda\) here, so the excess grows with the length of the window. It is the standard warning sign that counts are not Poisson.
Use it when you decide whether a Poisson model is adequate: compute the sample variance over the mean of real counts and compare it with 1 and with this prediction. A ratio well above 1 suggests a negative binomial or Cox model. It also shows how the choice of counting window changes the apparent noise.
model-validationclustering-of-events
def overdispersion(meanRate, rateSpread, t) ≔ coxSpread(meanRate, rateSpread, t) / nonzero(coxMean(meanRate, t))
Spatial points
Ripley’s K function counts, for a pattern of points in a square of side \(w\), how many neighbours lie within distance \(r\) of a typical point, scaled by the density. For complete spatial randomness \(K(r) = \pi r^2\); a larger value means clustering at that distance and a smaller one means regularity. This module scales by \(w^2/n^2\) rather than \(w^2/(n(n-1))\) and applies no edge correction, so compare patterns of similar size and expect some downward bias near the boundary.
neighbours
neighbours(xs, ys, r) · xs x coordinates vector n · ys y coordinates vector n · r distance scalar → count scalar
Let \(x, y \in \mathbb{R}^n\) (the arguments xs and ys) be the coordinates of \(n\) points, \(r \ge 0\) a distance, and \(D = \operatorname{apartness}(x, y)\). Then
The number of ordered pairs of points at distance at most r from each other. Each pair counts twice, and points at exactly the same location are not counted.
Intuition Draw a circle of radius \(r\) around every point and count the other points inside. Adding over all circles counts each close pair twice, once from each end.
Use it when you need the raw count behind ripleyK, or a direct closeness check: shops within walking distance of each other, disease cases within a radius, robots or aircraft inside a safety distance. Halve it for unordered pairs. Points at exactly the same location are not counted, so it cannot find duplicates. It builds an n×n distance matrix, which limits it to a few thousand points.
spatial-patternsproximity-detection
def neighbours(xs, ys, r) ≔ einsum("ij->", (apartness(xs, ys) ≤ r) · (apartness(xs, ys) > 0))
ripleyK
ripleyK(xs, ys, r, width) · xs x coordinates vector n · ys y coordinates vector n · r distance · width side of the square window → K(r) scalar
Let \(x, y \in \mathbb{R}^n\) (the arguments xs and ys) be the coordinates of \(n\) points in a square window of side \(w > 0\) (the argument width), and \(r \ge 0\) a distance. Then
Ripley’s K at distance r.
Intuition The average number of neighbours within \(r\) of a typical point, divided by the density of points. It asks how crowded a typical neighbourhood is compared with random scatter, which gives \(\pi r^2\).
Use it when you ask whether trees, crimes, cells or stores are clustered, random or evenly spaced, and at what distance: evaluate it over a range of \(r\) and compare with ripleyExpected. There is no edge correction, so at large \(r\) relative to width it is biased low. Dividing by \(n^2\) instead of \(n(n-1)\) also biases it low by a factor \((n-1)/n\), which matters only for small patterns.
spatial-patternsclustering-of-events
def ripleyK(xs, ys, r, width) ≔ width · width · neighbours(xs, ys, r) / nonzero(count(xs) · count(xs))
ripleyExpected
ripleyExpected(r) · r distance → πr² scalar
Let \(r \ge 0\) be a distance. Then
\[\operatorname{ripleyExpected}(r) = \pi r^2 .\]The value of K under complete spatial randomness.
Intuition Under pure randomness the expected number of neighbours within \(r\) is the density times the circle’s area, so after dividing by the density \(K\) is just the area \(\pi r^2\).
Use it when you need the benchmark for ripleyK: above it means clustering at that distance, below it means regular spacing. Plotting \(\sqrt{K/\pi} - r\) turns the benchmark into a flat zero line. Because ripleyK is biased low, random patterns land somewhat below this curve; for a firm verdict, simulate random patterns of the same size and compare with their range.
spatial-patternsmodel-validation
def ripleyExpected(r) ≔ 3.141592653589793 · r · r
Building blocks
apartness(xs, ys)— Let \(x, y \in \mathbb{R}^n\). Then \(\operatorname{apartness}(x, y)_{ij} = \sqrt{(x_j - x_i)^2 + (y_j - y_i)^2}\) for \(i, j = 1, \dots, n\). The n×n matrix of distances between every pair of points.
Fitting a point process
The maximum-likelihood rate of a Poisson process is the count over the horizon, with standard error \(\hat\lambda/\sqrt{n}\). Given their number, Poisson arrival times are uniform over the window, so their mean should be near half the horizon, which gives a cheap check of the assumption. For renewal data, the gap moments give the rate and the regime. A Hawkes process is fitted by moments with the decay known: the branching ratio is \(1 - \mu/\hat\lambda\), and \(\alpha\) follows.
fittedRate
fittedRate(times, horizon) · times arrivals vector n · horizon window length → rate scalar
Let \((t_1, \dots, t_n)\) (the argument times) be the arrival times observed in a window of length \(T > 0\) (the argument horizon). Then
The estimated rate.
Intuition Events per unit of observed time. For a Poisson process nothing else in the data carries information about the rate, so this simple ratio is the maximum-likelihood estimate.
Use it when you estimate an arrival rate from a log: requests per second, incidents per month, calls per hour. The result feeds poissonMean, simulations with timesOf, and the Hawkes moment fits. Use the length of the whole observation window, not the time of the last event, which would overstate the rate. Report it with rateError. When you only have the gaps between events, use gapRate.
parameter-estimationevent-arrivals
def fittedRate(times, horizon) ≔ count(times) / nonzero(horizon)
rateError
rateError(times, horizon) · times arrivals vector n · horizon window length → standard error scalar
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times in a window of length \(T > 0\) and \(\hat\lambda = \operatorname{fittedRate}(t_{1:n}, T)\). Then
The standard error of the estimated rate.
Intuition Precision depends on how many events you saw, not on how long you watched: the relative error is \(1/\sqrt{n}\), so 100 events pin the rate to about 10% and 10,000 to about 1%.
Use it when you put an interval on a rate (\(\hat\lambda \pm 2 \cdot\) error), compare two rates, or plan how long to observe (halving the error takes four times as many events). With no events the estimate and its error are both 0, which says nothing about precision; with very few events use an exact Poisson interval instead.
uncertainty-quantificationparameter-estimation
def rateError(times, horizon) ≔ fittedRate(times, horizon) / nonzero(sqrt(count(times)))
uniformity
uniformity(times, horizon) · times arrivals vector n · horizon window length → relative mean time scalar
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times in a window of length \(T > 0\) (the argument horizon). Then
The mean arrival time as a fraction of the horizon; about 0.5 for a Poisson process.
Intuition If the rate is constant, the events are spread evenly across the window, so their average time sits in the middle. A rising rate pulls the average later and a falling one earlier.
Use it when you want a cheap trend check on arrivals: demand growing through the window (above 0.5), activity fading after a launch (below 0.5). Under a constant rate its standard deviation is about \(1/\sqrt{12n}\), which tells you how far from 0.5 is surprising. It uses only the mean, so a pattern that is symmetric in time, such as bursts at both ends or a regular cycle, still gives 0.5; to look at bursts use gapRegime.
trend-detectionmodel-validation
def uniformity(times, horizon) ≔ mean(times) / nonzero(horizon)
clustered
clustered(times, horizon) · times arrivals vector n · horizon window length → flag 1 or 0
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times in a window of length \(T > 0\) and \(\rho = \operatorname{uniformity}(t_{1:n}, T)\). Then
1 when the mean arrival time is more than 0.1 of the horizon away from the middle, suggesting arrivals bunched early or late.
Intuition uniformity with a fixed tripwire: flag the window if the average event time is more than 10% of the window away from the middle.
Use it when you want a coarse dashboard flag for windows whose events are front- or back-loaded. The 0.1 threshold ignores sample size: with 5 events a constant rate crosses it often (the standard deviation is about 0.13), and with thousands of events a real but gentle trend slips under it. For a proper test, scale the distance from 0.5 by \(1/\sqrt{12n}\). Despite the name it does not detect bursts inside the window; use gapRegime for that.
trend-detectionclustering-of-events
def clustered(times, horizon) ≔ abs(uniformity(times, horizon) − 0.5) > 0.1
gapRate
gapRate(gaps) · gaps observed gaps vector n → rate scalar
Let \(g \in \mathbb{R}^n\) (the argument gaps) be observed gaps with mean \(\bar g = \frac{1}{n} \sum_i g_i\). Then
The renewal rate estimated as one over the mean gap.
Intuition If events are on average \(\bar g\) apart, they come at a rate of \(1/\bar g\). By the renewal theorem this holds whatever the shape of the gap distribution.
Use it when your data are times between events rather than a count in a window: mean time between failures turned into a failure rate, cycle times turned into throughput. It ignores the time before the first event and after the last, so on a short record it can differ from fittedRate, which uses the whole window.
parameter-estimationreliability-analysis
def gapRate(gaps) ≔ 1 / nonzero(mean(gaps))
gapVariation
gapVariation(gaps) · gaps observed gaps vector n → coefficient of variation scalar
Let \(g \in \mathbb{R}^n\) (the argument gaps), \(n \ge 2\), be observed gaps with mean \(\bar g\). Then
The sample standard deviation of the gaps over their mean.
Intuition The scatter of observed gaps relative to their mean, a unit-free burstiness score. Near 1 looks like random (Poisson) timing, well above 1 like bursts and lulls, well below 1 like a schedule.
Use it when you measure burstiness in real data: network traffic, user sessions, equipment failures. Match it to weibullVariation to pick a Weibull shape for simulation. It needs at least two gaps and is noisy with few; for a three-way label use gapRegime.
clustering-of-eventsparameter-estimation
def gapVariation(gaps) ≔ std(gaps) / nonzero(mean(gaps))
gapRegime
gapRegime(gaps) · gaps observed gaps vector n → regime +1, 0 or −1
Let \(g \in \mathbb{R}^n\) (the argument gaps) be observed gaps and \(c = \operatorname{gapVariation}(g)\). Then
\(+1\) when the gaps are bursty (coefficient of variation above 1.05), \(-1\) when they are regular (below 0.95), and 0 when they are consistent with Poisson.
Intuition gapVariation read as a sign, with a ±5% dead band around 1 so that a near-Poisson stream is not labelled either way.
Use it when you label streams to pick a planning method: bursty (\(+1\)) needs headroom for peaks, Poisson-like (0) suits standard queueing formulas, regular (\(-1\)) behaves like a schedule. The band is fixed and ignores sample size; with only tens of gaps the sampling noise in the coefficient of variation is wider than the band, so treat the label as indicative unless you have many gaps.
clustering-of-eventscapacity-planning
def gapRegime(gaps) ≔ (gapVariation(gaps) > 1.05) − (gapVariation(gaps) < 0.95)
branchingFrom
branchingFrom(times, horizon, background) · times arrivals vector n · horizon window length · background known μ → branching ratio scalar
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times in a window of length \(T > 0\) (the argument horizon), \(\mu \ge 0\) (the argument background) the known background rate, and \(\hat\lambda = n / \operatorname{nonzero}(T)\). Then
The moment estimate of the Hawkes branching ratio.
Intuition A stable Hawkes process runs at \(\mu/(1 - m)\), so any activity above the background is credited to triggering, and \(m = 1 - \mu/\hat\lambda\) follows. If the observed rate is below the background, the answer is clipped to 0.
Use it when the background rate is known from elsewhere (a quiet period, a count of external shocks) and you want a quick estimate of how much activity is self-generated: market reactions, social cascades, follow-on incidents. It uses only the event count, not the timing, and assumes a stable process observed over a long enough window. For a proper fit, use this as a starting value and maximise hawkesLikelihood.
parameter-estimationclustering-of-events
def branchingFrom(times, horizon, background) ≔ max(1 − background / nonzero(observedRate(times, horizon)), 0)
exciteFrom
exciteFrom(times, horizon, background, decay) · times arrivals vector n · horizon window length · background known μ · decay known β → α scalar
Let \((t_1, \dots, t_n)\) (the argument times) be arrival times in a window of length \(T > 0\), \(\mu \ge 0\) (the argument background) the known background rate, and \(\beta > 0\) (the argument decay) the known decay. Then
The moment estimate of the excitation \(\alpha\): the branching ratio times the decay.
Intuition The branching ratio is \(\alpha/\beta\), so once the decay is known, the jump each event adds to the rate is the ratio times the decay. A fast-fading fire must jump higher to trigger as many children.
Use it when you complete a parameter set \((\mu, \alpha, \beta)\) to pass to excitedAt, hawkesRates or hawkesLikelihood, or seed an optimiser. The decay has to come from outside the data given here, such as the typical time over which aftershocks or reactions fade. Any error in it scales \(\hat\alpha\) by the same factor.
parameter-estimationclustering-of-events
def exciteFrom(times, horizon, background, decay) ≔ branchingFrom(times, horizon, background) · decay
Building blocks
observedRate(times, horizon)— Let \((t_1, \dots, t_n)\) be arrival times in a window of length \(T > 0\). Then \(\operatorname{observedRate}(t_{1:n}, T) = \operatorname{fittedRate}(t_{1:n}, T) = n / \operatorname{nonzero}(T)\). The observed event rate, the same asfittedRate; used bybranchingFrom.
Examples
Arrival times from uniform draws
open Prelude.Points
fn arrivalTimes(u, rate) ↦ timesOf(u, rate)
fn arrivedBy(u, rate, t) ↦ countedBy(u, rate, t)
using Examples;
using MLambda.Hilbert.Runtime;
var u = new[] { 0.5, 0.25, 0.75 };
var times = Arrivals.ArrivalTimes((Tensor)u, 2.0);
var byOne = Arrivals.ArrivedBy((Tensor)u, 2.0, 1.0);
Output (verified, rounded to 6 decimal places):
times = [0.346574, 0.490415, 1.183562]
counted = 2
At rate 2 the gaps are \(\ln 2 / 2\), \(\ln(4/3)/2\) and \(\ln 4 / 2\); their running totals put two arrivals before time 1.
Evaluating a Hawkes process
Background rate 0.5, excitation 0.8, decay 1, with arrivals at times 1 and 2.
open Prelude.Points
fn longRunRate(background, excite, decay) ↦ hawkesRate(background, excite, decay)
fn intensityAt(times, background, excite, decay, t) ↦ excitedAt(times, background, excite, decay, t)
fn likelihoodOf(times, background, excite, decay, horizon) ↦ hawkesLikelihood(times, background, excite, decay, horizon)
using Examples;
using MLambda.Hilbert.Runtime;
var times = (Tensor)new[] { 1d, 2d };
var rate = Arrivals.LongRunRate(1.0, 0.5, 1.0);
var atThree = Arrivals.IntensityAt(times, 0.5, 0.8, 1.0, 3.0);
var logLik = Arrivals.LikelihoodOf(times, 0.5, 0.8, 1.0, 3.0);
Output (verified, rounded to 6 decimal places):
hawkesRate = 2
intensity = 0.902572
likelihood = -3.620865
With branching ratio ½ a background of 1 doubles to a long-run rate of 2. At time 3 the intensity is \(0.5 + 0.8(e^{-2} + e^{-1})\).
Renewal gaps: burstiness and the inspection paradox
open Prelude.Points
fn gapVariationOf(shape) ↦ weibullVariation(shape)
fn waitingLeft(shape, scale) ↦ residualLife(shape, scale)
using Examples;
var cvHalf = Arrivals.GapVariationOf(0.5);
var cvOne = Arrivals.GapVariationOf(1.0);
var poissonWait = Arrivals.WaitingLeft(1.0, 2.0);
var burstyWait = Arrivals.WaitingLeft(0.5, 2.0);
Output (verified, rounded to 6 decimal places):
cv(0.5) = 2.236068
cv(1) = 1
residual(1,2) = 2
residual(0.5,2) = 12
At shape ½, \(\Gamma(5)/\Gamma(3)^2 - 1 = 5\), so the coefficient of variation is \(\sqrt5\). Exponential gaps of mean 2 leave an observer a whole mean gap of 2; Weibull gaps with shape ½ and scale 2 have mean 4 but leave \(4 \cdot 6 / 2 = 12\).
Ripley’s K for three points
open Prelude.Points
fn closePairs(xs, ys, r) ↦ neighbours(xs, ys, r)
fn ripleyOf(xs, ys, r, width) ↦ ripleyK(xs, ys, r, width)
using Examples;
using MLambda.Hilbert.Runtime;
var xs = (Tensor)new[] { 0d, 1d, 0d };
var ys = (Tensor)new[] { 0d, 0d, 1d };
var pairs = Arrivals.ClosePairs(xs, ys, 1.0);
var k = Arrivals.RipleyOf(xs, ys, 1.0, 10.0);
Output (verified, rounded to 6 decimal places):
neighbours = 4
ripley = 44.444444
Points (1, 0) and (0, 1) are each at distance 1 from the origin and \(\sqrt2\) from each other, giving two pairs counted both ways. In a 10 × 10 window, \(K = 100 \cdot 4 / 9\), far above \(\pi \approx 3.14\): three points this close together are strongly clustered.
Fitting a rate and reading overdispersion
open Prelude.Points
fn rateOf(times, horizon) ↦ fittedRate(times, horizon)
fn rateSe(times, horizon) ↦ rateError(times, horizon)
fn regimeOf(gaps) ↦ gapRegime(gaps)
fn excessOf(meanRate, rateSpread, t) ↦ overdispersion(meanRate, rateSpread, t)
using Examples;
using MLambda.Hilbert.Runtime;
var times = (Tensor)new[] { 0.5, 1.5, 4.0, 6.5, 8.0, 9.5 };
var rate = Arrivals.RateOf(times, 10.0);
var se = Arrivals.RateSe(times, 10.0);
var regime = Arrivals.RegimeOf((Tensor)new[] { 1d, 1d, 1d, 1.2 });
var excess = Arrivals.ExcessOf(2.0, 1.0, 3.0);
Output (verified, rounded to 6 decimal places):
rate = 0.6
rateSe = 0.244949
regime = -1
overdispersion = 2.5
Six arrivals in 10 time units give \(0.6 \pm 0.6/\sqrt6\). Nearly equal gaps are regular (\(-1\)). A Cox process with mean rate 2 and rate spread 1 over 3 time units has variance \(6 + 9 = 15\) against a mean of 6.