Learn how machines learn from data: implement regression, trees, metrics, k-means, PCA and anomaly detection from scratch in NumPy, and practise honest evaluation.

Machine Learning Foundations

## Why this matters

Marketing teams want customer segments, biologists want groups of cells with similar gene expression, operations teams want to know which stores behave alike. None of these come with labels. **Clustering** finds groups of rows that are similar to each other and different from other groups. k-means is the simplest and most widely used clustering algorithm, and it is a good way to see how unsupervised learning works.

## The k-means idea

You choose the number of clusters `k`. Each cluster is represented by a **centroid**, the average point of its members. k-means looks for centroids that make the total squared distance from each row to its own centroid as small as possible. That total is called **inertia** (or within-cluster sum of squares).

## The algorithm

```text
1. pick k rows at random as the starting centroids
2. repeat until nothing changes:
     assign: give each row to its nearest centroid
     update: move each centroid to the mean of its rows
```

Each assign step and each update step can only lower the inertia, so the algorithm always settles, usually within a few dozen iterations.

### Worked example in one dimension

Points: `1, 2, 3, 10, 11, 12`, `k = 2`, starting centroids `1` and `2`.

```text
assign: {1} -> centroid 1;  {2, 3, 10, 11, 12} -> centroid 2
update: centroids become 1 and (2+3+10+11+12)/5 = 7.6
assign: {1, 2, 3} closer to 1;  {10, 11, 12} closer to 7.6
update: centroids become 2 and 11
assign: nothing changes -> done
```

Two clusters, centred at 2 and 11, exactly as you would draw them by eye.

### Vectorised assignment

With NumPy you compute all row-to-centroid distances at once by broadcasting:

```python
d = ((X[:, None, :] - centroids[None, :, :]) ** 2).sum(axis=2)   # shape (n_rows, k)
labels = d.argmin(axis=1)
```

## Scale your features first

k-means uses Euclidean distance, so a feature measured in large units dominates. If annual spend ranges over 0–120 (thousand dollars) and visits per month over 0–10, distances are almost entirely about spend. **Standardise** each column (subtract the mean, divide by the standard deviation) before clustering, then convert centroids back to original units to describe the segments.

## Initialisation and restarts

Random starting centroids can land badly, for example two centroids inside one true group, and k-means then gets stuck in a poor solution (a **local minimum**). Two standard fixes:

- **Restarts**: run k-means several times with different seeds and keep the result with the lowest inertia.
- **k-means++**: choose starting centroids spread far apart. It is scikit-learn's default.

## Choosing k

Inertia always falls as `k` grows (with `k = n` it is zero), so you cannot just minimise it. Common approaches:

- **Elbow plot**: plot inertia against `k` and look for the point where the curve bends and further clusters help little.
- **Silhouette score** (`sklearn.metrics.silhouette_score`): compares each row's distance to its own cluster with its distance to the next cluster; higher is better, from -1 to 1.
- **Usefulness**: can the business act on 4 segments? Would 12 be too many to name?

## Interpreting segments

After clustering, describe each cluster with its centroid (in original units) and size, and give it a plain name:

```text
segment   spend ($k)   visits/month   size   name
   A          20            2          60    occasional low spenders
   B          60            8          60    frequent regulars
   C         110            4          60    high-value planners
```

Then check that the segments are **stable** (similar on a new sample or another seed) and **different on things you did not cluster on**, such as churn rate. That is what makes them useful.

## In scikit-learn

```python
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans

X = StandardScaler().fit_transform(raw)
km = KMeans(n_clusters=3, n_init=10, random_state=0).fit(X)
print(km.cluster_centers_, km.inertia_)
labels = km.labels_
```

## Limits of k-means

k-means assumes clusters that are roughly round and of similar size. It struggles with long, curved or very unequal clusters, and every point is forced into some cluster, including outliers. Alternatives: **DBSCAN** (density-based, finds arbitrary shapes and labels noise), **Gaussian mixture models** (soft, elliptical clusters) and **hierarchical clustering** (a tree of merges you can cut at any level).

## Common mistakes

- **Not scaling features**, so one column decides everything.
- **A single run with one seed**, which can return a poor local minimum.
- **Treating clusters as ground truth.** They are a useful summary, not natural kinds.
- **Clustering on IDs or dates** encoded as numbers.

## Quick recap

- k-means alternates assign (nearest centroid) and update (cluster mean) to reduce inertia.
- Standardise first, use several restarts or k-means++, and pick `k` with an elbow plot, silhouette score and judgement.
- Describe each segment by its centroid in original units and check it is stable and useful.

## Try it

On the Do tab, write the assign and update steps and segment 180 customers into three groups.


Segment customers with k-means, choose the number of clusters, and describe each segment in plain terms.

Clustering

StatQuest: K-means clustering

## Why this matters

Real datasets often have dozens or hundreds of columns that move together: height, arm span and shoe size; the prices of related stocks; thousands of gene-expression levels driven by a few biological processes. **Dimensionality reduction** replaces many correlated columns with a few new ones that keep most of the information. That lets you plot high-dimensional data in 2D, speed up models, reduce noise, and see the hidden factors in the data. **Principal component analysis (PCA)** is the standard first tool.

## The idea

PCA finds new axes, called **principal components**, for the data:

1. The first component is the direction along which the data varies the most.
2. The second is the direction of most remaining variance, at a right angle (orthogonal) to the first.
3. And so on, up to the number of original features.

Projecting the data onto the first few components gives a compact version that keeps as much spread (variance) as possible. If five measurements are really driven by two hidden factors, the first two components will capture almost all the variance.

## The linear algebra

This is where your linear algebra pays off. For a data matrix `X` with `n` rows:

```text
1. centre:      Xc = X - column means
2. covariance:  C  = Xc.T @ Xc / (n - 1)          (features x features)
3. eigen:       C v = lambda v
4. sort eigenvalues from largest to smallest; the eigenvectors are the components
5. project:     Z = Xc @ V_k                      (n x k)
```

Each eigenvalue is the variance along its component, so the **explained variance ratio** of component `i` is `lambda_i / sum(lambda)`.

```python
import numpy as np

Xc = X - X.mean(axis=0)
cov = Xc.T @ Xc / (len(X) - 1)
values, vectors = np.linalg.eigh(cov)        # eigh: for symmetric matrices, ascending order
order = np.argsort(values)[::-1]
values, vectors = values[order], vectors[:, order]
ratio = values / values.sum()
```

Libraries usually compute the same thing with the **singular value decomposition** (SVD) of `Xc`, which is numerically more stable: the squared singular values divided by `n - 1` are the eigenvalues.

### Worked example: two correlated columns

Suppose height and arm span (in cm) have variances 100 and 110 and covariance 100, an almost perfect correlation of about 0.95. The covariance matrix's eigenvalues are about 205 and 5. The first component, roughly `(0.69, 0.73)`, is "overall body size" and explains `205 / 210 ≈ 97.6%` of the variance. The second, roughly `(-0.73, 0.69)`, is "arms long for your height" and explains the remaining 2.4%. Keeping one component loses very little.

## How many components?

- Plot the explained variance ratio for each component (a **scree plot**) and look for the elbow.
- Or keep enough components to reach a target such as 90% or 95% of cumulative variance.
- For visualisation, use 2 (or 3), whatever they explain, and report the percentage on the axes.

## Scale first (usually)

PCA chases variance, so columns in big units dominate, just as in k-means. If columns are in different units (dollars, years, counts), **standardise** them first; this is PCA on the correlation matrix. If they share units and their scale is meaningful (pixel intensities, expression levels already normalised), you may centre only.

## Reading the components

Each component is a weighted combination of the original features. The weights are called **loadings**. A component with large positive loadings on all the size measurements is a "size" factor. Eigenvector signs are arbitrary: `v` and `-v` are equally valid, so do not read meaning into the sign alone, and compare components across runs only up to sign.

## In scikit-learn

```python
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline

pipe = make_pipeline(StandardScaler(), PCA(n_components=0.95))   # keep 95% of variance
Z_train = pipe.fit_transform(X_train)
Z_test = pipe.transform(X_test)              # fitted on train only
print(pipe[-1].explained_variance_ratio_)
```

For plotting complex, non-linear structure, people also use **t-SNE** and **UMAP**. They make good pictures but distort distances, so use them for exploration, not as model features.

## Common mistakes

- **Forgetting to centre**, which makes the first component point at the mean instead of the spread.
- **Not scaling mixed-unit features**, so PCA just finds the largest-unit column.
- **Fitting PCA on all data before splitting**, which leaks test information (fit on train, transform test).
- **Expecting components to be causes.** They are directions of variance, which may or may not match real factors.
- **Discarding low-variance components that carry the signal.** PCA ignores the label; a rare but important pattern can live in a small component.

## Quick recap

- PCA finds orthogonal directions of maximum variance: the eigenvectors of the covariance matrix.
- Explained variance ratio = eigenvalue / total; keep enough components for your purpose.
- Centre always, standardise mixed units, fit on training data only.
- Loadings show what each component means; signs are arbitrary.

## Try it

On the Do tab, implement PCA with `np.linalg.eigh` and show that two components capture almost all the variance of five measurements.


Compress many correlated features into a few principal components and read what they mean.

Dimensionality reduction

StatQuest: PCA main ideas in only 5 minutes!!!

StatQuest: Principal Component Analysis (PCA), Step-by-Step

## Why this matters

Fraudulent payments, failing machines, hacked accounts and mislabelled lab samples have something in common: they are **rare and different**. You usually have few or no labelled examples of them, so supervised learning is hard. **Anomaly detection** scores how unusual each row is so that people can review the strangest ones first. Done well, it turns a million transactions into a short list worth looking at.

## Kinds of anomalies

- **Point anomalies**: a single row is extreme, such as a $40,000 purchase on a card that usually spends $50.
- **Contextual anomalies**: normal in one context, odd in another, such as 30 °C in January.
- **Collective anomalies**: a group of rows that is odd together, such as a burst of logins from one IP address.

This lesson focuses on point anomalies in tabular data.

## Simple baselines first

For one numeric column, classic rules work well:

- **z-score**: flag `|x - mean| / sd > 3`.
- **IQR rule**: flag values below `Q1 - 1.5*IQR` or above `Q3 + 1.5*IQR` (the box-plot rule, robust to the outliers themselves).

They break down with several features. A sensor reading can have a normal temperature and a normal vibration level individually, yet the **combination** is impossible for a healthy machine. You need a multivariate method.

## The isolation forest idea

Instead of modelling what normal looks like, an isolation forest asks: **how easy is this point to isolate?**

Build a tree by repeatedly picking a random feature and a random threshold between that feature's minimum and maximum, splitting the data until each point sits alone. A normal point lives in a dense crowd and needs many random cuts to separate. An outlier sits far from the crowd, so a random cut isolates it in just a few steps.

```text
dense normal point:  ~8-10 random splits to isolate
far-away outlier:    ~2-3 random splits to isolate
```

One random tree is noisy, so the forest builds many (often 100), each on a small random **subsample** of rows (256 by default), and averages each point's **path length** (the depth at which it gets isolated). Small subsamples actually help: outliers are less likely to be masked by other outliers.

## The anomaly score

Path lengths are normalised by `c(n)`, the average path length of an unsuccessful search in a binary search tree of `n` points:

```text
c(n) = 2 * (ln(n - 1) + 0.5772) - 2 * (n - 1) / n       (c(2) = 1, c(1) = 0)
score(x) = 2 ^ ( -mean_path(x) / c(sample_size) )
```

(0.5772 is the Euler–Mascheroni constant.) The score is between 0 and 1:

- close to **1**: isolated very quickly, a strong anomaly;
- around **0.5**: path length is average, nothing special;
- well **below 0.5**: deep inside the crowd, clearly normal.

`c` is also used inside the trees: when the depth limit stops a branch that still holds `m` points, the path gets `c(m)` added as an estimate of the remaining depth.

### Worked example

With a subsample of 256, `c(256) = 2 * (ln 255 + 0.5772) - 2 * 255/256 ≈ 10.24`. A point isolated after 3 splits on average scores `2^(-3/10.24) ≈ 0.82`, a strong anomaly. A point that needs 10.24 splits scores exactly 0.5.

## Turning scores into flags

Scores rank rows; you still need a cut-off. In practice you pick a **review budget**: "the analysts can check 2% of readings a day", so flag the top 2%. scikit-learn calls this `contamination`. If you do have a few confirmed anomalies, use them to check precision among the flagged rows (Lesson 6's metrics).

```python
from sklearn.ensemble import IsolationForest

iso = IsolationForest(n_estimators=200, contamination=0.02, random_state=0).fit(X_train)
flags = iso.predict(X_new) == -1        # -1 = anomaly
strength = -iso.score_samples(X_new)    # higher = more anomalous
```

Isolation forests need no feature scaling (cuts are random within each feature's range) and are fast. Alternatives include **Local Outlier Factor** (compares local density), **one-class SVM** and, for complex data, autoencoders.

## Common mistakes

- **Treating every flag as a confirmed problem.** An anomaly is unusual, not necessarily bad; a human or a follow-up rule decides.
- **Training on data full of the anomalies you want to find**, which makes them look normal. Train on a period you believe is mostly clean when you can.
- **Including ID or timestamp columns**, which makes every new row look unusual.
- **Ignoring context**, such as comparing night-time and daytime readings together.
- **Setting the contamination rate without thinking about review capacity.**

## Quick recap

- Anomaly detection ranks rows by how unusual they are, usually without labels.
- z-scores and IQR rules work for single columns; multivariate data needs methods like isolation forests.
- Isolation forests isolate points with random splits; short average paths mean anomalies.
- Score = `2^(-mean_path / c(n))`; near 1 is anomalous, about 0.5 is ordinary. Flag the top share your team can review.

## Try it

On the Do tab, write the normalising constant, the score and the flagging rule, then find the faulty sensor readings.


Score how unusual each row is with an isolation forest and flag the readings worth reviewing.

Anomaly detection

Isolation Forests: Identify Outliers in Data

Group similar rows with k-means, compress correlated features with PCA, and flag unusual records with an isolation forest.