Suppose the observed data are pairs $(x_i,y_i)$, collected as $\mathcal D=\{(x_i,y_i)\}_{i=1}^{n}$. The input $x_i$ might be a time or dosage, and the output $y_i$ might be a temperature or measured response. The goal is to learn the relationship $x\mapsto y$ well enough to predict the output at a new input. This prediction problem is called regression.

Begin with one fitted function $\hat f$: one deterministic input–output rule selected to explain the observations and predict new outputs. Several deterministic functions may explain the observed pairs while giving different outputs between them or beyond their range. Replacing the single commitment with a probability distribution over these plausible functions preserves that uncertainty.

Each plausible function is one complete deterministic mapping across the input domain. Once a particular function is selected, it assigns one definite value $f(x)$ to every input $x$. The observed input–output pairs change the probability of each complete function according to how well that one mapping explains the data as a whole. At a new input, the same collection of functions supplies the possible outputs.

At one chosen input $x_1$, every plausible function supplies one possible value $f(x_1)$. The probabilities assigned to the functions therefore produce a probability distribution over the possible values of $f(x_1)$. At two chosen inputs, every function supplies an ordered pair $(f(x_1),f(x_2))$. A probability distribution over these pairs describes the two function values together, so it is called their joint distribution.

These output vectors make a probability distribution over complete functions concrete. A Gaussian process assigns probabilities to functions so that, whenever inputs $x_1,\ldots,x_m$ are selected, their function values $(f(x_1),\ldots,f(x_m))$ have a multivariate Gaussian joint distribution.

The function itself covers the complete input domain, which may contain every real number in an interval. The word finite describes only how many inputs are selected in one question: perhaps two, ten, or one thousand. One selected input gives an ordinary Gaussian distribution, two selected inputs give a two-dimensional Gaussian distribution, and $m$ selected inputs give an $m$-dimensional Gaussian distribution. The two-input distribution also records whether the two function values tend to be high and low together.

A Gaussian process supplies one such joint distribution for every finite selection of inputs. Compatibility makes the distributions agree whenever their selected inputs overlap.

The original goal was to replace one fixed predictive function with a probability distribution over complete candidate functions. The mean function and covariance kernel now provide that distribution. The next task is to use measured input–output pairs to determine which candidates remain plausible and then predict at an input where no output has been measured.

The first panel summarizes the starting probabilities across functions. At each input, the cyan line marks their average value and the band contains the central 95% of their values. These probabilities vary according to the chosen mean and kernel: a zero mean function centers the candidates around zero, while the RBF kernel gives more probability to candidates whose values change smoothly across its length scale. This starting distribution is the prior distribution over functions.

The middle panel leaves that distribution unchanged and places the measured input–output pairs over it. The third panel shows the computed result after those measurements are used. Its mean moves toward the observed outputs, and its band narrows most near observed inputs. This change occurs because complete functions that explain all the measurements receive more relative probability than functions whose values lie farther away.

The RBF kernel is the covariance rule used to construct the joint Gaussian distribution of function values. A draw from that joint distribution supplies connected values belonging to one candidate function $f$. Its value $f(x)$ is the output that this candidate assigns to input $x$. The candidate function therefore comes from a Gaussian-process draw; it is the object being learned and used for prediction.

A measured pair $(x_i,y_i)$ compares one candidate function with reality. For example, a candidate may assign $f(2)=0.80$ while the measurement at that input is $y_i=0.92$. Their difference is $0.12$. This difference is the candidate’s residual at that input. Repeated measurements can vary around the underlying response because of sensor variation and influences outside the function. Represent that measurement variation by a random quantity $\varepsilon_i$ called observation noise. The measured output is:

The lowercase $\sigma_n$ is the standard deviation of measurement noise. Its subscript $n$ distinguishes it from the function scale $\sigma_f$. The RBF kernel uses $\sigma_f$ to describe the prior spread of underlying function values; the observation model uses $\sigma_n$ to describe how far measurements typically vary from those values. The noise draws are independent across measurements and independent of the candidate function.

For one complete candidate function, calculate a residual at every measured input. Small residuals across the entire dataset make the measurements more probable under that candidate; large residuals make them less probable. The resulting data-dependent multiplier is the likelihood. Write the measured inputs as $X=(x_1,\ldots,x_n)$ and their outputs as $\mathbf y=(y_1,\ldots,y_n)^{\mathsf T}$. Under the Gaussian observation model, the likelihood has the form:

The squared residuals appear in the exponent. A smaller $\sigma_n$ makes the likelihood fall rapidly as a candidate moves away from the measurements. A larger $\sigma_n$ makes that fall slower, so residual differences produce a smaller change in the candidates’ relative probabilities.

Multiply each candidate’s prior probability by its likelihood, then normalize the results so that the updated probabilities sum to one. This multiplication-and-normalization rule is Bayes’ rule. Using it to incorporate the observed data is called conditioning, and the updated distribution is the posterior distribution over functions:

Now select a new input $x_*$. Every posterior candidate supplies one value $f(x_*)$from the same complete function that was compared with all the measurements. The distribution across these values describes uncertainty about the underlying response at $x_*$. A future measured output also includes a new draw of observation noise:

The resulting distribution of possible future measurements $Y_*$ is the posterior predictive distribution. It is the probabilistic answer sought at the beginning: a center and a spread over plausible outputs at the selected input, informed by every observed pair.

Move the question location $x_*$ in the laboratory. Near several observations that support similar function values, the predictive distribution is narrow. In the large gap or near an edge, it widens. Increase the observation noise $\sigma_n$ and differences in residual error change the relative function probabilities less. Increase the RBF length scale $\ell$ and an observation changes predictions across a wider range of inputs. Increase the function scale $\sigma_f$ and the prior allows function values farther above and below the mean function.

With a zero mean function and the RBF kernel, the posterior mean function is a weighted sum of RBF kernel shapes centered at the observed inputs. The measured outputs determine the weights; the length scale determines how far each centered shape reaches.

The laboratory performs this Bayes update at every question location. Gaussian distributions make the update calculable without listing or scoring the candidate functions one by one. The same calculation can be organized in four steps:

1. 1
   Collect the measured pairs

   Write the observed inputs as $X=(x_1,\ldots,x_n)$ and their outputs as $\mathbf y=(y_1,\ldots,y_n)^{\mathsf T}$.

2. 2
   Relate every pair of observed inputs

   Apply the kernel to each pair and place the answers in $\mathbf K$, where $K_{ij}=k(x_i,x_j)$.

3. 3
   Connect the question location to the data

   Place $k(x_i,x_*)$ for every observed input in $k_*$, and compute its prior variance $k_{**}=k(x_*,x_*)$.

4. 4
   Condition the joint Gaussian distribution

   Combine $\mathbf y$, $\mathbf K$, $k_*$, $k_{**}$, and $\sigma_n$ to obtain the predictive center and spread at $x_*$.

Before the outputs are observed, the measured outputs $\mathbf y$ and future output $Y_*$ belong to one joint multivariate Gaussian distribution supplied by the process. Conditioning that distribution on $\mathbf y$ gives this posterior predictive distribution when the mean function is zero:

The first expression is the predictive mean. It combines observed values according to their kernel relationship with $x_*$. The second expression begins with the prior variance $k_{**}$ and removes the part resolved by kernel connections to the observed inputs. The final $\sigma_n^2$ adds the variation of a new measurement around its underlying function value.

One random probability measure $G$ drawn from the Dirichlet process places probability mass on particular component parameters. Several draws $\theta_1,\theta_2,\ldots$ from that same $G$ can therefore return exactly the same parameter. Observations whose parameters repeat belong to the same cluster; a newly drawn parameter begins a new cluster. A distribution that concentrates mass on individual points is called discrete, and a Dirichlet-process draw is discrete with probability one even when $H$ is continuous.

Picture each distinct component parameter as a table and each draw from $G$ as a customer. After $n$ draws, $K_n$ tables are occupied. Table $j$ represents one distinct parameter $\theta_j^*$, and its $n_j$ customers record how many times that parameter has repeated.

The next draw reuses the parameter at table $j$ with weight $n_j$. The concentration $\alpha$ supplies the weight for opening a new table whose parameter is drawn from $H$. Dividing by the total weight $\alpha+n$ turns these weights into probabilities. This table-and-customer predictive rule is called the Chinese restaurant process. The symbol $\delta_{\theta_j^*}$ represents a distribution whose entire mass sits on the existing parameter $\theta_j^*$:

This predictive rule averages across all possible values of $G$. The operation is called integrating out $G$, and it directly gives the probabilities of reusing an existing component parameter or drawing a new one from $H$.

The number of occupied tables adapts as observations arrive. The concentration $\alpha$, the base probability measure $H$, and the component density jointly shape how readily new clusters appear and how observations accumulate within existing clusters.

If the component parameters $\theta_{1:n}$ were known, they would also update the probability distribution over $G$ itself. Each known parameter adds probability mass at its location, while the total concentration grows from $\alpha$ to $\alpha+n$. The resulting updated distribution remains a Dirichlet process:

Each component parameter that receives positive probability mass from $G$ is called an atom. A second construction makes these atoms and their weights visible. Break a unit-length stick repeatedly, draw one atom from $H$ for each piece, and place that piece of probability on its atom:

The fraction $v_j$ breaks off part of the probability still available at step $j$. The resulting piece $\pi_j$ becomes the weight of component parameter $\theta_j$. Repeating the construction indefinitely produces the complete discrete probability measure $G$.

A Gaussian process supplies an uncertain function; a Dirichlet process supplies an uncertain probability measure. Now consider events recorded by a detector. We want to describe the entire collection of event times, so that we can ask how many events occur in any selected interval. A Poisson process provides this distribution over event configurations. We will first derive its count probabilities, then recover its waiting times, and finally let each event change the conditions governing the next one. That last extension describes the activation and deactivation events in a neural population.

One draw is an event configuration

Suppose we choose to record the detector for two seconds, starting at time zero. This choice fixes the observation interval; what remains unknown is how many events will occur and when they will occur within it. The Poisson process describes that uncertainty. A draw supplies an event configuration, whose events within our chosen interval might occur at 0.4, 1.1, and 1.8 seconds. Another draw can give different event times and a different total count in the same interval. We can choose a longer or shorter observation interval as needed.

For the configuration with events at 0.4, 1.1, and 1.8 seconds, count how many events have occurred by each time: zero initially, one after 0.4 seconds, two after 1.1, and three after 1.8. The event timeline and this staircase describe the same outcome. Call the staircase the counting function, $N(t)$, where $t$ is elapsed time and $N(0)=0$.

Select the interval (t₁, t₂], with t₁ = 0 seconds and t₂ = 2 seconds: it includes times after t₁ up to and including t₂. In the configuration shown above, this extracts a count of three events. Selecting (0, 1] instead extracts one event. Once a configuration is fixed, these counts are fixed, just as selecting inputs from a fixed Gaussian-process function extracts fixed function values. Before the configuration is drawn, the count in a selected interval is uncertain; its probabilities come from the process that generates configurations. For example, the probability of at least one event in (0, 1] is the probability of the set of configurations containing at least one event in that subinterval.

The endpoints specify where an interval lies; its duration specifies how long it lasts. If we extend the recording beyond two seconds, counting events between 5 and 10 seconds means taking N(10) − N(5): the count up to 10 seconds minus the count up to 5 seconds. N(5) itself counts events from 0 to 5 seconds. For the constant-rate Poisson process considered here, intervals of equal duration have the same count distribution. Thus (0, 5], (5, 10], and (15, 20] all have the same count distribution, even though their observed counts can differ. Their counts are also independent because these intervals are disjoint. Overlapping intervals share events and generally have dependent counts.

For that same configuration, assign every interval its event count. Counts add when disjoint intervals are combined: the one event in (0, 1] plus the two events in (1, 2] give three in (0, 2]. Because this equality holds for every configuration, counting the whole interval directly or adding its two half-interval counts gives the same probability distribution. This is an instance of compatibility: the joint count distributions agree when we select the same intervals from different larger collections. Deriving all these counts from one random configuration supplies the compatibility in our earlier definition of a stochastic process. For each fixed configuration, these counts determine a counting measure. It assigns whole-number event counts to intervals and, more generally, measurable sets of times.

From possible configurations to a count probability

The displayed configuration answers a counting question once its event times are known. To predict a count before observing the events, we need to consider the other possible configurations of the same two-second interval and their probabilities. Several configurations can give the same count, just as several Gaussian-process functions can give the same value at a selected input. We will build a finite approximation, assign probabilities to its configurations, and add the probabilities of configurations with the requested count. Making the approximation finer will give the continuous-time count distribution.

First divide the interval into eight equal time slots, each a quarter of a second long. In the approximation, a slot contains either one event or zero events. The original times 0.4, 1.1, and 1.8 seconds occupy slots 2, 5, and 8. Other choices of three occupied slots describe other configurations with the same count. Each occupied slot records that an event occurred somewhere inside it; its precise time is left unspecified. The figure shows three such choices.

The slot choices specify the possible configurations; a probability rule tells us how likely each choice is. The illustrated configurations each contain three events in two seconds, an observed rate of 1.5 events per second. Choose a model rate of 1.5 events per second to match this example; its expected count in two seconds is three. Individual draws can still have different counts.

Our approximation uses eight slots, each holding either zero or one event. Make each slot independently occupied with probability 3/8. Across eight slots, the expected number of occupied slots is 8 × (3/8) = 3, giving three events on average.

More generally, call the rate $\lambda$ and the duration $T$. Their product $\mu=\lambda T$ is the expected count. Write n for the number of slots. Let $p_n$ be the probability that a particular slot is occupied. With the same probability for every slot, the expected count is $n p_n$. To make that count equal to$\mu$, choose $p_n=\mu/n$. The number of slots must be large enough that this probability is at most one. Here $p_8=3/8$.

Configuration A's slot pattern has three occupied slots and five empty slots. Independence lets us multiply the eight probabilities: three factors of p₈ and five factors of 1 − p₈. Every other choice of three occupied slots has that same probability. Counting the choices gives the binomial coefficient, written “8 choose 3”.

For n slots and a requested count k, there are “n choose k” choices. Call the random occupied-slot count Nₙ(T). Adding the probabilities of these mutually exclusive choices gives its binomial distribution:

Eight slots are a coarse approximation: each permits at most one event. Refine to 16, then 32 slots, and onward, keeping T = 2 seconds and μ = 3 fixed. The slots become shorter, and their event probabilities fall to 3/16, then 3/32. Each finite approximation gives a consistent set of counting answers on its own grid; increasing n changes the approximate law. In the limit, events can occur at continuous times. Thus discretizing time constructs an approximation to the process; the process itself has no fixed time bins. A Poisson process supplies a Poisson count distribution for each selected interval, with these distributions satisfying compatibility across intervals.

To derive the Poisson distribution from the finite-slot approximation, let n increase without limit while keeping T, μ, and the requested count k fixed, where k is a nonnegative integer. The expression $\binom{n}{k}(\mu/n)^k$ tends to $\mu^k/k!$, and $(1-\mu/n)^{n-k}$ tends to $e^{-\mu}$. The first line rewrites its expression as a product of k factors, indexed by i from zero to k − 1; each factor $1-i/n$ tends to one. Both arrows show the limit as n tends to infinity.

The finite-slot probability is already the product of the two expressions whose limits we just calculated. By the product rule for limits, its limit is the product of their limits: replace$\binom{n}{k}(\mu/n)^k$ by $\mu^k/k!$ and$(1-\mu/n)^{n-k}$ by $e^{-\mu}$.

This gives the Poisson distribution for the count. Recall that μ = λT: the expected count equals the rate multiplied by the interval duration. Substituting λT for μ in both places gives$e^{-\mu}\mu^k/k!=e^{-\lambda T}(\lambda T)^k/k!$. Write N(T) for the count in the continuous-time limit. Its probability for each whole-number count is:

Applying this construction consistently across the time axis gives the homogeneous Poisson process. “Homogeneous” means the rate is constant: equal-length intervals have the same count distribution. For times s and t with s ≤ t, the number of events in the interval (s, t] is the difference of the two cumulative counts:

These differences are called increments. Consider N(1) and N(2), the cumulative counts up to 1 and 2 seconds. Their difference N(2) − N(1) counts events in (1, 2]. Learning that N(1) = 3 leaves the distribution of N(2) − N(1) unchanged. This is independence: learning one count leaves the other’s probabilities unchanged. In this constant-rate Poisson model, increments over disjoint intervals are independent.

Now compare the first-second count N(1) with the total two-second count N(2). Before observing anything, N(2) could be zero. Learning that N(1) = 3 tells us immediately that N(2) ≥ 3, because those three events are included in both counts. This is dependence: learning one count changes the probabilities for the other. With the positive rate used in our example, N(1) and N(2) are dependent.

• N(1) and N(2) − N(1): independent, because they count events in disjoint intervals.
• N(1) and N(2): dependent, because N(2) includes the events counted by N(1).

This distinction matters when calculating probabilities. For independent counts, multiplying their separate probabilities gives the probability of both outcomes occurring together. For dependent counts, we account for how knowing one changes the probabilities of the other.

The whole process specifies these relationships as well as each interval’s count distribution. Compatibility ensures that selecting the same counts from different larger collections gives the same joint distribution. Independence of counts in disjoint intervals is an additional property of this Poisson model. Together, these properties let us combine interval counts while preserving a consistent description of the same random configuration.

From event counts to the time of the first event

So far, we have selected an interval and asked how many events it contains. The same event configuration also answers another question: how much time passes before the first event occurs? In our illustrated configuration, the first event occurs at 0.4 seconds. This elapsed duration is called the waiting time.

Now consider a new random configuration, before its event times are known. Our model has a rate of 1.5 events per second. Select the interval from 0 to 0.5 seconds and ask: what is the probability of zero events in this interval, given the rate and interval duration?

To calculate this probability across new draws, use the Poisson count formula with k = 0, rate λ = 1.5 events per second, and duration T = 0.5 seconds. Here zero is the requested count; the expected count is λT = 1.5 × 0.5 = 0.75. The factor (λT)ᵏ/k! becomes 0.75⁰/0! = 1/1 = 1, leaving only the exponential factor:

Therefore our model says that the probability that the waiting time exceeds 0.5 seconds is approximately 47.2%.

Let τ be the random waiting time until the first event, and let u be a chosen elapsed duration, such as 0.5 seconds. Given the model rate λ, asking for the probability that τ exceeds u is equivalent to asking for the probability of zero events in (0, u]. Here N(u) is the count in that interval. This probability is called the survival probability. Substituting k = 0 and duration u into the Poisson count formula gives:

Since the probability that the waiting time exceeds u is $e^{-\lambda u}$, the probability that it is at most u is $1-e^{-\lambda u}$. These two possibilities exhaust the outcomes. In other words, $1-e^{-\lambda u}$ is the probability that the first event has occurred by time u. This is called the cumulative probability, written $F_\tau(u)$.

Plugging in u = 0.5 seconds and λ = 1.5 events per second gives$F_\tau(0.5)=1-e^{-0.75}\approx0.528$: about a 52.8% chance that the first event has occurred by then. As u grows, this probability accumulates toward one.

The two curves below answer complementary questions at the same rate of 1.5 events per second. At each elapsed time u, their probabilities add to one: either the first event has occurred by then, or it is still to come.

This waiting-time distribution is called the exponential distribution. Its mean waiting time is 1/λ for a positive rate. A zero rate produces no events.

Generate a configuration by sampling successive waits

We have calculated probabilities for possible configurations. To simulate a recording, we now need to produce one configuration with those probabilities. Start at time zero, draw a waiting time, and place the first event that far along the timeline. Draw another wait and add it to the first event’s time to place the next event. Repeating these steps produces the event configuration; stop when the next event would fall beyond the recording.

The cumulative probability tells us how to draw each wait. Choose a number U uniformly between zero and one, find that height on the cumulative probability curve, and read off its elapsed time. For example, a height near 0.528 gives a wait near 0.5 seconds at our rate of 1.5 events per second. Because 52.8% of uniform draws fall below that height, approximately 52.8% of the sampled waits are at most 0.5 seconds, as required. Solving for the time at height U gives this sampling rule:

Here U is a fresh uniform draw in (0, 1), λ is the fixed rate, and τ is the sampled delay in seconds. Choosing independent uniform draws gives independent exponential waits. For example, waits of 0.4, 0.7, and 0.7 seconds place events at 0.4, 1.1, and 1.8 seconds—the configuration used earlier.

At a constant rate, observing no event for a while leaves the same distribution for the remaining wait. This follows because counts in the subsequent interval are independent of earlier counts and depend only on its duration. This property is called memorylessness. Successive independent exponential waits generate the same Poisson process whose count probabilities we derived.

The cumulative probability supplies the sampling rule directly. To examine how likely the first event is within a short interval, the same distribution also has a density representation, developed below.

Sampling at a fixed rate now gives a complete event configuration. In a neural population, however, an event can change how many neurons are eligible for the next event. To simulate that feedback, we need to record what changed and update the rates before drawing the next wait.

Let events change a state

Consider ten neurons, with three currently active. An activation changes the active count from 3 to 4; a deactivation changes it from 3 to 2. This active count is the state, written $X(t)$. The simulation must now choose both when the next event occurs and which of these changes it makes. It will calculate the rates from the current state, sample a delay and event type, update the state, and repeat.

Suppose each inactive neuron activates at 0.2 per second and each active neuron deactivates at 0.1 per second. With three active neurons, seven can activate and three can deactivate. The combined activation rate is 7 × 0.2 = 1.4 per second; the combined deactivation rate is 3 × 0.1 = 0.3 per second. These replace the single fixed rate of 1.5 used for the earlier detector example.

While this state stays unchanged, model the two event types with independent exponential waits at their respective rates. For the next event to take longer than duration u, both types must have zero events during that duration. Multiplying those independent probabilities gives:

This is the survival probability of an exponential wait with total rate 1.7 per second. We can therefore use the same waiting-time sampling rule with λ = 1.7. The event is an activation with probability 1.4/1.7 and a deactivation with probability 0.3/1.7. To select its type, draw another independent uniform number V: choose activation if V is below 1.4/1.7, and deactivation otherwise.

Each event changes the rates for what follows

Suppose the sampled event is an activation. Advance time by the sampled delay and change the active count from 3 to 4. Now six neurons can activate and four can deactivate: the rates become 1.2 and 0.4 per second. The next wait uses total rate 1.6, and the next activation probability is 1.2/1.6. Recalculating after each event lets the generated configuration reflect its own effects on the population.

Write x for the current active count, with + labeling activation and − labeling deactivation.

At zero active neurons the deactivation rate is zero, and at ten the activation rate is zero. If the total rate is zero, the state stays fixed until an external change supplies a positive rate. In Wilson–Cowan, the per-neuron activation rate also depends on the shared input generated by the current population counts, giving further feedback.

Between events the state remains constant; each event makes a discrete change. This is a jump process. A draw now records both the event times and the resulting states, forming a state trajectory.

This explains the name continuous-time Markov jump process. “Continuous-time” describes the time axis; “jump” describes the discrete state changes. “Markov” says that the current state supplies all the information from the past needed to specify the future law. Here the rates depend only on x, and the state remains constant between events. A rate depending on additional unrecorded history would require a richer state to make that claim.

The ordinary Poisson counting process is the special case with one possible change, adding one to the count, at a constant rate. The Markov property is shared by both constructions; state dependence identifies the added feedback. “State-dependent Poisson event generation” describes the mechanism; “continuous-time Markov jump process” names the resulting state trajectory and its prediction rule.

The simulation procedure now follows directly: calculate the current rates, draw the next waiting time and event type, apply the state change, and repeat. This is the Gillespie direct method. Simulation time advances by each sampled wait; physical time is the same continuous variable used in the corresponding population equations.

For the explicit six-event neural example, continue to How the Gillespie Algorithm Works. It also distinguishes simulated activation counts from the illustrative neuron rows used in the demo raster.

Source: Anderson and Kurtz, Continuous Time Markov Chain Models for Chemical Reaction Networks, for the Poisson construction, state-dependent event rates, and random time-change representation.