ENH: add warm start for root finder in find_degrees_of_freedom

[dschmitz89]

As promised in #1385 I added a similar warm start for the other degrees of freedom finder. It is based on Hill's approximation for the t distribution quantile $t_{p, v}$ with probability $p$ and degrees of freedom $v$:

$$ t_{p, v} \approx z_p ( 1 + (1 + z_p^2)/(4v)) $$

where $z_p$ denotes the normal distribution quantile. This was again done with help from an LLM and iterated on manually for more robustness (step 5 in details). I did a few spot checks where I printed out the number of iterations of the root finder and found that it reduced by ca. 30%. I guess the root finder anyway spends most of its time fine tuning the solution to machine precision, so even a good guess does not give a large speedup.

Detailed derivation

Derivation of the Hill-Corrected Initial Guess for Student's $t$ Degrees of Freedom

To improve the numerical efficiency of finding degrees of freedom $\nu$ for a specified statistical power, we derive an analytical "warm start" using Hill's approximation.

1. The Power Equation

The required sample size $n = \nu + 1$ is defined by the relationship between the effect size $\delta$, standard deviation $\sigma$, and the quantiles of the $t$-distribution at significance level $\alpha$ and power $1-\beta$:

$$\delta = (t_{1-\alpha, \nu} + t_{1-\beta, \nu}) \cdot \frac{\sigma}{\sqrt{\nu+1}}$$

Squaring both sides and rearranging to solve for $\nu$:

$$\nu + 1 = \left( \frac{\sigma}{\delta} \right)^2 (t_{1-\alpha, \nu} + t_{1-\beta, \nu})^2$$

2. Hill's Approximation

We substitute the Student's $t$ quantiles with Hill’s approximation, which corrects the normal quantile $z_p$ for the heavy tails of the $t$-distribution:

$$t_{p, \nu} \approx z_p \left(1 + \frac{z_p^2 + 1}{4\nu}\right)$$

Let $z_a = \Phi^{-1}(1-\alpha)$ and $z_b = \Phi^{-1}(1-\beta)$. Substituting the approximation into the power equation:

$$\nu + 1 \approx \left( \frac{\sigma}{\delta} \right)^2 \left[ z_a \left(1 + \frac{z_a^2 + 1}{4\nu}\right) + z_b \left(1 + \frac{z_b^2 + 1}{4\nu}\right) \right]^2$$

3. First-Order Expansion

Let $S_1 = z_a + z_b$ and $S_2 = \frac{z_a(z_a^2+1) + z_b(z_b^2+1)}{4}$. The equation becomes:

$$\nu + 1 \approx \left( \frac{\sigma}{\delta} \right)^2 \left( S_1 + \frac{S_2}{\nu} \right)^2$$

Expanding the square and truncating terms of order $O(1/\nu^2)$ and higher:

$$\nu + 1 \approx \left( \frac{\sigma}{\delta} \right)^2 \left( S_1^2 + \frac{2 S_1 S_2}{\nu} \right)$$

Defining the zero-order normal approximation as $\nu_{norm} + 1 = \left( \frac{\sigma S_1}{\delta} \right)^2$:

$$\nu + 1 \approx (\nu_{norm} + 1) + (\nu_{norm} + 1) \frac{2 S_2}{S_1 \nu}$$

4. Solving for the Correction Term

Assuming $\frac{\nu+1}{\nu} \approx 1$ for the correction term, we isolate the shift from the normal approximation:

$$\Delta \nu = \nu - \nu_{norm} \approx \frac{2 S_2}{S_1}$$

Substituting the definitions of $S_1$ and $S_2$:

$$\Delta \nu \approx \frac{z_a(z_a^2+1) + z_b(z_b^2+1)}{2(z_a + z_b)} = \frac{z_a^3 + z_b^3 + z_a + z_b}{2(z_a + z_b)}$$

5. Algebraic Trick to Avoid Division

Using the sum of cubes identity $a^3 + b^3 = (a+b)(a^2 - ab + b^2)$, we can factor the numerator:

$$z_a^3 + z_b^3 + z_a + z_b = (z_a + z_b)(z_a^2 - z_a z_b + z_b^2) + (z_a + z_b)$$
$$= (z_a + z_b)(z_a^2 - z_a z_b + z_b^2 + 1)$$

The $(z_a + z_b)$ terms cancel out, leaving a robust, division-free correction:

$$\Delta \nu \approx \frac{z_a^2 - z_a z_b + z_b^2 + 1}{2}$$

Final Guess Formula

The implemented "warm start" is:

$$\mathbf{\nu_{guess} = \nu_{norm} + \frac{z_a^2 - z_a z_b + z_b^2 + 1}{2}}$$

[codecov]

Codecov Report

❌ Patch coverage is 96.46018% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 95.61%. Comparing base (febada3) to head (22da379).
⚠️ Report is 16 commits behind head on develop.

Files with missing lines	Patch %	Lines
include/boost/math/distributions/students_t.hpp	71.42%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff            @@
##           develop    #1403    +/-   ##
=========================================
  Coverage    95.61%   95.61%            
=========================================
  Files          825      825            
  Lines        68487    68596   +109     
=========================================
+ Hits         65482    65589   +107     
- Misses        3005     3007     +2     

Files with missing lines	Coverage Δ
include/boost/math/distributions/non_central_f.hpp	90.73% <100.00%> (+2.76%)	⬆️
test/test_nc_f.cpp	100.00% <100.00%> (ø)
include/boost/math/distributions/students_t.hpp	91.70% <71.42%> (-0.96%)	⬇️

... and 2 files with indirect coverage changes

---

Continue to review full report in Codecov by Harness.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update febada3...22da379. Read the comment docs.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[JacobHass8]

It seems like all the testing errors aren't related to this PR. The issue is in students_t_single_sample.cpp where there is an overflow error calling boost::math::students_t_quantile<double>(double,double,double).

If I understand the derivation, this approximation only seems valid for $\nu \gg 1$. Does find_degrees_of_freedom still return the correct value even for small degrees of freedom?

[jzmaddock]

@dschmitz89 : what seems to have happened here is that the new initial guess has exposed an old bug. Basically we're starting from a different location which results in a call to the quantile with df = 0.0012658573267757545 and p = 0.125. The true result of that is -INF at double precision, hence the overflow_error. Interestingly when evaluated without long double support, numeric instability results in a finite double result "by chance" and everything proceeds as normal since any large value will do as far as the root finder is concerned. Updated to catch the exceptions and return an appropriately large value.