[alejandroschuler]

Problem

LogitNormal is built with make_distribution without scipy_dist_cls or mean_fn, so it got the factory's Monte Carlo mean, np.mean(self.sample(1000), axis=0). As a result, predict() was random. Two calls on one fitted model (100 rows, 5 trees) gave values up to 0.05 apart. Pickle round-trip checks and other comparisons of predictions were also unreliable (#445 compares predicted parameters for this reason).

Change

LogitNormal now passes a deterministic mean_fn. The mean of expit(mu + sigma * Z) with Z ~ N(0, 1) has no closed form. The new function uses one of two quadrature rules:

• For sigma <= 1, 64-node Gauss-Hermite quadrature over Z (hermegauss).
• For sigma > 1, E[Phi((mu - L) / sigma)] with L standard logistic, integrated with the trapezoidal rule (step 0.5 on [-40, 40]).

Gauss-Hermite alone is not accurate enough for large sigma. With 64 nodes the error is 3e-3 at sigma = 10 and 3e-2 at sigma = 30, because the integrand is close to a step in Z. The logistic form becomes more accurate as sigma grows, so the two rules complement each other.

check	result
largest error against an mpmath reference, sigma 0.01 to 10, mu -10 to 10	6e-15
same, sigma up to 1000	below 1e-14
2e6 Monte Carlo draws per cell, 294 cells	largest |z| = 2.9
time for 1e6 rows	1.5 s

The function loops over the nodes and is vectorized over the rows, so its memory grows only with the number of rows.

No other built-in distribution uses the Monte Carlo fallback. Beta and BetaBinomial use scipy, BetaBinomial and BetaBinomialEstN pass mean_fn, and BetaBernoulli is a classifier. The fallback stays for user distributions. The make_distribution docstring and docs/sympy_distributions.md now say that this mean is random and that mean_fn makes it deterministic.

Tests

tests/test_logitnormal.py (15 tests):

• two predict() calls on a fitted model give equal results;
• the mean agrees with 1e6 Monte Carlo draws to within 5 standard errors, for sigma in {0.01, 0.1, 0.5, 1, 2, 5, 10} and 9 values of mu in [-10, 10];
• the mean agrees with scipy.integrate.quad to 1e-10 on the same grid, and a one-row slice gives the same mean as the full batch.

The old Monte Carlo mean fails all 15 tests. A Gauss-Hermite-only mean fails 3 (sigma = 5 and 10).

pytest tests --slow: 184 passed, 2 xfailed, 22 xpassed. black, isort, flake8 and pylint on ngboost pass. pylint tests fails with R0801 (duplicate-code) between tests/test_helpers.py and tests/test_pickling.py. Both files are unchanged here, so the failure is already on master.

🤖 Generated with Claude Code