Module multinomial_predictive 

Posterior-predictive class probabilities for the penalized multinomial logit, computed as a RATIO OF NORMALISING CONSTANTS.

Why this exists rather than integrating the Gaussian posterior

The published estimand is the posterior mean probability E[softmax(x'β) | data]. The obvious implementation — approximate the posterior of β by the Laplace Gaussian N(β̂, H⁻¹) and integrate softmax against it — is not a valid approximation of that estimand, and the failure is not small.

Write the posterior mean of any positive functional g(β) as a ratio of integrals. Laplace applied SEPARATELY to numerator and denominator (the “fully exponential” form of Tierney and Kadane) has its O(n⁻¹) errors cancel between the two, leaving O(n⁻²). Integrating g against the Gaussian instead keeps only the CURVATURE half of the O(n⁻¹) correction (½ tr(H⁻¹ ∇²g)) and silently drops the SKEWNESS half, which comes from the third derivative of the log-posterior. On a well-conditioned fit the two halves are both small and nobody notices. On a (quasi-)separated multinomial they are not small and they have opposite signs: the likelihood is flat toward more separation and steep away from it, so the true posterior is strongly skewed toward LARGER |η|, while the symmetric Gaussian puts half of its mass on the side the likelihood has already excluded. softmax is concave along the winning coordinate, so that misplaced mass converts directly into under-confidence: right argmax, flattened probabilities.

For g = p_c(x) = P(new row at x is class c | β) the ratio is not merely a device — it is exactly the posterior predictive, because the extra row’s likelihood factor IS the functional being averaged:

    E[p_c(x) | D]  =  Z(D ∪ {(x, c)}) / Z(D)

with Z the posterior normalising constant. Approximating each Z by Laplace at its own mode gives

    E[p_c(x)] ≈ exp( L⁺(β̂⁺) − L(β̂) ) · sqrt( det H / det H⁺ )

where L is the penalized log-posterior, β̂⁺ the mode with the extra row present, and H, H⁺ the corresponding negative Hessians. The (2π)^{d/2} factors cancel exactly (same dimension on both sides).

The identity Σ_c E[p_c(x)] = 1 is exact for the true integrals, so the deviation of the computed Σ_c from one is a MEASURED accuracy statement about this approximation, available at every prediction row and requiring no reference. MultinomialPredictiveModel refuses rather than publishing a row whose mass defect exceeds PREDICTIVE_MASS_DEFECT_TOLERANCE.

The same machinery supplies the second moments the standard-error surface consumes, with two extra rows instead of one:

    E[p_c(x) · p_d(x)]  =  Z(D ∪ {(x, c), (x, d)}) / Z(D)

Cost

One warm-started Newton solve per (row, class) — the augmented objective is strictly convex, so Newton with backtracking is unconditionally safe — plus K(K+1)/2 more per row when second moments are requested. Each Newton iteration is O(n·M²·P²) for the curvature (as M(M+1)/2 GEMMs) and O(d³) for the factorisation, so the whole predictive is O(R·K·iters·(n M² P² + d³)).

On the fixture this exists for that is a large improvement, not a cost: the Smolyak integrator it replaces spent ~930 s on one penguins prediction block, because its level requirement grows with exactly the posterior width that makes the Gaussian wrong in the first place, while the same block here is n = 228, P = 37, M = 2 — under a second. The scaling is different in kind, though, and worth stating plainly: this method’s cost grows with the TRAINING size, which the Gaussian route’s did not, because evaluating a posterior away from its mode is what the Gaussian route was avoiding by being wrong.

Structs

MultinomialPredictiveModel

The training data and penalty a saved multinomial model needs in order to evaluate its own log-posterior away from the mode.

MultinomialPredictiveMoments

Posterior-predictive moments at a block of prediction rows.

Constants

PREDICTIVE_MASS_DEFECT_TOLERANCE

Largest tolerated deviation of Σ_c E[p_c(x)] from one before a prediction row is refused.

Functions

predictive_standard_deviation

Per-class posterior standard deviation of the probability, from the moments above: sd(p_c) = sqrt(E[p_c²] − E[p_c]²).