ENH: add parameter finder for degrees or freedom for students_t distribution

[dschmitz89]

Towards #1305

This adds a parameter finder for the student t distribution with respect to the degrees or freedom.

Disclaimer: this PR is heavily LLM supported as C++ is not (yet?) my forte.

I verified the math with a simple Python program before. I also ran a simple sanity script to see that the function actually executes. I am not so confident about the tests though as I cannot really decipher the output of b2.

Approach: We use bracket_and_solve_root with an initial guess from a second order Edgeworth expansion of the CDF. The initial guess is very good for $df > 1$ where the t distribution behaves similar to the normal distribution but useless for low degrees of freedom (see the plot below). In this case or if it gives a relative error of the CDF worse than 10%, we fall back to an initial guess of 0.01.

Detailed approach for finding the initial guess

Given a quantile $x$ and a CDF value $p = P(T \leq x)$, we want to recover the degrees of freedom $\nu$.

Step 1 — Edgeworth warm start.

We use the 2nd-order Edgeworth expansion of the $t$ CDF in powers of $1/\nu$:

$$F(x;\nu) \approx \Phi(x) - \frac{\varphi(x)(x + x^3)}{4\nu} + \frac{\varphi(x)(3x + 5x^3 + 7x^5 - 3x^7)}{96\nu^2}$$
where $\Phi(x)$ is the standard normal CDF and $\phi(x)$ the standard normal PDF.
Setting $u = 1/\nu$ and $F(x;\nu) = p$ yields the quadratic

$$B u^2 - A u + (\Phi(x) - p) = 0$$

where

$$A = \frac{\varphi(x)(x + x^3)}{4}, \qquad B = \frac{\varphi(x)(3x + 5x^3 + 7x^5 - 3x^7)}{96}$$

The physically meaningful root (smallest positive $u$, i.e. largest $\nu$) gives a closed-form starting estimate $\hat\nu = 1/u$.

Step 2 — Validation.

We plug $\hat\nu$ into the exact $t$ CDF. If the relative residual $|F(x;\hat\nu) - p|/|p|$ exceeds 10%, we fall back to a safe low starting value $\nu_0 = 10^{-2}$.

[dschmitz89]

CC @JacobHass8 would you have time to help out with a cursory initial review?

[JacobHass8]

I see these test values are used in quite a few places. It would be nice to combine all these into a single function which checks every method. This ensures we only have to define the values for the spot check in one place. It also cuts down on the amount of code. Maybe this would be for a separate PR though?

[JacobHass8]

With all this in the function, it was initially hard for me to parse what was going on. Could this be put in a separate function called something like boost::math::detail::edgeworth_approximation? Then you could replace all this with hint = boost::math::detail::edgeworth_approximation(x_abs, p).