ENH: add iv_ratioinv

[dschmitz89]

Reference issue

Towards scipy/scipy#20253

What does this implement/fix?

This PR implements the new function iv_ratioinv which solves $r=I_v(x)/I_{v-1}(x)$ for x given $r$ and $v$. This equation occurs in MLE of distributions such as von Mises-Fisher. Currently, an adhoc solution exists in SciPy.

Additional information

The function is a simple derivative free root finding procedure. Tight brackets for the root were taken from
this paper. In case those bounds fail, we fall back to the bracketing routine for monotonic functions. Once, bounds are found we run Chandrupatla's algorithm to find the exact root.

The SciPy issue contains a gradient based approach but I found this one easier to implement and assume that it is more robust as the objective can be near flat in the tails.

AI Generation Disclosure

LLMs helped inverting the bounds equations.

[dschmitz89]

Cc @fancidev in case you want to take a look

[dschmitz89]

This function should be renamed as it works with any monotonic function, not just CDFs.

[dschmitz89]

@fbourgey Would you have time to review?

[fbourgey]

When $v=\frac12$, it appears we have $I_{1/2}(x)/I_{-1/2}(x)=\tanh(x) = r$, hence $x = \text{atanh}(r)$.

Could you add some tests for $v=0.5$, ideally for the extreme tails?

I wonder if std::nextafter(1.0, 0.0) and std::nextafter(0.0, 1.0) pass?

[dschmitz89]

When v = 1 2 , it appears we have I 1 / 2 ( x ) / I − 1 / 2 ( x ) = tanh ⁡ ( x ) = r , hence x = atanh ( r ) .

Could you add some tests for v = 0.5 , ideally for the extreme tails?

I wonder if std::nextafter(1.0, 0.0) and std::nextafter(0.0, 1.0) pass?

They do not unfortunately pass a roundtrip test. I think it is not possible in double precision. Could you take one last look?

Edit: the iv_ratio function itself does not use the analytical method yet, will try that out.

[fbourgey]

Suggested change

	if (v == 0.5) {
	// Closed-form solution for v = 0.5: iv_ratio(0.5, x) = tanh(x)
	// Since tanh(x) = 2*expit(2*x) - 1 = r, we have expit(2*x) = (1+r)/2
	// For r > 0.5, use logit((1+r)/2) for better stability near r=1
	// For r <= 0.5, atanh is sufficiently stable
	if (r > 0.5) {
	return 0.5 * logit((1.0 + r) * 0.5);
	} else {
	return std::atanh(r);
	}
	if (v == 0.5) {
	return std::atanh(r);
	}

would not be enough?

On my machine os-arm64, iv_ratioinv(0.5, nextafter(1.0, 0.0)) gives inf while it should be 18.714973875118524. I think this is because (1+r)*0.5 rounds to 1 in that case.