# twaarden.py - de kritieke t-waarden zelf berekenen, en controleren met scipy
import math
from scipy import stats


def t_dichtheid(t, nu):
    """Kansdichtheid van Student-t met nu vrijheidsgraden."""
    c = math.exp(math.lgamma((nu + 1) / 2) - math.lgamma(nu / 2)) / math.sqrt(nu * math.pi)
    return c * (1 + t * t / nu) ** (-(nu + 1) / 2)


def kans_tot(t, nu, n=2000):
    """P(T <= t) voor t >= 0: de helft plus de oppervlakte van 0 tot t (regel van Simpson, n even)."""
    h = t / n
    som = t_dichtheid(0, nu) + t_dichtheid(t, nu)
    for i in range(1, n):
        som += (4 if i % 2 else 2) * t_dichtheid(i * h, nu)
    return 0.5 + som * h / 3


def t_kritiek(p, nu):
    """Zoek t met P(T <= t) = p door bisectie (p > 0,5)."""
    a, b = 0.0, 100.0
    for _ in range(60):
        m = (a + b) / 2
        if kans_tot(m, nu) < p:
            a = m
        else:
            b = m
    return (a + b) / 2


print(" nu   t(95 %)  scipy     verschil")
for nu in (1, 2, 5, 10, 30):
    eigen = t_kritiek(0.975, nu)
    print("%3d   %7.4f  %7.4f   %.1e" % (nu, eigen, stats.t.ppf(0.975, nu), abs(eigen - stats.t.ppf(0.975, nu))))