{"spec_id":"biplot-pca","library":"d3","language":"javascript","code":"// anyplot.ai\n// biplot-pca: PCA Biplot with Scores and Loading Vectors\n// Library: d3 7.9.0 | JavaScript 22.23.2\n// Quality: 90/100 | Created: 2026-09-01\n//# anyplot-orientation: square\n\nconst t = window.ANYPLOT_TOKENS;\nconst { width, height } = window.ANYPLOT_SIZE;\n\n// --- Data: synthetic plant morphology survey, 3 growth-habit groups -------\nfunction lcg(seed) {\n  let state = seed >>> 0;\n  return () => {\n    state = (state * 1664525 + 1013904223) >>> 0;\n    return state / 4294967296;\n  };\n}\nconst rand = lcg(42);\nfunction gaussian() {\n  let u1 = 0;\n  while (u1 === 0) u1 = rand();\n  const u2 = rand();\n  return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);\n}\n\n// Each feature loads on the two latent factors (size, shape) at its own angle,\n// so the six loading arrows fan out around the circle instead of bunching up.\nconst featureSpecs = [\n  { name: \"Sepal Length\", mean: 5.0, spread: 0.7, angleDeg: 10 },\n  { name: \"Petal Length\", mean: 3.5, spread: 0.9, angleDeg: 55 },\n  { name: \"Petal Width\", mean: 1.2, spread: 0.35, angleDeg: 95 },\n  { name: \"Sepal Width\", mean: 3.0, spread: 0.4, angleDeg: 150 },\n  { name: \"Leaf Area\", mean: 15, spread: 3.5, angleDeg: 205 },\n  { name: \"Stem Height\", mean: 40, spread: 7, angleDeg: 320 },\n];\nconst featureNames = featureSpecs.map((f) => f.name);\nconst groups = [\n  { name: \"Compact\", sizeMean: -1.3 },\n  { name: \"Standard\", sizeMean: 0 },\n  { name: \"Vigorous\", sizeMean: 1.3 },\n];\n\nconst observations = [];\nfor (const grp of groups) {\n  for (let i = 0; i < 30; i++) {\n    const size = grp.sizeMean + gaussian() * 0.5;\n    const shape = gaussian() * 0.5;\n    observations.push({\n      group: grp.name,\n      values: featureSpecs.map((f) => {\n        const angleRad = (f.angleDeg * Math.PI) / 180;\n        const signal = Math.cos(angleRad) * size + Math.sin(angleRad) * shape;\n        return f.mean + f.spread * 0.9 * signal + gaussian() * f.spread * 0.45;\n      }),\n    });\n  }\n}\n\n// --- Standardize columns, then PCA on the correlation matrix ---------------\nconst n = observations.length;\nconst p = featureNames.length;\n\nconst means = Array(p).fill(0);\nfor (const obs of observations) obs.values.forEach((v, j) => (means[j] += v / n));\n\nconst stds = Array(p).fill(0);\nfor (const obs of observations) obs.values.forEach((v, j) => (stds[j] += (v - means[j]) ** 2 / (n - 1)));\nstds.forEach((s, j) => (stds[j] = Math.sqrt(s)));\n\nconst standardized = observations.map((obs) => obs.values.map((v, j) => (v - means[j]) / stds[j]));\n\nconst corr = Array.from({ length: p }, () => Array(p).fill(0));\nfor (let a = 0; a < p; a++) {\n  for (let b = 0; b < p; b++) {\n    let sum = 0;\n    for (let i = 0; i < n; i++) sum += standardized[i][a] * standardized[i][b];\n    corr[a][b] = sum / (n - 1);\n  }\n}\n\n// Cyclic Jacobi eigenvalue algorithm for the symmetric correlation matrix\nfunction jacobiEigen(matrix) {\n  const dim = matrix.length;\n  const a = matrix.map((row) => row.slice());\n  const v = Array.from({ length: dim }, (_, i) => Array.from({ length: dim }, (_, j) => (i === j ? 1 : 0)));\n  for (let sweep = 0; sweep < 100; sweep++) {\n    let offDiag = 0;\n    for (let i = 0; i < dim; i++) for (let j = i + 1; j < dim; j++) offDiag += a[i][j] * a[i][j];\n    if (offDiag < 1e-12) break;\n    for (let pi = 0; pi < dim; pi++) {\n      for (let qi = pi + 1; qi < dim; qi++) {\n        if (Math.abs(a[pi][qi]) < 1e-14) continue;\n        const theta = (a[qi][qi] - a[pi][pi]) / (2 * a[pi][qi]);\n        const tt = Math.sign(theta || 1) / (Math.abs(theta) + Math.sqrt(theta * theta + 1));\n        const c = 1 / Math.sqrt(tt * tt + 1);\n        const s = tt * c;\n        const app = a[pi][pi];\n        const aqq = a[qi][qi];\n        const apq = a[pi][qi];\n        a[pi][pi] = c * c * app - 2 * s * c * apq + s * s * aqq;\n        a[qi][qi] = s * s * app + 2 * s * c * apq + c * c * aqq;\n        a[pi][qi] = 0;\n        a[qi][pi] = 0;\n        for (let i = 0; i < dim; i++) {\n          if (i === pi || i === qi) continue;\n          const aip = a[i][pi];\n          const aiq = a[i][qi];\n          a[i][pi] = c * aip - s * aiq;\n          a[pi][i] = a[i][pi];\n          a[i][qi] = s * aip + c * aiq;\n          a[qi][i] = a[i][qi];\n        }\n        for (let i = 0; i < dim; i++) {\n          const vip = v[i][pi];\n          const viq = v[i][qi];\n          v[i][pi] = c * vip - s * viq;\n          v[i][qi] = s * vip + c * viq;\n        }\n      }\n    }\n  }\n  return { eigenvalues: Array.from({ length: dim }, (_, i) => a[i][i]), eigenvectors: v };\n}\n\nconst { eigenvalues, eigenvectors } = jacobiEigen(corr);\nconst order = eigenvalues.map((_, idx) => idx).sort((a, b) => eigenvalues[b] - eigenvalues[a]);\nconst [pc1Idx, pc2Idx] = order;\nconst totalVariance = eigenvalues.reduce((sum, val) => sum + val, 0);\nconst pc1Pct = (eigenvalues[pc1Idx] / totalVariance) * 100;\nconst pc2Pct = (eigenvalues[pc2Idx] / totalVariance) * 100;\n\nconst scores = observations.map((obs, i) => {\n  let pc1 = 0;\n  let pc2 = 0;\n  for (let j = 0; j < p; j++) {\n    pc1 += standardized[i][j] * eigenvectors[j][pc1Idx];\n    pc2 += standardized[i][j] * eigenvectors[j][pc2Idx];\n  }\n  return { group: obs.group, pc1, pc2 };\n});\n\n// Correlation-biplot loadings: variable-PC correlation, magnitude <= 1\nconst loadings = featureNames.map((name, j) => ({\n  name,\n  pc1: eigenvectors[j][pc1Idx] * Math.sqrt(eigenvalues[pc1Idx]),\n  pc2: eigenvectors[j][pc2Idx] * Math.sqrt(eigenvalues[pc2Idx]),\n}));\n\n// The feature whose loading correlates most strongly with PC1 is the single\n// best explanation for the group separation visible along that axis — call\n// it out visually instead of leaving every arrow at the same default weight.\nconst topDriver = loadings.reduce((best, d) => (Math.abs(d.pc1) > Math.abs(best.pc1) ? d : best), loadings[0]);\n\nconst maxScoreRadius = d3.max(scores, (d) => Math.hypot(d.pc1, d.pc2));\nconst arrowScale = 0.85 * maxScoreRadius;\nconst extent = maxScoreRadius * 1.3;\n\n// --- Layout: equal-aspect plot area so loading-vector angles read true -----\nconst margin = { top: 160, right: 190, bottom: 110, left: 130 };\nconst availW = width - margin.left - margin.right;\nconst availH = height - margin.top - margin.bottom;\nconst plotSize = Math.min(availW, availH);\nconst xOffset = margin.left + (availW - plotSize) / 2;\nconst yOffset = margin.top + (availH - plotSize) / 2;\n\nconst x = d3.scaleLinear().domain([-extent, extent]).range([xOffset, xOffset + plotSize]);\nconst y = d3.scaleLinear().domain([-extent, extent]).range([yOffset + plotSize, yOffset]);\nconst color = d3.scaleOrdinal().domain(groups.map((g) => g.name)).range(t.palette.slice(0, groups.length));\n\n// --- SVG mount ---------------------------------------------------------------\nconst svg = d3.select(\"#container\").append(\"svg\").attr(\"width\", width).attr(\"height\", height);\n\nsvg\n  .append(\"defs\")\n  .append(\"marker\")\n  .attr(\"id\", \"loading-arrowhead\")\n  .attr(\"viewBox\", \"0 0 10 10\")\n  .attr(\"refX\", 8)\n  .attr(\"refY\", 5)\n  .attr(\"markerWidth\", 7)\n  .attr(\"markerHeight\", 7)\n  .attr(\"orient\", \"auto-start-reverse\")\n  .append(\"path\")\n  .attr(\"d\", \"M 0 0 L 10 5 L 0 10 z\")\n  .attr(\"fill\", t.ink);\n\n// --- Unit circle: reference for correlation-loading magnitude --------------\nsvg\n  .append(\"circle\")\n  .attr(\"cx\", x(0))\n  .attr(\"cy\", y(0))\n  .attr(\"r\", x(arrowScale) - x(0))\n  .attr(\"fill\", \"none\")\n  .attr(\"stroke\", t.grid)\n  .attr(\"stroke-width\", 1.5)\n  .attr(\"stroke-dasharray\", \"6,5\");\n\n// --- Axes ---------------------------------------------------------------------\nconst xAxisG = svg.append(\"g\").attr(\"transform\", `translate(0,${yOffset + plotSize})`).call(d3.axisBottom(x).ticks(6));\nconst yAxisG = svg.append(\"g\").attr(\"transform\", `translate(${xOffset},0)`).call(d3.axisLeft(y).ticks(6));\nfor (const axisG of [xAxisG, yAxisG]) {\n  axisG.selectAll(\"text\").attr(\"fill\", t.inkSoft).style(\"font-size\", \"14px\");\n  axisG.selectAll(\"line\").attr(\"stroke\", t.grid);\n  axisG.select(\".domain\").attr(\"stroke\", t.inkSoft);\n}\n\nsvg\n  .append(\"text\")\n  .attr(\"x\", xOffset + plotSize / 2)\n  .attr(\"y\", yOffset + plotSize + 70)\n  .attr(\"text-anchor\", \"middle\")\n  .attr(\"fill\", t.ink)\n  .style(\"font-size\", \"17px\")\n  .text(`PC1 (${pc1Pct.toFixed(1)}%)`);\n\nsvg\n  .append(\"text\")\n  .attr(\"transform\", `translate(${xOffset - 80},${yOffset + plotSize / 2}) rotate(-90)`)\n  .attr(\"text-anchor\", \"middle\")\n  .attr(\"fill\", t.ink)\n  .style(\"font-size\", \"17px\")\n  .text(`PC2 (${pc2Pct.toFixed(1)}%)`);\n\n// --- Observation scores (drawn first so loading arrows/labels stay on top) -----\nsvg\n  .selectAll(\".score-point\")\n  .data(scores)\n  .join(\"circle\")\n  .attr(\"cx\", (d) => x(d.pc1))\n  .attr(\"cy\", (d) => y(d.pc2))\n  .attr(\"r\", 6)\n  .attr(\"fill\", (d) => color(d.group))\n  .attr(\"fill-opacity\", 0.65)\n  .attr(\"stroke\", t.pageBg)\n  .attr(\"stroke-width\", 1);\n\n// --- Loading arrows + labels ---------------------------------------------------\nsvg\n  .selectAll(\".loading-arrow\")\n  .data(loadings)\n  .join(\"line\")\n  .attr(\"x1\", x(0))\n  .attr(\"y1\", y(0))\n  .attr(\"x2\", (d) => x(d.pc1 * arrowScale))\n  .attr(\"y2\", (d) => y(d.pc2 * arrowScale))\n  .attr(\"stroke\", t.ink)\n  .attr(\"stroke-width\", (d) => (d.name === topDriver.name ? 3.5 : 2.5))\n  .attr(\"marker-end\", \"url(#loading-arrowhead)\");\n\nsvg\n  .selectAll(\".loading-label\")\n  .data(loadings)\n  .join(\"text\")\n  .attr(\"x\", (d) => x(d.pc1 * arrowScale * 1.2))\n  .attr(\"y\", (d) => y(d.pc2 * arrowScale * 1.2))\n  .attr(\"text-anchor\", (d) => (d.pc1 >= 0 ? \"start\" : \"end\"))\n  .attr(\"dominant-baseline\", (d) => (d.pc2 >= 0 ? \"auto\" : \"hanging\"))\n  .attr(\"fill\", t.ink)\n  .style(\"font-size\", (d) => (d.name === topDriver.name ? \"17px\" : \"15px\"))\n  .style(\"font-weight\", (d) => (d.name === topDriver.name ? \"700\" : \"600\"))\n  .style(\"paint-order\", \"stroke\")\n  .attr(\"stroke\", t.pageBg)\n  .attr(\"stroke-width\", 6.5)\n  .text((d) => d.name);\n\n// --- Legend -----------------------------------------------------------------------\nconst legendX = xOffset + plotSize + 30;\nconst legendY = yOffset + plotSize / 2 - (groups.length * 44) / 2;\nsvg\n  .append(\"text\")\n  .attr(\"x\", legendX)\n  .attr(\"y\", legendY - 22)\n  .attr(\"fill\", t.inkSoft)\n  .style(\"font-size\", \"14px\")\n  .style(\"font-weight\", \"600\")\n  .text(\"Growth habit\");\ngroups.forEach((grp, i) => {\n  const rowY = legendY + i * 44;\n  svg.append(\"rect\").attr(\"x\", legendX).attr(\"y\", rowY).attr(\"width\", 22).attr(\"height\", 22).attr(\"fill\", color(grp.name));\n  svg\n    .append(\"text\")\n    .attr(\"x\", legendX + 32)\n    .attr(\"y\", rowY + 16)\n    .attr(\"fill\", t.ink)\n    .style(\"font-size\", \"16px\")\n    .text(grp.name);\n});\n\n// --- Title --------------------------------------------------------------------------\nsvg\n  .append(\"text\")\n  .attr(\"x\", width / 2)\n  .attr(\"y\", 60)\n  .attr(\"text-anchor\", \"middle\")\n  .attr(\"fill\", t.ink)\n  .style(\"font-size\", \"26px\")\n  .style(\"font-weight\", \"600\")\n  .text(\"biplot-pca · javascript · d3 · anyplot.ai\");\n\n// Subtitle: call out the strongest driver of the PC1 separation seen in the\n// scores below, so the story is more than \"points happen to cluster by color\".\nsvg\n  .append(\"text\")\n  .attr(\"x\", width / 2)\n  .attr(\"y\", 92)\n  .attr(\"text-anchor\", \"middle\")\n  .attr(\"fill\", t.inkSoft)\n  .style(\"font-size\", \"15px\")\n  .style(\"font-style\", \"italic\")\n  .text(`${topDriver.name} loads most strongly on PC1 — the axis separating growth-habit groups`);\n"}