Vectorizing the Survivor-Density Update

Why Vectorize a Survivor-Density Update?

A recursive boundary calculation repeatedly performs the same operation: take the distribution of states that have survived so far, propagate it through a Gaussian transition, and discard states that cross a stopping boundary. Written as an integral, this is conceptually simple. Written as nested loops over thousands of grid points at every analysis, it can be slow and vulnerable to indexing errors. Vectorization expresses the same calculation as a weighted matrix–vector product, so numerical libraries can perform the repeated arithmetic efficiently.

Vectorization does not change the statistical model or make a coarse grid accurate. Its value is computational organization: the transition, integration, and continuation restriction are explicit array operations, and probability accounting can be checked alongside them.

Model and Discretized Update

Use the canonical score scale \(Y_k\), for which the score at information \(I_k\) has mean \(\theta I_k\) and variance \(I_k\). Here \(\theta\) is the drift under the parameter value being evaluated. Under the independent-increment Gaussian model, conditional on \(Y_k=x\), the next state satisfies

At look \(k+1\), let \(a_{k+1}\) and \(b_{k+1}\) be the lower and upper continuation boundaries on this score scale, with \(a_{k+1}<b_{k+1}\). A path continues if its state lies between them. The boundaries can be obtained by expressing the study's \(Z\)-scale boundaries on the score scale; this tutorial focuses on propagating the density once the boundaries are specified.

Choose source-grid nodes \(x_1,\ldots,x_n\) at look \(k\) and destination-grid nodes \(y_1,\ldots,y_m\) at look \(k+1\). Let \(w_i\) and \(v_j\) be the corresponding quadrature weights. The source probability represented by node \(i\) is approximately \(q_i=w_i f_k(x_i)\). The Gaussian transition density evaluated at each source-destination pair is

where \(\phi\) is the standard normal density. In an array convention with sources in rows and destinations in columns, the propagated density values are

The destination continuation mask is \(c_j=1\) when \(a_{k+1}<y_j<b_{k+1}\), and \(c_j=0\) otherwise. The updated survivor density and its quadrature mass are

The multiplication by \(w_i\) belongs on the source side because the integral is over source states. The destination weights \(v_j\) are used afterward to integrate the updated density. Keeping these roles distinct is central to implementing the update correctly.

Derivation: From an Integral to an Array Operation

Conditional on survival through look \(k\), the density at the next state before applying its boundary is the integral of the previous survivor density times the transition density. Since \(f_k\) is unnormalized, the integral also remains unconditional with respect to survival.

Applying a quadrature rule to this integral at each destination node gives \(g(y_j)\approx\sum_i w_i f_k(x_i)K_{ij}\). Defining \(q_i=w_if_k(x_i)\) turns those repeated sums into \(\mathbf K^{\mathsf T}\mathbf q\). Finally, multiplying each destination value by \(c_j\) enforces the continuation condition. Thus the vectorized update is the discretized integral itself, not a different approximation method.

The array shape provides a useful orientation check: \(\mathbf K\) has \(n\) rows and \(m\) columns, \(\mathbf q\) has length \(n\), and \(\mathbf K^{\mathsf T}\mathbf q\) has length \(m\). If source and destination grids have the same number of nodes, a mistaken transpose can still produce a conformable result, so shape alone is not enough; the source-versus-destination definition of \(K_{ij}\) must also be checked.

Probability Accounting Alongside the Update

A useful independent check comes from the conditional normal distribution. For source node \(x_i\), define \(\mu_i=x_i+\theta\Delta I_k\) and \(\sigma_\Delta=\sqrt{\Delta I_k}\). The probabilities of crossing below the lower boundary, continuing, and crossing above the upper boundary are

For every source state, \(p_{L,i}+p_{C,i}+p_{U,i}=1\), up to numerical rounding. Weighting by the source masses gives the corresponding event probabilities for paths that survived through look \(k\): \(P_L\approx\sum_iq_ip_{L,i}\), \(P_C\approx\sum_iq_ip_{C,i}\), and \(P_U\approx\sum_iq_ip_{U,i}\). Their sum should equal the source survivor mass \(\sum_iq_i\), subject to the accuracy of the source quadrature and arithmetic.

This CDF calculation is a probability-accounting check, not a replacement for evaluating the propagated density at destination nodes. Agreement between \(P_C\) and the destination-grid estimate \(S_{k+1}\) depends on the destination quadrature and grid coverage. The CDF values integrate the normal transition over each boundary-defined region directly; the density update samples the transition across the destination grid.

Implementation, Diagnostics, and Alternatives

In a vectorized implementation, construct the \(n\)-by-\(m\) transition array, multiply its transpose by the source-mass vector, apply the destination mask, and integrate with destination weights. For multiple parameter values, form the appropriate transition array for each drift; the drift changes the conditional means and therefore the kernel. Vectorization removes interpreted loops over individual source-destination pairs, but the calculation still requires order \(nm\) kernel evaluations and multiplications. An array with \(n=m=1200\) has \(1{,}440{,}000\) entries; in double precision this is about \(11.5\) megabytes before temporary arrays are counted.

A diagnostic discrepancy between the destination-grid continuation mass and the CDF-based continuation probability is not automatically a model failure. It can indicate coarse spacing, insufficient coverage of the state space, inaccurate quadrature, or a boundary mask applied at the wrong nodes. Refine or extend the grid and check whether the estimates stabilize. Also verify that \(\Delta I_k>0\), the boundaries are ordered, the normal density includes its \(1/\sqrt{\Delta I_k}\) factor, and the mask retains the intended interior.

The main alternative is a scalar implementation that loops over source and destination nodes and accumulates the same weighted sums. It is easier to inspect on a tiny grid and can use less peak memory if values are generated as needed, but interpreted loops may be slower. A matrix implementation is often faster and its algebra is compact, while requiring attention to array orientation and memory use. Gaussian quadrature or other quadrature rules concern how integrals are approximated; they can be used with either vectorized or scalar arithmetic. Vectorization itself does not improve the quadrature rule or eliminate truncation error.

For very large grids, storing the full kernel may use substantial memory. Processing destination or source nodes in blocks preserves the same sums without holding every matrix entry at once. More specialized methods can exploit additional structure in equally spaced grids, but those methods require their own conditions and numerical checks. The direct matrix update remains a transparent reference calculation.

Common Pitfall: A Density Is Not a Probability Mass

A frequent error is to sum the entries of the propagated density vector and call that the survivor probability. The entries are density values, not probability masses; a sum only approximates an integral when it includes the appropriate destination weights. The correct grid estimate is \(\sum_jv_jc_jg_j\). Conversely, source values must be weighted before propagation, because the source integral is \(\sum_iw_if_k(x_i)K_{ij}\).

Another misleading check is to expect the sampled transition density to integrate to exactly one over a finite grid. A Gaussian transition may put probability outside the grid, and finite-grid quadrature is approximate. Use the CDF partition at the actual stopping boundaries to check event accounting, while remembering that it does not remove approximation in the source survivor density. Probability accounting is meaningful only relative to the mass that reached the transition.