Practical Lab

12 experiments, each set out as 1. Question, 2. Aim, 3. Steps, 4. Programme, 5. Execution and Results.

Code lives in labs/course-12a-ml/.

pip install -r tools/requirements.txt
python3 tools/data-science/run_ml_labs.py

Two cross-course checks run as part of the suite

If the two courses ever disagree, one of them is wrong, and the suite says so. That is the point of reusing the datasets rather than inventing new ones.

The three datasets

Experiment 1 — Import and export data with pandas

1. Question

Import and export data in CSV, Excel, JSON and Parquet with pandas.

2. Aim

Round-trip a frame through each format, and catch what each one silently changes.

3. Steps

1. Write a CSV and read it back.
2. See what the default index does.
3. See which types a CSV loses.
4. Write and read Excel, with several sheets.
5. Compare the JSON orientations.
6. Round-trip through Parquet.
7. Read part of a large file.

4. Programme

"""Experiment 1 — Import and export data with pandas.

The first practical, and the one that decides whether the rest of the course is
frustrating. Every format has a way of silently changing your data on the round
trip, and this script demonstrates each one failing and then fixing it.
"""
import io
import json
import pathlib
import tempfile

import pandas as pd

from fixtures import STUDY, churn

TMP = pathlib.Path(tempfile.mkdtemp(prefix="ml_lab1_"))


def csv_round_trip():
    path = TMP / "study.csv"
    STUDY.to_csv(path, index=False)
    back = pd.read_csv(path)

    assert back.shape == STUDY.shape == (10, 2)
    assert list(back.columns) == ["hours", "score"]
    assert back["score"].sum() == STUDY["score"].sum() == 710
    assert back.equals(STUDY)

    print(f"  CSV round trip: {back.shape}, score total {back['score'].sum()} -- identical")


def index_false_matters():
    """to_csv(index=True) adds a phantom column that read_csv turns into data."""
    with_index = TMP / "with_index.csv"
    STUDY.to_csv(with_index)                      # index=True is the DEFAULT
    back = pd.read_csv(with_index)

    assert back.shape == (10, 3), back.shape
    assert back.columns[0] == "Unnamed: 0", list(back.columns)

    fixed = pd.read_csv(with_index, index_col=0)
    assert fixed.shape == (10, 2)
    assert fixed.equals(STUDY)

    print(f"  to_csv() with the default index -> read back as {back.shape}, "
          f"first column '{back.columns[0]}'")
    print("       'Unnamed: 0' is the commonest junk column in data science.")
    print("       Write with index=False, or read with index_col=0")


def dtypes_are_not_preserved_by_csv():
    """CSV is text. Types are GUESSED on the way back in."""
    df = pd.DataFrame({"id": ["001", "002", "010"],
                       "flag": [True, False, True],
                       "when": pd.to_datetime(["2026-01-15", "2026-02-10", "2026-04-05"])})
    path = TMP / "types.csv"
    df.to_csv(path, index=False)
    back = pd.read_csv(path)

    # The zero-padded id became an integer and lost its padding.
    assert str(df["id"].dtype) == "str" or df["id"].dtype == object
    assert back["id"].tolist() == [1, 2, 10], back["id"].tolist()
    assert back["id"].tolist() != df["id"].tolist()

    # The date became a plain string. (In pandas 3 the string dtype prints as
    # "str" rather than the old "object" -- assert the property, not the name.)
    assert not pd.api.types.is_datetime64_any_dtype(back["when"])
    assert isinstance(back["when"].iloc[0], str), back["when"].iloc[0]

    fixed = pd.read_csv(path, dtype={"id": str}, parse_dates=["when"])
    assert fixed["id"].tolist() == ["001", "002", "010"]
    assert pd.api.types.is_datetime64_any_dtype(fixed["when"])

    print(f"  zero-padded id '001','002','010' -> read back as {back['id'].tolist()}")
    print(f"  dates -> dtype {back['when'].dtype} (text, not a date)")
    print("  with dtype={'id': str}, parse_dates=['when'] -> both correct")
    print("       CSV stores TEXT. Every type is re-guessed on read, and a")
    print("       pin code, phone number or account id loses its leading zeros")


def excel_and_multiple_sheets():
    path = TMP / "data.xlsx"
    with pd.ExcelWriter(path) as xl:
        STUDY.to_excel(xl, sheet_name="study", index=False)
        churn().head(20).to_excel(xl, sheet_name="churn", index=False)

    sheets = pd.read_excel(path, sheet_name=None)     # None -> a dict of ALL
    assert set(sheets) == {"study", "churn"}
    assert sheets["study"].shape == (10, 2)
    assert sheets["churn"].shape == (20, 4)

    one = pd.read_excel(path, sheet_name="study")
    assert one.equals(STUDY)

    print(f"  Excel: {len(sheets)} sheets -> {[f'{k} {v.shape}' for k, v in sheets.items()]}")
    print("       sheet_name=None returns a DICT of every sheet. Unlike CSV,")
    print("       Excel preserves dtypes, which is its one real advantage")


def json_orientations():
    """The 'orient' argument changes the file completely."""
    small = STUDY.head(3)
    shapes = {}
    for orient in ("records", "columns", "index", "split", "values"):
        text = small.to_json(orient=orient)
        shapes[orient] = len(text)
        if orient in ("records", "columns", "index", "split"):
            back = pd.read_json(io.StringIO(text), orient=orient)
            assert back.shape == small.shape, (orient, back.shape)

    assert json.loads(small.to_json(orient="records"))[0] == {"hours": 2, "score": 52}
    assert list(json.loads(small.to_json(orient="columns"))) == ["hours", "score"]

    print("  the same 3 rows, five JSON orientations:")
    for orient, size in shapes.items():
        print(f"    orient='{orient}':{'':<{10 - len(orient)}} {size} chars")
    print("       'records' is the one APIs use -- a list of objects. You must")
    print("       pass the SAME orient to read_json, or the frame comes back")
    print("       transposed or empty")


def parquet_preserves_everything():
    """The format to use between Python programs."""
    df = pd.DataFrame({"id": ["001", "002"],
                       "when": pd.to_datetime(["2026-01-15", "2026-02-10"]),
                       "flag": [True, False],
                       "value": [1.5, 2.5]})
    path = TMP / "data.parquet"
    try:
        df.to_parquet(path, index=False)
    except (ImportError, ValueError) as exc:
        print(f"  parquet unavailable ({type(exc).__name__}) -- skipped")
        return

    back = pd.read_parquet(path)
    assert back["id"].tolist() == ["001", "002"], "leading zeros SURVIVE"
    assert pd.api.types.is_datetime64_any_dtype(back["when"]), "dates survive"
    assert back["flag"].dtype == bool
    assert back.equals(df)

    csv_size = len(df.to_csv(index=False).encode())
    parquet_size = path.stat().st_size
    print(f"  Parquet: every dtype preserved, frame identical on round trip")
    print(f"    csv {csv_size} bytes, parquet {parquet_size} bytes "
          f"(parquet wins only on LARGE data -- it has a header)")
    print("       columnar, typed, compressed. Use it between programs;")
    print("       use CSV only when a human or a foreign tool must read it")


def reading_a_subset_of_a_large_file():
    """usecols and nrows -- what you do when the file will not fit."""
    big = churn(n=400)
    path = TMP / "big.csv"
    big.to_csv(path, index=False)

    head = pd.read_csv(path, nrows=5)
    assert head.shape == (5, 4)

    two_cols = pd.read_csv(path, usecols=["tenure_months", "churned"])
    assert two_cols.shape == (400, 2), two_cols.shape

    chunks = list(pd.read_csv(path, chunksize=150))
    assert [len(c) for c in chunks] == [150, 150, 100]
    assert sum(len(c) for c in chunks) == 400

    print(f"  nrows=5              -> {head.shape}")
    print(f"  usecols=[2 columns]  -> {two_cols.shape}")
    print(f"  chunksize=150        -> {[len(c) for c in chunks]} rows per chunk")
    print("       read the head FIRST to see the columns, then usecols to load")
    print("       only what you need. chunksize when it will not fit in memory")


def main():
    print("Experiment 1 -- Importing and exporting data with pandas")
    # Step 1: Write a CSV and read it back
    csv_round_trip()
    # Step 2: See what the default index does
    index_false_matters()
    # Step 3: See which types a CSV loses
    dtypes_are_not_preserved_by_csv()
    # Step 4: Write and read Excel, with several sheets
    excel_and_multiple_sheets()
    # Step 5: Compare the JSON orientations
    json_orientations()
    # Step 6: Round-trip through Parquet
    parquet_preserves_everything()
    # Step 7: Read part of a large file
    reading_a_subset_of_a_large_file()


if __name__ == "__main__":
    main()

5. Execution and Results

Experiment 1 -- Importing and exporting data with pandas
  CSV round trip: (10, 2), score total 710 -- identical
  to_csv() with the default index -> read back as (10, 3), first column 'Unnamed: 0'
       'Unnamed: 0' is the commonest junk column in data science.
       Write with index=False, or read with index_col=0
  zero-padded id '001','002','010' -> read back as [1, 2, 10]
  dates -> dtype str (text, not a date)
  with dtype={'id': str}, parse_dates=['when'] -> both correct
       CSV stores TEXT. Every type is re-guessed on read, and a
       pin code, phone number or account id loses its leading zeros
  Excel: 2 sheets -> ['study (10, 2)', 'churn (20, 4)']
       sheet_name=None returns a DICT of every sheet. Unlike CSV,
       Excel preserves dtypes, which is its one real advantage
  the same 3 rows, five JSON orientations:
    orient='records':    70 chars
    orient='columns':    60 chars
    orient='index':      82 chars
    orient='split':      75 chars
    orient='values':     22 chars
       'records' is the one APIs use -- a list of objects. You must
       pass the SAME orient to read_json, or the frame comes back
       transposed or empty
  Parquet: every dtype preserved, frame identical on round trip
    csv 68 bytes, parquet 2661 bytes (parquet wins only on LARGE data -- it has a header)
       columnar, typed, compressed. Use it between programs;
       use CSV only when a human or a foreign tool must read it
  nrows=5              -> (5, 4)
  usecols=[2 columns]  -> (400, 2)
  chunksize=150        -> [150, 150, 100] rows per chunk
       read the head FIRST to see the columns, then usecols to load
       only what you need. chunksize when it will not fit in memory

Experiment 2 — Data pre-processing techniques

1. Question

Apply data pre-processing techniques: scaling, encoding and imputation.

2. Aim

Scale, encode and impute, and measure what each choice does to the data and to a model.

3. Steps

1. Compare the scalers, with an outlier.
2. Encode the categories.
3. Meet the dummy variable trap.
4. See mean imputation shrink the variance.
5. Keep a was-missing flag.
6. Scale after splitting, in a pipeline.

4. Programme

"""Experiment 2 — Data pre-processing techniques.

Unit 2 section 2.1 claims the ORDER of these steps matters and that getting it
wrong corrupts your evaluation silently. This script measures that, rather than
asserting it in prose.
"""
import numpy as np
import pandas as pd
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import (MinMaxScaler, OneHotEncoder, OrdinalEncoder,
                                   RobustScaler, StandardScaler)

from fixtures import RANDOM_STATE, churn


def scaling_compared():
    """Standard, min-max and robust, on data with one deliberate outlier."""
    x = np.array([10.0, 12.0, 11.0, 13.0, 12.0, 11.0, 250.0]).reshape(-1, 1)

    std = StandardScaler().fit_transform(x).ravel()
    mm = MinMaxScaler().fit_transform(x).ravel()
    rob = RobustScaler().fit_transform(x).ravel()

    assert round(std.mean(), 10) == 0.0
    assert round(std.std(), 10) == 1.0
    assert round(mm.min(), 10) == 0.0 and round(mm.max(), 10) == 1.0

    # The outlier crushes min-max: the six normal points land in a tiny band.
    normal_spread_mm = mm[:6].max() - mm[:6].min()
    normal_spread_rob = rob[:6].max() - rob[:6].min()
    # min-max: the six normal points span 3 of a total range of 240 -> 0.0125.
    # robust: they span 3 divided by the IQR of 1.5 -> 2.0. A factor of 160.
    assert round(normal_spread_mm, 4) == 0.0125, round(normal_spread_mm, 4)
    assert round(normal_spread_rob, 4) == 2.0, round(normal_spread_rob, 4)
    assert round(normal_spread_rob / normal_spread_mm, 6) == 160.0

    print("  x = [10, 12, 11, 13, 12, 11, 250]   (one outlier)")
    print(f"    {'standardised':16s} {np.round(std, 3)}")
    print(f"    {'min-max':16s} {np.round(mm, 4)}")
    print(f"    {'robust':16s} {np.round(rob, 3)}")
    print(f"  spread of the six NORMAL points after scaling:")
    print(f"    min-max {normal_spread_mm:.4f}   robust {normal_spread_rob:.4f}")
    print(f"    robust keeps them {normal_spread_rob / normal_spread_mm:.0f}x further apart")
    print("       min-max squeezed six distinct values into a 0.0125-wide band,")
    print("       because the outlier owns the whole range. Robust scaling uses")
    print("       the median and IQR, so the outlier does not set the scale")


def encoding_categoricals():
    """One-hot for nominal, ordinal for ordinal -- and why swapping them hurts."""
    colours = pd.DataFrame({"colour": ["Red", "Green", "Blue", "Red"]})
    sizes = pd.DataFrame({"size": ["S", "M", "L", "M"]})

    onehot = OneHotEncoder(sparse_output=False).fit_transform(colours)
    assert onehot.shape == (4, 3), onehot.shape
    assert onehot.sum(axis=1).tolist() == [1.0, 1.0, 1.0, 1.0], "exactly one 1 per row"

    ordinal = OrdinalEncoder(categories=[["S", "M", "L"]]).fit_transform(sizes)
    assert ordinal.ravel().tolist() == [0.0, 1.0, 2.0, 1.0]

    # The mistake: label-encoding a NOMINAL feature invents distances.
    bad = OrdinalEncoder().fit_transform(colours).ravel()
    order = dict(zip(colours["colour"], bad))
    d_red_green = abs(order["Red"] - order["Green"])
    d_red_blue = abs(order["Red"] - order["Blue"])
    assert d_red_green != d_red_blue, "the encoding has invented a distance"

    print(f"  nominal 'colour' one-hot -> shape {onehot.shape}, one 1 per row")
    print(f"  ordinal 'size' S<M<L     -> {ordinal.ravel().tolist()}  (order kept)")
    print(f"  colour LABEL-encoded     -> {bad.tolist()}")
    print(f"    distance Red-Green {d_red_green:.0f}, Red-Blue {d_red_blue:.0f}")
    print("       the encoder just told every distance-based model that Red is")
    print("       nearer one colour than another. That is not a fact about")
    print("       colours; it is a fact about alphabetical order")


def the_dummy_variable_trap():
    """k one-hot columns are perfectly collinear. Drop one for linear models."""
    colours = pd.DataFrame({"colour": ["Red", "Green", "Blue", "Red", "Green"]})

    full = OneHotEncoder(sparse_output=False).fit_transform(colours)
    dropped = OneHotEncoder(sparse_output=False, drop="first").fit_transform(colours)

    assert full.shape == (5, 3) and dropped.shape == (5, 2)
    # Every row of the full encoding sums to 1 -- perfect collinearity with the
    # intercept column a linear model adds.
    assert np.allclose(full.sum(axis=1), 1.0)
    # Dropping one column loses NO information: the dropped category is the
    # one where all remaining columns are zero.
    baseline_rows = (dropped.sum(axis=1) == 0).sum()
    assert baseline_rows == 1, baseline_rows

    print(f"  3 categories -> full one-hot {full.shape}, drop='first' {dropped.shape}")
    print(f"    every full row sums to 1.0 -> collinear with the intercept")
    print(f"    dropped encoding: {baseline_rows} row is all-zero = the baseline category")
    print("       drop one for LINEAR and LOGISTIC regression. Trees and")
    print("       regularised models do not care")


def imputation_shrinks_variance():
    """Mean imputation is not free -- it distorts the distribution."""
    rng = np.random.default_rng(RANDOM_STATE)
    full = rng.normal(50, 10, 200)
    holed = full.copy()
    holed[rng.choice(200, 60, replace=False)] = np.nan

    imputed = SimpleImputer(strategy="mean").fit_transform(
        holed.reshape(-1, 1)).ravel()

    assert np.isnan(holed).sum() == 60
    assert not np.isnan(imputed).any()
    assert round(abs(imputed.mean() - np.nanmean(holed)), 10) == 0.0, \
        "mean imputation preserves the MEAN exactly"
    assert imputed.std() < np.nanstd(holed), "but it SHRINKS the spread"

    print(f"  200 values, 60 made missing (30%):")
    print(f"    observed      mean {np.nanmean(holed):7.4f}  sd {np.nanstd(holed):7.4f}")
    print(f"    mean-imputed  mean {imputed.mean():7.4f}  sd {imputed.std():7.4f}")
    print(f"    sd fell by {(1 - imputed.std() / np.nanstd(holed)) * 100:.1f}%")
    print("       the mean is preserved exactly and the SPREAD is not. Sixty")
    print("       identical values were inserted at the centre. Every")
    print("       correlation involving this column is now weaker than it was")


def a_was_missing_flag_keeps_the_signal():
    """When missingness is informative, the flag is worth more than the value."""
    rng = np.random.default_rng(RANDOM_STATE)
    n = 400
    y = rng.integers(0, 2, n)
    # Income is missing FAR more often for class 1 -- missing-not-at-random.
    income = rng.normal(50000, 12000, n)
    missing = rng.random(n) < np.where(y == 1, 0.6, 0.1)
    income[missing] = np.nan

    imputed = pd.Series(income).fillna(np.nanmean(income)).to_numpy()
    flag = missing.astype(int)

    Xa = imputed.reshape(-1, 1)
    Xb = np.c_[imputed, flag]
    cv_a = cross_val_score(LogisticRegression(), Xa, y, cv=5).mean()
    cv_b = cross_val_score(LogisticRegression(), Xb, y, cv=5).mean()

    assert cv_b > cv_a + 0.10, (cv_a, cv_b)

    print(f"  income missing 60% of the time for class 1, 10% for class 0:")
    print(f"    imputed value only        CV accuracy {cv_a:.4f}")
    print(f"    imputed value + was_missing flag       {cv_b:.4f}")
    print(f"    the flag alone is worth {cv_b - cv_a:+.4f}")
    print("       WHY a value is missing can matter more than what it was.")
    print("       Imputing erases that; one extra column keeps it")


def scaling_before_splitting_leaks():
    """The rule from unit-2.md 2.1, measured on the churn data."""
    df = churn()
    X = df[["tenure_months", "support_calls", "irrelevant"]].to_numpy()
    y = df["churned"].to_numpy()

    # WRONG: scale the whole dataset, then split.
    X_all_scaled = StandardScaler().fit_transform(X)
    Xtr_w, Xte_w, ytr_w, yte_w = train_test_split(
        X_all_scaled, y, test_size=0.25, random_state=RANDOM_STATE, stratify=y)
    leaked = LogisticRegression(random_state=RANDOM_STATE).fit(
        Xtr_w, ytr_w).score(Xte_w, yte_w)

    # RIGHT: split first, and let a Pipeline fit the scaler on train only.
    Xtr, Xte, ytr, yte = train_test_split(
        X, y, test_size=0.25, random_state=RANDOM_STATE, stratify=y)
    clean = make_pipeline(
        StandardScaler(),
        LogisticRegression(random_state=RANDOM_STATE)).fit(Xtr, ytr).score(Xte, yte)

    # The scaler's own parameters differ, which is the mechanism of the leak.
    mean_all = StandardScaler().fit(X).mean_
    mean_train = StandardScaler().fit(Xtr).mean_
    assert not np.allclose(mean_all, mean_train), \
        "the two scalers learned DIFFERENT means -- that difference is the leak"

    print(f"  scale-then-split (leaky) test accuracy : {leaked:.4f}")
    print(f"  split-then-scale (Pipeline)            : {clean:.4f}")
    print(f"  scaler mean fitted on ALL data   : {np.round(mean_all, 4)}")
    print(f"  scaler mean fitted on TRAIN only : {np.round(mean_train, 4)}")
    print("       on this data the two scores happen to be close -- the point is")
    print("       NOT that leakage always inflates the number, but that the")
    print("       leaky score is not an estimate of anything, because the")
    print("       transformer saw the test rows. Use a Pipeline and the question")
    print("       cannot arise")


def main():
    print("Experiment 2 -- Data pre-processing techniques")
    # Step 1: Compare the scalers, with an outlier
    scaling_compared()
    # Step 2: Encode the categories
    encoding_categoricals()
    # Step 3: Meet the dummy variable trap
    the_dummy_variable_trap()
    # Step 4: See mean imputation shrink the variance
    imputation_shrinks_variance()
    # Step 5: Keep a was-missing flag
    a_was_missing_flag_keeps_the_signal()
    # Step 6: Scale after splitting, in a pipeline
    scaling_before_splitting_leaks()


if __name__ == "__main__":
    main()

5. Execution and Results

Experiment 2 -- Data pre-processing techniques
  x = [10, 12, 11, 13, 12, 11, 250]   (one outlier)
    standardised     [-0.426 -0.402 -0.414 -0.39  -0.402 -0.414  2.449]
    min-max          [0.     0.0083 0.0042 0.0125 0.0083 0.0042 1.    ]
    robust           [ -1.333   0.     -0.667   0.667   0.     -0.667 158.667]
  spread of the six NORMAL points after scaling:
    min-max 0.0125   robust 2.0000
    robust keeps them 160x further apart
       min-max squeezed six distinct values into a 0.0125-wide band,
       because the outlier owns the whole range. Robust scaling uses
       the median and IQR, so the outlier does not set the scale
  nominal 'colour' one-hot -> shape (4, 3), one 1 per row
  ordinal 'size' S<M<L     -> [0.0, 1.0, 2.0, 1.0]  (order kept)
  colour LABEL-encoded     -> [2.0, 1.0, 0.0, 2.0]
    distance Red-Green 1, Red-Blue 2
       the encoder just told every distance-based model that Red is
       nearer one colour than another. That is not a fact about
       colours; it is a fact about alphabetical order
  3 categories -> full one-hot (5, 3), drop='first' (5, 2)
    every full row sums to 1.0 -> collinear with the intercept
    dropped encoding: 1 row is all-zero = the baseline category
       drop one for LINEAR and LOGISTIC regression. Trees and
       regularised models do not care
  200 values, 60 made missing (30%):
    observed      mean 49.9646  sd  8.7813
    mean-imputed  mean 49.9646  sd  7.3470
    sd fell by 16.3%
       the mean is preserved exactly and the SPREAD is not. Sixty
       identical values were inserted at the centre. Every
       correlation involving this column is now weaker than it was
  income missing 60% of the time for class 1, 10% for class 0:
    imputed value only        CV accuracy 0.4550
    imputed value + was_missing flag       0.7575
    the flag alone is worth +0.3025
       WHY a value is missing can matter more than what it was.
       Imputing erases that; one extra column keeps it
  scale-then-split (leaky) test accuracy : 0.9400
  split-then-scale (Pipeline)            : 0.9400
  scaler mean fitted on ALL data   : [ 2.14935e+01  2.59040e+00 -9.10000e-03]
  scaler mean fitted on TRAIN only : [21.263   2.6028 -0.0415]
       on this data the two scores happen to be close -- the point is
       NOT that leakage always inflates the number, but that the
       leaky score is not an estimate of anything, because the
       transformer saw the test rows. Use a Pipeline and the question
       cannot arise

Experiment 3 — Dimensionality reduction with PCA

1. Question

Reduce the dimensions of the iris data with Principal Component Analysis.

2. Aim

Compute the components and their explained variance, choose how many to keep, and see what PCA cannot do.

3. Steps

1. Compute the explained variance.
2. See why scikit-learn's variances sum to more than 4.
3. Choose the number of components.
4. Compare PCA with LDA.
5. Standardise before PCA.
6. Measure the reconstruction error.

4. Programme

"""Experiment 3 — Dimensionality reduction with PCA.

Reproduces every figure in unit-2.md section 2.10, including the check that the
correlation-matrix eigenvalues sum to exactly p, and the disagreement between
the Kaiser criterion and the 95%-variance rule.
"""
import numpy as np
from sklearn.decomposition import PCA
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from fixtures import RANDOM_STATE, iris_frame

_, DATA = iris_frame()
X = DATA.data.to_numpy()
y = DATA.target


def explained_variance():
    """The table in unit-2.md 2.10."""
    Z = StandardScaler().fit_transform(X)
    pca = PCA().fit(Z)
    ratio = pca.explained_variance_ratio_
    cumulative = np.cumsum(ratio)

    assert [round(v, 4) for v in ratio] == [0.7296, 0.2285, 0.0367, 0.0052]
    assert round(cumulative[1], 4) == 0.9581, round(cumulative[1], 4)
    assert round(cumulative[-1], 10) == 1.0

    # Eigenvalues of the CORRELATION matrix, which is what Kaiser refers to.
    eig = np.sort(np.linalg.eigvalsh(np.corrcoef(X, rowvar=False)))[::-1]
    assert [round(v, 4) for v in eig] == [2.9185, 0.9140, 0.1468, 0.0207]
    assert round(eig.sum(), 10) == 4.0, "eigenvalues sum to p -- a free check"
    assert np.allclose(eig / eig.sum(), ratio, atol=1e-10), \
        "the RATIOS are identical however you compute the eigenvalues"

    print("  component  eigenvalue  explained  cumulative")
    for i, (e, r, c) in enumerate(zip(eig, ratio, cumulative), 1):
        print(f"    PC{i}       {e:9.4f}   {r * 100:7.2f}%   {c * 100:8.2f}%")
    print(f"  eigenvalues sum to {eig.sum():.4f} = p = {X.shape[1]}  <- check")
    print(f"  TWO components carry {cumulative[1] * 100:.2f}% of the variance")


def sklearn_uses_n_minus_one():
    """Why sklearn's explained_variance_ differs slightly, stated precisely."""
    Z = StandardScaler().fit_transform(X)
    pca = PCA().fit(Z)
    eig_corr = np.sort(np.linalg.eigvalsh(np.corrcoef(X, rowvar=False)))[::-1]

    assert [round(v, 4) for v in pca.explained_variance_] == \
        [2.9381, 0.9202, 0.1477, 0.0209]
    assert round(pca.explained_variance_.sum(), 4) == 4.0268
    assert round(eig_corr.sum(), 4) == 4.0

    n = len(X)
    assert np.allclose(pca.explained_variance_, eig_corr * n / (n - 1), atol=1e-9), \
        "sklearn divides by n-1; StandardScaler divided by n"

    print(f"  sklearn explained_variance_ : "
          f"{[round(v, 4) for v in pca.explained_variance_]}  sum "
          f"{pca.explained_variance_.sum():.4f}")
    print(f"  correlation eigenvalues     : "
          f"{[round(v, 4) for v in eig_corr]}  sum {eig_corr.sum():.4f}")
    print(f"  the factor is exactly n/(n-1) = {n}/{n - 1} = {n / (n - 1):.6f}")
    print("       the RATIOS are identical, which is why the ratio is what you")
    print("       report and the raw eigenvalue is not")


def how_many_components():
    """Kaiser says 1, the 95% rule says 2. They disagree, and that is normal."""
    eig = np.sort(np.linalg.eigvalsh(np.corrcoef(X, rowvar=False)))[::-1]
    cumulative = np.cumsum(eig / eig.sum())

    kaiser = int((eig > 1).sum())
    rule_90 = int(np.searchsorted(cumulative, 0.90) + 1)
    rule_95 = int(np.searchsorted(cumulative, 0.95) + 1)

    assert kaiser == 1, kaiser
    assert rule_90 == 2 and rule_95 == 2

    print(f"  Kaiser (eigenvalue > 1)      -> {kaiser} component")
    print(f"  cumulative variance >= 90%   -> {rule_90} components")
    print(f"  cumulative variance >= 95%   -> {rule_95} components")
    print("       THE RULES DISAGREE, and that is normal -- Kaiser is known to")
    print("       under-select when p is small. k is a hyperparameter: state")
    print("       which rule you used, or validate it downstream")


def pca_is_unsupervised_and_it_shows():
    """PCA maximises variance, not separation. Compare with LDA."""
    from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

    Z = StandardScaler().fit_transform(X)
    pcs = PCA(n_components=2).fit_transform(Z)
    lds = LinearDiscriminantAnalysis(n_components=2).fit_transform(Z, y)

    # Both give 2-D; judge them by how separable the classes become.
    pca_cv = cross_val_score(LogisticRegression(max_iter=1000), pcs, y, cv=5).mean()
    lda_cv = cross_val_score(LogisticRegression(max_iter=1000), lds, y, cv=5).mean()
    raw_cv = cross_val_score(
        make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000)),
        X, y, cv=5).mean()

    assert pcs.shape == lds.shape == (150, 2)
    assert lda_cv >= pca_cv, (lda_cv, pca_cv)

    print(f"  logistic regression, 5-fold CV accuracy on 2 dimensions:")
    print(f"    PCA (unsupervised, maximises VARIANCE)   {pca_cv:.4f}")
    print(f"    LDA (supervised, maximises SEPARATION)   {lda_cv:.4f}")
    print(f"    all 4 original features                  {raw_cv:.4f}")
    print("       PCA never looks at y. When the goal is separating classes,")
    print("       LDA is the supervised alternative -- and here it matches or")
    print("       beats PCA on the same number of dimensions")


def pca_needs_standardising():
    """Without scaling, PCA reports whichever feature has the largest units."""
    # Blow up one feature's units: sepal length in micrometres.
    X_mixed = X.copy()
    X_mixed[:, 0] *= 10000

    raw = PCA().fit(X_mixed)
    scaled = PCA().fit(StandardScaler().fit_transform(X_mixed))

    # Unscaled, PC1 is essentially the inflated column alone.
    loading = abs(raw.components_[0])
    assert loading.argmax() == 0
    assert round(loading[0], 6) == 1.0, round(loading[0], 6)
    assert round(raw.explained_variance_ratio_[0], 6) == 1.0

    assert round(scaled.explained_variance_ratio_[0], 4) == 0.7296, \
        "standardised, the units cancel and the real answer returns"

    print("  sepal length rescaled to micrometres (x10,000):")
    print(f"    UNSCALED PCA: PC1 explains "
          f"{raw.explained_variance_ratio_[0] * 100:.4f}% and its loading on the")
    print(f"      inflated column is {loading[0]:.4f} -- PC1 IS that column")
    print(f"    STANDARDISED: PC1 explains "
          f"{scaled.explained_variance_ratio_[0] * 100:.2f}%, unchanged from before")
    print("       PCA maximises variance and variance has units. Without")
    print("       standardising it reports your choice of measurement scale")


def reconstruction_error():
    """What 'losing 4.19% of the information' actually means."""
    Z = StandardScaler().fit_transform(X)
    for k in (1, 2, 3, 4):
        pca = PCA(n_components=k).fit(Z)
        back = pca.inverse_transform(pca.transform(Z))
        mse = float(((Z - back) ** 2).mean())
        kept = pca.explained_variance_ratio_.sum()
        if k == 2:
            assert round(kept, 4) == 0.9581
            assert round(mse, 4) == 0.0419, round(mse, 4)
        if k == 4:
            assert round(mse, 10) == 0.0, "keeping every component loses nothing"
        print(f"    k={k}: variance kept {kept * 100:6.2f}%   "
              f"reconstruction MSE {mse:.4f}")
    print("       at k=2 the MSE is 0.0419 -- exactly 1 minus the 0.9581 kept.")
    print("       'Losing 4.19% of the variance' is literally the squared error")
    print("       you would make rebuilding the original four columns")


def main():
    print("Experiment 3 -- Principal Component Analysis")
    # Step 1: Compute the explained variance
    explained_variance()
    # Step 2: See why scikit-learn's variances sum to more than 4
    sklearn_uses_n_minus_one()
    # Step 3: Choose the number of components
    how_many_components()
    # Step 4: Compare PCA with LDA
    pca_is_unsupervised_and_it_shows()
    # Step 5: Standardise before PCA
    pca_needs_standardising()
    # Step 6: Measure the reconstruction error
    print("  reconstruction error by number of components:")
    reconstruction_error()


if __name__ == "__main__":
    main()

5. Execution and Results

Experiment 3 -- Principal Component Analysis
  component  eigenvalue  explained  cumulative
    PC1          2.9185     72.96%      72.96%
    PC2          0.9140     22.85%      95.81%
    PC3          0.1468      3.67%      99.48%
    PC4          0.0207      0.52%     100.00%
  eigenvalues sum to 4.0000 = p = 4  <- check
  TWO components carry 95.81% of the variance
  sklearn explained_variance_ : [np.float64(2.9381), np.float64(0.9202), np.float64(0.1477), np.float64(0.0209)]  sum 4.0268
  correlation eigenvalues     : [np.float64(2.9185), np.float64(0.914), np.float64(0.1468), np.float64(0.0207)]  sum 4.0000
  the factor is exactly n/(n-1) = 150/149 = 1.006711
       the RATIOS are identical, which is why the ratio is what you
       report and the raw eigenvalue is not
  Kaiser (eigenvalue > 1)      -> 1 component
  cumulative variance >= 90%   -> 2 components
  cumulative variance >= 95%   -> 2 components
       THE RULES DISAGREE, and that is normal -- Kaiser is known to
       under-select when p is small. k is a hyperparameter: state
       which rule you used, or validate it downstream
  logistic regression, 5-fold CV accuracy on 2 dimensions:
    PCA (unsupervised, maximises VARIANCE)   0.9133
    LDA (supervised, maximises SEPARATION)   0.9800
    all 4 original features                  0.9600
       PCA never looks at y. When the goal is separating classes,
       LDA is the supervised alternative -- and here it matches or
       beats PCA on the same number of dimensions
  sepal length rescaled to micrometres (x10,000):
    UNSCALED PCA: PC1 explains 100.0000% and its loading on the
      inflated column is 1.0000 -- PC1 IS that column
    STANDARDISED: PC1 explains 72.96%, unchanged from before
       PCA maximises variance and variance has units. Without
       standardising it reports your choice of measurement scale
  reconstruction error by number of components:
    k=1: variance kept  72.96%   reconstruction MSE 0.2704
    k=2: variance kept  95.81%   reconstruction MSE 0.0419
    k=3: variance kept  99.48%   reconstruction MSE 0.0052
    k=4: variance kept 100.00%   reconstruction MSE 0.0000
       at k=2 the MSE is 0.0419 -- exactly 1 minus the 0.9581 kept.
       'Losing 4.19% of the variance' is literally the squared error
       you would make rebuilding the original four columns

Also asserted:

• Reconstruction MSE equals 1 − variance kept, exactly. At k=2 the MSE is 0.0419 against 0.9581 kept. "Losing 4.19% of the variance" is literally the squared error of rebuilding the four original columns.

• Kaiser says 1 component, the 95% rule says 2. They disagree, and that is normal.

• PCA is unsupervised. On two dimensions, LDA scores 0.9800 against PCA's 0.9133 — and all four raw features score 0.9600.

• PCA without standardising reports your units. Rescale sepal length to micrometres and PC1 explains 100.0000% with a loading of 1.0000 on that one column.

• sklearn's explained_variance_ sums to 4.0268 rather than 4.0 because it divides by n−1 — exactly a factor of 150/149. The ratios are identical.

Experiment 4 — Data visualization techniques

1. Question

Apply data visualization techniques to explore a dataset.

2. Aim

Compute what each chart shows and check it, starting with Anscombe's quartet.

3. Steps

1. Plot Anscombe's quartet.
2. Vary the histogram's bin count.
3. Find outliers with a boxplot.
4. Read the scatter matrix and the correlations.
5. Plot the class balance first.
6. Sort x before a line plot.

4. Programme

"""Experiment 4 — Data visualization techniques.

A plot cannot be asserted, but the numbers behind it can, and the point of a
plot is nearly always a number you could have missed. Each function here
computes what its chart shows and checks it; the images are written to plots/.

The centrepiece is Anscombe's quartet, which is the single best argument for
plotting before modelling.
"""
import pathlib

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from fixtures import STUDY, churn, iris_frame

# [Changed: the charts went to output/, which now holds what each program printed;
# they are drawn into plots/, and the lab page shows them.]
OUT = pathlib.Path(__file__).parent / "plots"
OUT.mkdir(exist_ok=True)

# Anscombe's quartet: four datasets, identical summary statistics.
ANSCOMBE = {
    "I":   ([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
            [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68]),
    "II":  ([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
            [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74]),
    "III": ([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5],
            [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73]),
    "IV":  ([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8],
            [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89]),
}


def anscombes_quartet():
    """Four datasets. Same mean, variance, correlation and regression line."""
    stats = {}
    for name, (xs, ys) in ANSCOMBE.items():
        x, y = np.array(xs, dtype=float), np.array(ys, dtype=float)
        slope, intercept = np.polyfit(x, y, 1)
        stats[name] = (x.mean(), y.mean(), y.var(ddof=1),
                       np.corrcoef(x, y)[0, 1], slope, intercept)

    # The agreement is close but not exact, and the precision differs per
    # statistic. Assert exactly what is true, not the folklore:
    #   mean x   exact       9.0000 in all four
    #   slope    3 decimals  0.50009, 0.50000, 0.49973, 0.49991
    #   mean y   2 decimals  7.5009, 7.5009, 7.5000, 7.5009
    #   icept    2 decimals  3.00009, 3.00091, 3.00245, 3.00173
    #   r        2 decimals  0.81642, 0.81624, 0.81629, 0.81652
    #   var y    1 decimal   4.1273, 4.1276, 4.1226, 4.1232
    for name, s in stats.items():
        assert round(s[0], 4) == 9.0000, (name, s[0])        # mean x -- exact
        assert round(s[4], 3) == 0.500, (name, s[4])         # slope
        assert round(s[1], 2) == 7.50, (name, s[1])          # mean y
        assert round(s[5], 2) == 3.00, (name, s[5])          # intercept
        assert round(s[3], 2) == 0.82, (name, s[3])          # correlation
        assert round(s[2], 1) == 4.1, (name, s[2])           # var y

    fig, axes = plt.subplots(2, 2, figsize=(7, 6))
    for ax, (name, (xs, ys)) in zip(axes.ravel(), ANSCOMBE.items()):
        ax.scatter(xs, ys, color="#0f4c81")
        xx = np.linspace(3, 20, 2)
        ax.plot(xx, 3.0 + 0.5 * xx, color="#dc2626", linewidth=1)
        ax.set_title(f"Anscombe {name}")
        ax.set_xlim(2, 20); ax.set_ylim(2, 14)
    fig.suptitle("Identical statistics, four different datasets")
    fig.tight_layout()
    fig.savefig(OUT / "04_anscombe.png", dpi=110)
    plt.close(fig)

    print("  set  mean x  mean y   var y     r    slope  intercept")
    for name, s in stats.items():
        print(f"   {name:3s} {s[0]:7.2f} {s[1]:7.2f} {s[2]:7.2f} {s[3]:6.3f} "
              f"{s[4]:7.3f} {s[5]:9.2f}")
    print("       mean x is IDENTICAL, the slope agrees to three decimals, and")
    print("       mean y, intercept and r to two -- yet the four")
    print("       datasets look nothing alike: I is a genuine linear")
    print("       relationship, II is a parabola, III is a line plus one")
    print("       outlier, IV is a vertical stack plus one leverage point.")
    print("       PLOT BEFORE YOU MODEL. Summary statistics cannot see shape")


def histogram_and_the_bin_count():
    """The same data, four bin counts, four different impressions."""
    rng = np.random.default_rng(0)
    data = np.concatenate([rng.normal(20, 3, 200), rng.normal(35, 3, 200)])

    def interior_peaks(bins):
        hist, _ = np.histogram(data, bins=bins)
        return sum(1 for i in range(1, len(hist) - 1)
                   if hist[i] > hist[i - 1] and hist[i] > hist[i + 1])

    found = {bins: interior_peaks(bins) for bins in (2, 5, 15, 30, 60)}

    assert len(data) == 400
    assert found[2] == 0, "two bins cannot show ANY interior structure"
    assert found[5] == 2 and found[15] == 2, "the two real modes"
    assert found[30] == 4 and found[60] == 15, found
    assert found[60] > 7 * found[15], "60 bins invents peaks that are not there"

    fig, axes = plt.subplots(1, 4, figsize=(12, 3))
    for ax, bins in zip(axes, (2, 5, 15, 60)):
        ax.hist(data, bins=bins, color="#0f4c81")
        ax.set_title(f"{bins} bins")
    fig.tight_layout()
    fig.savefig(OUT / "04_histogram_bins.png", dpi=110)
    plt.close(fig)

    print("  a genuinely BIMODAL sample (peaks at 20 and 35), by bin count:")
    for bins, peaks in found.items():
        verdict = ("no structure visible" if peaks == 0 else
                   "correct" if peaks == 2 else "SPURIOUS peaks from noise")
        print(f"    {bins:2d} bins -> {peaks:2d} interior peak(s)   {verdict}")
    print("       2 bins shows nothing; 5 and 15 find both real modes; 60 finds")
    print("       FIFTEEN, which are sampling noise. The bin count is a CHOICE")
    print("       that changes the conclusion in both directions, which is why")
    print("       a histogram should never be the only thing you look at")


def boxplot_finds_outliers():
    """The IQR rule, computed -- Course 4's method, applied."""
    x = pd.Series([10, 12, 11, 13, 12, 11, 14, 12, 250])
    q1, q3 = x.quantile(0.25), x.quantile(0.75)
    iqr = q3 - q1
    lo, hi = q1 - 1.5 * iqr, q3 + 1.5 * iqr
    outliers = x[(x < lo) | (x > hi)]

    assert q1 == 11.0 and q3 == 13.0 and iqr == 2.0
    assert lo == 8.0 and hi == 16.0
    assert outliers.tolist() == [250]

    # The mean is destroyed by the outlier; the median is not.
    assert round(x.mean(), 4) == 38.3333
    assert x.median() == 12.0

    print(f"  Q1 {q1}, Q3 {q3}, IQR {iqr}  ->  fences [{lo}, {hi}]")
    print(f"  outliers: {outliers.tolist()}")
    print(f"  mean {x.mean():.4f} vs median {x.median():.1f}")
    print("       one value moved the mean by 26 and the median by 0. A box")
    print("       plot shows both the outlier and that robustness at once")


def scatter_matrix_and_correlation():
    """Which iris features are redundant -- the question a pair plot answers."""
    df, data = iris_frame()
    features = list(data.feature_names)
    corr = data.data.corr()

    petal = corr.loc["petal length (cm)", "petal width (cm)"]
    sepal_petal = corr.loc["sepal length (cm)", "petal length (cm)"]
    sepal_w = corr.loc["sepal length (cm)", "sepal width (cm)"]

    assert round(petal, 4) == 0.9629, round(petal, 4)
    assert round(sepal_petal, 4) == 0.8718, round(sepal_petal, 4)
    assert round(sepal_w, 4) == -0.1176, round(sepal_w, 4)
    assert petal > sepal_petal > abs(sepal_w)

    print("  iris feature correlations:")
    print(f"    petal length ~ petal width  {petal:+.4f}   <- nearly redundant")
    print(f"    sepal length ~ petal length {sepal_petal:+.4f}")
    print(f"    sepal length ~ sepal width  {sepal_w:+.4f}   <- almost none")
    print("       petal length and width carry nearly the same information,")
    print("       which is exactly why PCA compressed four columns to two")


def class_balance_should_be_plotted_first():
    """The first plot to make on any classification problem."""
    df = churn()
    counts = df["churned"].value_counts().sort_index()
    rate = df["churned"].mean()

    assert counts.tolist() == [340, 60]
    assert round(rate, 4) == 0.15
    assert round(counts[0] / len(df), 4) == 0.85, "the majority-class baseline"

    print(f"  class balance: {counts[0]} negative, {counts[1]} positive "
          f"({rate * 100:.0f}% positive)")
    print(f"  a model predicting 'never' scores {counts[0] / len(df) * 100:.0f}% accuracy")
    print("       plot this FIRST. It tells you what accuracy will mean before")
    print("       you fit anything, and it is one line of code")


def line_plot_needs_ordered_x():
    df = STUDY.sort_values("hours")
    assert df["hours"].is_monotonic_increasing
    assert df["score"].is_monotonic_increasing

    fig, ax = plt.subplots(figsize=(5, 3))
    ax.scatter(df["hours"], df["score"], color="#0f4c81")
    m, c = np.polyfit(df["hours"], df["score"], 1)
    ax.plot(df["hours"], m * df["hours"] + c, color="#dc2626")
    ax.set_xlabel("hours studied"); ax.set_ylabel("exam score")
    ax.set_title(f"score = {c:.2f} + {m:.2f} x hours")
    fig.tight_layout()
    fig.savefig(OUT / "04_scatter_fit.png", dpi=110)
    plt.close(fig)

    assert round(m, 4) == 4.3030 and round(c, 4) == 43.0303

    print(f"  study data: fitted line score = {c:.4f} + {m:.4f} x hours")
    print("       scatter plot for a RELATIONSHIP, with the fit drawn over it.")
    print("       Points only -- never join scattered points with lines")


def main():
    print("Experiment 4 -- Data visualization techniques")
    # Step 1: Plot Anscombe's quartet
    anscombes_quartet()
    # Step 2: Vary the histogram's bin count
    histogram_and_the_bin_count()
    # Step 3: Find outliers with a boxplot
    boxplot_finds_outliers()
    # Step 4: Read the scatter matrix and the correlations
    scatter_matrix_and_correlation()
    # Step 5: Plot the class balance first
    class_balance_should_be_plotted_first()
    # Step 6: Sort x before a line plot
    line_plot_needs_ordered_x()
    print(f"  charts written to {OUT.name}/")


if __name__ == "__main__":
    main()

5. Execution and Results

Experiment 4 -- Data visualization techniques
  set  mean x  mean y   var y     r    slope  intercept
   I      9.00    7.50    4.13  0.816   0.500      3.00
   II     9.00    7.50    4.13  0.816   0.500      3.00
   III    9.00    7.50    4.12  0.816   0.500      3.00
   IV     9.00    7.50    4.12  0.817   0.500      3.00
       mean x is IDENTICAL, the slope agrees to three decimals, and
       mean y, intercept and r to two -- yet the four
       datasets look nothing alike: I is a genuine linear
       relationship, II is a parabola, III is a line plus one
       outlier, IV is a vertical stack plus one leverage point.
       PLOT BEFORE YOU MODEL. Summary statistics cannot see shape
  a genuinely BIMODAL sample (peaks at 20 and 35), by bin count:
     2 bins ->  0 interior peak(s)   no structure visible
     5 bins ->  2 interior peak(s)   correct
    15 bins ->  2 interior peak(s)   correct
    30 bins ->  4 interior peak(s)   SPURIOUS peaks from noise
    60 bins -> 15 interior peak(s)   SPURIOUS peaks from noise
       2 bins shows nothing; 5 and 15 find both real modes; 60 finds
       FIFTEEN, which are sampling noise. The bin count is a CHOICE
       that changes the conclusion in both directions, which is why
       a histogram should never be the only thing you look at
  Q1 11.0, Q3 13.0, IQR 2.0  ->  fences [8.0, 16.0]
  outliers: [250]
  mean 38.3333 vs median 12.0
       one value moved the mean by 26 and the median by 0. A box
       plot shows both the outlier and that robustness at once
  iris feature correlations:
    petal length ~ petal width  +0.9629   <- nearly redundant
    sepal length ~ petal length +0.8718
    sepal length ~ sepal width  -0.1176   <- almost none
       petal length and width carry nearly the same information,
       which is exactly why PCA compressed four columns to two
  class balance: 340 negative, 60 positive (15% positive)
  a model predicting 'never' scores 85% accuracy
       plot this FIRST. It tells you what accuracy will mean before
       you fit anything, and it is one line of code
  study data: fitted line score = 43.0303 + 4.3030 x hours
       scatter plot for a RELATIONSHIP, with the fit drawn over it.
       Points only -- never join scattered points with lines
  charts written to plots/

Also measured: the histogram bin count changes the conclusion in both directions — 2 bins show no structure, 5 and 15 find both real modes, and 60 bins find fifteen, which are sampling noise. Plus the IQR outlier rule, iris correlations (petal length ~ petal width 0.9629), and the class balance plot that should be the first thing you make.

The three charts above are Anscombe's quartet, the histogram at four bin counts, and the scatter with its fitted line. Changed: they were written to output/, which now holds what each program printed; the program writes them to plots/.

Experiment 5 — Maximum likelihood estimation

1. Question

Estimate parameters by maximum likelihood.

2. Aim

Find the MLE of a coin's p two ways, and connect least squares to likelihood.

3. Steps

1. Estimate a coin's p by grid search and by calculus.
2. See why logs are taken.
3. Show least squares is the MLE under normal errors.
4. Estimate a normal mean and variance.

4. Programme

"""Experiment 5 — Maximum likelihood estimation.

unit-3.md section 3.7 derives the MLE for a coin by hand and claims two things
this script checks numerically:

  * the MLE of a Bernoulli parameter is the sample proportion (7 heads in 10
    gives p-hat = 0.7)
  * least squares IS the MLE for linear regression under normal errors

The second is the connection worth knowing and the one rarely demonstrated.
"""
import numpy as np
from scipy import optimize, stats

from fixtures import STUDY, study_xy


def coin_mle_by_grid_and_by_calculus():
    """Ten flips, seven heads. The hand derivation says p-hat = 0.7."""
    heads, n = 7, 10

    def log_likelihood(p):
        return heads * np.log(p) + (n - heads) * np.log(1 - p)

    grid = np.linspace(0.001, 0.999, 999)
    numeric = grid[np.argmax([log_likelihood(p) for p in grid])]
    analytic = heads / n

    assert round(analytic, 4) == 0.7
    assert abs(numeric - analytic) < 0.002, (numeric, analytic)

    # The likelihood at the MLE beats every alternative -- the definition.
    assert log_likelihood(0.7) > log_likelihood(0.5)
    assert log_likelihood(0.7) > log_likelihood(0.9)

    print(f"  {heads} heads in {n} flips")
    print(f"    analytic MLE  (calculus)   p-hat = {analytic:.4f}")
    print(f"    numeric  MLE  (grid search) p-hat = {numeric:.4f}")
    print(f"    log-likelihood at p=0.5 {log_likelihood(0.5):8.4f}")
    print(f"                   at p=0.7 {log_likelihood(0.7):8.4f}  <- maximum")
    print(f"                   at p=0.9 {log_likelihood(0.9):8.4f}")
    print("       the MLE is the SAMPLE PROPORTION, which is what intuition")
    print("       already suggested -- MLE gives it a principled derivation")


def why_take_logs():
    """Products of many small probabilities underflow. Sums of logs do not."""
    p = 0.5
    n = 2000

    product = p ** n
    log_sum = n * np.log(p)

    assert product == 0.0, "underflowed to EXACTLY zero"
    assert log_sum < -1000 and np.isfinite(log_sum)

    print(f"  0.5 ** {n} as a float  = {product}   <- underflowed to zero")
    print(f"  {n} * ln(0.5)          = {log_sum:.4f}   <- finite and usable")
    print("       a likelihood is a PRODUCT of n probabilities each below 1.")
    print("       Beyond a few hundred observations it is exactly 0.0 in")
    print("       floating point, and every candidate parameter ties at zero.")
    print("       The log is monotonic, so the maximum is unchanged")


def least_squares_is_the_mle_under_normal_errors():
    """unit-3.md's claim, verified: two different derivations, one answer."""
    X, y = study_xy()
    x = X.ravel()

    # (a) Ordinary least squares, in closed form.
    x_bar, y_bar = x.mean(), y.mean()
    ls_slope = ((x - x_bar) * (y - y_bar)).sum() / ((x - x_bar) ** 2).sum()
    ls_intercept = y_bar - ls_slope * x_bar

    # (b) Maximum likelihood, assuming y ~ Normal(b0 + b1 x, sigma).
    def negative_log_likelihood(theta):
        b0, b1, log_sigma = theta
        sigma = np.exp(log_sigma)
        residual = y - (b0 + b1 * x)
        return -np.sum(stats.norm.logpdf(residual, loc=0, scale=sigma))

    fit = optimize.minimize(negative_log_likelihood, x0=[0.0, 1.0, 0.0],
                            method="Nelder-Mead",
                            options={"xatol": 1e-10, "fatol": 1e-10,
                                     "maxiter": 20000, "maxfev": 20000})
    mle_intercept, mle_slope, log_sigma = fit.x

    assert fit.success, fit.message
    assert round(ls_slope, 4) == 4.3030
    assert round(ls_intercept, 4) == 43.0303
    assert abs(mle_slope - ls_slope) < 1e-3, (mle_slope, ls_slope)
    assert abs(mle_intercept - ls_intercept) < 1e-2, (mle_intercept, ls_intercept)

    # And the MLE of sigma is the ROOT MEAN squared residual (divide by n),
    # not the unbiased estimate (divide by n-2).
    residual = y - (ls_intercept + ls_slope * x)
    sigma_mle = np.sqrt((residual ** 2).mean())
    assert abs(np.exp(log_sigma) - sigma_mle) < 1e-2

    print("                     slope    intercept")
    print(f"    least squares  {ls_slope:9.4f}  {ls_intercept:11.4f}")
    print(f"    maximum likelihood {mle_slope:5.4f}  {mle_intercept:11.4f}")
    print(f"    agree to within {abs(mle_slope - ls_slope):.2e} and "
          f"{abs(mle_intercept - ls_intercept):.2e}")
    print(f"    MLE of sigma = RMS residual = {sigma_mle:.4f}")
    print("       TWO DIFFERENT DERIVATIONS, ONE ANSWER. Minimising squared")
    print("       error and maximising likelihood are the same thing when the")
    print("       errors are normal -- which is why least squares is not an")
    print("       arbitrary choice but a consequence of an assumption")


def mle_for_a_normal_mean_and_variance():
    """The other standard derivation: mu-hat is the sample mean."""
    rng = np.random.default_rng(0)
    sample = rng.normal(50, 8, 500)

    mu_hat = sample.mean()
    sigma_hat = np.sqrt(((sample - mu_hat) ** 2).mean())    # divides by n
    sigma_unbiased = sample.std(ddof=1)                     # divides by n-1

    assert abs(mu_hat - 50) < 1.0
    assert sigma_hat < sigma_unbiased, "the MLE of sigma is BIASED downward"
    ratio = sigma_unbiased / sigma_hat
    assert abs(ratio - np.sqrt(500 / 499)) < 1e-12

    print(f"  500 draws from Normal(50, 8):")
    print(f"    mu-hat    (MLE)          {mu_hat:.4f}   = the sample mean")
    print(f"    sigma-hat (MLE, /n)      {sigma_hat:.4f}")
    print(f"    sigma     (unbiased, /n-1) {sigma_unbiased:.4f}")
    print(f"    ratio = sqrt(n/(n-1)) = {ratio:.6f}")
    print("       the MLE of the mean is the sample mean, but the MLE of the")
    print("       VARIANCE divides by n and is biased low. Course 4 divided by")
    print("       n-1 for exactly this reason -- MLE is not automatically")
    print("       unbiased, and knowing that is worth a mark")


def main():
    print("Experiment 5 -- Maximum likelihood estimation")
    # Step 1: Estimate a coin's p by grid search and by calculus
    coin_mle_by_grid_and_by_calculus()
    # Step 2: See why logs are taken
    why_take_logs()
    # Step 3: Show least squares is the MLE under normal errors
    least_squares_is_the_mle_under_normal_errors()
    # Step 4: Estimate a normal mean and variance
    mle_for_a_normal_mean_and_variance()


if __name__ == "__main__":
    main()

5. Execution and Results

Experiment 5 -- Maximum likelihood estimation
  7 heads in 10 flips
    analytic MLE  (calculus)   p-hat = 0.7000
    numeric  MLE  (grid search) p-hat = 0.7000
    log-likelihood at p=0.5  -6.9315
                   at p=0.7  -6.1086  <- maximum
                   at p=0.9  -7.6453
       the MLE is the SAMPLE PROPORTION, which is what intuition
       already suggested -- MLE gives it a principled derivation
  0.5 ** 2000 as a float  = 0.0   <- underflowed to zero
  2000 * ln(0.5)          = -1386.2944   <- finite and usable
       a likelihood is a PRODUCT of n probabilities each below 1.
       Beyond a few hundred observations it is exactly 0.0 in
       floating point, and every candidate parameter ties at zero.
       The log is monotonic, so the maximum is unchanged
                     slope    intercept
    least squares     4.3030      43.0303
    maximum likelihood 4.3030      43.0303
    agree to within 6.66e-09 and 3.13e-08
    MLE of sigma = RMS residual = 0.8015
       TWO DIFFERENT DERIVATIONS, ONE ANSWER. Minimising squared
       error and maximising likelihood are the same thing when the
       errors are normal -- which is why least squares is not an
       arbitrary choice but a consequence of an assumption
  500 draws from Normal(50, 8):
    mu-hat    (MLE)          49.7849   = the sample mean
    sigma-hat (MLE, /n)      8.1088
    sigma     (unbiased, /n-1) 8.1169
    ratio = sqrt(n/(n-1)) = 1.001002
       the MLE of the mean is the sample mean, but the MLE of the
       VARIANCE divides by n and is biased low. Course 4 divided by
       n-1 for exactly this reason -- MLE is not automatically
       unbiased, and knowing that is worth a mark

Experiment 6 — Simple and multiple linear regression

1. Question

Fit simple and multiple linear regression models.

2. Aim

Reproduce Statistical Foundations' regression, then see extrapolation, multicollinearity and overfitting go wrong.

3. Steps

1. Fit the simple regression, and check it against Statistical Foundations.
2. Predict, and see the extrapolation trap.
3. Fit a multiple regression with collinear predictors.
4. Compare R² with adjusted R².
5. Raise the polynomial degree.
6. Compare Ridge with Lasso.

4. Programme

"""Experiment 6 — Simple and multiple linear regression.

Refits COURSE 4's worked example with scikit-learn and asserts it reproduces
the hand-computed slope 4.3030, intercept 43.0303 and R^2 0.9958. If the two
courses ever disagree, this test says so.

Then multiple regression, multicollinearity, and the polynomial degree /
conditioning result from unit-3.md section 3.5.
"""
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression, Lasso, Ridge
from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

from fixtures import RANDOM_STATE, STUDY, study_xy


def simple_linear_regression_matches_course_4():
    X, y = study_xy()
    model = LinearRegression().fit(X, y)
    pred = model.predict(X)

    slope = float(model.coef_[0])
    intercept = float(model.intercept_)

    assert round(slope, 4) == 4.3030, slope
    assert round(intercept, 4) == 43.0303, intercept
    assert round(model.score(X, y), 4) == 0.9958
    assert round(float(np.sqrt(mean_squared_error(y, pred))), 4) == 0.8015
    assert round(float(mean_absolute_error(y, pred)), 4) == 0.7394

    # Two checks Course 4 gave, verified rather than asserted in prose.
    assert abs((y - pred).sum()) < 1e-10, "least-squares residuals sum to zero"
    r = float(np.corrcoef(X.ravel(), y)[0, 1])
    assert round(r ** 2, 4) == round(model.score(X, y), 4), "R^2 = r^2"

    print(f"  slope     {slope:.4f}   (Course 4 computed 4.3030)")
    print(f"  intercept {intercept:.4f}  (Course 4 computed 43.0303)")
    print(f"  R^2       {model.score(X, y):.4f}   (Course 4 computed 0.9958)")
    print(f"  RMSE {np.sqrt(mean_squared_error(y, pred)):.4f}   "
          f"MAE {mean_absolute_error(y, pred):.4f}")
    print(f"  residuals sum to {(y - pred).sum():.1e}   r^2 = {r ** 2:.4f} = R^2")
    print("       scikit-learn reproduces Course 4's hand arithmetic exactly")
    return model


def prediction_and_the_extrapolation_trap(model):
    inside = float(model.predict([[7.5]])[0])
    outside = float(model.predict([[50]])[0])

    assert round(inside, 4) == 75.3030, inside
    assert round(outside, 4) == 258.1818, outside
    assert outside > 100, "a percentage score above 100 -- impossible"

    print(f"  x = 7.5 (inside the observed range 2-11) -> {inside:.4f} marks")
    print(f"  x = 50  (far outside it)                 -> {outside:.4f} marks")
    print("       258 marks out of 100. Nothing in the code warns you; the")
    print("       model has no evidence beyond the range it was fitted on")


def multiple_regression_and_multicollinearity():
    """A duplicate-ish predictor destabilises coefficients, not predictions."""
    rng = np.random.default_rng(RANDOM_STATE)
    n = 200
    x1 = rng.normal(0, 1, n)
    x2 = x1 * 0.98 + rng.normal(0, 0.05, n)      # almost a copy of x1
    x3 = rng.normal(0, 1, n)
    y = 3 * x1 + 2 * x3 + rng.normal(0, 0.5, n)

    clean = LinearRegression().fit(np.c_[x1, x3], y)
    collinear = LinearRegression().fit(np.c_[x1, x2, x3], y)

    corr = float(np.corrcoef(x1, x2)[0, 1])
    assert corr > 0.99, corr

    # Predictions barely change; the COEFFICIENTS do.
    r2_clean = clean.score(np.c_[x1, x3], y)
    r2_coll = collinear.score(np.c_[x1, x2, x3], y)
    assert abs(r2_clean - r2_coll) < 0.01, (r2_clean, r2_coll)
    assert abs(collinear.coef_[0] - 3.0) > abs(clean.coef_[0] - 3.0), \
        "the x1 coefficient moved away from its true value of 3"

    # Ridge is the standard remedy.
    ridge = make_pipeline(StandardScaler(), Ridge(alpha=1.0)).fit(
        np.c_[x1, x2, x3], y)

    print(f"  x1 and x2 correlate {corr:.4f}; the truth is y = 3*x1 + 2*x3")
    print(f"    without x2: coefficients {np.round(clean.coef_, 4)}  R^2 {r2_clean:.4f}")
    print(f"    with    x2: coefficients {np.round(collinear.coef_, 4)}  R^2 {r2_coll:.4f}")
    print(f"    ridge      : coefficients {np.round(ridge[-1].coef_, 4)}")
    print("       R^2 barely moved -- PREDICTION is fine. The coefficients split")
    print("       the shared effect arbitrarily between x1 and x2, so they can")
    print("       no longer be interpreted. Multicollinearity breaks explanation,")
    print("       not prediction, and Ridge is the standard remedy")


def r2_never_decreases_but_adjusted_r2_does():
    """Add PURE NOISE as predictors, over many random draws.

    One draw proves nothing here, because adjusted R^2 can rise by chance.
    The claim that is actually true is about the DISTRIBUTION: R^2 rises every
    single time, while adjusted R^2 rises only about as often as it falls and
    has no upward drift. That is what this measures.
    """
    n, trials, extra = 100, 300, 5

    def adjusted(r2, p):
        return 1 - (1 - r2) * (n - 1) / (n - p - 1)

    rose_r2 = rose_adj = 0
    deltas_r2, deltas_adj = [], []
    for seed in range(trials):
        rng = np.random.default_rng(seed)
        x = rng.normal(0, 1, n)
        y = 3 * x + rng.normal(0, 1, n)

        X1 = x.reshape(-1, 1)
        r2_one = LinearRegression().fit(X1, y).score(X1, y)

        Xk = np.c_[x, rng.normal(0, 1, (n, extra))]     # pure noise columns
        r2_many = LinearRegression().fit(Xk, y).score(Xk, y)

        rose_r2 += r2_many > r2_one
        rose_adj += adjusted(r2_many, 1 + extra) > adjusted(r2_one, 1)
        deltas_r2.append(r2_many - r2_one)
        deltas_adj.append(adjusted(r2_many, 1 + extra) - adjusted(r2_one, 1))

    assert rose_r2 == trials, (rose_r2, trials)
    assert rose_adj < trials // 2, rose_adj
    assert np.mean(deltas_r2) > 0.004
    assert abs(np.mean(deltas_adj)) < 0.001, np.mean(deltas_adj)
    assert min(deltas_r2) > 0, "R^2 never fell, not once"

    print(f"  {trials} random datasets, adding {extra} columns of PURE NOISE:")
    print(f"    R^2 rose in          {rose_r2:3d}/{trials}  "
          f"({rose_r2 / trials * 100:5.1f}%)   mean change {np.mean(deltas_r2):+.6f}")
    print(f"    adjusted R^2 rose in {rose_adj:3d}/{trials}  "
          f"({rose_adj / trials * 100:5.1f}%)   mean change {np.mean(deltas_adj):+.6f}")
    print(f"    smallest R^2 change seen: {min(deltas_r2):+.6f}  <- still positive")
    print("       R^2 rose EVERY time, because adding a column can only give")
    print("       least squares more freedom. Adjusted R^2 penalises the extra")
    print("       parameters, rises about as often as it falls, and has no")
    print("       upward drift -- which is why it can compare models of")
    print("       different sizes and R^2 cannot")


def polynomial_degree_and_conditioning():
    """unit-3.md 3.5: degree 9 interpolates 10 points -- if you standardise."""
    X, y = study_xy()

    raw, scaled, conditions = {}, {}, {}
    for d in (1, 2, 3, 5, 9):
        raw[d] = make_pipeline(PolynomialFeatures(d),
                               LinearRegression()).fit(X, y).score(X, y)
        scaled[d] = make_pipeline(StandardScaler(), PolynomialFeatures(d),
                                  LinearRegression()).fit(X, y).score(X, y)
        conditions[d] = float(np.linalg.cond(PolynomialFeatures(d).fit_transform(X)))

    assert round(raw[1], 6) == 0.995812
    assert round(scaled[9], 6) == 1.000000, "10 points, degree 9 -> exact fit"
    assert round(raw[9], 6) == 0.981546, "unscaled, it fails NUMERICALLY"
    assert raw[9] < raw[1], "and ends up worse than a straight line"
    assert conditions[9] > 1e13 and conditions[1] < 1e2

    print(f"    {'degree':>7} {'R2 unscaled':>13} {'R2 scaled':>11} {'condition no.':>15}")
    for d in (1, 2, 3, 5, 9):
        print(f"    {d:>7} {raw[d]:13.6f} {scaled[d]:11.6f} {conditions[d]:15.2e}")
    print("       standardised, degree 9 fits 10 points EXACTLY (R^2 = 1.0) --")
    print("       9 coefficients plus an intercept for 10 points. That is")
    print("       overfitting made visible.")
    print("       UNSTANDARDISED it scores 0.9815, worse than a straight line,")
    print("       because the condition number reaches 2.7e13 and the solve")
    print("       loses its precision. Standardise before PolynomialFeatures")


def ridge_and_lasso():
    """L2 shrinks; L1 shrinks to exactly zero and therefore selects."""
    rng = np.random.default_rng(RANDOM_STATE)
    n, p = 100, 12
    X = rng.normal(0, 1, (n, p))
    true = np.zeros(p)
    true[:3] = [4.0, -3.0, 2.0]                  # only 3 of 12 matter
    y = X @ true + rng.normal(0, 0.5, n)

    ols = LinearRegression().fit(X, y)
    ridge = Ridge(alpha=10.0).fit(X, y)
    lasso = Lasso(alpha=0.3).fit(X, y)

    ols_zero = int(np.sum(np.abs(ols.coef_) < 1e-8))
    ridge_zero = int(np.sum(np.abs(ridge.coef_) < 1e-8))
    lasso_zero = int(np.sum(np.abs(lasso.coef_) < 1e-8))

    assert ols_zero == 0 and ridge_zero == 0
    assert lasso_zero >= 7, lasso_zero
    assert all(abs(lasso.coef_[i]) > 1.0 for i in range(3)), \
        "Lasso kept the three real predictors"

    print(f"  12 predictors, only 3 of which matter:")
    print(f"    OLS   : {ols_zero:2d} coefficients exactly zero")
    print(f"    Ridge : {ridge_zero:2d} coefficients exactly zero  (shrunk, not removed)")
    print(f"    Lasso : {lasso_zero:2d} coefficients exactly zero  <- it SELECTED")
    print(f"    Lasso's first three: {np.round(lasso.coef_[:3], 4)}  (true 4, -3, 2)")
    print("       that is the difference: L2 shrinks toward zero, L1 shrinks TO")
    print("       zero, which makes Lasso an embedded feature-selection method")


def main():
    print("Experiment 6 -- Simple and multiple linear regression")
    # Step 1: Fit the simple regression, and check it against Statistical Foundations
    model = simple_linear_regression_matches_course_4()
    # Step 2: Predict, and see the extrapolation trap
    prediction_and_the_extrapolation_trap(model)
    # Step 3: Fit a multiple regression with collinear predictors
    multiple_regression_and_multicollinearity()
    # Step 4: Compare R² with adjusted R²
    r2_never_decreases_but_adjusted_r2_does()
    # Step 5: Raise the polynomial degree
    print("  polynomial degree, with and without standardising:")
    polynomial_degree_and_conditioning()
    # Step 6: Compare Ridge with Lasso
    ridge_and_lasso()


if __name__ == "__main__":
    main()