BUG: compute Probit loglike, score and hessian in log space - #10229

• closes BUG: Probit loglike is clipped flat and its Hessian is nan for extreme linear predictors #10228
• tests added / passed.
• code/documentation is well formatted.
• properly formatted commit message. See
  NumPy's guide.

AI Disclosure

Contributions must comply with the statsmodels AI Policy.

• No AI tools were used to develop this pull request.
• AI tools were used.
  Tool(s): Claude (Claude Code).
  Used for: locating the bug (probing Probit against scipy references after the Logit BUG/ENH: overflow in Logit #3923 fix), drafting the implementation and the two regression tests, running the verification (scipy references, complex-step and finite-difference derivative checks, min-version scipy 1.8.1 check), and drafting this description and the linked issue.
  I have personally read, understood, and can explain every line of this diff, and
  have verified the statistical/numerical correctness of the change.

Summary

Probit evaluated its log-likelihood as log(clip(cdf(z), FLOAT_EPS, 1)), which is flat at log(eps) = -36.04 for q * XB below about -8 while the true log(Phi(z)) keeps falling (-53 at z = -10, -805 at z = -40). The score used the same clipped cdf in the Mills ratio pdf / cdf, where pdf underflows first, so extreme observations contributed ~0 score instead of a value growing like |z|; hessian and hessian_factor did not clip at all and returned nan below z = -38. With a single extreme misclassified observation in the sample this gave a wrong llf, a wrong-signed score, a nan Hessian, and fit(method="bfgs") stopping at a spurious optimum with nan bse (details and repro in #10228). This is the Probit counterpart of the Logit overflow fixed in gh-3923.

The fix computes everything in log space with the exact scipy pieces:

• _norm_logcdf(z) = log_ndtr(z) replaces log(clip(cdf)) in loglike and loglikeobs.
• _norm_mills(z) = exp(norm._logpdf(z) - log_ndtr(z)) replaces the clipped pdf / cdf in score, score_obs, score_factor, hessian and hessian_factor, so the inverse Mills ratio stays finite and accurate where cdf underflows (it behaves like -z there).
• _deriv_score_obs_dendog uses exp(logpdf - log_ndtr(z) - log_ndtr(-z)) for pdf / (cdf * (1 - cdf)), via 1 - Phi(z) = Phi(-z).

log_ndtr accepts complex input (checked on scipy 1.8.1, the minimum supported version), so complex-step differentiation of the new loglike still works and agrees with the analytic score. In the normal range the new values agree with the old ones to ~1e-14, and all existing Probit tests pass unchanged.

Tests

test_probit_loglikeobs_extreme_linpred mirrors the gh-3923 Logit test (z from -1e4 to 1e4: finite loglikeobs, score_obs, hessian; agreement with scipy.stats.norm.logcdf at 1e-14 and with the Mills ratio reference at 1e-12). test_probit_extreme_observation_fit uses the outlier sample: exact loglike/score/hessian against scipy references, complex-step and finite-difference derivative checks, and both newton and bfgs converging to the same optimum as a direct scipy.optimize.minimize of the exact log-likelihood with finite bse. Both tests fail on main and pass with the fix. The reference values are scipy's log_ndtr/norm.logcdf and norm.logpdf, which are the same quantities R's pnorm(log.p = TRUE) and dnorm(log = TRUE) compute.

pytest statsmodels/discrete/tests and statsmodels/tools/tests/test_numdiff.py pass; ruff check is clean. No release note added since the recent BUG merges did not touch docs/source/release/; happy to add one wherever you want it.