fix: replace divergent fixed-point iteration with bisection in fixed_point_approx

[shayan0373n]

Problem

fixed_point_approx uses direct fixed-point iteration to solve for kappa in the Gaussian approximation of Poisson observations (Discrete model). When r * sigma² is large — which happens naturally as sigma shrinks during inference with score differences >= 3 — the fixed-point map becomes non-contractive: its derivative exceeds 1 in magnitude, causing the iteration to oscillate between two values that never converge.

This leads to either:

• ValueError: math domain error (log of a negative number)
• Wildly incorrect posterior estimates (the iteration ends at a non-fixed-point after hitting max_iter)

Reproduction

from trueskillthroughtime import fixed_point_approx

# Diverges: r*sigma² = 6*0.0625 = 0.375
fixed_point_approx(r=6, mu=0.5, sigma=0.25)  # ValueError or garbage

# Even moderate cases diverge:
fixed_point_approx(r=3, mu=1.0, sigma=0.5)   # Oscillates after ~8 iterations

Root Cause

The fixed-point equation is:

$$\kappa = \log\left(\frac{\mu + r\sigma^2 - 1 - \kappa + \sqrt{(\kappa - \mu - r\sigma^2 - 1)^2 + 2\sigma^2}}{2\sigma^2}\right)$$

The map f(kappa) has a unique fixed point, but |f'(kappa*)| > 1 for many practical parameter combinations, making direct iteration divergent. The fixed point still exists — the iteration just can't reach it.

Fix

Replace direct iteration with bisection on g(kappa) = f(kappa) - kappa:

• g is monotonically decreasing (proven by the structure of the equation)
• The root always exists in [-20, 20] for practical inputs
• Bisection is guaranteed to converge
• 50 steps reach ~1e-10 tolerance (microseconds of compute)
• Zero new dependencies — only uses math.sqrt and math.log

Verified against scipy.optimize.brentq

r	mu	sigma	Bisection mu_new	Bisection sigma_new	Error vs brentq
3	1.0	0.50	1.011399	0.379202	2e-11
6	0.5	0.25	0.740305	0.234693	3e-12
6	-0.5	0.25	-0.178838	0.243530	9e-13
6	0.0	0.20	0.190666	0.195242	1e-12
6	2.0	0.15	1.971799	0.139080	3e-12