NTN: fix the lambda-floor cancellation and the sub-NHN fit (#55)

[davidhbernstein]

Fixes #55. Both defects in that issue, on one branch, because shipping the first
alone would ship a measured regression — see "Why both halves" below.

1. The log-likelihood cancelled at the lambda floor

l4 + l5 was log Phi(aa) - log Phi(bb) with

aa = (mu/lam - eps*lam)/sig      bb = (mu/sig)*sqrt(1 + lam^-2)

Both arguments carry the same mu/lam. With mu < 0 and lam at its floor
each term reaches about -1e13 while the difference is O(1), so the reported
value was the rounding error of two enormous numbers of opposite sign. Nothing
is wrong with either pnorm() call.

The signature is unmistakable. On seed 1011401 / y_pcs_st, where the fit
reached mu = -32578, sig = 213, the per-observation error came out as an
exact multiple of 512 — the floating-point spacing of aa^2/2 at 2.3e18 —
and the total over 200 observations was exactly 200 * 512 = 102400.

aa^2 - bb^2 = (-2*mu*eps + eps^2*lam^2 - mu^2)/sig^2 analytically, with
mu^2/lam^2 cancelling exactly; against l3 the remaining mu^2 cancels too,
leaving

l3 + l4 + l5 = -eps^2 (1 + lam^2)/(2 sig^2) + tilt(aa) - tilt(bb)

Since sig_v^2 = sig^2/(1 + lam^2) that leading term is -eps^2/(2 sig_v^2) —
structurally the residual #28 left behind in TSL, which is why the two fixes
look alike. Taken through .log_phi_tilt(), as "NE", "NGE" and "TSL"
already do.

Only when both arguments are negative. With mu > 0 both go to +Inf,
log Phi of each goes to 0, and the original expression has no cancellation —
there tilt(aa) - tilt(bb) would be a difference of two numbers near 5e15 and
would manufacture the problem. That branch keeps the direct path, and a test
pins it to the pre-fix expression to 1e-12.

Verified two independent ways, neither reusing the package's algebra

1. The analytic lam -> 0 limit. As lam -> 0, bb/aa -> 1, so
   tilt(aa) - tilt(bb) -> 0 and the density collapses to a normal: the limit
   is -n log sig - n log(2pi)/2 - sum(eps^2)/(2 sig^2), in closed form. The
   approach rate is the real test — the gap must fall like lam^2:

   lam	1e-3	1e-4	1e-5	1e-6	1e-7	1e-8
   before	2.4e-6	3.3e-8	7.1e-7	1.0e-4	5.2e-3	6.3e-1
   after	2.4e-6	2.4e-8	2.4e-10	2.4e-12	2.3e-14	2.2e-16

   The old form tracks the limit to 1e-4 and then turns around.

2. Level-space log(Phi(aa)/Phi(bb)), valid while neither Phi underflows.
   On 116 grid points in (lam, mu, sig) where it applies, old and new both
   match it to better than 1e-9 — which is what rules out the rewrite having
   broken the benign regime.

The four samples in the issue reproduce exactly, including its independently
computed Mills-ratio values (-372.7305, -374.8290, -371.1794).

2. The fit could end far below the nested NHN

At mu = 0 the NTN likelihood is NHN's with the same (lambda, sigma), and
mu = 0 is interior — start_cs() gives NTN
lower_bob <- c(rep(.Machine$double.eps, 2), rep(-Inf, n_x_vars + 1)), so only
lambda and sigma are bounded. So a maximum below NHN's is a failure, not a
limit that is merely approached, unlike TSL's OLS supremum in #28 or the
sigma_v boundary in #56.

Fixed on the pattern #45 built for NG and #30 for NNAK: fit NTN's likelihood
with mu held at 0 from NHN's own starting vector plus five variance
splits, then check the finished fit against that point and polish from it. NTN
is added to the existing c("NG", "NNAK") list and .nhn_ref to the references
the loop already walks.

It is an end check, not a start, for the reason that block already records
for NG: as a candidate it moves where the stages go on samples that were never
in trouble.

Measured on 2700 fits

N = 200, seeds 1000001–1000150, all 18 y_pcs_* columns, against main
1edf823:

	main	this branch
reported above the true value at its own optimum	40 (max 1.2978)	0
on the lam floor and above OLS — impossible	9	0
below the nested NHN by >1e-3	200 (worst 764)	13 (worst 0.39)

Scoring both builds' optima under one correct likelihood: this branch is
higher on 278, lower on 27 (worst 3.78), equal on 2395 — +12123 in total
log-likelihood.

The 13 that remain below NHN are all y_pcs_ez with mu = 0 and lambda
between 1320 and 6140 — the flat sigma_v boundary of #56/A59, where one polish
stops short on a ridge. Reported rather than tuned away.

A scoping note, since it is easy to get wrong: 2090 of the 2700 fits are above
OLS at some lambda, which is entirely legitimate and measures nothing. Only
"on the floor and above OLS" is a defect.

Why both halves are in one PR

The cancellation fix alone is a clear net improvement (higher on 124, lower on
73, +4818 total) but it perturbs the optimizer's path, and 73 fits ended lower
— worst -30.58, and verified genuine by scoring main's own optimum under
the corrected likelihood, not an illusion being removed. The nested check
absorbs that tail (73 → 27, worst -3.78). Shipping part 1 alone would knowingly
ship that regression.

mu -> -Inf is a supremum — an extra fix measured and then reverted

Worth recording because it looked right. The mu = 0 reference sits at exactly
0, so lower.start(differ = 0.5) boxes mu into [-0.5, 0.5] and it stops on
the edge; three of the four issue samples came back at exactly mu = -0.5000.
Re-centring gained 0.33–1.25 more. mu then stopped at exactly -2.5 (five
passes of 0.5). Doubling the window gave -7.5, -15.5, -31.5 — every one an
exact window boundary.

A continuation profile holding mu fixed at each rung settles it:

       mu     profile ll      lambda       sigma    gain over previous
     -2.5    -369.802659      2.0784      2.6420       +0.650575
     -7.5    -369.041415      2.8787      3.4986       +0.761244
    -31.5    -368.554971      5.2054      6.1128       +0.486444
     -100    -368.438577      8.9765     10.4371       +0.116394
    -1e+03   -368.391070     27.9802     32.3861       +0.012217
    -1e+04   -368.386385     88.3168    102.1952       +0.001220

Gains fall geometrically to a finite limit near -368.3862 with lambda
and sigma growing in proportion to |mu|: a ray to infinity along which the
supremum is approached and never attained — the same kind of object as #56's
boundary. So a widening window is the wrong instrument: there is no point to
walk to, and the fit would stop wherever the cap fell, reporting arbitrarily
large |mu| with no standard errors. Reverted; what ships is the one-pass check.

(An intermediate build with a worse likelihood had wandered onto that same ray
and returned -368.386084 on this sample from mu = -23530. The shipped
likelihood reproduces that value exactly at the same parameter vector — a wrong
likelihood found a right point by accident, and it is the number the profile
converges to.)

Also in this PR

• NTN now also gets the A25 "never return worse than your own start" check,
  because it joins the loop that performs it. Monotone in the objective.
• tests/testthat/test-issue55-ntn-lambda-floor.R: 23 assertions fail on
  main, all pass here. Covers the limit and its O(lam^2) rate, the
  cannot-beat-OLS bound near the floor, the mu > 0 path being unchanged, the
  NTN >= NHN floor, the mu = 0 identity against a real NHN fit, mu = 0
  being interior, and every NTN parameter layout (uhet, muhet, both, no
  intercept, numeric start_val) since .nhn_ref inserts into a
  variable-length vector.
• test-vcov-bhhh-advice.R reworked — see the section below.

Recorded outside the repository, in the research workspace (so neither the
tarball nor this PR carries them): the long-form derivation in
notes/code_history/sfm.md, three evidence scripts in horserace/, and gap
A61 done in horserace/FUNCTIONALITY_GAPS.md — whose status counts were
stale and whose own documented recount command was wrong, since
^### [A-Z][0-9]+\. misses A56a, A52a, G2 / I3 and H6 / I5. Command
widened and the numbers rerun.

Cost

An NTN fit goes from 0.084 to 0.264 s (3.2x), because .nhn_ref costs six
bounded optim() calls. NG pays four for the same pattern, so it is in
precedent. The full suite is 35 min against a 39 min baseline, i.e. the
cost is not visible there: FAIL 0 | WARN 43 | SKIP 8 | PASS 5368.

The seed budget was measured rather than assumed, and the first measurement was
misleading. Scoring .nhn_ref in isolation, the NHN start alone is already best
68% of the time and the five splits add a median of 0.0000 and a max of 0.1466 —
which reads as "drop them". But what matters is the FINISHED fit. Re-fitting the
200 samples that were below NHN on main with a 4-seed build (NHN start +
w = 0.3, 0.5, 0.7, matching NG's budget):

	below NHN	worst	vs the other build
6 seeds (ships)	13 of 200	-0.3938	higher on 26
4 seeds	13 of 200	-0.3976	higher on 12

Identical on 162. So the splits do earn their place on the subset the check
exists for, and they stay.

One existing test had to be adapted, and why that was its fault not this PR's

test-vcov-bhhh-advice.R (from #54) failed here. Worth reading carefully,
because the invariant it guards was never in danger.

It triggers the "bread is undefined" path with a real degenerate NTN fit on
rand = 2. On this branch that fit slides from lambda = 9.7e5 to 4.4e6 --
a log-likelihood move of 9e-5 along the flat ridge it already sat on.
solve(hessian) still fails, every standard error is still NA, and
vcov(type = "bhhh") still fails. But the sandwich bread became invertible,
so vcov(type = "sandwich") returned a matrix instead of erroring and the
expect_s3_class(e, "error") line failed. The trigger stopped triggering; the
agreement it then asserts was fine.

That makes the fixture too brittle for the job -- any correct change to NTN can
nudge it. The test now searches rand = c(2, 4, 5, 12) for a fit whose
sandwich path actually errors, asserts the invariant on every one that does, and
ends with expect_gt(exercised, 0) so a search that found nothing cannot make
the assertions vacuous. It also still requires all four candidates to be the
degenerate kind (lambda > 1e4, all SEs NA).

Verified both ways: 34 assertions pass on this branch, 40 on main 1edf823
(more candidates reach the error path there). The redundant first clause of
claims_available was dropped -- it was ||-ed with the broader
grepl("is defined here") that follows it.

Left open deliberately

• The 13 residual below-NHN fits on the sigma_v boundary (NTN: warn and flag when sigma_v collapses to the boundary #56/A59 territory).
• .ne_ref holds NG's shape at exactly 1, so NG's end-check box is centred on a
  constrained coordinate the same way NTN's was, and may be pinned at 0.5 or
  1.5. Unmeasured; not touched here.
• Issue NTN: log-likelihood cancels at the lambda floor (reports above OLS), and the fit can end far below the nested NHN #55 names its third sample y_pcs_tht. There is no such column in
  data_gen_cs(); the one reproducing its figures is y_pcs_st.

🤖 Generated with Claude Code