BUG: Probit loglike is clipped flat and its Hessian is nan for extreme linear predictors #10228

Describe the bug

Probit computes its log-likelihood as log(clip(cdf(z), FLOAT_EPS, 1)), so every observation with q * XB below about -8 contributes the constant log(eps) = -36.04 instead of log(Phi(z)), which keeps falling (-53.2 at z = -10, -804.6 at z = -40). The score uses the same clipped cdf in the inverse Mills ratio pdf / cdf, and pdf underflows to 0 first, so the score of those observations is ~0 where the true value grows like |z|. hessian and hessian_factor do not clip at all, so pdf / cdf becomes 0 / 0 = nan below about z = -38.

Practical consequences with one extreme, misclassified observation in the sample:

• reported llf (and therefore AIC/BIC) is wrong away from the optimum: -155.9 where the exact value is -1924.9 in the example below,
• the score has the wrong sign on the slope, so gradient-based optimizers can stop at spurious points: fit(method="bfgs") stops at params [-0.036, 1.72] with llf = -155.1 and bse = [nan, nan], while the exact optimum is [-0.129, 0.077] with llf = -202.7,
• the Hessian is nan at any params where some q * XB < -38.

Logit had the analogous overflow fixed via logaddexp (gh-3923, and the recent _log_1pexp follow-up); the clipping in Probit was noted in that thread as something that "could be put into the cdf method directly" but was never replaced. scipy provides the exact, stable pieces: log_ndtr for log(Phi(z)) and exp(logpdf - log_ndtr) for the Mills ratio, and log_ndtr also accepts complex input so complex-step differentiation keeps working.

Code Sample, a copy-pastable example

import numpy as np
from scipy import stats
import statsmodels.api as sm

rng = np.random.default_rng(1)
n = 300
x = rng.standard_normal(n)
y = (1.5 * x + rng.standard_normal(n) > 0).astype(int)
x = np.append(x, 40.0)   # one extreme, misclassified observation
y = np.append(y, 0)
X = np.column_stack([np.ones(n + 1), x])

model = sm.Probit(y, X)
params = np.array([0.0, 1.5])
q = 2 * y - 1
z = q * (X @ params)

print(model.loglike(params), stats.norm.logcdf(z).sum())
# -155.934  vs exact  -1924.904
print(model.score(params))
# [-4.874  7.393]  vs exact  [-64.891  -2393.273]
print(np.isnan(model.hessian(params)).any())
# True
res = model.fit(method="bfgs", disp=0)
print(res.params, res.bse, res.llf)
# [-0.036  1.72]  [nan nan]  -155.067   (exact optimum: [-0.129  0.077], llf -202.674)

Expected Output

loglike, score and hessian equal the exact quantities (log_ndtr based), and fit converges to the same optimum as a direct minimization of the exact log-likelihood with finite standard errors, as Logit does on the same data.

Actual Output

See the inline comments in the code sample: loglike -155.934 vs exact -1924.904, score [-4.874, 7.393] vs exact [-64.891, -2393.273], nan Hessian, and fit(method="bfgs") returning params [-0.036, 1.72], bse = [nan, nan], llf = -155.067 where the exact optimum is [-0.129, 0.077] with llf = -202.674.

Versions

statsmodels 0.15.0 and current main (2fedaa9), numpy 2.x, scipy 1.x, Python 3.12, macOS arm64.

I have a fix ready that replaces the clipped cdf in loglike, loglikeobs, score, score_obs, score_factor, hessian, hessian_factor and _deriv_score_obs_dendog with log_ndtr based computations; PR to follow.

Checklist

• I have searched the issue tracker for a similar issue and did not find one.
• I have confirmed this bug exists on the latest released version of statsmodels, or on the main branch.

AI Disclosure

AI tools were used (Claude / Claude Code) to locate this bug by probing Probit against scipy references, to produce the repro, and to draft this report. The numbers above were run and checked by me against scipy.stats.norm.logcdf and a direct scipy.optimize.minimize of the exact log-likelihood.