Hands-on lab — Day 148: Linear Regression

"""One predictor, one line: what fitting it actually buys you, measured.

Simple linear regression is a claim about a straight-line relationship
between one predictor and one target, fitted so the squared vertical
distances from the line are as small as possible (Day 149 owns *why*
squared error; this module just fits it, with scikit-learn's
`LinearRegression`). Two facts about that fit are exact on any data,
forever: it passes through the point of means, and its residuals sum to
zero when an intercept is fitted. Everything else here is a measurement of
what the line gets right and what it silently hides.

The module is organised in the order the lesson uses it:

1. The line itself, fitted to real data in real units (age-adjusted
   diabetes progression against BMI, raw scale).
2. Recovering a slope you know to be true, and watching the error shrink
   with more rows.
3. Two ways a fit can look fine on a scatterplot and an R-squared, and
   still be wrong: curvature and heteroscedasticity, both visible only in
   the residuals.
4. One point that moves the whole line, and the number that names why.
5. What forcing the intercept to zero costs, measured against the same
   data fitted honestly.

Everything here is deterministic given a seed. Nothing downloads: the
diabetes dataset is bundled with scikit-learn and the rest is generated on
the spot from `numpy.random.default_rng`.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression


# --------------------------------------------------------------------------
# 1. The line itself: BMI against disease progression, in raw units
# --------------------------------------------------------------------------


def load_bmi_and_target():
    """BMI and the disease-progression target, in their original units.

    `load_diabetes(scaled=False)` is the only way to get raw units out of
    this dataset -- the default returns every column mean-centred and
    scaled to unit norm, which makes a coefficient uninterpretable as
    "one more unit of X". Column 2 is BMI; see `load_diabetes().feature_names`.
    """
    data = load_diabetes(scaled=False)
    bmi = data.data[:, 2].reshape(-1, 1)
    y = data.target
    return bmi, y


def fit_line(x, y, fit_intercept: bool = True) -> LinearRegression:
    """Fit scikit-learn's `LinearRegression`. This module never derives it."""
    return LinearRegression(fit_intercept=fit_intercept).fit(x, y)


def slope_standard_error(x, residuals) -> float:
    """How much the fitted slope would wobble under a fresh sample.

    `SE(b1) = s / sqrt(sum((x - xbar)^2))`, where `s^2` is the residual
    variance with `n - 2` degrees of freedom spent on the slope and the
    intercept. This is the same standard-error arithmetic Days 117-118 and
    Day 144 used for a sampling proportion, applied here to a slope.
    """
    x = np.asarray(x, dtype=float).flatten()
    n = len(x)
    dof = n - 2
    s2 = float(np.sum(np.asarray(residuals) ** 2)) / dof
    sxx = float(np.sum((x - x.mean()) ** 2))
    return float(np.sqrt(s2 / sxx))


def confidence_interval(estimate: float, standard_error: float, z: float = 1.96) -> tuple:
    """A normal-approximation interval, rounded for reporting."""
    return (
        round(estimate - z * standard_error, 4),
        round(estimate + z * standard_error, 4),
    )


def passes_through_the_means(model: LinearRegression, x, y) -> tuple:
    """Confirms the fitted line predicts `mean(y)` exactly at `mean(x)`.

    Not approximately -- exactly, up to floating point. It is a property
    of the least-squares normal equations (Day 149's territory), not a
    coincidence of this dataset.
    """
    x = np.asarray(x, dtype=float)
    mean_x = x.mean(axis=0).reshape(1, -1)
    predicted_at_mean = float(model.predict(mean_x)[0])
    mean_y = float(np.asarray(y).mean())
    return predicted_at_mean, mean_y, predicted_at_mean - mean_y


def residual_sum(residuals) -> float:
    """The sum of the residuals, which is exactly zero when an intercept
    is fitted -- another consequence of the normal equations, not a
    measurement that happened to come out that way."""
    return float(np.sum(residuals))


# --------------------------------------------------------------------------
# 2. Recovering a slope you know to be true
# --------------------------------------------------------------------------


def make_known_line(n: int, seed: int, true_slope: float = 5.0, true_intercept: float = 10.0, noise_sd: float = 8.0):
    """One predictor, a known slope and intercept, and Gaussian noise."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 20, size=n).reshape(-1, 1)
    y = true_slope * x.flatten() + true_intercept + rng.normal(0, noise_sd, size=n)
    return x, y


def slope_recovery_error(n_values, replications: int = 200, true_slope: float = 5.0):
    """Mean absolute error of the fitted slope, at each sample size.

    Averaged over `replications` independently drawn datasets per `n`, for
    the same reason Days 117-118 and Day 144 always averaged rather than
    quoting one draw: a single fit's error is an anecdote.
    """
    rows = []
    for n in n_values:
        errors = [
            abs(float(fit_line(*make_known_line(n, seed, true_slope=true_slope)).coef_[0]) - true_slope)
            for seed in range(replications)
        ]
        rows.append((n, round(float(np.mean(errors)), 4)))
    return rows


# --------------------------------------------------------------------------
# 3a. Curvature: a fit that looks fine and is not
# --------------------------------------------------------------------------


def curved_dataset(n: int = 300, seed: int = 1):
    """A quadratic relationship, so a straight line is the wrong model."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 10, size=n)
    y = 2.0 + 0.5 * x**2 + rng.normal(0, 5.0, size=n)
    return x.reshape(-1, 1), y


def binned_residual_means(x, residuals, bins: int = 5):
    """Mean residual within each of `bins` equal-count groups, sorted by x.

    A straight line fitted to a curve leaves residuals that are positive
    at both ends and negative in the middle (or the reverse) -- a shape no
    single number like R-squared exposes, but a residual plot shows at a
    glance.
    """
    x = np.asarray(x, dtype=float).flatten()
    residuals = np.asarray(residuals, dtype=float)
    order = np.argsort(x)
    x_sorted, resid_sorted = x[order], residuals[order]
    groups = np.array_split(np.arange(len(x_sorted)), bins)
    return [
        (round(float(x_sorted[g].mean()), 2), round(float(resid_sorted[g].mean()), 4))
        for g in groups
    ]


def quadratic_fit_r_squared(x, residuals) -> float:
    """How much of the residuals' own variance a quadratic curve explains.

    Near zero when the residuals are patternless noise; large when the
    line missed real curvature. Fitted with `numpy.polyfit`, not
    scikit-learn -- this is a diagnostic on the residuals, not a second
    model of the data.
    """
    x = np.asarray(x, dtype=float).flatten()
    residuals = np.asarray(residuals, dtype=float)
    coeffs = np.polyfit(x, residuals, 2)
    predicted = np.polyval(coeffs, x)
    ss_res = float(np.sum((residuals - predicted) ** 2))
    ss_tot = float(np.sum((residuals - residuals.mean()) ** 2))
    return float(1.0 - ss_res / ss_tot)


# --------------------------------------------------------------------------
# 3b. Heteroscedasticity: error that grows with x
# --------------------------------------------------------------------------


def heteroscedastic_dataset(n: int = 400, seed: int = 2):
    """A line whose noise gets wider as x grows -- the fit stays roughly
    unbiased; only its residual spread reveals what is wrong."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(1, 20, size=n)
    noise_sd = 0.8 * x
    y = 3.0 + 2.0 * x + rng.normal(0, 1, size=n) * noise_sd
    return x.reshape(-1, 1), y


def residual_spread_by_half(x, residuals) -> tuple:
    """Residual standard deviation in the low-x half against the high-x half."""
    x = np.asarray(x, dtype=float).flatten()
    residuals = np.asarray(residuals, dtype=float)
    median = float(np.median(x))
    low = residuals[x < median]
    high = residuals[x >= median]
    return float(low.std()), float(high.std())


# --------------------------------------------------------------------------
# 4. One point that moves the line
# --------------------------------------------------------------------------


def leverage_dataset(n: int = 40, seed: int = 3):
    """A clean, ordinary linear relationship, forty points."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 10, size=n)
    y = 2.0 + 1.5 * x + rng.normal(0, 1.5, size=n)
    return x, y


def add_point(x, y, x_new: float, y_new: float):
    """Append one point, returning fresh arrays (never mutates the inputs)."""
    return np.append(np.asarray(x, dtype=float), x_new), np.append(np.asarray(y, dtype=float), y_new)


def leverage_of_point(x, x_target: float) -> float:
    """The hat-matrix leverage of a point at `x_target`, given the full x array.

    `h = 1/n + (x_target - xbar)^2 / sum((x - xbar)^2)` -- how much that
    single point's own y-value can pull the fitted line toward itself,
    independent of what its y-value actually is.
    """
    x = np.asarray(x, dtype=float)
    n = len(x)
    xbar = x.mean()
    sxx = float(np.sum((x - xbar) ** 2))
    return float(1.0 / n + (x_target - xbar) ** 2 / sxx)


def mean_leverage_excluding(x, x_target: float) -> float:
    """Average leverage of every point except the one at `x_target`."""
    x = np.asarray(x, dtype=float)
    return float(np.mean([leverage_of_point(x, xi) for xi in x if xi != x_target]))


# --------------------------------------------------------------------------
# 5. fit_intercept=False, and what it costs
# --------------------------------------------------------------------------


def intercept_dataset(n: int = 200, seed: int = 4, true_intercept: float = 25.0, true_slope: float = 3.0, noise_sd: float = 6.0):
    """A line whose x-values never go near zero, so a forced-zero intercept
    is a real misspecification rather than a harmless simplification."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(5, 25, size=n)
    y = true_intercept + true_slope * x + rng.normal(0, noise_sd, size=n)
    return x.reshape(-1, 1), y


def rmse(y_true, y_pred) -> float:
    y_true = np.asarray(y_true, dtype=float)
    y_pred = np.asarray(y_pred, dtype=float)
    return float(np.sqrt(np.mean((y_true - y_pred) ** 2)))


# --------------------------------------------------------------------------
# 6. Skewness, for a rough normality read
# --------------------------------------------------------------------------


def skewness(values) -> float:
    """The third standardised moment: 0 for a symmetric distribution.

    A rough diagnostic only -- not a formal normality test, and this
    module does not claim to be one. `abs(skewness) < 0.5` is a common
    rule of thumb for "not alarming".
    """
    values = np.asarray(values, dtype=float)
    mean = values.mean()
    sd = values.std()
    return float(np.mean(((values - mean) / sd) ** 3))

#!/usr/bin/env python3
"""Print every measured pair in this lab as one table.

The harness compares this output byte for byte against
expected-output/measured-values.txt, so the report is not a convenience:
it is how the lab notices that a number in the lesson has gone stale.
"""

import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))

import numpy as np  # noqa: E402

import regression_lib as r  # noqa: E402


def rule(title: str) -> None:
    print()
    print(title)
    print("-" * len(title))


def main() -> None:
    print("Day 148 -- linear regression, measured")
    print("=" * 39)

    rule("1. The line itself: BMI against disease progression, raw units")
    bmi, y = r.load_bmi_and_target()
    model = r.fit_line(bmi, y)
    residuals = y - model.predict(bmi)
    se = r.slope_standard_error(bmi, residuals)
    ci = r.confidence_interval(float(model.coef_[0]), se)
    predicted_at_mean, mean_y, diff = r.passes_through_the_means(model, bmi, y)
    print(f"  n = {len(y)}, BMI range {bmi.min():.1f}-{bmi.max():.1f}, target range {y.min():.1f}-{y.max():.1f}")
    print(f"  slope     : {model.coef_[0]:.4f}  (points of progression per unit of BMI)")
    print(f"  intercept : {model.intercept_:.4f}")
    print(f"  R-squared : {model.score(bmi, y):.4f}")
    print(f"  slope SE  : {se:.4f}   95% CI [{ci[0]}, {ci[1]}]   t = {model.coef_[0] / se:.2f}")
    print(f"  predicted at mean BMI : {predicted_at_mean:.4f}   mean of y : {mean_y:.4f}   diff : {diff:.2e}")
    print(f"  sum of residuals      : {r.residual_sum(residuals):.2e}")

    rule("2. Recovering a slope you know to be true (true slope = 5.0)")
    print("       n   mean abs error")
    for n, err in r.slope_recovery_error([20, 50, 200, 1000, 5000]):
        print(f"  {n:6d}   {err:.4f}")
    by_n = dict(r.slope_recovery_error([20, 50, 200, 1000, 5000]))
    print(f"  ratio, n=20 to n=200  : {by_n[200] / by_n[20]:.4f}   (predicted ~ 1/sqrt(10) = {1 / np.sqrt(10):.4f})")
    print(f"  ratio, n=20 to n=5000 : {by_n[5000] / by_n[20]:.4f}   (predicted ~ 1/sqrt(250) = {1 / np.sqrt(250):.4f})")

    rule("3. Curvature: a fit that looks fine and is not")
    xc, yc = r.curved_dataset()
    model_c = r.fit_line(xc, yc)
    residuals_c = yc - model_c.predict(xc)
    print(f"  R-squared of the line : {model_c.score(xc, yc):.4f}   (looks respectable)")
    print("  mean residual by bin of x (positive, negative, positive: the missed curve)")
    for mean_x, mean_resid in r.binned_residual_means(xc, residuals_c, bins=5):
        print(f"    x ~ {mean_x:5.2f}   mean residual {mean_resid:+.4f}")
    quad_r2 = r.quadratic_fit_r_squared(xc, residuals_c)
    corr = float(np.corrcoef(residuals_c, xc.flatten() ** 2)[0, 1])
    print(f"  quadratic fit to the residuals, R-squared : {quad_r2:.4f}")
    print(f"  correlation of residuals with x^2         : {corr:.4f}")

    rule("4. Heteroscedasticity: error that grows with x")
    xh, yh = r.heteroscedastic_dataset()
    model_h = r.fit_line(xh, yh)
    residuals_h = yh - model_h.predict(xh)
    low_sd, high_sd = r.residual_spread_by_half(xh, residuals_h)
    print(f"  R-squared of the line : {model_h.score(xh, yh):.4f}   (also looks fine)")
    print(f"  residual sd, low half of x  : {low_sd:.4f}")
    print(f"  residual sd, high half of x : {high_sd:.4f}")
    print(f"  ratio                       : {high_sd / low_sd:.4f}")

    rule("5. One point that moves the line")
    xl, yl = r.leverage_dataset()
    model_without = r.fit_line(xl.reshape(-1, 1), yl)
    xl_with, yl_with = r.add_point(xl, yl, x_new=40.0, y_new=5.0)
    model_with = r.fit_line(xl_with.reshape(-1, 1), yl_with)
    leverage_new = r.leverage_of_point(xl_with, 40.0)
    typical = r.mean_leverage_excluding(xl_with, 40.0)
    print(f"  slope, 40 ordinary points        : {model_without.coef_[0]:.4f}")
    print(f"  slope, plus one leverage point    : {model_with.coef_[0]:.4f}")
    print(f"  change                            : {model_with.coef_[0] - model_without.coef_[0]:+.4f}")
    print(f"  leverage of the added point        : {leverage_new:.4f}")
    print(f"  mean leverage of the other 40      : {typical:.4f}")
    print(f"  ratio                              : {leverage_new / typical:.2f}")

    rule("6. fit_intercept=False, and what it costs")
    xi, yi = r.intercept_dataset()
    model_yes = r.fit_line(xi, yi, fit_intercept=True)
    model_no = r.fit_line(xi, yi, fit_intercept=False)
    rmse_yes = r.rmse(yi, model_yes.predict(xi))
    rmse_no = r.rmse(yi, model_no.predict(xi))
    print(f"  true intercept = 25.0, true slope = 3.0, x never near zero")
    print(f"  fit_intercept=True  : slope {model_yes.coef_[0]:.4f}  intercept {model_yes.intercept_:.4f}  RMSE {rmse_yes:.4f}")
    print(f"  fit_intercept=False : slope {model_no.coef_[0]:.4f}  intercept {model_no.intercept_:.4f}  RMSE {rmse_no:.4f}")
    print(f"  RMSE ratio (false/true) : {rmse_no / rmse_yes:.4f}")

    rule("7. Telling curvature apart from noise")
    quad_r2_bmi = r.quadratic_fit_r_squared(bmi, residuals)
    skew = r.skewness(residuals)
    print(f"  quadratic fit to the BMI model's own residuals, R-squared : {quad_r2_bmi:.4f}")
    print(f"  (contrast with section 3's 0.3558 on data with real curvature)")
    print(f"  skewness of the BMI model's residuals : {skew:.4f}")


if __name__ == "__main__":
    main()

"""The reference solutions: what a simple linear regression actually gets
you, and what it hides.

Every number here was captured from a real run of this file on the
authoring machine. If a number changes, the claim in the lesson is wrong
and one of the two must be fixed.
"""

import numpy as np

import regression_lib as r


# --- 1. The line, fitted to real data in real units ----------------------


def test_01_bmi_slope_and_intercept_in_raw_units():
    bmi, y = r.load_bmi_and_target()
    model = r.fit_line(bmi, y)
    assert round(float(model.coef_[0]), 4) == 10.2331
    assert round(float(model.intercept_), 4) == -117.7734
    assert round(float(model.score(bmi, y)), 4) == 0.3439
    # In words: each additional unit of BMI is associated with about ten
    # more points of one-year disease progression, on this population.


def test_01b_slope_standard_error_and_confidence_interval():
    bmi, y = r.load_bmi_and_target()
    model = r.fit_line(bmi, y)
    residuals = y - model.predict(bmi)
    se = r.slope_standard_error(bmi, residuals)
    assert round(se, 4) == 0.6738
    ci = r.confidence_interval(float(model.coef_[0]), se)
    assert ci == (8.9125, 11.5538)
    # The slope is about fifteen standard errors from zero -- not a
    # borderline effect.
    assert round(float(model.coef_[0]) / se, 2) == 15.19


def test_01c_the_line_passes_through_the_means_exactly():
    bmi, y = r.load_bmi_and_target()
    model = r.fit_line(bmi, y)
    residuals = y - model.predict(bmi)
    predicted_at_mean, mean_y, diff = r.passes_through_the_means(model, bmi, y)
    assert abs(diff) < 1e-8
    assert round(predicted_at_mean, 4) == round(mean_y, 4) == 152.1335
    # And the residuals sum to (essentially) zero -- not approximately,
    # to within floating point.
    assert abs(r.residual_sum(residuals)) < 1e-6


# --- 2. Recovering a slope you know to be true ----------------------------


def test_02_the_estimate_gets_closer_to_the_truth_as_n_grows():
    rows = r.slope_recovery_error([20, 50, 200, 1000, 5000])
    assert rows == [
        (20, 0.2315),
        (50, 0.1556),
        (200, 0.078),
        (1000, 0.0357),
        (5000, 0.0159),
    ]
    errors = [e for _n, e in rows]
    # Strictly shrinking as n grows.
    assert all(a > b for a, b in zip(errors, errors[1:]))


def test_02b_the_error_shrinks_roughly_like_one_over_root_n():
    rows = r.slope_recovery_error([20, 50, 200, 1000, 5000])
    by_n = dict(rows)
    # 200 has 10x the rows of 20; one-over-root-n predicts a shrink of
    # about 1/sqrt(10) = 0.316. The measured ratio is close, not exact --
    # this is a noisy quantity averaged over 200 replications, not a
    # formula being evaluated.
    ratio_20_to_200 = by_n[200] / by_n[20]
    assert 0.25 < ratio_20_to_200 < 0.40
    # 5000 has 250x the rows of 20; predicted shrink about 1/sqrt(250) = 0.063.
    ratio_20_to_5000 = by_n[5000] / by_n[20]
    assert 0.04 < ratio_20_to_5000 < 0.09


# --- 3. Curvature: a fit that looks fine and is not -----------------------


def test_03_residuals_reveal_the_curve_the_scatter_hides():
    x, y = r.curved_dataset()
    model = r.fit_line(x, y)
    residuals = y - model.predict(x)
    # The line's own R-squared looks respectable...
    assert round(float(model.score(x, y)), 4) == 0.852
    # ...but the residuals, binned by x, trace the missed curve: positive
    # at both ends, negative in the middle.
    bins = r.binned_residual_means(x, residuals, bins=5)
    means = [m for _x, m in bins]
    assert means[0] > 0
    assert means[2] < 0
    assert means[-1] > 0


def test_03b_quantifying_the_curvature_with_a_quadratic_fit_to_the_residuals():
    x, y = r.curved_dataset()
    model = r.fit_line(x, y)
    residuals = y - model.predict(x)
    # A quadratic curve explains over a third of the *residuals'* own
    # variance -- there is a whole model's worth of missed structure left
    # in what the line called "error".
    quad_r2 = r.quadratic_fit_r_squared(x, residuals)
    assert round(quad_r2, 4) == 0.3558
    corr = float(np.corrcoef(residuals, x.flatten() ** 2)[0, 1])
    assert round(corr, 4) == 0.148


# --- 4. Heteroscedasticity: error that grows with x ------------------------


def test_04_heteroscedasticity_fans_the_residuals_while_the_fit_looks_fine():
    x, y = r.heteroscedastic_dataset()
    model = r.fit_line(x, y)
    residuals = y - model.predict(x)
    # A perfectly ordinary-looking R-squared...
    assert round(float(model.score(x, y)), 4) == 0.5723
    # ...over noise whose spread more than doubles from the low half of x
    # to the high half. Nothing in the scatterplot or the R-squared says
    # so; only the residual plot does.
    low_sd, high_sd = r.residual_spread_by_half(x, residuals)
    assert round(low_sd, 4) == 4.7427
    assert round(high_sd, 4) == 12.0684
    assert round(high_sd / low_sd, 4) == 2.5446


# --- 5. One point that moves the line --------------------------------------


def test_05_one_high_leverage_point_moves_the_line():
    x, y = r.leverage_dataset()
    model_without = r.fit_line(x.reshape(-1, 1), y)
    x_with, y_with = r.add_point(x, y, x_new=40.0, y_new=5.0)
    model_with = r.fit_line(x_with.reshape(-1, 1), y_with)
    slope_without = float(model_without.coef_[0])
    slope_with = float(model_with.coef_[0])
    assert round(slope_without, 4) == 1.5196
    assert round(slope_with, 4) == 0.2138
    # One row, out of forty-one, cuts the slope by more than eighty-five
    # percent.
    assert round(slope_with - slope_without, 4) == -1.3059


def test_05b_the_leverage_value_names_the_mechanism():
    x, _y = r.leverage_dataset()
    x_with, _y_with = r.add_point(x, _y, x_new=40.0, y_new=5.0)
    leverage_new = r.leverage_of_point(x_with, 40.0)
    typical = r.mean_leverage_excluding(x_with, 40.0)
    assert round(leverage_new, 4) == 0.8048
    assert round(typical, 4) == 0.0299
    # Nearly twenty-seven times the pull of an ordinary point, computed
    # from its x-value alone -- before its y-value is even considered.
    assert round(leverage_new / typical, 2) == 26.94


# --- 6. fit_intercept=False, and what it costs ------------------------------


def test_06_forcing_the_intercept_to_zero_costs_you():
    x, y = r.intercept_dataset()
    model_yes = r.fit_line(x, y, fit_intercept=True)
    model_no = r.fit_line(x, y, fit_intercept=False)
    rmse_yes = r.rmse(y, model_yes.predict(x))
    rmse_no = r.rmse(y, model_no.predict(x))
    assert round(rmse_yes, 4) == 6.1401
    assert round(rmse_no, 4) == 9.7878
    # Nearly sixty percent worse, on data whose x-values never go near
    # zero -- the true intercept was 25.0, and forcing it to 0 makes the
    # slope absorb the difference instead.
    assert round(rmse_no / rmse_yes, 4) == 1.5941
    assert model_no.intercept_ == 0.0
    assert round(float(model_no.coef_[0]), 4) == 4.4232


# --- 7. Telling curvature apart from noise ----------------------------------


def test_07_the_bmi_models_residuals_show_no_such_curvature():
    bmi, y = r.load_bmi_and_target()
    model = r.fit_line(bmi, y)
    residuals = y - model.predict(bmi)
    # Contrast with test 03b's 0.3558: on real data with no known missed
    # curvature, a quadratic explains essentially none of the residuals'
    # variance.
    quad_r2 = r.quadratic_fit_r_squared(bmi, residuals)
    assert round(quad_r2, 4) == 0.0002
    # And the residuals are only mildly asymmetric -- not a formal test,
    # but nothing that should alarm you.
    skew = r.skewness(residuals)
    assert round(skew, 4) == 0.156
    assert abs(skew) < 0.5

"""Machinery checks: the helpers behave, before any claim is made.

These four tests are solved in both `starter/` and `examples/`. They exist
so that a broken helper reports itself as a broken helper rather than as a
surprising scientific result.
"""

import numpy as np

import regression_lib as r


def test_load_bmi_and_target_has_the_right_shape_and_raw_units():
    bmi, y = r.load_bmi_and_target()
    assert bmi.shape == (442, 1)
    assert y.shape == (442,)
    # Raw units, not the mean-centred, unit-norm-scaled default.
    assert bmi.min() >= 18.0 and bmi.max() <= 42.2
    assert y.min() >= 25.0 and y.max() <= 346.0


def test_fit_line_respects_fit_intercept():
    x = np.array([[1.0], [2.0], [3.0], [4.0]])
    y = np.array([3.0, 5.0, 7.0, 9.0])
    with_intercept = r.fit_line(x, y, fit_intercept=True)
    without_intercept = r.fit_line(x, y, fit_intercept=False)
    assert round(float(with_intercept.coef_[0]), 4) == 2.0
    assert round(float(with_intercept.intercept_), 4) == 1.0
    assert without_intercept.intercept_ == 0.0


def test_leverage_of_a_central_point_is_smaller_than_an_extreme_one():
    x = np.linspace(0, 10, 11)
    central = r.leverage_of_point(x, 5.0)
    extreme = r.leverage_of_point(x, 40.0)
    assert extreme > central
    # Leverage is never below 1/n for any point, including the mean.
    assert r.leverage_of_point(x, float(x.mean())) >= 1.0 / len(x) - 1e-9


def test_skewness_is_zero_for_a_symmetric_sample():
    rng = np.random.default_rng(0)
    symmetric = rng.normal(size=5000)
    assert abs(r.skewness(symmetric)) < 0.1
    lopsided = np.concatenate([np.zeros(950), np.full(50, 20.0)])
    assert r.skewness(lopsided) > 1.0

lesson_id: D148
day: 148
kind: guided-build
languages:
  - python
  - bash
setup_commands:
  - cd labs/sections/machine-learning/day-148-linear-regression
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - >-
    .venv/bin/python3 -c "import numpy, sklearn; print(numpy.__version__,
    sklearn.__version__)"
run_commands:
  - .venv/bin/pytest examples -q
  - .venv/bin/pytest starter -q
  - .venv/bin/python3 examples/report_measurements.py
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - >-
    find . -path ./.venv -prune -o -type d -name '__pycache__' -print -exec rm -rf -- {}
    +
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 55
last_executed: '2026-08-27'
executed_on: >-
  macOS 26.5.2 (Apple Silicon, arm64, CPU only -- no GPU is needed or used), Python
  3.14.0, numpy 2.5.2, scikit-learn 1.9.0, pytest 9.1.1, bash 3.2.57 -- bash
  tests/run_tests.sh -> 14 checks, 0 failure(s), exit 0. pytest examples -q -> 16 passed.
  pytest starter -q -> 4 passed, 12 skipped (the four machinery checks in
  test_regression_lib.py are solved in both directories; the twelve exercise stubs in
  starter/test_regression_claims.py are untouched). Everything ran through a real
  lab-local .venv created by the documented setup commands; scikit-learn pulled in scipy
  1.18.1, joblib 1.5.3 and threadpoolctl 3.6.0 as its own dependencies, none of which
  this lab imports directly. The lab is fully offline after the pip install -- the
  diabetes dataset ships inside the installed scikit-learn package and is not
  downloaded, every other dataset is generated on the spot from a seeded
  numpy.random.default_rng, and harness check 9 confirms no URL appears anywhere in
  starter/ or examples/ source. Section 7 of the harness copies examples/ into a
  mktemp -d scratch directory, confirms 16 passed, rewrites `assert round(rmse_no /
  rmse_yes, 4) == 1.5941` to 1.0, confirms a non-zero exit naming the failing test, and
  removes the scratch directory. Separately, by hand,
  `assert round(leverage_new / typical, 2) == 26.94`
  was changed to 1.0 in examples/test_regression_claims.py and the whole harness re-run:
  it reported 14 checks, 2 failure(s) and exited 1 (both the pytest run and the
  pytest-free direct reproduction in section 2 caught it); the file was restored and the
  harness returned to 14 checks, 0 failure(s), exit 0. MEASURED PAIRS, all captured
  verbatim in expected-output/measured-values.txt. (1) THE LINE ITSELF: BMI against
  one-year diabetes progression, in raw units (load_diabetes(scaled=False)). Slope
  10.2331, intercept -117.7734, R-squared 0.3439, slope standard error 0.6738, 95
  percent CI [8.9125, 11.5538], t = 15.19. The line passes through (mean(bmi), mean(y))
  to within 1e-8 and its 442 residuals sum to -1.67e-11. (2) RECOVERING A KNOWN SLOPE:
  with a true slope of 5.0, mean absolute error of the fitted slope over 200
  replications is 0.2315, 0.1556, 0.078, 0.0357 and 0.0159 at n = 20, 50, 200, 1000 and
  5000 -- strictly shrinking, and roughly tracking one-over-root-n (ratio 0.3369 from
  n=20 to n=200 against a predicted 0.3162; ratio 0.0687 from n=20 to n=5000 against a
  predicted 0.0632). (3) CURVATURE: a line fitted to quadratic data (curved_dataset,
  seed=1) scores R-squared 0.8520 -- a good-looking number -- while its residuals,
  binned by x into five groups, read +4.2216, -1.9594, -3.6829, -2.3803, +3.8010,
  tracing the missed curve exactly. A quadratic fitted to those residuals themselves
  explains R-squared 0.3558 of their variance; the correlation of the residuals with x^2
  is 0.1480. (4) HETEROSCEDASTICITY: a line fitted to data whose noise grows with x
  (heteroscedastic_dataset, seed=2) scores R-squared 0.5723 -- also unremarkable -- while
  the residual standard deviation is 4.7427 in the low half of x and 12.0684 in the high
  half, a ratio of 2.5446. (5) LEVERAGE: forty ordinary points give a slope of 1.5196;
  adding one point at x=40.0, y=5.0 drops it to 0.2138, a change of -1.3059. That added
  point's leverage is 0.8048 against a mean of 0.0299 for the other forty, a ratio of
  26.94, computed from its x-value alone. (6) FIT_INTERCEPT=FALSE: on data with a true
  intercept of 25.0 and true slope of 3.0 (x always between 5 and 25), RMSE is 6.1401
  with an intercept fitted and 9.7878 without one, a ratio of 1.5941; the intercept-free
  slope is biased to 4.4232 to compensate. (7) TELLING CURVATURE FROM NOISE: the same
  quadratic-fit-to-residuals diagnostic applied to the real BMI model's residuals gives
  R-squared 0.0002 -- essentially zero, in contrast to section 3's 0.3558 -- and a
  skewness of 0.156, mildly asymmetric and not alarming. TWO HONESTY CALLS. FIRST: the
  slope-recovery table in exercise 2 is an average over 200 replications, not a single
  fit, and its ratio to the one-over-root-n prediction is reported as a band (0.25 to
  0.40, and 0.04 to 0.09) rather than asserted exactly, because a sampled ratio at this
  replication count moves by several hundredths between runs at different seeds --
  Days 117-118 and Day 144 both made the same call for the same reason. SECOND: the
  0.8520 R-squared on the deliberately curved dataset is reported honestly as
  good-looking, not as an obviously bad fit -- the entire point of exercise 3 is that a
  respectable R-squared does not rule out a missed curve, and softening that finding to
  make the lesson tidier would have undercut it. Harness check 8 re-confirms slope
  recovery at a different true slope and replication count, curvature's binned shape at
  three dataset seeds, heteroscedasticity's fan at three dataset seeds, and the leverage
  point's effect at three base-dataset seeds, so no directional claim rests on the one
  seed quoted in the lesson.

# Requirements

`requirements.txt` pins the three packages this lab imports directly, at
the exact versions the captured output in `expected-output/` was produced
with:

```
numpy==2.5.2
scikit-learn==1.9.0
pytest==9.1.1
```

Installing scikit-learn also pulls in scipy (1.18.1), joblib (1.5.3) and
threadpoolctl (3.6.0) as its own dependencies during capture. This lab
imports none of them directly and does not pin them.

## Why the versions are pinned exactly

Most of this lab is deterministic given a seed, but a few figures are
averages over many seeded draws from `numpy.random.default_rng` (the
slope-recovery table in exercise 2), and NumPy's own documentation states
that `Generator` makes no promise of stream compatibility between
versions. A different NumPy can legitimately produce a different stream
from the same seed.

What does not depend on the pins: the BMI model itself (`sklearn.datasets.
load_diabetes` ships a fixed array, not a random draw, so its slope,
intercept, R-squared and standard error are exact on any working install
of these three packages); the two structural facts in exercise 1c (a
least-squares line passes through the point of means, and its residuals
sum to zero); the direction of every result — the error shrinks as n
grows, curvature shows up in binned residuals, heteroscedasticity fans the
residual spread, a leverage point moves the slope, and forcing the
intercept to zero costs accuracy.

`expected-output/FIELDS.md` separates the two categories in full.

## Installing

From the lab directory:

```bash
python3 -m venv .venv
.venv/bin/pip install -r requirements/requirements.txt
```

The install step needs the network. Everything after it is offline: the
diabetes dataset ships inside scikit-learn's own package data, and every
other dataset in this lab is generated on the spot from a seeded
generator.

## Free and open-source status

All three packages are free and open source — NumPy and scikit-learn under
the BSD 3-Clause licence, pytest under the MIT licence. There is no paid
tier, no account and no API key anywhere in this lab.

numpy==2.5.2
scikit-learn==1.9.0
pytest==9.1.1

# Day 148 lab brief — One Line, Measured

Everybody has seen a fitted line drawn through a scatterplot. Rather fewer
people can read the slope off in real units, say what the four assumptions
behind it are, or tell from a residual plot when the line is quietly
wrong.

This lab measures all of that on one predictor, one target, one line.

## The claim you are here to measure

> A slope is a real-unit statement, and it is only interpretable in real units.

Exercise 1 fits BMI against one-year diabetes-progression score, in the
dataset's *raw* units — not the mean-centred, unit-scaled default scikit-learn
ships, which throws away exactly the interpretability this lesson is
about.

| quantity | measured |
| --- | --- |
| slope | 10.2331 |
| intercept | −117.7734 |
| R-squared | 0.3439 |
| slope standard error | 0.6738 |
| 95% confidence interval | [8.9125, 11.5538] |

In one sentence a clinician could read: **each additional unit of BMI is
associated with about ten more points of one-year disease progression**,
give or take about 0.67 either way — and that slope sits roughly fifteen
standard errors from zero, which is not a borderline effect.

Two facts about that line hold on *any* dataset, exactly, forever, not
just this one: it passes through the point `(mean(x), mean(y))`, and its
residuals sum to zero. Exercise 1c measures both to floating-point
precision.

## The part that should worry you a little

A line can pass every glance-test — a clean scatterplot, a respectable
R-squared — and still be wrong in a way only the residuals reveal.

Exercise 3 fits a line to genuinely curved data. **The R-squared is
0.852.** That looks like a good fit. But bin the residuals by x and they
trace the missed curve exactly: positive at both ends, negative in the
middle. A quadratic curve explains over a third of the residuals' own
variance — a whole model's worth of structure the line called "noise".

Exercise 4 does the same with **heteroscedasticity**: noise whose spread
grows with x. The R-squared is 0.5723, also perfectly ordinary-looking.
The residual standard deviation more than doubles from the low half of x
to the high half — a fan shape invisible in the scatterplot and in the
single number, visible immediately in a residual plot.

Exercise 7 closes the loop: run the same quadratic-fit-to-residuals check
on the *real* BMI model, and it comes back at 0.0002 — essentially zero.
That contrast is how you tell "there is a missed curve" from "there is
just noise".

## The four assumptions, and what breaks each one

| Assumption | What breaks it, in this lab |
| --- | --- |
| Linearity | Exercise 3: a quadratic relationship fitted with a straight line |
| Constant variance (homoscedasticity) | Exercise 4: noise that grows with x |
| No point dominates the fit | Exercise 5: one point out of forty-one changes the slope by more than 85 percent |
| Roughly normal residuals | Exercise 7: a skewness check, and what "roughly" means in practice |

Exercise 5 is the one that should alarm you most. Forty ordinary points,
an unremarkable linear relationship — add one point far out on x with a
y-value that does not follow the trend, and the slope drops from 1.5196
to 0.2138. The mechanism has a name and a number: **leverage**, computed
from the point's x-value alone, before its y-value is even considered.
The added point's leverage is almost 27 times the average of the other
forty.

## How to work

1. Build the environment (see the lab `README.md`).
2. Run `.venv/bin/pytest starter -q`. You will see four passes (the
   machinery checks in `test_regression_lib.py`) and twelve skips.
3. Replace one `pytest.skip(...)` at a time with real code. The skip text
   names the exact helper and the exact value to assert.
4. Print the measured pair in every exercise. A number you did not print
   is a number you did not look at.
5. When you want the whole measured table at once, run
   `.venv/bin/python3 examples/report_measurements.py`.

Do not run `pytest starter examples` in one invocation. Both directories
define `regression_lib.py`, `test_regression_lib.py` and
`test_regression_claims.py`; pytest aborts on the module-name collision.
Run them separately, always.

## One scope note

This lab fits exactly one predictor at a time with scikit-learn's
`LinearRegression`. It never derives *why* squared error is the thing to
minimize (Day 149), never adds a second predictor (Day 150), and never
implements ordinary least squares from first principles (Day 153). The
question here is narrower and comes first: given a fitted line, what does
it actually tell you, and how do you catch it lying to you.

"""One predictor, one line: what fitting it actually buys you, measured.

Simple linear regression is a claim about a straight-line relationship
between one predictor and one target, fitted so the squared vertical
distances from the line are as small as possible (Day 149 owns *why*
squared error; this module just fits it, with scikit-learn's
`LinearRegression`). Two facts about that fit are exact on any data,
forever: it passes through the point of means, and its residuals sum to
zero when an intercept is fitted. Everything else here is a measurement of
what the line gets right and what it silently hides.

The module is organised in the order the lesson uses it:

1. The line itself, fitted to real data in real units (age-adjusted
   diabetes progression against BMI, raw scale).
2. Recovering a slope you know to be true, and watching the error shrink
   with more rows.
3. Two ways a fit can look fine on a scatterplot and an R-squared, and
   still be wrong: curvature and heteroscedasticity, both visible only in
   the residuals.
4. One point that moves the whole line, and the number that names why.
5. What forcing the intercept to zero costs, measured against the same
   data fitted honestly.

Everything here is deterministic given a seed. Nothing downloads: the
diabetes dataset is bundled with scikit-learn and the rest is generated on
the spot from `numpy.random.default_rng`.
"""

from __future__ import annotations

import numpy as np

from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression


# --------------------------------------------------------------------------
# 1. The line itself: BMI against disease progression, in raw units
# --------------------------------------------------------------------------


def load_bmi_and_target():
    """BMI and the disease-progression target, in their original units.

    `load_diabetes(scaled=False)` is the only way to get raw units out of
    this dataset -- the default returns every column mean-centred and
    scaled to unit norm, which makes a coefficient uninterpretable as
    "one more unit of X". Column 2 is BMI; see `load_diabetes().feature_names`.
    """
    data = load_diabetes(scaled=False)
    bmi = data.data[:, 2].reshape(-1, 1)
    y = data.target
    return bmi, y


def fit_line(x, y, fit_intercept: bool = True) -> LinearRegression:
    """Fit scikit-learn's `LinearRegression`. This module never derives it."""
    return LinearRegression(fit_intercept=fit_intercept).fit(x, y)


def slope_standard_error(x, residuals) -> float:
    """How much the fitted slope would wobble under a fresh sample.

    `SE(b1) = s / sqrt(sum((x - xbar)^2))`, where `s^2` is the residual
    variance with `n - 2` degrees of freedom spent on the slope and the
    intercept. This is the same standard-error arithmetic Days 117-118 and
    Day 144 used for a sampling proportion, applied here to a slope.
    """
    x = np.asarray(x, dtype=float).flatten()
    n = len(x)
    dof = n - 2
    s2 = float(np.sum(np.asarray(residuals) ** 2)) / dof
    sxx = float(np.sum((x - x.mean()) ** 2))
    return float(np.sqrt(s2 / sxx))


def confidence_interval(estimate: float, standard_error: float, z: float = 1.96) -> tuple:
    """A normal-approximation interval, rounded for reporting."""
    return (
        round(estimate - z * standard_error, 4),
        round(estimate + z * standard_error, 4),
    )


def passes_through_the_means(model: LinearRegression, x, y) -> tuple:
    """Confirms the fitted line predicts `mean(y)` exactly at `mean(x)`.

    Not approximately -- exactly, up to floating point. It is a property
    of the least-squares normal equations (Day 149's territory), not a
    coincidence of this dataset.
    """
    x = np.asarray(x, dtype=float)
    mean_x = x.mean(axis=0).reshape(1, -1)
    predicted_at_mean = float(model.predict(mean_x)[0])
    mean_y = float(np.asarray(y).mean())
    return predicted_at_mean, mean_y, predicted_at_mean - mean_y


def residual_sum(residuals) -> float:
    """The sum of the residuals, which is exactly zero when an intercept
    is fitted -- another consequence of the normal equations, not a
    measurement that happened to come out that way."""
    return float(np.sum(residuals))


# --------------------------------------------------------------------------
# 2. Recovering a slope you know to be true
# --------------------------------------------------------------------------


def make_known_line(n: int, seed: int, true_slope: float = 5.0, true_intercept: float = 10.0, noise_sd: float = 8.0):
    """One predictor, a known slope and intercept, and Gaussian noise."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 20, size=n).reshape(-1, 1)
    y = true_slope * x.flatten() + true_intercept + rng.normal(0, noise_sd, size=n)
    return x, y


def slope_recovery_error(n_values, replications: int = 200, true_slope: float = 5.0):
    """Mean absolute error of the fitted slope, at each sample size.

    Averaged over `replications` independently drawn datasets per `n`, for
    the same reason Days 117-118 and Day 144 always averaged rather than
    quoting one draw: a single fit's error is an anecdote.
    """
    rows = []
    for n in n_values:
        errors = [
            abs(float(fit_line(*make_known_line(n, seed, true_slope=true_slope)).coef_[0]) - true_slope)
            for seed in range(replications)
        ]
        rows.append((n, round(float(np.mean(errors)), 4)))
    return rows


# --------------------------------------------------------------------------
# 3a. Curvature: a fit that looks fine and is not
# --------------------------------------------------------------------------


def curved_dataset(n: int = 300, seed: int = 1):
    """A quadratic relationship, so a straight line is the wrong model."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 10, size=n)
    y = 2.0 + 0.5 * x**2 + rng.normal(0, 5.0, size=n)
    return x.reshape(-1, 1), y


def binned_residual_means(x, residuals, bins: int = 5):
    """Mean residual within each of `bins` equal-count groups, sorted by x.

    A straight line fitted to a curve leaves residuals that are positive
    at both ends and negative in the middle (or the reverse) -- a shape no
    single number like R-squared exposes, but a residual plot shows at a
    glance.
    """
    x = np.asarray(x, dtype=float).flatten()
    residuals = np.asarray(residuals, dtype=float)
    order = np.argsort(x)
    x_sorted, resid_sorted = x[order], residuals[order]
    groups = np.array_split(np.arange(len(x_sorted)), bins)
    return [
        (round(float(x_sorted[g].mean()), 2), round(float(resid_sorted[g].mean()), 4))
        for g in groups
    ]


def quadratic_fit_r_squared(x, residuals) -> float:
    """How much of the residuals' own variance a quadratic curve explains.

    Near zero when the residuals are patternless noise; large when the
    line missed real curvature. Fitted with `numpy.polyfit`, not
    scikit-learn -- this is a diagnostic on the residuals, not a second
    model of the data.
    """
    x = np.asarray(x, dtype=float).flatten()
    residuals = np.asarray(residuals, dtype=float)
    coeffs = np.polyfit(x, residuals, 2)
    predicted = np.polyval(coeffs, x)
    ss_res = float(np.sum((residuals - predicted) ** 2))
    ss_tot = float(np.sum((residuals - residuals.mean()) ** 2))
    return float(1.0 - ss_res / ss_tot)


# --------------------------------------------------------------------------
# 3b. Heteroscedasticity: error that grows with x
# --------------------------------------------------------------------------


def heteroscedastic_dataset(n: int = 400, seed: int = 2):
    """A line whose noise gets wider as x grows -- the fit stays roughly
    unbiased; only its residual spread reveals what is wrong."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(1, 20, size=n)
    noise_sd = 0.8 * x
    y = 3.0 + 2.0 * x + rng.normal(0, 1, size=n) * noise_sd
    return x.reshape(-1, 1), y


def residual_spread_by_half(x, residuals) -> tuple:
    """Residual standard deviation in the low-x half against the high-x half."""
    x = np.asarray(x, dtype=float).flatten()
    residuals = np.asarray(residuals, dtype=float)
    median = float(np.median(x))
    low = residuals[x < median]
    high = residuals[x >= median]
    return float(low.std()), float(high.std())


# --------------------------------------------------------------------------
# 4. One point that moves the line
# --------------------------------------------------------------------------


def leverage_dataset(n: int = 40, seed: int = 3):
    """A clean, ordinary linear relationship, forty points."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(0, 10, size=n)
    y = 2.0 + 1.5 * x + rng.normal(0, 1.5, size=n)
    return x, y


def add_point(x, y, x_new: float, y_new: float):
    """Append one point, returning fresh arrays (never mutates the inputs)."""
    return np.append(np.asarray(x, dtype=float), x_new), np.append(np.asarray(y, dtype=float), y_new)


def leverage_of_point(x, x_target: float) -> float:
    """The hat-matrix leverage of a point at `x_target`, given the full x array.

    `h = 1/n + (x_target - xbar)^2 / sum((x - xbar)^2)` -- how much that
    single point's own y-value can pull the fitted line toward itself,
    independent of what its y-value actually is.
    """
    x = np.asarray(x, dtype=float)
    n = len(x)
    xbar = x.mean()
    sxx = float(np.sum((x - xbar) ** 2))
    return float(1.0 / n + (x_target - xbar) ** 2 / sxx)


def mean_leverage_excluding(x, x_target: float) -> float:
    """Average leverage of every point except the one at `x_target`."""
    x = np.asarray(x, dtype=float)
    return float(np.mean([leverage_of_point(x, xi) for xi in x if xi != x_target]))


# --------------------------------------------------------------------------
# 5. fit_intercept=False, and what it costs
# --------------------------------------------------------------------------


def intercept_dataset(n: int = 200, seed: int = 4, true_intercept: float = 25.0, true_slope: float = 3.0, noise_sd: float = 6.0):
    """A line whose x-values never go near zero, so a forced-zero intercept
    is a real misspecification rather than a harmless simplification."""
    rng = np.random.default_rng(seed)
    x = rng.uniform(5, 25, size=n)
    y = true_intercept + true_slope * x + rng.normal(0, noise_sd, size=n)
    return x.reshape(-1, 1), y


def rmse(y_true, y_pred) -> float:
    y_true = np.asarray(y_true, dtype=float)
    y_pred = np.asarray(y_pred, dtype=float)
    return float(np.sqrt(np.mean((y_true - y_pred) ** 2)))


# --------------------------------------------------------------------------
# 6. Skewness, for a rough normality read
# --------------------------------------------------------------------------


def skewness(values) -> float:
    """The third standardised moment: 0 for a symmetric distribution.

    A rough diagnostic only -- not a formal normality test, and this
    module does not claim to be one. `abs(skewness) < 0.5` is a common
    rule of thumb for "not alarming".
    """
    values = np.asarray(values, dtype=float)
    mean = values.mean()
    sd = values.std()
    return float(np.mean(((values - mean) / sd) ** 3))

"""Twelve exercises in what a simple linear regression actually gets you,
and what it hides.

Read `00_brief.md` first. Each function below is a `pytest.skip` naming
exactly what to build and what to assert; replace the skip with real code.
`regression_lib.py` is complete -- it is the machinery, not the exercise.

Run this suite on its own:

    .venv/bin/pytest starter -q

Never run `pytest starter examples` in one invocation: both directories
define modules with the same names and pytest aborts on the collision.
"""

import numpy as np  # noqa: F401  (you will need it)
import pytest

import regression_lib as r  # noqa: F401  (you will need it)


# --- 1. The line, fitted to real data in real units ----------------------


def test_01_bmi_slope_and_intercept_in_raw_units():
    pytest.skip(
        "Load bmi, y with r.load_bmi_and_target() and fit r.fit_line(bmi, y). "
        "Assert model.coef_[0] rounds to 10.2331, model.intercept_ rounds to "
        "-117.7734, and model.score(bmi, y) rounds to 0.3439. That slope "
        "means one more unit of BMI is associated with about ten more points "
        "of one-year disease progression, on this population."
    )


def test_01b_slope_standard_error_and_confidence_interval():
    pytest.skip(
        "Fit the same BMI model, compute its residuals, and pass them to "
        "r.slope_standard_error(bmi, residuals). Assert it rounds to 0.6738. "
        "Then assert r.confidence_interval(slope, se) equals (8.9125, "
        "11.5538), and that slope / se rounds to 15.19 -- about fifteen "
        "standard errors from zero."
    )


def test_01c_the_line_passes_through_the_means_exactly():
    pytest.skip(
        "Call r.passes_through_the_means(model, bmi, y) and assert the "
        "difference is smaller than 1e-8, with both the predicted value at "
        "mean(bmi) and mean(y) rounding to 152.1335. Then assert "
        "r.residual_sum(residuals) is smaller than 1e-6 in absolute value -- "
        "not approximately zero, to within floating point."
    )


# --- 2. Recovering a slope you know to be true ----------------------------


def test_02_the_estimate_gets_closer_to_the_truth_as_n_grows():
    pytest.skip(
        "Call r.slope_recovery_error([20, 50, 200, 1000, 5000]) and assert it "
        "equals [(20, 0.2315), (50, 0.1556), (200, 0.078), (1000, 0.0357), "
        "(5000, 0.0159)]. Then assert the errors are strictly decreasing as n "
        "grows. The true slope is 5.0 and was never told to the fit."
    )


def test_02b_the_error_shrinks_roughly_like_one_over_root_n():
    pytest.skip(
        "Using the same rows, compute the ratio of the error at n=200 to the "
        "error at n=20 and assert it falls between 0.25 and 0.40 (predicted "
        "by one-over-root-n: 1/sqrt(10) = 0.316). Then assert the ratio from "
        "n=20 to n=5000 falls between 0.04 and 0.09 (predicted 1/sqrt(250) = "
        "0.063). This is a noisy quantity, so assert a range, not a formula."
    )


# --- 3. Curvature: a fit that looks fine and is not ------------------------


def test_03_residuals_reveal_the_curve_the_scatter_hides():
    pytest.skip(
        "Fit r.curved_dataset() with r.fit_line and assert the R-squared "
        "rounds to 0.852 -- a respectable-looking number. Then call "
        "r.binned_residual_means(x, residuals, bins=5) and assert the first "
        "bin's mean residual is positive, the middle bin's is negative, and "
        "the last bin's is positive: the missed curve, visible only in the "
        "residuals."
    )


def test_03b_quantifying_the_curvature_with_a_quadratic_fit_to_the_residuals():
    pytest.skip(
        "Assert r.quadratic_fit_r_squared(x, residuals) on the curved "
        "dataset rounds to 0.3558 -- a quadratic explains over a third of "
        "the RESIDUALS' own variance. Then assert "
        "np.corrcoef(residuals, x.flatten() ** 2)[0, 1] rounds to 0.148."
    )


# --- 4. Heteroscedasticity: error that grows with x -------------------------


def test_04_heteroscedasticity_fans_the_residuals_while_the_fit_looks_fine():
    pytest.skip(
        "Fit r.heteroscedastic_dataset() and assert the R-squared rounds to "
        "0.5723. Then call r.residual_spread_by_half(x, residuals) and "
        "assert the low-x-half standard deviation rounds to 4.7427, the "
        "high-x-half rounds to 12.0684, and their ratio rounds to 2.5446. "
        "Nothing in the scatterplot or the R-squared shows this; only the "
        "residual plot does."
    )


# --- 5. One point that moves the line ----------------------------------------


def test_05_one_high_leverage_point_moves_the_line():
    pytest.skip(
        "Fit r.leverage_dataset() with and without one extra point at "
        "x_new=40.0, y_new=5.0 added via r.add_point. Assert the slope "
        "without it rounds to 1.5196, the slope with it rounds to 0.2138, "
        "and the change rounds to -1.3059 -- one row out of forty-one "
        "cutting the slope by more than eighty-five percent."
    )


def test_05b_the_leverage_value_names_the_mechanism():
    pytest.skip(
        "With the extra point added, call r.leverage_of_point(x_with, 40.0) "
        "and r.mean_leverage_excluding(x_with, 40.0). Assert the first "
        "rounds to 0.8048, the second to 0.0299, and their ratio to 26.94 -- "
        "computed from the point's x-value alone, before its y-value is even "
        "considered."
    )


# --- 6. fit_intercept=False, and what it costs -------------------------------


def test_06_forcing_the_intercept_to_zero_costs_you():
    pytest.skip(
        "Fit r.intercept_dataset() with fit_intercept=True and False. Assert "
        "the RMSE with an intercept rounds to 6.1401 and without it rounds "
        "to 9.7878, a ratio of 1.5941. Then assert the intercept-free "
        "model's intercept_ is exactly 0.0 and its slope rounds to 4.4232 -- "
        "the true slope was 3.0 and the true intercept was 25.0; forcing the "
        "intercept to zero makes the slope absorb the difference."
    )


# --- 7. Telling curvature apart from noise ------------------------------------


def test_07_the_bmi_models_residuals_show_no_such_curvature():
    pytest.skip(
        "On the BMI model's own residuals, assert "
        "r.quadratic_fit_r_squared(bmi, residuals) rounds to 0.0002 -- "
        "contrast with exercise 3b's 0.3558 on data with real curvature. "
        "Then assert r.skewness(residuals) rounds to 0.156 and its absolute "
        "value is below 0.5: mildly asymmetric, nothing alarming."
    )

"""Machinery checks: the helpers behave, before any claim is made.

These four tests are solved in both `starter/` and `examples/`. They exist
so that a broken helper reports itself as a broken helper rather than as a
surprising scientific result.
"""

import numpy as np

import regression_lib as r


def test_load_bmi_and_target_has_the_right_shape_and_raw_units():
    bmi, y = r.load_bmi_and_target()
    assert bmi.shape == (442, 1)
    assert y.shape == (442,)
    # Raw units, not the mean-centred, unit-norm-scaled default.
    assert bmi.min() >= 18.0 and bmi.max() <= 42.2
    assert y.min() >= 25.0 and y.max() <= 346.0


def test_fit_line_respects_fit_intercept():
    x = np.array([[1.0], [2.0], [3.0], [4.0]])
    y = np.array([3.0, 5.0, 7.0, 9.0])
    with_intercept = r.fit_line(x, y, fit_intercept=True)
    without_intercept = r.fit_line(x, y, fit_intercept=False)
    assert round(float(with_intercept.coef_[0]), 4) == 2.0
    assert round(float(with_intercept.intercept_), 4) == 1.0
    assert without_intercept.intercept_ == 0.0


def test_leverage_of_a_central_point_is_smaller_than_an_extreme_one():
    x = np.linspace(0, 10, 11)
    central = r.leverage_of_point(x, 5.0)
    extreme = r.leverage_of_point(x, 40.0)
    assert extreme > central
    # Leverage is never below 1/n for any point, including the mean.
    assert r.leverage_of_point(x, float(x.mean())) >= 1.0 / len(x) - 1e-9


def test_skewness_is_zero_for_a_symmetric_sample():
    rng = np.random.default_rng(0)
    symmetric = rng.normal(size=5000)
    assert abs(r.skewness(symmetric)) < 0.1
    lopsided = np.concatenate([np.zeros(950), np.full(50, 20.0)])
    assert r.skewness(lopsided) > 1.0

#!/usr/bin/env bash
# Day 148 lab harness: "One Line, Measured"
#
# Prints "N checks, M failure(s)" and exits 0 only when M is zero.
set -u

LAB_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"
cd "$LAB_DIR"

PYTHON="${PYTHON:-.venv/bin/python3}"
PYTEST="${PYTEST:-.venv/bin/pytest}"

# Clear caches at the START so the final cleanliness check measures what
# THIS run left behind, not what a previous `pytest starter -q` left.
find . -path ./.venv -prune -o -type d -name '__pycache__' -exec rm -rf -- {} + 2>/dev/null
rm -rf .pytest_cache

CHECKS=0
FAILURES=0

ok() {
  CHECKS=$((CHECKS + 1))
  echo "  ok: $1"
}

fail() {
  CHECKS=$((CHECKS + 1))
  FAILURES=$((FAILURES + 1))
  echo "  FAIL: $1"
}

if [ ! -x "$PYTHON" ]; then
  echo "No lab .venv found at $PYTHON."
  echo "Run: python3 -m venv .venv && .venv/bin/pip install -r requirements/requirements.txt"
  exit 2
fi

echo "1. Installed versions match requirements/requirements.txt"
VERSION_CHECK=$("$PYTHON" - <<'PYEOF'
import numpy, sklearn, pytest
print("numpy", numpy.__version__)
print("scikit-learn", sklearn.__version__)
print("pytest", pytest.__version__)
PYEOF
)
echo "$VERSION_CHECK" | sed 's/^/    /'
while read -r pkg pin; do
  pin_version="${pin#*==}"
  installed=$(echo "$VERSION_CHECK" | awk -v p="$pkg" '$1==p {print $2}')
  if [ "$installed" = "$pin_version" ]; then
    ok "$pkg $installed matches the pin"
  else
    fail "$pkg installed=$installed pinned=$pin_version"
  fi
done < <(sed 's/==/ ==/' requirements/requirements.txt)

echo ""
echo "2. Every published claim, reproduced directly (no pytest involved)"
DIRECT_CHECK=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")

import numpy as np

import regression_lib as r

errors = []


def expect(label, got, want):
    if got != want:
        errors.append(f"{label}: expected {want}, got {got}")


# 1. The BMI line
bmi, y = r.load_bmi_and_target()
model = r.fit_line(bmi, y)
residuals = y - model.predict(bmi)
expect("bmi slope", round(float(model.coef_[0]), 4), 10.2331)
expect("bmi intercept", round(float(model.intercept_), 4), -117.7734)
expect("bmi r2", round(float(model.score(bmi, y)), 4), 0.3439)
se = r.slope_standard_error(bmi, residuals)
expect("bmi slope se", round(se, 4), 0.6738)
expect("bmi slope ci", r.confidence_interval(float(model.coef_[0]), se), (8.9125, 11.5538))
predicted_at_mean, mean_y, diff = r.passes_through_the_means(model, bmi, y)
if abs(diff) >= 1e-8:
    errors.append(f"line does not pass through the means: diff={diff}")
if abs(r.residual_sum(residuals)) >= 1e-6:
    errors.append(f"residuals do not sum to zero: {r.residual_sum(residuals)}")

# 2. Recovering a known slope
rows = r.slope_recovery_error([20, 50, 200, 1000, 5000])
expect(
    "slope recovery",
    rows,
    [(20, 0.2315), (50, 0.1556), (200, 0.078), (1000, 0.0357), (5000, 0.0159)],
)
errs = [e for _n, e in rows]
if not all(a > b for a, b in zip(errs, errs[1:])):
    errors.append("slope recovery error did not strictly decrease with n")

# 3. Curvature
xc, yc = r.curved_dataset()
model_c = r.fit_line(xc, yc)
residuals_c = yc - model_c.predict(xc)
expect("curved r2", round(float(model_c.score(xc, yc)), 4), 0.852)
bins = r.binned_residual_means(xc, residuals_c, bins=5)
means = [m for _x, m in bins]
if not (means[0] > 0 and means[2] < 0 and means[-1] > 0):
    errors.append(f"curvature bins did not show positive/negative/positive: {means}")
expect("curvature quadratic r2", round(r.quadratic_fit_r_squared(xc, residuals_c), 4), 0.3558)

# 4. Heteroscedasticity
xh, yh = r.heteroscedastic_dataset()
model_h = r.fit_line(xh, yh)
residuals_h = yh - model_h.predict(xh)
expect("hetero r2", round(float(model_h.score(xh, yh)), 4), 0.5723)
low_sd, high_sd = r.residual_spread_by_half(xh, residuals_h)
expect("hetero low sd", round(low_sd, 4), 4.7427)
expect("hetero high sd", round(high_sd, 4), 12.0684)
expect("hetero ratio", round(high_sd / low_sd, 4), 2.5446)

# 5. Leverage
xl, yl = r.leverage_dataset()
model_without = r.fit_line(xl.reshape(-1, 1), yl)
xl_with, yl_with = r.add_point(xl, yl, x_new=40.0, y_new=5.0)
model_with = r.fit_line(xl_with.reshape(-1, 1), yl_with)
expect("leverage slope without", round(float(model_without.coef_[0]), 4), 1.5196)
expect("leverage slope with", round(float(model_with.coef_[0]), 4), 0.2138)
expect(
    "leverage slope change",
    round(float(model_with.coef_[0]) - float(model_without.coef_[0]), 4),
    -1.3059,
)
leverage_new = r.leverage_of_point(xl_with, 40.0)
typical = r.mean_leverage_excluding(xl_with, 40.0)
expect("leverage value", round(leverage_new, 4), 0.8048)
expect("leverage typical", round(typical, 4), 0.0299)
expect("leverage ratio", round(leverage_new / typical, 2), 26.94)

# 6. fit_intercept=False
xi, yi = r.intercept_dataset()
model_yes = r.fit_line(xi, yi, fit_intercept=True)
model_no = r.fit_line(xi, yi, fit_intercept=False)
rmse_yes = r.rmse(yi, model_yes.predict(xi))
rmse_no = r.rmse(yi, model_no.predict(xi))
expect("rmse with intercept", round(rmse_yes, 4), 6.1401)
expect("rmse without intercept", round(rmse_no, 4), 9.7878)
expect("rmse ratio", round(rmse_no / rmse_yes, 4), 1.5941)

# 7. Telling curvature apart from noise
expect("bmi quadratic r2", round(r.quadratic_fit_r_squared(bmi, residuals), 4), 0.0002)
expect("bmi residual skew", round(r.skewness(residuals), 4), 0.156)

if errors:
    for e in errors:
        print("ERROR:", e)
    sys.exit(1)
print("all direct checks passed")
PYEOF
)
if echo "$DIRECT_CHECK" | grep -q "all direct checks passed"; then
  ok "exercises 1-7 reproduced directly against regression_lib, no pytest involved"
else
  fail "direct library checks failed"
  echo "$DIRECT_CHECK" | sed 's/^/    /'
fi

echo ""
echo "3. examples/ passes in full"
EXAMPLES_OUT=$("$PYTEST" examples -q 2>&1)
if echo "$EXAMPLES_OUT" | tail -1 | grep -qE "^16 passed"; then
  ok "pytest examples -q -> 16 passed"
else
  fail "pytest examples -q did not report 16 passed"
  echo "$EXAMPLES_OUT" | tail -20 | sed 's/^/    /'
fi

echo ""
echo "4. starter/ is an untouched skeleton"
STARTER_OUT=$("$PYTEST" starter -q 2>&1)
if echo "$STARTER_OUT" | tail -1 | grep -qE "4 passed, 12 skipped"; then
  ok "pytest starter -q -> 4 passed, 12 skipped (the machinery checks pass; the twelve exercises are stubs)"
else
  fail "pytest starter -q did not report 4 passed, 12 skipped"
  echo "$STARTER_OUT" | tail -20 | sed 's/^/    /'
fi

echo ""
echo "5. pytest examples starter (one invocation) aborts on the module-name collision"
COMBINED_OUT=$("$PYTEST" examples starter 2>&1)
if echo "$COMBINED_OUT" | grep -q "import file mismatch"; then
  ok "combined invocation reports import file mismatch, as documented -- never run starter and examples together"
else
  fail "combined invocation did not fail with import file mismatch as expected"
fi

echo ""
echo "6. The report reproduces the captured table exactly"
REPORT_OUT=$("$PYTHON" examples/report_measurements.py 2>&1)
if [ "$REPORT_OUT" = "$(cat expected-output/measured-values.txt)" ]; then
  ok "report_measurements.py output is byte-identical to expected-output/measured-values.txt"
else
  fail "report_measurements.py drifted from expected-output/measured-values.txt"
  echo "$REPORT_OUT" | diff - expected-output/measured-values.txt | head -20 | sed 's/^/    /'
fi

echo ""
echo "7. Proof the harness can fail"
SCRATCH=$(mktemp -d "${TMPDIR:-/tmp}/d148-scratch.XXXXXX")
cp examples/*.py "$SCRATCH"/
SCRATCH_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
if echo "$SCRATCH_OUT" | tail -1 | grep -qE "^16 passed"; then
  ok "scratch copy of examples/ passes before it is broken"
else
  fail "scratch copy did not pass before being broken: $(echo "$SCRATCH_OUT" | tail -3)"
fi
"$PYTHON" - "$SCRATCH/test_regression_claims.py" <<'PYEOF'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert round(rmse_no / rmse_yes, 4) == 1.5941"
replacement = "assert round(rmse_no / rmse_yes, 4) == 1.0"
assert needle in text, "could not find the assertion to break"
open(path, "w").write(text.replace(needle, replacement, 1))
PYEOF
BROKEN_OUT=$("$PYTEST" "$SCRATCH" -q 2>&1)
BROKEN_STATUS=$?
if [ "$BROKEN_STATUS" -ne 0 ] && echo "$BROKEN_OUT" | grep -q "test_06_forcing_the_intercept_to_zero_costs_you"; then
  ok "breaking exercise 6's assertion produces a non-zero exit and names the failing test"
else
  fail "broken copy did not fail as expected (exit=$BROKEN_STATUS)"
fi
rm -rf "$SCRATCH"

echo ""
echo "8. The direction of every result holds beyond the quoted seed"
DIRECTION=$("$PYTHON" - <<'PYEOF'
import sys
sys.path.insert(0, "examples")
import numpy as np
import regression_lib as r

problems = []

# Slope recovery error keeps shrinking at a different true slope and
# replication count.
rows = r.slope_recovery_error([20, 200, 2000], replications=80, true_slope=-3.0)
errs = [e for _n, e in rows]
if not all(a > b for a, b in zip(errs, errs[1:])):
    problems.append(f"slope recovery error did not shrink at true_slope=-3.0: {rows}")

# Curvature bins keep the same shape at other dataset seeds.
for seed in (2, 3, 4):
    x, y = r.curved_dataset(seed=seed)
    model = r.fit_line(x, y)
    resid = y - model.predict(x)
    means = [m for _x, m in r.binned_residual_means(x, resid, bins=5)]
    if not (means[0] > 0 and means[2] < 0 and means[-1] > 0):
        problems.append(f"seed {seed}: curvature bins lost their shape: {means}")

# Heteroscedasticity fans at other seeds.
for seed in (3, 4, 5):
    x, y = r.heteroscedastic_dataset(seed=seed)
    model = r.fit_line(x, y)
    resid = y - model.predict(x)
    low, high = r.residual_spread_by_half(x, resid)
    if high <= low * 1.5:
        problems.append(f"seed {seed}: heteroscedasticity ratio too small: {high / low:.4f}")

# The leverage point moves the slope substantially at other base seeds.
for seed in (10, 11, 12):
    x, y = r.leverage_dataset(seed=seed)
    m0 = r.fit_line(x.reshape(-1, 1), y)
    x2, y2 = r.add_point(x, y, 40.0, 5.0)
    m1 = r.fit_line(x2.reshape(-1, 1), y2)
    change = float(m1.coef_[0]) - float(m0.coef_[0])
    if change > -0.5:
        problems.append(f"seed {seed}: leverage point did not move the slope enough: {change:.4f}")

if problems:
    for p in problems:
        print("ERROR:", p)
else:
    print("every direction held")
PYEOF
)
if [ "$DIRECTION" = "every direction held" ]; then
  ok "slope recovery, curvature, heteroscedasticity and leverage all hold at seeds the lesson does not quote"
else
  fail "a direction failed beyond the quoted seed"
  echo "$DIRECTION" | sed 's/^/    /'
fi

echo ""
echo "9. Offline, and nothing left behind"
if ! grep -rInE "https?://" examples/*.py starter/*.py > /dev/null 2>&1; then
  ok "no URLs inside examples/ or starter/ source -- this lab reaches no network"
else
  fail "found a URL inside examples/ or starter/"
fi
if [ -z "$(find . -path ./.venv -prune -o -type d -name '__pycache__' -print 2>/dev/null)" ]; then
  ok "no __pycache__ left behind"
else
  find . -path ./.venv -prune -o -type d -name '__pycache__' -exec rm -rf -- {} + 2>/dev/null
  ok "no __pycache__ left behind (cleaned during this run)"
fi
if [ ! -d .pytest_cache ]; then
  ok "no .pytest_cache left behind"
else
  rm -rf .pytest_cache
  ok "no .pytest_cache left behind (cleaned during this run)"
fi

echo ""
echo "---------------------------------------------------------------"
echo "$CHECKS checks, $FAILURES failure(s)"
if [ "$FAILURES" -ne 0 ]; then
  exit 1
fi
exit 0