Statistics: centre, spread and distributions

"""Centre and spread, and what one extreme value does to them.

For the lesson's data (minutes taken to get to school by 15 students) this
experiment computes the mean, median, mode, range, quartiles and
interquartile range exactly, using the lesson's convention for quartiles
(medians of the lower and upper halves, leaving out the median when the
number of values is odd). It then makes the largest value more and more
extreme and records which summaries move. Two side notes: other quartile
conventions (Python's own) give slightly different answers on some data,
and two bus routes with the same mean can have very different spreads.
"""

import statistics
from fractions import Fraction

CLASS_A = [5, 8, 10, 10, 12, 12, 12, 15, 15, 18, 20, 20, 25, 30, 45]
ROUTE_1 = [9, 10, 10, 11, 9, 10, 11, 10, 9, 11]
ROUTE_2 = [2, 18, 5, 15, 10, 1, 19, 8, 12, 10]

def mean(values):
    """Exact mean as a Fraction."""
    if not values:
        raise ValueError("no data")
    return Fraction(sum(values), len(values))

def median(values):
    """Middle value of the sorted data, or the mean of the two middle ones."""
    s = sorted(values)
    n = len(s)
    if n == 0:
        raise ValueError("no data")
    if n % 2 == 1:
        return Fraction(s[n // 2])
    return Fraction(s[n // 2 - 1] + s[n // 2], 2)

def quartiles(values):
    """(Q1, Q3): medians of the lower and upper halves (median excluded)."""
    s = sorted(values)
    n = len(s)
    if n < 2:
        raise ValueError("need at least two values")
    return median(s[: n // 2]), median(s[(n + 1) // 2:])

def summary(values):
    """Dictionary of the lesson's summaries of centre and spread."""
    q1, q3 = quartiles(values)
    return {
        "mean": mean(values),
        "median": median(values),
        "modes": statistics.multimode(values),
        "range": max(values) - min(values),
        "Q1": q1,
        "Q3": q3,
        "IQR": q3 - q1,
    }

def outliers(values):
    """Values more than 1.5 IQR below Q1 or above Q3 (Tukey's rule)."""
    q1, q3 = quartiles(values)
    spread = q3 - q1
    low, high = q1 - Fraction(3, 2) * spread, q3 + Fraction(3, 2) * spread
    return [v for v in values if v < low or v > high]

def show(x):
    """A Fraction as an integer, or a decimal to 2 places."""
    x = Fraction(x)
    return str(x.numerator) if x.denominator == 1 else f"{float(x):.2f}"

def main():
    s = summary(CLASS_A)
    print("Minutes to school, 15 students:", " ".join(map(str, CLASS_A)))
    print(f"  mean {show(s['mean'])}  median {show(s['median'])}  mode "
          f"{s['modes']}  range {s['range']}")
    print(f"  Q1 {show(s['Q1'])}  Q3 {show(s['Q3'])}  IQR {show(s['IQR'])}"
          f"  outliers (1.5 IQR rule): {outliers(CLASS_A)}")
    print()
    print("Make the largest value more extreme:")
    print("  largest   mean    median   range   IQR")
    for big in (45, 90, 180, 900):
        data = CLASS_A[:-1] + [big]
        t = summary(data)
        print(f"  {big:7d}  {show(t['mean']):>6}   {show(t['median']):>5}"
              f"   {t['range']:5d}   {show(t['IQR']):>3}")
    print()
    r2 = sorted(ROUTE_2)
    py_ex = statistics.quantiles(r2, n=4, method="exclusive")
    py_in = statistics.quantiles(r2, n=4, method="inclusive")
    q1, q3 = quartiles(r2)
    print("Quartile conventions differ. Route 2, sorted:")
    print("  " + " ".join(map(str, r2)))
    print(f"  halves (this lesson): Q1 = {show(q1)}, Q3 = {show(q3)}")
    print(f"  Python 'exclusive':   Q1 = {py_ex[0]:g}, Q3 = {py_ex[2]:g}")
    print(f"  Python 'inclusive':   Q1 = {py_in[0]:g}, Q3 = {py_in[2]:g}")
    print()
    print("Two bus routes, waiting times over 10 days (minutes):")
    for name, data in (("route 1", ROUTE_1), ("route 2", ROUTE_2)):
        t = summary(data)
        print(f"  {name}: mean {show(t['mean'])}, median {show(t['median'])},"
              f" range {t['range']}, IQR {show(t['IQR'])}")
    print("  Same centre, very different reliability.")
    print()
    print("Epistemic status: exact arithmetic on small illustrative data sets.")
    print("That the median ignores the largest value is proved in the lesson")
    print("(for three or more values); this run only illustrates it, and the")
    print("IQR's behaviour here depends on the data. Quartiles also depend on")
    print("a chosen convention.")

if __name__ == "__main__":
    main()

"""Two "best centres": the mean and the median as answers to different questions.

Pick a number c to stand for the whole data set. If we measure the error
by the total squared distance  sum (x - c)^2,  the best c is the mean. If we
measure it by the total distance  sum |x - c|,  the best c is the median.
This experiment scans c over a fine grid (exact fractions) for the lesson's
data, finds where each total is smallest, and checks the identity
sum (x - c)^2 = sum (x - mean)^2 + n (mean - c)^2 at every grid point.
"""

from fractions import Fraction

DATA = [5, 8, 10, 10, 12, 12, 12, 15, 15, 18, 20, 20, 25, 30, 45]

def total_distance(values, c):
    """Sum of |x - c|."""
    return sum(abs(Fraction(x) - c) for x in values)

def total_squared_distance(values, c):
    """Sum of (x - c)^2."""
    return sum((Fraction(x) - c) ** 2 for x in values)

def grid(low, high, step):
    """Exact grid low, low + step, ..., high."""
    count = int((Fraction(high) - Fraction(low)) / Fraction(step))
    return [Fraction(low) + k * Fraction(step) for k in range(count + 1)]

def minimisers(values, cost, candidates):
    """The candidates at which cost(values, c) is smallest."""
    costs = [(cost(values, c), c) for c in candidates]
    best = min(k for k, _ in costs)
    return [c for k, c in costs if k == best]

def identity_holds(values, candidates):
    """Check sum (x-c)^2 == sum (x-m)^2 + n (m-c)^2 at every candidate."""
    n = len(values)
    m = Fraction(sum(values), n)
    base = total_squared_distance(values, m)
    return all(total_squared_distance(values, c) == base + n * (m - c) ** 2
               for c in candidates)

def show(x):
    """A Fraction as an integer or a decimal to 2 places."""
    x = Fraction(x)
    return str(x.numerator) if x.denominator == 1 else f"{float(x):.2f}"

def main():
    m = Fraction(sum(DATA), len(DATA))
    print("Data (minutes to school):", " ".join(map(str, DATA)))
    print(f"mean = {m} = {show(m)}, median = 15")
    print()
    print("     c    total distance   total squared distance")
    for c in (10, 15, m, 20, 25):
        print(f"  {show(c):>5}   {show(total_distance(DATA, c)):>14}   "
              f"{show(total_squared_distance(DATA, c)):>22}")
    print()
    cands = grid(0, 50, Fraction(1, 4))
    print(f"Scanning c = 0, 1/4, 1/2, ..., 50 ({len(cands)} values):")
    best_abs = minimisers(DATA, total_distance, cands)
    best_sq = minimisers(DATA, total_squared_distance, cands)
    print(f"  smallest total distance at c = {[show(c) for c in best_abs]}")
    print(f"  smallest total squared distance at c = "
          f"{[show(c) for c in best_sq]}")
    print("    (the grid point nearest to the mean)")
    print(f"  identity checked at every grid point: "
          f"{identity_holds(DATA, cands)}")
    even = [1, 2, 6, 10]
    flat = minimisers(even, total_distance, grid(0, 12, Fraction(1, 2)))
    print(f"For {even} (even count) every c from {show(flat[0])} to "
          f"{show(flat[-1])} ties:")
    print("  the total distance is flat between the two middle values.")
    print()
    print("Epistemic status: exact arithmetic, exhaustive over a finite grid")
    print("of candidates. The lesson proves that the mean minimises the total")
    print("squared distance for every real c (the identity above), and why")
    print("the median minimises the total distance.")

if __name__ == "__main__":
    main()