Part XVII — Statistical Computing

Part XVI closed by naming what it had been avoiding: "every page here has either assumed the required integral was available in closed form or computed it on a grid small enough to enumerate, and both routes fail at the dimension a real hierarchical or state-space model occupies." This part is the machinery that removes the restriction — the optimizers, the quadrature rules, the expectation-maximization recursion and the chains whose stationary distribution is the posterior — and it delivers on that promise completely. A chain reproduces the posterior mean Part XVI computed on a \(601\times601\) grid, \(5.5229\) basis points against \(5.5687\), a miss of \(1.13\) Monte Carlo standard errors, and where importance sampling's effective sample size per draw collapses from \(0.692041\) to \(0.000052\) as dimension runs \(2\) to \(50\), the chain's falls only to \(0.007023\). What arrives with the machinery is a new kind of failure, and it is the same one in six costumes: every method on these pages reports a diagnostic computed from its own output, and every diagnostic is silent about exactly what the method missed.

The dependencies run in file order and split into two halves that meet at the end. Numerical Optimization and Numerical Integration are the two operations underneath everything else — maximizing and integrating — and they are a matched pair, the first bought against the curvature of a surface and the second against the smoothness of an integrand. The EM Algorithm is what happens when the function to be maximized is itself an integral over variables nobody observed, and it is the last page that returns a point. The three that follow return distributions instead. Markov Chain Monte Carlo carries the contract — invariance, discrete-time detailed balance, the ergodic average, integrated autocorrelation time and what \(\hat R\) can and cannot detect — while calling a kernel as a black box; Metropolis–Hastings builds the general kernel and prices the one parameter its correctness proof declines to constrain; Gibbs Sampling builds the kernel that never rejects and shows that the parameter did not disappear but moved to the moment the model was written. Where this part stops is worth stating: it builds no Hamiltonian sampler and no gradient-based proposal, derives variational inference no further than naming it as the E-step relaxed, runs no sequential Monte Carlo, which is Particle and Kalman Filters; it re-derives neither the dimensional cost of a grid nor the collapse of an envelope, which are Monte Carlo Simulation and Rejection Sampling; it specializes no recursion to a regime model, which is Hidden Markov Models; and it computes no risk number, prices no instrument and sizes no position, which is Part XVIII.

The failure running through the part is that every method here self-assesses, and the self-assessment is built from the same assumption as the answer, so the two fail together rather than independently. An optimizer reports success on \(0.7000\) of the starts that collapsed onto a zero-width spike against \(0.3172\) of the starts that found the honest optimum, and the highest likelihood over four hundred restarts belongs to a fit whose narrow component has standard deviation \(0.000000\) times the gap between two adjacent days. An adaptive integrator asked for a strategy's evidence returns \(0.000000\) with a reported error of \(7.27\times10^{-18}\), and then with a reported error of exactly zero, against a true value of \(0.833333\). Expectation-maximization records \(0\) decreases across \(2{,}291{,}200\) sweeps while \(0.0434\) of its runs climb toward a likelihood that is unbounded, so keeping the best of thirty-two restarts selects a degenerate fit on \(0.4925\) of datasets and puts the probability of a \(-3\%\) day at \(0.2949\%\) against a truth of \(1.2119\%\) — searching harder makes the answer worse. Four chains started from one mode of a bimodal posterior report \(\hat R=1.0001\) and an effective sample size of \(25{,}490\), eighty times the value from the harder case where mixing genuinely occurred, while the posterior mean is wrong by \(563.8\) Monte Carlo standard errors. An independence sampler accepts \(0.8744\) of proposals — a rate that reads as healthy — while one run sits frozen at a single value for \(9{,}708\) consecutive draws, its error bar \(3.74\) times too small and its nominal interval covering \(0.5433\). And a componentwise Gibbs sampler on a fifty-strategy hierarchy runs \(59\) times slower than an algebraically identical rewrite while \(\hat R\) separates them by four thousandths.

Topics

Topic	Focus
Numerical Optimization	A convergence flag as a statement about which stopping rule fired, the first-order condition being necessary rather than sufficient once concavity is dropped; the condition number as a property of the coordinates rather than the problem, a volatility span of \(3.1\) to \(87.5\) driving \(\kappa\) from \(11\) to \(9{,}363\) and gradient descent from \(85\) to \(77{,}561\) iterations while Newton takes \(2\) in every row and a diagonal rescaling requiring no new information restores \(\kappa\) to \(3.5\) and the count to \(30\); a two-component mixture from \(400\) starts reaching the honest optimum on \(0.9300\) at log-likelihood \(1584.4864\) and a \(1\%\) daily loss of \(-3.1147\%\) against a truth of \(-2.8774\%\), while \(0.0250\) collapse onto zero-width spikes carrying a higher likelihood of \(1592.1922\) and a loss of \(-2.4013\%\), with success reported on \(0.7000\) of the spikes against \(0.3172\) of the honest fits; and the U-shaped finite-difference error floors of \(8.410\times10^{-08}\) and \(1.477\times10^{-11}\) degrading to a relative error of \(1.248\times10^{+06}\) on a simulated objective, where BFGS never takes a step and forfeits \(11.5932\%\) of annualized growth until common random numbers freeze the draw and the shortfall falls to \(0.1985\%\)
Numerical Integration	Order as a claim about the integrand's derivatives rather than a property of the rule, an \(n\)-point Gauss rule buying degree \(2n-1\) against a \(2n\)-th derivative a payoff kink does not possess; the trapezoid rule reaching \(6.537\times10^{-13}\) at \(32\) panels on a smooth expectation against Simpson's \(5.511\times10^{-08}\) because every Euler–Maclaurin boundary term vanishes, while one kink at the strike collapses all three rules to observed order \(2.00\) and inverts the ranking, Gauss–Legendre's \(1.664\times10^{-01}\) trailing the trapezoid rule's \(7.113\times10^{-03}\) by a factor of \(23\), and a jump takes every rule to order \(1.00\); splitting at the strike restoring Simpson to order \(4.00\) and \(4.669\times10^{-10}\) while \(32\) split Gauss nodes beat \(2048\) unsplit ones by four orders of magnitude, against a split misplaced by \(0.1280\) of a panel width surrendering \(0.6565\) of the benefit; and scipy.integrate.quad returning \(0.828159\) with a reported error understating the truth by \(1.50\times10^{+06}\), then \(0.000000\) with reported errors of \(7.27\times10^{-18}\) and exactly zero, against a true evidence of \(0.833333\)
The EM Algorithm	The ELBO/KL decomposition making ascent a property of the bound rather than of the model, so a merely-increasing M-step inherits it intact; linear convergence at the spectral radius of \(I_{\mathrm{mis}}I_{\mathrm{com}}^{-1}\), the rate read off the error decay matching the Jacobian of the EM map at \(0.0656\) against \(0.0652\) and then exactly at \(0.2858\), \(0.6021\), \(0.7949\) and \(0.8662\), with sweep counts of \(8\), \(17\), \(43\), \(94\) and \(147\) against predictions of \(6.3\), \(15.3\), \(40.4\), \(91.5\) and \(147.5\); Louis' identity making the complete-data standard error too small by \(\sqrt{1-r}\), verified two independent ways to four decimals, the naive figure flat at \(5.3666\) basis points where the truth reaches \(6.9912\) and coverage falls to \(0.8686\) against the observed information's \(0.9496\); and monotone ascent holding across \(2{,}291{,}200\) sweeps with \(0\) decreases while \(0.0434\) of starts climb an unbounded likelihood, degenerate selection rising \(0.0400\), \(0.0600\), \(0.1325\), \(0.2275\), \(0.4925\) with the restart count and the probability of a \(-3\%\) day falling from \(0.6449\%\) to \(0.2949\%\) against a truth of \(1.2119\%\), until a fifteen-basis-point variance floor removes it entirely
Markov Chain Monte Carlo	Invariance as the whole contract and detailed balance a sufficient rather than necessary route to it, so a correctness proof transfers no information about a finite run; a chain reproducing Part XVI's independently published grid answer at \(5.5229\) basis points against \(5.5687\), a miss of \(1.13\) Monte Carlo standard errors, with \(\hat R=1.0001\) and an integrated autocorrelation time of \(6.87\) converting \(120{,}000\) draws into \(17{,}474\); the trade forced by dimension, importance sampling's effective sample size per draw collapsing \(0.692041\), \(0.398052\), \(0.158109\), \(0.009939\), \(0.000052\) as \(d\) runs \(2\) to \(50\) with its tail estimate decaying to \(0.0285\) against a truth of \(0.0500\), where the chain's falls only to \(0.007023\), the estimate holds at \(0.0464\), and acceptance settles at \(0.2406\) unbidden; and four chains started from one mode of a bimodal posterior returning \(\hat R=1.0001\) and an effective sample size of \(25{,}490\) — eighty times the value from the harder row where mixing genuinely occurred — with a posterior mean wrong by \(563.8\) standard errors that dispersed starts expose at \(\hat R=3.0048\) without correcting
Metropolis–Hastings	Detailed balance holding for any proposal with adequate support and the target entering only through a ratio, so the normalizing constant Part XVI called the only expensive step is never computed; correctness uniform across tuning where cost is not, seven scales spanning five orders of magnitude all landing within \(2.37\) Monte Carlo standard errors of the grid answer while effective sample size per draw ran \(0.00021\) to \(0.14598\), a factor of \(768\), the two worst settings failing oppositely at acceptance rates of \(0.9948\) and \(0.0000\); optimal scaling reproduced, acceptance at the best step falling \(0.4429\), \(0.3595\), \(0.2903\), \(0.2611\), \(0.2456\), \(0.2398\) toward \(0.234\), and an isotropic proposal on a target of condition number \(1{,}000\) delivering \(0.000638\) against a preconditioned \(0.028584\), a factor of \(44.8\) recovered by the change of variables page 1 used on an optimizer; and a proposal with lighter tails than the target accepting \(0.9328\) and \(0.8744\) while freezing for \(1{,}765\) and \(9{,}708\) consecutive draws, error bars \(1.27\) and \(3.74\) times too small and coverage of \(0.8633\) and \(0.5433\), the violated condition being the \(\sup p/q<\infty\) of importance sampling arriving as a statement about time
Gibbs Sampling	A full-conditional update as a Metropolis step whose acceptance probability is identically one, and a composition of invariant kernels being invariant while a composition of reversible ones need not be, so the systematic scan is correct and not reversible; on the semi-conjugate model Part XVI predicted, acceptance of \(1.0000\) and \(0.9796\) effective draws per draw against a tuned random walk's \(0.1209\), because the posterior correlation is \(-0.0206\); the cost being that correlation exactly, \(\tau=(1+\rho^{2})/(1-\rho^{2})\) measured at \(1.03\), \(1.69\), \(9.73\), \(89.87\) and \(733.92\) against \(1.00\), \(1.68\), \(9.45\), \(107.37\) and \(1067.15\), so two factors correlated at \(0.999\) need \(25{,}292{,}185{,}213\) draws for a hundredth of a basis point where blocking needs \(34{,}847{,}581\); and an algebraically identical rewrite of a fifty-strategy hierarchy taking the effective sample size from \(644\) to \(37{,}735\) while \(\hat R\) moves only from \(1.0042\) to \(1.0000\), with every full conditional capable of being proper while the joint distribution does not exist at all