Practical — Statistical Methods using Python (STS-105)

# tails.py -- the four tail areas every test below needs, written out because
# no statistical package may be used.  Save this file once and import from it.
import math

# Step 1: The continued fraction for the incomplete beta
def betacf(a, b, x):
    """Continued fraction for the incomplete beta function, by Lentz's method."""
    tiny = 1e-300
    qab, qap, qam = a + b, a + 1, a - 1
    c, d = 1.0, 1 - qab * x / qap
    d = 1 / (tiny if abs(d) < tiny else d)
    h = d
    for m in range(1, 300):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1 + aa * d
        c = 1 + aa / c
        d = 1 / (tiny if abs(d) < tiny else d)
        h *= d * (tiny if abs(c) < tiny else c)
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1 + aa * d
        c = 1 + aa / c
        d = 1 / (tiny if abs(d) < tiny else d)
        de = d * c
        h *= de
        if abs(de - 1) < 3e-16:
            break
    return h

# Step 2: The regularised incomplete beta
def betainc(a, b, x):
    """The regularised incomplete beta function I_x(a, b)."""
    if x <= 0:
        return 0.0
    if x >= 1:
        return 1.0
    lf = math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
    if x < (a + 1) / (a + b + 2):
        return math.exp(lf + a * math.log(x) + b * math.log(1 - x)) * betacf(a, b, x) / a
    return 1 - math.exp(lf + b * math.log(1 - x) + a * math.log(x)) * betacf(b, a, 1 - x) / b

# Step 3: The t and F tails from it
def t_two_sided(t, df):
    """P(|T_df| > |t|)."""
    return betainc(df / 2, 0.5, df / (df + t * t))

def F_upper(f, d1, d2):
    """P(F_{d1,d2} > f)."""
    return 1.0 if f <= 0 else betainc(d2 / 2, d1 / 2, d2 / (d2 + d1 * f))

# Step 4: The chi-square tail from the incomplete gamma
def chi2_upper(q, df):
    """P(chi^2_df > q): the series for the lower gamma when q is small, and the
    continued fraction for the upper gamma when it is large, so that a tiny tail
    is computed directly instead of as 1 minus a number close to 1."""
    a, xq = df / 2, q / 2
    if xq <= 0:
        return 1.0
    front = math.exp(-xq + a * math.log(xq) - math.lgamma(a))
    if xq < a + 1:
        term, s, k = 1.0 / a, 1.0 / a, 1
        while abs(term) > 1e-15 * abs(s) and k < 10000:
            term *= xq / (a + k)
            s += term
            k += 1
        return 1 - s * front
    tiny = 1e-300                          # Lentz's method, as in betacf
    b = xq + 1 - a
    c, d = 1 / tiny, 1 / b
    h = d
    for i in range(1, 10000):
        an = -i * (i - a)
        b += 2
        d = an * d + b
        d = 1 / (tiny if abs(d) < tiny else d)
        c = b + an / c
        c = tiny if abs(c) < tiny else c
        de = d * c
        h *= de
        if abs(de - 1) < 3e-16:
            break
    return front * h

# Step 5: The normal tail
def normal_two_sided(z):
    """P(|Z| > |z|) for the standard normal."""
    return math.erfc(abs(z) / math.sqrt(2))

# Step 6: Check all four against the tables
if __name__ == "__main__":
    # three published values, as a check that the functions are right
    print(f"P(|t_10| > 2.228)      = {t_two_sided(2.228, 10):.4f}   (table: 0.0500)")
    print(f"P(F_(5,10) > 3.326)    = {F_upper(3.326, 5, 10):.4f}   (table: 0.0500)")
    print(f"P(chi^2_8 > 15.507)    = {chi2_upper(15.507, 8):.4f}   (table: 0.0500)")
    print(f"P(|Z| > 1.959964)      = {normal_two_sided(1.959964):.4f}   (table: 0.0500)")

P(|t_10| > 2.228)      = 0.0500   (table: 0.0500)
P(F_(5,10) > 3.326)    = 0.0500   (table: 0.0500)
P(chi^2_8 > 15.507)    = 0.0500   (table: 0.0500)
P(|Z| > 1.959964)      = 0.0500   (table: 0.0500)

# Practical 1 -- sum and product of two matrices, no packages
# Step 1: Check the shapes
def shape(M):
    return len(M), len(M[0])

# Step 2: Add elementwise
def add(A, B):
    if shape(A) != shape(B):
        raise ValueError("addition needs matrices of the same shape")
    return [[A[i][j] + B[i][j] for j in range(len(A[0]))] for i in range(len(A))]

# Step 3: Multiply by inner products
def multiply(A, B):
    ra, ca = shape(A)
    rb, cb = shape(B)
    if ca != rb:
        raise ValueError("columns of A must equal rows of B")
    return [[sum(A[i][k] * B[k][j] for k in range(ca)) for j in range(cb)]
            for i in range(ra)]

# Step 4: Print a matrix
def show(name, M):
    print(name)
    for row in M:
        print("  " + "  ".join(f"{v:6.2f}" for v in row))

# Step 5: Enter the two matrices
A = [[2, 1, 1], [1, 3, 2], [1, 0, 4]]
B = [[1, 0, 2], [2, 1, 0], [0, 3, 1]]

# Step 6: Form A + B, AB and BA, and compare AB with BA
show("A + B", add(A, B))
show("A x B", multiply(A, B))
show("B x A", multiply(B, A))
print("A x B equals B x A?", multiply(A, B) == multiply(B, A))

A + B
    3.00    1.00    3.00
    3.00    4.00    2.00
    1.00    3.00    5.00
A x B
    4.00    4.00    5.00
    7.00    9.00    4.00
    1.00   12.00    6.00
B x A
    4.00    1.00    9.00
    5.00    5.00    4.00
    4.00    9.00   10.00
A x B equals B x A? False

# Practical 2 -- determinant by cofactor expansion and inverse by the adjoint
from fractions import Fraction

# Step 1: Delete a row and a column
def minor(M, i, j):
    return [[M[r][c] for c in range(len(M)) if c != j]
            for r in range(len(M)) if r != i]

# Step 2: Expand the determinant along row 1
def det(M):
    n = len(M)
    if n == 1:
        return M[0][0]
    if n == 2:
        return M[0][0] * M[1][1] - M[0][1] * M[1][0]
    return sum((-1) ** j * M[0][j] * det(minor(M, 0, j)) for j in range(n))

# Step 3: Invert by the transposed cofactors
def inverse(M):
    d = det(M)
    if d == 0:
        raise ValueError("matrix is singular -- no inverse")
    n = len(M)
    cof = [[(-1) ** (i + j) * det(minor(M, i, j)) for j in range(n)] for i in range(n)]
    # the adjoint is the TRANSPOSE of the cofactor matrix
    return [[Fraction(cof[j][i], d) for j in range(n)] for i in range(n)]

# Step 4: Print a matrix of fractions
def show(name, M):
    print(name)
    for row in M:
        print("  " + "  ".join(f"{str(v):>8}" for v in row))

# Step 5: Find det(A) and the inverse
A = [[2, 1, 1], [1, 3, 2], [1, 0, 4]]
print("det(A) =", det(A))
Ainv = inverse(A)
show("A inverse", Ainv)

# Step 6: Check that A times its inverse is I
# check: A times its inverse must be the identity, exactly
prod = [[sum(Fraction(A[i][k]) * Ainv[k][j] for k in range(3)) for j in range(3)]
        for i in range(3)]
show("A x A inverse", prod)

# Step 7: Try a singular matrix
S = [[1, 2, 3], [2, 4, 6], [1, 0, 1]]
print("det(S) =", det(S), "-- row 2 is twice row 1, so S is singular")

det(A) = 19
A inverse
     12/19     -4/19     -1/19
     -2/19      7/19     -3/19
     -3/19      1/19      5/19
A x A inverse
         1         0         0
         0         1         0
         0         0         1
det(S) = 0 -- row 2 is twice row 1, so S is singular

# Practical 3 -- four sorts and two searches, each counting its own work
# Step 1: Bubble sort, counting comparisons
def bubble(a):
    a = a[:]
    comps = 0
    for i in range(len(a) - 1):
        swapped = False
        for j in range(len(a) - 1 - i):
            comps += 1
            if a[j] > a[j + 1]:
                a[j], a[j + 1] = a[j + 1], a[j]
                swapped = True
        if not swapped:
            break
    return a, comps

# Step 2: Insertion sort, counting comparisons
def insertion(a):
    a = a[:]
    comps = 0
    for i in range(1, len(a)):
        key = a[i]
        j = i - 1
        while j >= 0:
            comps += 1
            if a[j] <= key:
                break
            a[j + 1] = a[j]
            j -= 1
        a[j + 1] = key
    return a, comps

# Step 3: Merge sort
def merge_sort(a):
    if len(a) <= 1:
        return a
    mid = len(a) // 2
    left, right = merge_sort(a[:mid]), merge_sort(a[mid:])
    out, i, j = [], 0, 0
    while i < len(left) and j < len(right):
        if left[i] <= right[j]:
            out.append(left[i]); i += 1
        else:
            out.append(right[j]); j += 1
    return out + left[i:] + right[j:]

# Step 4: Quick sort
def quick_sort(a):
    if len(a) <= 1:
        return a
    pivot = a[len(a) // 2]
    lo = [v for v in a if v < pivot]
    eq = [v for v in a if v == pivot]
    hi = [v for v in a if v > pivot]
    return quick_sort(lo) + eq + quick_sort(hi)

# Step 5: Linear and binary search, counting comparisons
def linear_search(a, key):
    for i, v in enumerate(a):
        if v == key:
            return i, i + 1          # position, comparisons used
    return -1, len(a)

def binary_search(a, key):
    lo, hi, comps = 0, len(a) - 1, 0
    while lo <= hi:
        mid = (lo + hi) // 2
        comps += 1
        if a[mid] == key:
            return mid, comps
        if a[mid] < key:
            lo = mid + 1
        else:
            hi = mid - 1
    return -1, comps

# Step 6: Sort the data four ways and compare
data = [42, 17, 93, 8, 55, 23, 71, 4, 66, 30]
b, cb = bubble(data)
i_, ci = insertion(data)
print("bubble   ", b, "comparisons:", cb)
print("insertion", i_, "comparisons:", ci)
print("merge    ", merge_sort(data))
print("quick    ", quick_sort(data))
print("all four agree?", b == i_ == merge_sort(data) == quick_sort(data))

# Step 7: Search the sorted list
srt = b
print("linear search for 66:", linear_search(srt, 66))
print("binary search for 66:", binary_search(srt, 66))
print("binary search for 50:", binary_search(srt, 50))

bubble    [4, 8, 17, 23, 30, 42, 55, 66, 71, 93] comparisons: 44
insertion [4, 8, 17, 23, 30, 42, 55, 66, 71, 93] comparisons: 29
merge     [4, 8, 17, 23, 30, 42, 55, 66, 71, 93]
quick     [4, 8, 17, 23, 30, 42, 55, 66, 71, 93]
all four agree? True
linear search for 66: (7, 8)
binary search for 66: (7, 2)
binary search for 50: (-1, 4)

# Practical 4 -- median and mode of an array, written out
# Step 1: Median: sort, then take the middle
def median(a):
    s = sorted(a)
    n = len(s)
    mid = n // 2
    if n % 2:
        return s[mid]
    return (s[mid - 1] + s[mid]) / 2

# Step 2: Mode: count, then keep every value at the top count
def mode(a):
    freq = {}
    for v in a:
        freq[v] = freq.get(v, 0) + 1
    top = max(freq.values())
    modes = sorted(k for k, c in freq.items() if c == top)
    return modes, top, freq

# Step 3: Enter the test arrays
odd = [7, 3, 9, 3, 5, 8, 3, 9, 4]
even = [12, 4, 9, 4, 15, 7]
bimodal = [2, 5, 2, 7, 5, 9]

# Step 4: Medians for odd and even n
print("odd  n =", len(odd), " median =", median(odd))
print("even n =", len(even), " median =", median(even))
# Step 5: Modes: one, two and none
m, c, f = mode(odd)
print("odd  mode(s) =", m, "occurring", c, "times; frequencies", dict(sorted(f.items())))
m, c, f = mode(bimodal)
print("bimodal mode(s) =", m, "occurring", c, "times -- two modes, so the mode is not unique")
allsame = [4, 6, 9, 11]
m, c, _ = mode(allsame)
print("no repeats:", m, "each", c, "time -- the mode is undefined for this data, not 'all of them'")

odd  n = 9  median = 5
even n = 6  median = 8.0
odd  mode(s) = [3] occurring 3 times; frequencies {3: 3, 4: 1, 5: 1, 7: 1, 8: 1, 9: 2}
bimodal mode(s) = [2, 5] occurring 2 times -- two modes, so the mode is not unique
no repeats: [4, 6, 9, 11] each 1 time -- the mode is undefined for this data, not 'all of them'

# Practical 5 -- frequency table, then mean, median, mode, variance, SD
# both from the raw values and from the grouped table, so the grouping error shows
import math

# Step 1: Enter the raw data
raw = [23, 45, 12, 37, 52, 28, 41, 19, 33, 47,
       55, 26, 38, 44, 31, 17, 49, 35, 22, 40,
       29, 51, 36, 43, 25, 48, 32, 39, 27, 46,
       34, 21, 42, 30, 53, 24, 37, 45, 33, 41]

# Step 2: Count the classes
def freq_table(values, low, width, k):
    edges = [low + i * width for i in range(k + 1)]
    counts = [0] * k
    for v in values:
        idx = int((v - low) // width)
        if idx == k:            # the top edge belongs to the last class
            idx = k - 1
        counts[idx] += 1
    return edges, counts

# Step 3: Grouped mean
def grouped_mean(edges, counts):
    mids = [(edges[i] + edges[i + 1]) / 2 for i in range(len(counts))]
    n = sum(counts)
    return sum(m * f for m, f in zip(mids, counts)) / n, mids, n

# Step 4: Grouped median
def grouped_median(edges, counts):
    n = sum(counts)
    cum = 0
    for i, f in enumerate(counts):
        if cum + f >= n / 2:
            L = edges[i]
            h = edges[i + 1] - edges[i]
            return L + (n / 2 - cum) / f * h
        cum += f

# Step 5: Grouped mode
def grouped_mode(edges, counts):
    i = counts.index(max(counts))
    f1 = counts[i]
    f0 = counts[i - 1] if i > 0 else 0
    f2 = counts[i + 1] if i + 1 < len(counts) else 0
    L = edges[i]
    h = edges[i + 1] - edges[i]
    return L + (f1 - f0) / (2 * f1 - f0 - f2) * h

# Step 6: Print the table and the grouped summaries
edges, counts = freq_table(raw, 10, 10, 5)
print("class        frequency")
for i, f in enumerate(counts):
    print(f"  {edges[i]:2d} - {edges[i+1]:2d}      {f:2d}")
print("  total        ", sum(counts))

gm, mids, n = grouped_mean(edges, counts)
gmed = grouped_median(edges, counts)
gmode = grouped_mode(edges, counts)
gvar = sum(f * (m - gm) ** 2 for m, f in zip(mids, counts)) / n
print(f"grouped:  mean {gm:.4f}  median {gmed:.4f}  mode {gmode:.4f} "
      f" variance {gvar:.4f}  sd {math.sqrt(gvar):.4f}")

# Step 7: The same summaries from the raw values
rm = sum(raw) / len(raw)
s = sorted(raw)
rmed = (s[19] + s[20]) / 2
rvar = sum((v - rm) ** 2 for v in raw) / len(raw)
print(f"raw:      mean {rm:.4f}  median {rmed:.4f} "
      f"                variance {rvar:.4f}  sd {math.sqrt(rvar):.4f}")
print(f"grouping error in the mean: {gm - rm:+.4f}")

# Step 8: The tied modal class
# the two middle classes tie at 12, so the modal class is not unique --
# apply the formula to the SECOND of them and see what happens
i = 3
f1, f0, f2 = counts[i], counts[i-1], counts[i+1]
alt = edges[i] + (f1 - f0) / (2*f1 - f0 - f2) * 10
print(f"mode from the other tied class: {alt:.4f} -- the same boundary, 40")

class        frequency
  10 - 20       3
  20 - 30       9
  30 - 40      12
  40 - 50      12
  50 - 60       4
  total         40
grouped:  mean 36.2500  median 36.6667  mode 40.0000  variance 120.9375  sd 10.9972
raw:      mean 35.7500  median 36.5000                 variance 113.2375  sd 10.6413
grouping error in the mean: +0.5000
mode from the other tied class: 40.0000 -- the same boundary, 40

# Practical 6 -- the first four raw and central moments, skewness and kurtosis
import math

# Step 1: Raw moments
def raw_moments(x, k=4):
    n = len(x)
    return [sum(v ** r for v in x) / n for r in range(1, k + 1)]

# Step 2: Central moments from the raw moments
def central_from_raw(m):
    m1, m2, m3, m4 = m
    mu2 = m2 - m1 ** 2
    mu3 = m3 - 3 * m1 * m2 + 2 * m1 ** 3
    mu4 = m4 - 4 * m1 * m3 + 6 * m1 ** 2 * m2 - 3 * m1 ** 4
    return [0.0, mu2, mu3, mu4]

# Step 3: Central moments directly, as a check
def central_direct(x):
    n = len(x)
    xb = sum(x) / n
    return [0.0] + [sum((v - xb) ** r for v in x) / n for r in (2, 3, 4)]

# Step 4: Compute both ways and compare
x = [12, 15, 11, 18, 22, 14, 16, 19, 13, 25, 17, 20, 14, 16, 21]
m = raw_moments(x)
print("raw moments      m1..m4 :", "  ".join(f"{v:.4f}" for v in m))
c = central_from_raw(m)
d = central_direct(x)
print("central, from raw moments:", "  ".join(f"{v:.4f}" for v in c[1:]))
print("central, computed direct :", "  ".join(f"{v:.4f}" for v in d[1:]))
print("agree to within 1e-8      :", max(abs(a - b) for a, b in zip(c, d)) < 1e-8)

# Step 5: Skewness and kurtosis
mu2, mu3, mu4 = c[1], c[2], c[3]
b1 = mu3 ** 2 / mu2 ** 3
g1 = mu3 / mu2 ** 1.5
b2 = mu4 / mu2 ** 2
g2 = b2 - 3
print(f"beta1 {b1:.6f}   gamma1 {g1:+.6f}   beta2 {b2:.6f}   gamma2 {g2:+.6f}")
print("skew  :", "positive -- the long tail is to the right" if g1 > 0 else "negative")
print("kurt  :", "platykurtic -- flatter than normal" if g2 < 0 else "leptokurtic")
print("sd    :", f"{math.sqrt(mu2):.6f}")

raw moments      m1..m4 : 16.8667  299.1333  5563.2667  108043.1333
central, from raw moments: 14.6489  23.7286  506.5591
central, computed direct : 14.6489  23.7286  506.5591
agree to within 1e-8      : True
beta1 0.179114   gamma1 +0.423219   beta2 2.360591   gamma2 -0.639409
skew  : positive -- the long tail is to the right
kurt  : platykurtic -- flatter than normal
sd    : 3.827387

# Practical 7 -- random numbers from five distributions, built from an LCG upward
import math

# Step 1: The base generator
class LCG:
    """Lehmer's minimal standard generator: x <- 16807 x mod (2**31 - 1)."""
    def __init__(self, seed=12345):
        self.x = seed
    def uniform(self):
        self.x = (16807 * self.x) % 2147483647
        return self.x / 2147483647

rng = LCG(12345)

# Step 2: Uniform, Bernoulli and binomial
def uniform(a, b):                      # inverse transform on U(0,1)
    return a + (b - a) * rng.uniform()

def bernoulli(p):
    return 1 if rng.uniform() < p else 0

def binomial(n, p):                     # sum of n Bernoulli trials
    return sum(bernoulli(p) for _ in range(n))

# Step 3: Poisson by Knuth's method
def poisson(lam):                       # Knuth: multiply uniforms until below e^-lam
    L, k, prod = math.exp(-lam), 0, 1.0
    while True:
        prod *= rng.uniform()
        if prod <= L:
            return k
        k += 1

# Step 4: Normal by Box-Muller
def normal(mu=0.0, sigma=1.0):          # Box-Muller, one of the pair kept
    u1, u2 = rng.uniform(), rng.uniform()
    z = math.sqrt(-2 * math.log(u1)) * math.cos(2 * math.pi * u2)
    return mu + sigma * z

# Step 5: Exponential by inverse transform
def exponential(lam):                   # inverse transform: -ln(U)/lambda
    return -math.log(rng.uniform()) / lam

# Step 6: Compare a sample with the theory
def summarise(name, sample, mean_th, var_th):
    n = len(sample)
    m = sum(sample) / n
    v = sum((x - m) ** 2 for x in sample) / n
    print(f"{name:12s} n={n}  mean {m:8.4f} (theory {mean_th:7.4f})"
          f"   var {v:9.4f} (theory {var_th:8.4f})")

# Step 7: Draw 20,000 of each and compare
N = 20000
summarise("Uniform(2,8)", [uniform(2, 8) for _ in range(N)], 5.0, 36 / 12)
summarise("Bin(10,0.3)", [binomial(10, 0.3) for _ in range(N)], 3.0, 10 * 0.3 * 0.7)
summarise("Poisson(4)", [poisson(4) for _ in range(N)], 4.0, 4.0)
summarise("Normal(5,2)", [normal(5, 2) for _ in range(N)], 5.0, 4.0)
summarise("Exp(0.5)", [exponential(0.5) for _ in range(N)], 2.0, 4.0)

Uniform(2,8) n=20000  mean   5.0051 (theory  5.0000)   var    3.0352 (theory   3.0000)
Bin(10,0.3)  n=20000  mean   3.0009 (theory  3.0000)   var    2.0996 (theory   2.1000)
Poisson(4)   n=20000  mean   4.0315 (theory  4.0000)   var    4.0750 (theory   4.0000)
Normal(5,2)  n=20000  mean   5.0204 (theory  5.0000)   var    4.0090 (theory   4.0000)
Exp(0.5)     n=20000  mean   1.9800 (theory  2.0000)   var    3.9398 (theory   4.0000)

# Practical 8 -- fit Binomial, Poisson and Negative Binomial to ONE data set
# and let the goodness-of-fit test choose between them
import math
from tails import chi2_upper

# Step 1: Enter the data; the mean, the variance and their ratio
x = [0, 1, 2, 3, 4, 5]
f = [447, 132, 42, 21, 3, 2]
n = sum(f)
mean = sum(xi * fi for xi, fi in zip(x, f)) / n
var = sum(fi * (xi - mean) ** 2 for xi, fi in zip(x, f)) / n
print(f"n = {n}   mean = {mean:.6f}   variance = {var:.6f}   var/mean = {var/mean:.6f}")

# Step 2: The three probability functions
def C(n, r):
    return math.factorial(n) // (math.factorial(r) * math.factorial(n - r))

def binomial_pmf(N, p, k):
    return C(N, k) * p ** k * (1 - p) ** (N - k)

def poisson_pmf(lam, k):
    return math.exp(-lam) * lam ** k / math.factorial(k)

def nbinom_pmf(r, p, k):
    # P(X = k) = Gamma(r+k)/(k! Gamma(r)) p^r (1-p)^k
    coef = 1.0
    for i in range(k):
        coef *= (r + i) / (i + 1)
    return coef * p ** r * (1 - p) ** k

# Step 3: Estimate the parameters
N = 5
p_bin = mean / N
lam = mean
p_nb = mean / var
r_nb = mean ** 2 / (var - mean)
print(f"Binomial  N = {N}, p = {p_bin:.6f}     (p estimated by the method of moments)")
print(f"Poisson   lambda = {lam:.6f}")
print(f"NegBin    r = {r_nb:.6f}, p = {p_nb:.6f}")

# Step 4: Expected frequencies
fits = {
    "Binomial": [n * binomial_pmf(N, p_bin, k) for k in x],
    "Poisson": [n * poisson_pmf(lam, k) for k in x],
    "NegBin": [n * nbinom_pmf(r_nb, p_nb, k) for k in x],
}
# the tail beyond k = 5 belongs to the last class for the two unbounded fits
for name in ("Poisson", "NegBin"):
    fits[name][-1] = n - sum(fits[name][:-1])

print()
print("  k   observed   Binomial    Poisson     NegBin")
for i, k in enumerate(x):
    print(f"  {k}   {f[i]:8d}   {fits['Binomial'][i]:8.2f}   "
          f"{fits['Poisson'][i]:8.2f}   {fits['NegBin'][i]:8.2f}")

# Step 5: Pool the tail, then chi-square
def pooled_chi2(obs, exp, params):
    """Pool the upper tail until every expected frequency is at least 5."""
    o, e = list(obs), list(exp)
    while len(e) > 2 and e[-1] < 5:
        last_e, last_o = e.pop(), o.pop()
        e[-1] += last_e
        o[-1] += last_o
    chi2 = sum((oi - ei) ** 2 / ei for oi, ei in zip(o, e))
    df = len(e) - 1 - params
    return chi2, df, len(e)

print()
# Step 6: The test for each fit
for name, params in (("Binomial", 1), ("Poisson", 1), ("NegBin", 2)):
    chi2, df, k = pooled_chi2(f, fits[name], params)
    print(f"{name:9s} chi-square {chi2:9.4f} on {df} df "
          f"({k} classes after pooling, {params} parameter(s) estimated)")

for name, params in (("Binomial", 1), ("Poisson", 1), ("NegBin", 2)):
    chi2, df, _ = pooled_chi2(f, fits[name], params)
    print(f"{name:9s} p-value {chi2_upper(chi2, df):.4g}")

# Step 7: Draw the comparison with characters
# the "curve plot", drawn with characters since no plotting package is allowed
print()
print("observed (#) against the negative binomial fit (o), one column per count")
scale = 460 / 46
for i, k in enumerate(x):
    bar = "#" * max(1, round(f[i] / scale))
    fit = round(fits["NegBin"][i] / scale)
    print(f"  {k}  {bar:<48s}| fit {'o' * max(1, fit)}")

n = 647   mean = 0.465224   variance = 0.690831   var/mean = 1.484942
Binomial  N = 5, p = 0.093045     (p estimated by the method of moments)
Poisson   lambda = 0.465224
NegBin    r = 0.959340, p = 0.673427

  k   observed   Binomial    Poisson     NegBin
  0        447     397.04     406.31     442.77
  1        132     203.66     189.03     138.72
  2         42      41.79      43.97      44.38
  3         21       4.29       6.82      14.30
  4          3       0.22       0.79       4.62
  5          2       0.00       0.08       2.22

Binomial  chi-square   41.6739 on 1 df (3 classes after pooling, 1 parameter(s) estimated)
Poisson   chi-square   64.9467 on 2 df (4 classes after pooling, 1 parameter(s) estimated)
NegBin    chi-square    4.1299 on 2 df (5 classes after pooling, 2 parameter(s) estimated)
Binomial  p-value 1.078e-10
Poisson   p-value 7.889e-15
NegBin    p-value 0.1268

observed (#) against the negative binomial fit (o), one column per count
  0  #############################################   | fit oooooooooooooooooooooooooooooooooooooooooooo
  1  #############                                   | fit oooooooooooooo
  2  ####                                            | fit oooo
  3  ##                                              | fit o
  4  #                                               | fit o
  5  #                                               | fit o

# Practical 9 -- fit Normal, Exponential and Cauchy to the same grouped data
import math
from tails import chi2_upper

# Step 1: Enter the grouped data; the mean and the sd
edges = [10, 20, 30, 40, 50, 60, 70, 80]
obs = [9, 15, 37, 55, 36, 17, 6]
n = sum(obs)
mids = [(edges[i] + edges[i + 1]) / 2 for i in range(len(obs))]

mean = sum(m * f for m, f in zip(mids, obs)) / n
var = sum(f * (m - mean) ** 2 for m, f in zip(mids, obs)) / n
sd = math.sqrt(var)
print(f"n = {n}   mean = {mean:.6f}   sd = {sd:.6f}")

# Step 2: Quartiles by interpolation, for the Cauchy
def quantile(q):
    """The q-quantile of the grouped distribution, by linear interpolation."""
    target, cum = q * n, 0
    for i, f in enumerate(obs):
        if cum + f >= target:
            return edges[i] + (target - cum) / f * (edges[i + 1] - edges[i])
        cum += f

med = quantile(0.5)
q1, q3 = quantile(0.25), quantile(0.75)
gamma = (q3 - q1) / 2
print(f"median = {med:.6f}   Q1 = {q1:.6f}   Q3 = {q3:.6f}   half-IQR = {gamma:.6f}")

# Step 3: The three distribution functions
def normal_cdf(x, mu, s):
    return 0.5 * (1 + math.erf((x - mu) / (s * math.sqrt(2))))

def expon_cdf(x, lam):
    return 0.0 if x <= 0 else 1 - math.exp(-lam * x)

def cauchy_cdf(x, x0, g):
    return 0.5 + math.atan((x - x0) / g) / math.pi

# Step 4: The estimates
lam = 1 / mean
print(f"Normal    mu = {mean:.6f}, sigma = {sd:.6f}   (moment estimates)")
print(f"Exponent  lambda = 1/mean = {lam:.6f}")
print(f"Cauchy    x0 = median = {med:.6f}, gamma = half-IQR = {gamma:.6f}")
print("  -- the Cauchy has no mean or variance, so it cannot be fitted by moments;")
print("     the median and the half-interquartile range are the natural substitutes.")

# Step 5: Expected frequencies from the distribution function
def expected(cdf, *args):
    """Expected frequencies, with the two outer classes opened to +-infinity
    so that the probabilities sum to exactly 1."""
    e = []
    for i in range(len(obs)):
        lo = 0.0 if i == 0 else cdf(edges[i], *args)
        hi = 1.0 if i == len(obs) - 1 else cdf(edges[i + 1], *args)
        e.append(n * (hi - lo))
    return e

fits = {
    "Normal": expected(normal_cdf, mean, sd),
    "Exponential": expected(expon_cdf, lam),
    "Cauchy": expected(cauchy_cdf, med, gamma),
}

print()
print("  class      observed    Normal   Exponential    Cauchy")
for i in range(len(obs)):
    print(f"  {edges[i]:2d} - {edges[i+1]:2d}   {obs[i]:8d}  {fits['Normal'][i]:9.2f}  "
          f"{fits['Exponential'][i]:11.2f}  {fits['Cauchy'][i]:8.2f}")

print()
# Step 6: Chi-square for each fit
for name, params in (("Normal", 2), ("Exponential", 1), ("Cauchy", 2)):
    e = fits[name]
    chi2 = sum((o - ei) ** 2 / ei for o, ei in zip(obs, e))
    df = len(obs) - 1 - params
    print(f"{name:12s} chi-square {chi2:10.4f} on {df} df   p = {chi2_upper(chi2, df):.4g}")

n = 175   mean = 44.657143   sd = 13.852164
median = 44.818182   Q1 = 35.337838   Q3 = 54.236111   half-IQR = 9.449137
Normal    mu = 44.657143, sigma = 13.852164   (moment estimates)
Exponent  lambda = 1/mean = 0.022393
Cauchy    x0 = median = 44.818182, gamma = half-IQR = 9.449137
  -- the Cauchy has no mean or variance, so it cannot be fitted by moments;
     the median and the half-interquartile range are the natural substitutes.

  class      observed    Normal   Exponential    Cauchy
  10 - 20          9       6.57        63.18     20.26
  20 - 30         15      18.81        22.44     11.36
  30 - 40         37      39.09        17.93     29.61
  40 - 50         55      49.31        14.34     54.21
  50 - 60         36      37.77        11.46     28.55
  60 - 70         17      17.56         9.16     11.02
  70 - 80          6       5.89        36.50     20.00

Normal       chi-square     2.5410 on 4 df   p = 0.6373
Exponential  chi-square   269.2814 on 5 df   p = 3.992e-56
Cauchy       chi-square    24.2799 on 4 df   p = 7.019e-05

# Practical 10 -- correlation and BOTH regression lines
import math

# Step 1: Enter the data; the sums of squares and products
x = [65, 63, 67, 64, 68, 62, 70, 66, 68, 67]
y = [68, 66, 68, 65, 69, 66, 68, 65, 71, 67]
n = len(x)

xb, yb = sum(x) / n, sum(y) / n
Sxx = sum((v - xb) ** 2 for v in x)
Syy = sum((v - yb) ** 2 for v in y)
Sxy = sum((u - xb) * (v - yb) for u, v in zip(x, y))
print(f"n = {n}   xbar = {xb:.4f}   ybar = {yb:.4f}")
print(f"Sxx = {Sxx:.4f}   Syy = {Syy:.4f}   Sxy = {Sxy:.4f}")

# Step 2: The correlation
r = Sxy / math.sqrt(Sxx * Syy)
print(f"r = {r:.6f}   r^2 = {r*r:.6f}")

# Step 3: Both regression lines, and the checks
b_yx = Sxy / Sxx
a_yx = yb - b_yx * xb
b_xy = Sxy / Syy
a_xy = xb - b_xy * yb
print(f"y on x :  y = {a_yx:.6f} + {b_yx:.6f} x")
print(f"x on y :  x = {a_xy:.6f} + {b_xy:.6f} y")
print(f"check   :  b_yx * b_xy = {b_yx*b_xy:.6f}  and  r^2 = {r*r:.6f}")
print(f"the two lines meet at the point (xbar, ybar) = ({xb:.4f}, {yb:.4f})")
print(f"  y on x at xbar : {a_yx + b_yx*xb:.4f}")
print(f"  x on y at ybar : {a_xy + b_xy*yb:.4f}")

# Step 4: The angle between the lines
# the angle between the two lines, which is zero only when |r| = 1
m1, m2 = b_yx, 1 / b_xy
theta = math.degrees(math.atan(abs((m2 - m1) / (1 + m1 * m2))))
print(f"angle between the lines = {theta:.4f} degrees "
      f"(it would be 0 if r were +-1, and 90 if r were 0)")

# Step 5: The residual sd and a prediction
# residual standard deviation and the prediction at a new x
see = math.sqrt((Syy - b_yx * Sxy) / (n - 2))
print(f"residual sd (n-2 divisor) = {see:.6f}")
print(f"predicted y at x = 69 : {a_yx + b_yx*69:.6f}")

n = 10   xbar = 66.0000   ybar = 67.3000
Sxx = 56.0000   Syy = 32.1000   Sxy = 27.0000
r = 0.636821   r^2 = 0.405541
y on x :  y = 35.478571 + 0.482143 x
x on y :  x = 9.392523 + 0.841121 y
check   :  b_yx * b_xy = 0.405541  and  r^2 = 0.405541
the two lines meet at the point (xbar, ybar) = (66.0000, 67.3000)
  y on x at xbar : 67.3000
  x on y at ybar : 66.0000
angle between the lines = 24.1914 degrees (it would be 0 if r were +-1, and 90 if r were 0)
residual sd (n-2 divisor) = 1.544431
predicted y at x = 69 : 68.746429

# Practical 11 -- tests for means, variances and correlations,
# with every tail area taken from tails.py rather than from a table
import math
from tails import t_two_sided, F_upper, chi2_upper, normal_two_sided

# Step 1: Size, mean and variance (n - 1 divisor)
def stats(a):
    n = len(a)
    m = sum(a) / n
    s2 = sum((v - m) ** 2 for v in a) / (n - 1)
    return n, m, s2

# Step 2: One-sample t
sample = [48, 52, 55, 49, 53, 51, 47, 54, 50, 52, 56, 50]
n1, m1, s21 = stats(sample)
mu0 = 50
t1 = (m1 - mu0) / math.sqrt(s21 / n1)
print(f"one-sample t : n {n1}  mean {m1:.6f}  s^2 {s21:.6f}")
print(f"   H0: mu = {mu0}    t = {t1:.6f} on {n1-1} df   p = {t_two_sided(t1, n1-1):.6f}")

# Step 3: Two independent samples, pooled variance
A = [23, 27, 25, 29, 24, 26, 28]
B = [31, 28, 33, 30, 35, 29]
na, ma, s2a = stats(A)
nb, mb, s2b = stats(B)
sp2 = ((na - 1) * s2a + (nb - 1) * s2b) / (na + nb - 2)
t2 = (ma - mb) / math.sqrt(sp2 * (1 / na + 1 / nb))
print(f"two-sample t : means {ma:.6f} and {mb:.6f}   pooled s^2 {sp2:.6f}")
print(f"   t = {t2:.6f} on {na+nb-2} df   p = {t_two_sided(t2, na+nb-2):.6f}")

# Step 4: Paired t
before = [72, 68, 75, 70, 66, 74, 69, 71]
after = [70, 65, 71, 68, 65, 70, 67, 69]
d = [b - a for b, a in zip(before, after)]
nd, md, s2d = stats(d)
t3 = md / math.sqrt(s2d / nd)
print(f"paired t     : differences {d}")
print(f"   mean d {md:.6f}   t = {t3:.6f} on {nd-1} df   p = {t_two_sided(t3, nd-1):.8f}")

# Step 5: One variance against a stated value
sigma0_sq = 6.0
chi = (n1 - 1) * s21 / sigma0_sq
p_lo, p_hi = 1 - chi2_upper(chi, n1 - 1), chi2_upper(chi, n1 - 1)
print(f"one variance : H0 sigma^2 = {sigma0_sq}   chi-square = {chi:.6f} on {n1-1} df")
print(f"   two-sided p = {2*min(p_lo, p_hi):.6f}")

# Step 6: Two variances
F = max(s2a, s2b) / min(s2a, s2b)
df1, df2 = (na - 1, nb - 1) if s2a > s2b else (nb - 1, na - 1)
print(f"two variances: s^2 {s2a:.6f} and {s2b:.6f}   F = {F:.6f} on ({df1}, {df2}) df")
print(f"   two-sided p = {2*F_upper(F, df1, df2):.6f}  -- pooling was justified")

# Step 7: A correlation against zero
x = [65, 63, 67, 64, 68, 62, 70, 66, 68, 67]
y = [68, 66, 68, 65, 69, 66, 68, 65, 71, 67]
nx = len(x)
xb, yb = sum(x) / nx, sum(y) / nx
Sxx = sum((v - xb) ** 2 for v in x)
Syy = sum((v - yb) ** 2 for v in y)
Sxy = sum((u - xb) * (v - yb) for u, v in zip(x, y))
r = Sxy / math.sqrt(Sxx * Syy)
tr = r * math.sqrt((nx - 2) / (1 - r * r))
print(f"correlation  : r = {r:.6f}   t = {tr:.6f} on {nx-2} df   "
      f"p = {t_two_sided(tr, nx-2):.6f}")

# Step 8: A correlation against a stated non-zero value, by Fisher's z
rho0 = 0.8
z = 0.5 * math.log((1 + r) / (1 - r))
z0 = 0.5 * math.log((1 + rho0) / (1 - rho0))
zstat = (z - z0) * math.sqrt(nx - 3)
print(f"   Fisher z = {z:.6f}, z0 = {z0:.6f}   test statistic {zstat:.6f}   "
      f"p = {normal_two_sided(zstat):.6f}")

one-sample t : n 12  mean 51.416667  s^2 7.719697
   H0: mu = 50    t = 1.766274 on 11 df   p = 0.105050
two-sample t : means 26.000000 and 31.000000   pooled s^2 5.636364
   t = -3.785502 on 11 df   p = 0.003018
paired t     : differences [2, 3, 4, 2, 1, 4, 2, 2]
   mean d 2.500000   t = 6.614378 on 7 df   p = 0.00030027
one variance : H0 sigma^2 = 6.0   chi-square = 14.152778 on 11 df
   two-sided p = 0.449309
two variances: s^2 4.666667 and 6.800000   F = 1.457143 on (5, 6) df
   two-sided p = 0.654305  -- pooling was justified
correlation  : r = 0.636821   t = 2.336152 on 8 df   p = 0.047701
   Fisher z = 0.752807, z0 = 1.098612   test statistic -0.914914   p = 0.360237

# Practical 12 -- one-way and two-way analysis of variance, built from the
# correction factor upward, with the F tail area taken from tails.py
import math
from tails import F_upper

# Step 1: Print an analysis of variance table
def anova_table(rows, source_names):
    print(f"  {'Source':<14}{'SS':>12}{'df':>5}{'MS':>12}{'F':>10}{'p':>12}")
    err_ss, err_df = rows[-1][1], rows[-1][2]
    err_ms = err_ss / err_df
    for name, ss, df in rows[:-1]:
        ms = ss / df
        F = ms / err_ms
        print(f"  {name:<14}{ss:12.4f}{df:5d}{ms:12.4f}{F:10.4f}{F_upper(F, df, err_df):>12.4g}")
    print(f"  {'Error':<14}{err_ss:12.4f}{err_df:5d}{err_ms:12.4f}")
    tot_ss = sum(r[1] for r in rows)
    tot_df = sum(r[2] for r in rows)
    print(f"  {'Total':<14}{tot_ss:12.4f}{tot_df:5d}")

# Step 2: One-way: correction factor, then the sums of squares
groups = {"A": [20, 22, 19, 24, 25],
          "B": [27, 25, 30, 28, 26],
          "C": [23, 21, 24, 22, 25]}
allv = [v for g in groups.values() for v in g]
N = len(allv)
G = sum(allv)
CF = G * G / N
SST = sum(v * v for v in allv) - CF
SSTr = sum(sum(g) ** 2 / len(g) for g in groups.values()) - CF
SSE = SST - SSTr
print(f"one-way:  N = {N}  G = {G}  CF = {CF:.4f}")
anova_table([("Treatments", SSTr, len(groups) - 1),
             (None, SSE, N - len(groups))], None)

# Step 3: Two-way, one observation per cell
# rows = varieties, columns = blocks
tab = [[18, 22, 20, 16],
       [23, 25, 26, 21],
       [15, 19, 18, 14]]
r, c = len(tab), len(tab[0])
vals = [v for row in tab for v in row]
N2 = r * c
G2 = sum(vals)
CF2 = G2 * G2 / N2
SST2 = sum(v * v for v in vals) - CF2
SSR = sum(sum(row) ** 2 for row in tab) / c - CF2
SSC = sum(sum(tab[i][j] for i in range(r)) ** 2 for j in range(c)) / r - CF2
SSE2 = SST2 - SSR - SSC
print()
print(f"two-way:  N = {N2}  G = {G2}  CF = {CF2:.4f}")
anova_table([("Varieties", SSR, r - 1),
             ("Blocks", SSC, c - 1),
             (None, SSE2, (r - 1) * (c - 1))], None)
# Step 4: Check that the parts add back to the total
print(f"  check: SSR + SSC + SSE = {SSR + SSC + SSE2:.4f}  and  SST = {SST2:.4f}")

one-way:  N = 15  G = 361  CF = 8688.0667
  Source                  SS   df          MS         F           p
  Treatments         76.1333    2     38.0667    8.9921    0.004109
  Error              50.8000   12      4.2333
  Total             126.9333   14

two-way:  N = 12  G = 237  CF = 4680.7500
  Source                  SS   df          MS         F           p
  Varieties         108.5000    2     54.2500  114.8824   1.648e-05
  Blocks             48.9167    3     16.3056   34.5294   0.0003516
  Error               2.8333    6      0.4722
  Total             160.2500   11
  check: SSR + SSC + SSE = 160.2500  and  SST = 160.2500