graduate_code/solutions/evidence.py

"""G04: explicit statistical procedures; validity still needs a sampling model."""
import math
import random


def wilson(successes, trials, z=1.96):
    if type(trials) is not int or type(successes) is not int or not 0 <= successes <= trials or trials < 1:
        raise ValueError("integer counts with 0 <= successes <= trials required")
    if not math.isfinite(z) or z <= 0:
        raise ValueError("finite positive critical value required")
    p = successes / trials
    denominator = 1 + z*z/trials
    center = (p + z*z/(2*trials)) / denominator
    half = z*math.sqrt(p*(1-p)/trials + z*z/(4*trials*trials)) / denominator
    return max(0., center-half), min(1., center+half)


def entropy(probabilities):
    p = list(probabilities)
    if not p or any(not math.isfinite(x) or x < 0 or x > 1 for x in p):
        raise ValueError("nonempty finite probability vector required")
    if not math.isclose(math.fsum(p), 1., rel_tol=0., abs_tol=1e-12):
        raise ValueError("probabilities must sum to one")
    return -math.fsum(x*math.log2(x) for x in p if x > 0)


def quantile(values, probability):
    values = sorted(values)
    if not values or not math.isfinite(probability) or not 0 <= probability <= 1:
        raise ValueError("nonempty data and probability in [0,1] required")
    if any(not math.isfinite(x) for x in values):
        raise ValueError("finite observations required")
    position = (len(values)-1)*probability
    lo = math.floor(position)
    hi = min(lo+1, len(values)-1)
    fraction = position-lo
    return (1-fraction)*values[lo]+fraction*values[hi]


def paired_bootstrap(old, new, repeats=2000, seed=21, confidence=.95):
    if len(old) != len(new) or not 2 <= len(old) <= 10000:
        raise ValueError("2..10000 matched pairs required")
    if type(repeats) is not int or not 100 <= repeats <= 10000:
        raise ValueError("100..10000 resamples required")
    if not math.isfinite(confidence) or not 0 < confidence < 1:
        raise ValueError("confidence must lie in (0,1)")
    if any(not math.isfinite(v) for v in [*old, *new]):
        raise ValueError("finite measurements required")
    differences = [b-a for a,b in zip(old,new)]
    if any(not math.isfinite(d) for d in differences):
        raise ValueError("difference overflow")
    rng = random.Random(seed)
    n = len(differences)
    means = [math.fsum(rng.choice(differences) for _ in range(n))/n
             for _ in range(repeats)]
    tail = (1-confidence)/2
    return {"mean_difference": math.fsum(differences)/n,
            "interval": (quantile(means,tail),quantile(means,1-tail)),
            "seed":seed,"pairs":n,"resamples":repeats}


if __name__ == "__main__":
    print("Wilson:",wilson(18,20))
    print("Entropy:",entropy([.5,.25,.25]))
    print(paired_bootstrap([10,20,30],[8,18,28]))