[alejandroschuler]

Fixes the crash reported in #275, which was closed in 2021 without a fix.

Problem

NGBRegressor(Dist=T) with default settings fails with ValueError: Input y contains NaN after a few dozen iterations:

import numpy as np
from ngboost import NGBRegressor
from ngboost.distns import T

rng = np.random.default_rng(0)
X = rng.uniform(-1, 1, size=(500, 3))
y = np.sin(np.pi * X[:, 0]) + X[:, 1] + 0.5 * rng.standard_t(5, size=500)
NGBRegressor(Dist=T).fit(X, y)  # ValueError: Input y contains NaN

It failed in 15 of 15 fits that I tried (5 seeds, Y in units of 1, 100 and 1e4), at iterations 26 to 60. Scaling Y does not help. It is not a regression: v0.4.0 and v0.5.2 fail the same way. tests/test_distns.py marked this case as a known failure.

Cause

The gradients in TLogScore.d_score are correct. The problem is how log df behaves as df grows.

• The Fisher information for log df goes to zero as df grows. After loc and log scale are accounted for, it is about 1.5 / df². So the natural gradient for log df on one row has a size of about 0.8 df.
• The loss is almost flat toward large df, because a t with a large df is a normal. The line search stops a step that sends df toward 0 on one row, but not a step that sends it up.
• The default tree can put one row in a leaf of its own. A row whose df has drifted up gets a step of about learning_rate × df. That raises df, which makes the next step larger. In one trace, the largest log df on a row went 4.5, 5.0, 6.4, 10.4 and 288 in five iterations.
• When log df passes 709, np.exp(params[2]) overflows. The gradient is then NaN, and the tree fit fails.

The 100-draw Monte Carlo metric that T used made this worse. For each row it underestimated the information for log df, by a factor of 1.4 to 2.8 at the median.

Changes

• Cap on df. T.max_df = 1000. T.__init__ clips log df before exp. A t with df 1000 is very close to a normal. A row at the cap gets no push further up, and its log df is not coupled to log scale in the metric, so the fit can still bring df down. Without this, the raw parameter of capped rows kept rising (to about 1400 in 500 trees) and they stayed stuck at the cap. On Cauchy noise that cost 0.43 nats of test NLL (table below).
• Exact metric. TLogScore.metric is the Fisher information of (loc, log scale, log df) from Lange, Little and Taylor (1989). It agrees with quadrature of d_score to 1e-6 or better. Fits are about three times faster (2.7 s against 8.2 s for 500 trees on 500 rows). Fits with the same random_state are now identical. Before, the metric drew its samples from the global numpy random state.
• Start of df, unchanged. T.fit still starts every row at df = 3, now through T.init_df with a comment. I tried fitting df to Y with scipy.stats.t.fit, and it is worse. Y mixes the signal with the noise, so its tails say little about the noise. On these data scipy returned df = 1.6e10, every row started at the cap, and the fitted df stayed far above the truth (median 182 against 5 on t₅ noise). It also gave scipy RuntimeWarnings inside t.fit. See "fitted start" in the second table.
• Docstrings. The docstrings of T, TFixedDf and TFixedDfFixedVar called loc and scale the mean and standard deviation. They now say that loc is the mean (when df > 1) and give the standard deviation.
• Release notes.

Results

Data:

• Synthetic: X ~ U(-1, 1)³, mean sin(πx₀) + x₁, noise 0.5 × {t₅, N(0, 1), Cauchy}, 500 training rows and 2000 test rows.
• California housing: the first 1000 rows, with an 80/20 split by seed, which is the setting of tests/test_distns.py.

Every model uses the library defaults, because the bug is in the default fit. There are 10 seeds for each cell. The finish counts pool Y in units of 1 and 1e4 for the synthetic data. Test NLL is the mean over seeds at units of 1. The numbers in brackets are the mean paired difference from this PR, ± its Monte Carlo SE (positive means worse than this PR).

data	fits that finish, master / this PR	median fitted df (true)	test NLL, this PR	TFixedDf	Normal	true model
t₅ noise	0/20 / 20/20	3.6 (5)	1.185	1.121 [−0.063 ± 0.037]	1.387 [+0.20 ± 0.05]	0.938
normal noise	0/20 / 20/20	3.5 (∞)	0.985	0.938 [−0.047 ± 0.021]	1.106 [+0.12 ± 0.04]	0.727
Cauchy noise	0/20 / 20/20	2.0 (1)	2.133	2.231 [+0.098 ± 0.033]	679	1.851
California	1/10 / 10/10	3.0	0.627	0.447 [−0.18 ± 0.07]	0.917 [+0.29 ± 0.09]

Every fit now finishes. This PR beats Normal on all four data sets. It beats TFixedDf on Cauchy noise. It is a little worse than TFixedDf on t₅ and normal noise, and worse on California (see "Limits" below).

The design choices, with the same data and seeds:

variant	t₅	normal	Cauchy	California
this PR	1.185	0.985	2.133	0.627
cap only, rows above the cap still pushed up	1.122 [−0.063 ± 0.044]	0.951 [−0.034 ± 0.044]	2.558 [+0.43 ± 0.15]	0.637 [+0.01 ± 0.07]
fitted start (df from scipy.stats.t.fit)	1.428 [+0.24 ± 0.04]	1.346 [+0.36 ± 0.03]	2.223 [+0.09 ± 0.11]	0.807 [+0.18 ± 0.07]

Limits

T still learns more slowly than TFixedDf. On seed 0 of the t₅, Cauchy and California data, the median line-search scale is 0.016 to 0.5 for T and 1.0 for TFixedDf, and 27% to 87% of the T iterations take a scale below 0.5. The cause is the same scaling as above. The line search uses one scale for all parameters. A candidate step that is large in log df on one row makes the loss worse, so the step shrinks for loc and scale too. In 500 trees df moves only a little from its start at 3 (median 3.5 on normal noise). A fix is to put the tail on a scale with bounded information near the normal limit, such as 1/df, for which the information tends to 1.5. That changes the meaning of the third parameter, so I left it for a separate PR.

Trial steps that the line search rejects can still give a few scipy RuntimeWarnings (divide by zero encountered in log), when a trial step sends df close to 0. The fit does not change. #434 silences warnings from rejected trial steps.

Interaction with other open PRs

I ran trial merges of this branch with each of them.

• Mark only the T case as xfail in the logscore distribution test #437 merges cleanly, but the result fails CI. It marks the T case of test_dists_runs_on_examples_logscore as a strict xfail. This PR makes that case pass, so a strict xfail fails. If Mark only the T case as xfail in the logscore distribution test #437 merges first, the T marker has to go when this PR is merged. If this PR merges first, Mark only the T case as xfail in the logscore distribution test #437 should drop its T marker (its fix of the always-true condition still applies).
• Deprecate the var attribute of the Normal and T distributions #439 renames self.var to self._var. It conflicts with this PR on adjacent lines of T.__init__. Keep both changes. The new metric() reads self.scale**2, so it does not depend on the name.
• Add df to the params of the T distributions #440 adds df to params of the T classes. This PR leaves params alone so that the two do not overlap. The conflicts are the import lines of tests/test_distns.py and the release notes. Keep both sides.
• Stop overflow warnings from rejected line-search steps #434 conflicts in the release notes only.
• Add mu and sigma to LogNormal and raise AttributeError for unknown names #433 merges cleanly.

Tests

• test_t_df_stays_finite fits the reproducer above with 100 trees, at Y units 1 and 1e4. It fails on master with Input y contains NaN.
• test_t_df_is_capped: a log df of 800 gives df = 1000, a finite score and a finite natural gradient.
• test_t_no_push_above_cap: at the cap, a row that favors a larger df gets a log-df step of 0, and an outlier still gets a step down.
• test_t_metric checks the metric against quadrature of E[s sᵀ] at df 1.5, 3, 30 and 500.
• T moves from the gradient-only list to TEST_METRIC in tests/test_score.py.
• The xfail on test_dists_runs_on_examples_logscore is gone. T passes for both learners.
• All 8 new test cases fail on master. With --slow, the logscore distribution and metric tests pass (39 passed). The full suite without --slow passes (143 passed, 59 skipped).
• black, isort, flake8 and pylint are clean on the changed files, with the options from .pre-commit-config.yaml.

🤖 Generated with Claude Code

[ryan-wolbeck]

@alejandroschuler I'll try to get a release PR created in the next couple days. There are few breaking changes merged that will need to be documented as part of the release (to 0.5.0) then move to merge the rest you have created.

Good with that?
Ryan

[alejandroschuler]

yes, thanks @ryan-wolbeck! I admin-merged a couple of things that seemed unobjectionable, hope that's ok!

[ryan-wolbeck]

yes, thanks @ryan-wolbeck! I admin-merged a couple of things that seemed unobjectionable, hope that's ok!

No problem at all! Just want to make the releases small so if we need to roll things back it's easier to parse what happened and since there are things like sign flips now merged users will need to update code that may be running in prod that they install from pypi.