Two workhorses of data analysis, from the exercises in Theoretical Neuroscience. Everything on this page runs in your browser; nothing to install. Change a slider while a simulation is running and it reacts immediately.
exFont = "16px " + getComputedStyle(document.body).fontFamily
exFg = () => getComputedStyle(document.body).color
mkCanvas = (W, H) => {
const c = DOM.context2d(W, H);
c.canvas.style.width = "100%"; c.canvas.style.maxWidth = W + "px"; c.canvas.style.height = "auto";
return c;
}
axisBox = (c, x, y, w, h) => { c.strokeStyle = "#999"; c.lineWidth = 1; c.strokeRect(x, y, w, h); }Fisher’s iris data: 150 flowers, four measurements each. PCA finds the directions of largest variance, the eigenvectors of the covariance matrix.
\[C = \frac{1}{n} X^{\mathsf T} X, \qquad C\, e_k = \lambda_k e_k\]
irisRaw = FileAttachment("iris.json").json()
viewof pcStd = Inputs.toggle({label: "Standardise each variable (z-score)", value: true})
pcAxisOpts = new Map([["PC 1", "pc0"], ["PC 2", "pc1"], ["PC 3", "pc2"], ["PC 4", "pc3"],
...irisRaw.features.map((f, i) => [f, "v" + i])])
viewof pcX = Inputs.select(pcAxisOpts, {value: "pc0", label: "Horizontal axis"})
viewof pcY = Inputs.select(pcAxisOpts, {value: "pc1", label: "Vertical axis"})
viewof pcColour = Inputs.toggle({label: "Colour by species", value: true})symEig = A => { // Jacobi eigenvalue iteration for a small symmetric matrix
const n = A.length, a = A.map(r => r.slice()), v = Array.from({length: n}, (_, i) => Array.from({length: n}, (_, j) => +(i === j)));
for (let sweep = 0; sweep < 60; sweep++) {
let off = 0; for (let i = 0; i < n; i++) for (let j = i + 1; j < n; j++) off += a[i][j] ** 2;
if (off < 1e-18) break;
for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) {
if (Math.abs(a[p][q]) < 1e-15) continue;
const th = (a[q][q] - a[p][p]) / (2 * a[p][q]), t = Math.sign(th || 1) / (Math.abs(th) + Math.sqrt(th * th + 1));
const cs = 1 / Math.sqrt(t * t + 1), sn = t * cs;
for (let k = 0; k < n; k++) { const akp = a[k][p], akq = a[k][q]; a[k][p] = cs * akp - sn * akq; a[k][q] = sn * akp + cs * akq; }
for (let k = 0; k < n; k++) { const apk = a[p][k], aqk = a[q][k]; a[p][k] = cs * apk - sn * aqk; a[q][k] = sn * apk + cs * aqk; }
for (let k = 0; k < n; k++) { const vkp = v[k][p], vkq = v[k][q]; v[k][p] = cs * vkp - sn * vkq; v[k][q] = sn * vkp + cs * vkq; }
}
}
const order = [...Array(n).keys()].sort((i, j) => a[j][j] - a[i][i]);
return {values: order.map(i => a[i][i]), vectors: order.map(i => { const col = v.map(r => r[i]); const big = col.reduce((b, x) => Math.abs(x) > Math.abs(b) ? x : b, 0); return big < 0 ? col.map(x => -x) : col; })};
}
pca = {
const rows = irisRaw.rows, d = 4, n = rows.length;
const mean = Array.from({length: d}, (_, j) => rows.reduce((s, r) => s + r[j], 0) / n);
const sd = Array.from({length: d}, (_, j) => Math.sqrt(rows.reduce((s, r) => s + (r[j] - mean[j]) ** 2, 0) / n));
const X = rows.map(r => Array.from({length: d}, (_, j) => (r[j] - mean[j]) / (pcStd ? sd[j] : 1)));
const C = Array.from({length: d}, (_, i) => Array.from({length: d}, (_, j) => X.reduce((s, x) => s + x[i] * x[j], 0) / n));
const {values, vectors} = symEig(C);
const scores = X.map(x => vectors.map(e => e.reduce((s, ei, j) => s + ei * x[j], 0)));
const total = values.reduce((s, l) => s + l, 0);
return {X, C, values, vectors, scores, total, species: rows.map(r => r[4]), mean, sd};
}
pcCanvas = mkCanvas(760, 580){
const c = pcCanvas, fg = exFg(), W = 760, H = 580, cols = ["#1f77b4", "#ff7f0e", "#2ca02c"];
const {X, values, vectors, scores, total, species} = pca, feats = irisRaw.features;
c.clearRect(0, 0, W, H); c.font = exFont;
const get = (key, i) => key.startsWith("pc") ? scores[i][+key[2]] : X[i][+key[1]];
const label = key => key.startsWith("pc") ? `PC ${+key[2] + 1}` : feats[+key[1]] + (pcStd ? " (z-scored)" : " (cm, centred)");
const xs = X.map((_, i) => get(pcX, i)), ys = X.map((_, i) => get(pcY, i));
const lim = Math.max(...xs.map(Math.abs), ...ys.map(Math.abs)) * 1.1;
const sP = {x0: 90, y0: 12, s: 360};
const px = v => sP.x0 + sP.s * (v + lim) / (2 * lim), py = v => sP.y0 + sP.s * (1 - (v + lim) / (2 * lim));
axisBox(c, sP.x0, sP.y0, sP.s, sP.s);
c.strokeStyle = "#ccc"; c.lineWidth = 1; c.beginPath(); c.moveTo(px(0), sP.y0); c.lineTo(px(0), sP.y0 + sP.s); c.moveTo(sP.x0, py(0)); c.lineTo(sP.x0 + sP.s, py(0)); c.stroke();
// principal axes drawn in the plane of two raw variables
if (!pcX.startsWith("pc") && !pcY.startsWith("pc")) {
const i = +pcX[1], j = +pcY[1];
for (const [k, col] of [[0, "#d62728"], [1, "#9467bd"]]) {
const e = vectors[k], L = 2 * Math.sqrt(values[k]);
c.strokeStyle = col; c.lineWidth = 3; c.beginPath(); c.moveTo(px(-L * e[i]), py(-L * e[j])); c.lineTo(px(L * e[i]), py(L * e[j])); c.stroke();
c.fillStyle = col; c.textAlign = "left"; c.fillText(`PC ${k + 1}`, px(L * e[i]) + 6, py(L * e[j]) + 5);
}
}
X.forEach((_, i) => { c.fillStyle = pcColour ? cols[species[i]] : "#555"; c.globalAlpha = 0.6; c.beginPath(); c.arc(px(xs[i]), py(ys[i]), 4.5, 0, 2 * Math.PI); c.fill(); });
c.globalAlpha = 1; c.fillStyle = fg; c.textAlign = "center";
c.fillText(label(pcX), sP.x0 + sP.s / 2, sP.y0 + sP.s + 24);
c.fillText((-lim).toFixed(1), sP.x0, sP.y0 + sP.s + 44); c.fillText(lim.toFixed(1), sP.x0 + sP.s, sP.y0 + sP.s + 44);
c.save(); c.translate(sP.x0 - 60, sP.y0 + sP.s / 2); c.rotate(-Math.PI / 2); c.fillText(label(pcY), 0, 0); c.restore();
c.textAlign = "right"; c.fillText(lim.toFixed(1), sP.x0 - 8, sP.y0 + 6); c.fillText((-lim).toFixed(1), sP.x0 - 8, sP.y0 + sP.s + 2);
if (pcColour) irisRaw.species.forEach((s, k) => { c.fillStyle = cols[k]; c.textAlign = "left"; c.fillText(s, sP.x0 + 10, sP.y0 + 22 + 22 * k); });
// explained variance bars, right
const bP = {x0: 520, y0: 12, w: 220, h: 150};
axisBox(c, bP.x0, bP.y0, bP.w, bP.h);
values.forEach((l, k) => {
const frac = l / total, bw = bP.w / 4 * 0.7, bx = bP.x0 + bP.w * (k + 0.15) / 4;
c.fillStyle = "#1f77b4"; c.fillRect(bx, bP.y0 + bP.h * (1 - frac), bw, bP.h * frac);
c.fillStyle = fg; c.textAlign = "center"; c.fillText(`${(100 * frac).toFixed(0)}%`, bx + bw / 2, bP.y0 + bP.h * (1 - frac) - 6);
c.fillText(`PC ${k + 1}`, bx + bw / 2, bP.y0 + bP.h + 22);
});
c.fillStyle = fg; c.textAlign = "center"; c.fillText("explained variance", bP.x0 + bP.w / 2, bP.y0 + bP.h + 44);
// loadings table, full width below the scatter
const tY = 470, colX = j => 250 + 125 * j;
c.fillStyle = fg; c.textAlign = "left"; c.fillText("loadings: components of each eₖ", 90, tY);
c.font = "14px " + getComputedStyle(document.body).fontFamily;
feats.forEach((f, j) => { c.textAlign = "center"; c.fillText(f, colX(j), tY + 24); });
values.forEach((_, k) => {
c.textAlign = "left"; c.fillText(`PC ${k + 1}`, 90, tY + 46 + 20 * k);
vectors[k].forEach((e, j) => { c.textAlign = "center"; c.fillText(e.toFixed(2), colX(j), tY + 46 + 20 * k); });
});
return c.canvas;
}What you see. Left: each flower as a dot in the chosen coordinates. With species colouring on, blue is setosa, orange versicolor, green virginica; off, every dot is grey, which is what PCA itself sees, since it never uses the labels. When both axes are raw measurements, the first two principal components are drawn as lines through the origin, PC 1 in red and PC 2 in purple, each scaled to two standard deviations of the data along it. Right, top: the share of the total variance carried by each component, as blue bars. Right, bottom: the loadings, how much of each measurement goes into each component.
Try: switch standardisation off and PC 1 is dominated by petal length, simply because it has the largest spread in centimetres. Look at the raw-variable planes to see where PC 1 points. From the dimensionality-reduction exercise in Theoretical Neuroscience.
A signal made of two sine waves is sampled at \(N\) points spaced \(\Delta x\) apart over an interval of length \(T\). The discrete Fourier transform (DFT) decomposes the samples into frequencies \(f_m = m / T\):
\[y(x) = a_1 \sin(2\pi f_1 x) + a_2 \sin(2\pi f_2 x), \qquad Y_m = \sum_{n=0}^{N-1} y_n\, e^{-2\pi i\, m n / N}\]
viewof ftF1 = Inputs.range([0, 12], {value: 1, step: 0.1, label: "f₁ (cycles / unit)"})
viewof ftA1 = Inputs.range([0, 2], {value: 1, step: 0.1, label: "a₁"})
viewof ftF2 = Inputs.range([0, 12], {value: 3.5, step: 0.1, label: "f₂ (cycles / unit)"})
viewof ftA2 = Inputs.range([0, 2], {value: 0.5, step: 0.1, label: "a₂"})
viewof ftDx = Inputs.range([0.02, 0.5], {value: 0.1, step: 0.01, label: "Sampling interval Δx"})
viewof ftT = Inputs.range([1, 8], {value: 4, step: 0.5, label: "Interval length T"})
viewof ftK = Inputs.range([0, 10], {value: 2, step: 1, label: "K components kept for reconstruction"})ft = {
const T = ftT, dx = ftDx, N = Math.max(4, Math.floor(T / dx)), x0 = -T / 2;
const sig = x => ftA1 * Math.sin(2 * Math.PI * ftF1 * x) + ftA2 * Math.sin(2 * Math.PI * ftF2 * x);
const xs = Array.from({length: N}, (_, n) => x0 + n * dx), ys = xs.map(sig);
const M = Math.floor(N / 2), re = new Float64Array(M + 1), im = new Float64Array(M + 1);
for (let m = 0; m <= M; m++) for (let n = 0; n < N; n++) { const a = -2 * Math.PI * m * n / N; re[m] += ys[n] * Math.cos(a); im[m] += ys[n] * Math.sin(a); }
const amp = Array.from({length: M + 1}, (_, m) => ({m, f: m / T, A: (m === 0 ? 1 : 2) * Math.hypot(re[m], im[m]) / N, ph: Math.atan2(im[m], re[m])}));
const kept = amp.slice(1).sort((a, b) => b.A - a.A).slice(0, ftK);
const rec = x => amp[0].A * Math.cos(amp[0].ph) + kept.reduce((s, c) => s + c.A * Math.cos(2 * Math.PI * c.f * (x - x0) + c.ph), 0);
const fine = Array.from({length: 600}, (_, i) => x0 + T * i / 599);
const fNy = 1 / (2 * dx);
const alias = f => { const fs = 2 * fNy; let g = Math.abs(f - Math.round(f / fs) * fs); return g; };
return {T, dx, N, x0, xs, ys, fine, sig, rec, amp, kept, fNy, alias};
}
ftCanvas = mkCanvas(760, 450){
const c = ftCanvas, fg = exFg(), W = 760, H = 450, blue = "#1f77b4", orange = "#ff7f0e";
const {T, N, x0, xs, ys, fine, sig, rec, amp, kept, fNy, alias} = ft;
c.clearRect(0, 0, W, H); c.font = exFont;
// signal panel
const sP = {x0: 80, y0: 12, w: W - 100, h: 170};
const ymax = Math.max(0.5, ftA1 + ftA2) * 1.15;
const sx = x => sP.x0 + sP.w * (x - x0) / T, sy = y => sP.y0 + sP.h * (1 - (y + ymax) / (2 * ymax));
axisBox(c, sP.x0, sP.y0, sP.w, sP.h);
c.strokeStyle = "#888"; c.lineWidth = 1; c.beginPath(); c.moveTo(sP.x0, sy(0)); c.lineTo(sP.x0 + sP.w, sy(0)); c.stroke();
c.strokeStyle = "#aaa"; c.lineWidth = 1.5; c.beginPath(); fine.forEach((x, i) => i ? c.lineTo(sx(x), sy(sig(x))) : c.moveTo(sx(x), sy(sig(x)))); c.stroke();
c.strokeStyle = orange; c.lineWidth = 2; c.beginPath(); fine.forEach((x, i) => i ? c.lineTo(sx(x), sy(rec(x))) : c.moveTo(sx(x), sy(rec(x)))); c.stroke();
c.fillStyle = blue; xs.forEach((x, n) => { c.beginPath(); c.arc(sx(x), sy(ys[n]), 3.5, 0, 2 * Math.PI); c.fill(); });
c.fillStyle = fg; c.textAlign = "center"; c.fillText(`x · N = ${N} samples`, sP.x0 + sP.w / 2, sP.y0 + sP.h + 22);
c.textAlign = "right"; c.fillText("y", sP.x0 - 10, sP.y0 + sP.h / 2 + 5);
c.fillText(ymax.toFixed(1), sP.x0 - 8, sP.y0 + 6); c.fillText((-ymax).toFixed(1), sP.x0 - 8, sP.y0 + sP.h + 2);
// spectrum panel
const pP = {x0: 80, y0: 240, w: W - 100, h: 140};
const fmax = Math.max(fNy * 1.05, 0.5), Amax = Math.max(0.5, ...amp.slice(1).map(a => a.A)) * 1.15;
const fx = f => pP.x0 + pP.w * f / fmax, fy = A => pP.y0 + pP.h * (1 - A / Amax);
axisBox(c, pP.x0, pP.y0, pP.w, pP.h);
c.setLineDash([4, 4]); c.strokeStyle = "red"; c.beginPath(); c.moveTo(fx(fNy), pP.y0); c.lineTo(fx(fNy), pP.y0 + pP.h); c.stroke();
c.strokeStyle = "#2ca02c";
for (const [f, a] of [[ftF1, ftA1], [ftF2, ftA2]]) if (a > 0) { const g = alias(f); c.beginPath(); c.moveTo(fx(g), pP.y0); c.lineTo(fx(g), pP.y0 + pP.h); c.stroke(); }
c.setLineDash([]);
for (const a of amp.slice(1)) {
const isKept = kept.includes(a);
c.strokeStyle = isKept ? orange : blue; c.lineWidth = isKept ? 3 : 1.5;
c.beginPath(); c.moveTo(fx(a.f), fy(0)); c.lineTo(fx(a.f), fy(a.A)); c.stroke();
c.fillStyle = isKept ? orange : blue; c.beginPath(); c.arc(fx(a.f), fy(a.A), 3, 0, 2 * Math.PI); c.fill();
}
c.fillStyle = fg; c.textAlign = "center";
c.fillText(`frequency (cycles / unit) · Nyquist ${fNy.toFixed(2)}`, pP.x0 + pP.w / 2, pP.y0 + pP.h + 22);
c.fillText("0", fx(0), pP.y0 + pP.h + 42); c.fillText(fmax.toFixed(1), fx(fmax), pP.y0 + pP.h + 42);
c.textAlign = "right"; c.fillText("|Y|·2/N", pP.x0 - 10, pP.y0 + pP.h / 2 + 5);
c.fillText(Amax.toFixed(1), pP.x0 - 8, pP.y0 + 6); c.fillText("0", pP.x0 - 8, pP.y0 + pP.h + 2);
return c.canvas;
}What you see. Top: the true continuous signal in grey, its samples as blue dots, and in orange the signal rebuilt from the \(K\) strongest DFT components. Bottom: the amplitude spectrum, one stem per DFT frequency \(f_m = m/T\), up to the Nyquist frequency marked by the red dashed line. Stems used in the reconstruction are orange, the rest blue. The green dashed lines mark where the two true frequencies \(f_1\) and \(f_2\) land in the spectrum; when a frequency lies above Nyquist, its green line sits at the aliased frequency instead.
Try: set \(f_2\) above the Nyquist frequency and watch its energy appear at the wrong place. Make \(T\) a non-integer multiple of the period and the peak smears over neighbouring frequencies (leakage). With \(K = 2\) and well-resolved peaks, the orange reconstruction sits exactly on the grey signal. From the Fourier exercises in Theoretical Neuroscience.