{"spec_id":"biplot-pca","library":"highcharts","language":"javascript","code":"// anyplot.ai\n// biplot-pca: PCA Biplot with Scores and Loading Vectors\n// Library: highcharts 12.6.0 | JavaScript 22.23.2\n// Quality: 93/100 | Created: 2026-09-01\n\n//# anyplot-orientation: square\nconst t = window.ANYPLOT_TOKENS;\n\n// --- Reproducible PRNG (LCG + Box-Muller) -----------------------------------\nlet seed = 42;\nfunction lcg() {\n  seed = (seed * 1103515245 + 12345) & 0x7fffffff;\n  return seed / 0x7fffffff;\n}\nfunction randNormal() {\n  const u1 = Math.max(lcg(), 1e-9);\n  const u2 = lcg();\n  return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);\n}\n\n// --- Data: synthetic wine-cultivar physicochemical measurements ------------\n// Six correlated features driven by two latent factors, three cultivar groups.\nconst featureNames = [\"Alcohol\", \"Malic Acid\", \"Ash\", \"Alkalinity\", \"Phenols\", \"Flavanoids\"];\nconst featureCoefs = [\n  { base: 13.0, c1: 0.9, c2: 0.15, noise: 0.35 },\n  { base: 2.4, c1: -0.7, c2: 0.25, noise: 0.4 },\n  { base: 2.3, c1: 0.05, c2: 0.35, noise: 0.25 },\n  { base: 19.5, c1: -0.55, c2: 0.6, noise: 1.2 },\n  { base: 2.3, c1: 0.8, c2: -0.15, noise: 0.3 },\n  { base: 2.0, c1: 0.85, c2: -0.2, noise: 0.3 },\n];\nconst groups = [\n  { name: \"Cultivar A\", latent1: 1.3, latent2: 0.4 },\n  { name: \"Cultivar B\", latent1: 0.0, latent2: -0.7 },\n  { name: \"Cultivar C\", latent1: -1.3, latent2: 0.3 },\n];\nconst nPerGroup = 20;\nconst nFeatures = featureNames.length;\n\nconst rawRows = [];\nconst groupIndex = [];\ngroups.forEach((g, gi) => {\n  for (let i = 0; i < nPerGroup; i++) {\n    const latent1 = g.latent1 + randNormal() * 0.5;\n    const latent2 = g.latent2 + randNormal() * 0.5;\n    const row = featureCoefs.map(\n      (f) => f.base + f.c1 * latent1 + f.c2 * latent2 + randNormal() * f.noise\n    );\n    rawRows.push(row);\n    groupIndex.push(gi);\n  }\n});\nconst nObs = rawRows.length;\n\n// --- Standardize (z-score) each feature — correlation-based PCA ------------\nconst means = featureCoefs.map((_, k) => rawRows.reduce((s, r) => s + r[k], 0) / nObs);\nconst stds = featureCoefs.map((_, k) => {\n  const variance = rawRows.reduce((s, r) => s + (r[k] - means[k]) ** 2, 0) / (nObs - 1);\n  return Math.sqrt(variance);\n});\nconst standardized = rawRows.map((row) => row.map((v, k) => (v - means[k]) / stds[k]));\n\n// --- Correlation matrix ------------------------------------------------------\nconst corr = Array.from({ length: nFeatures }, (_, i) =>\n  Array.from({ length: nFeatures }, (_, j) => {\n    let s = 0;\n    for (let o = 0; o < nObs; o++) s += standardized[o][i] * standardized[o][j];\n    return s / (nObs - 1);\n  })\n);\n\n// --- Jacobi eigenvalue decomposition (symmetric matrix) ---------------------\nfunction jacobiEigen(matrix, n) {\n  const A = matrix.map((row) => row.slice());\n  const V = Array.from({ length: n }, (_, i) =>\n    Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))\n  );\n  for (let sweep = 0; sweep < 100; sweep++) {\n    let off = 0;\n    for (let p = 0; p < n; p++) for (let q = p + 1; q < n; q++) off += A[p][q] * A[p][q];\n    if (off < 1e-12) break;\n    for (let p = 0; p < n; p++) {\n      for (let q = p + 1; q < n; q++) {\n        if (Math.abs(A[p][q]) < 1e-14) continue;\n        const theta = (A[q][q] - A[p][p]) / (2 * A[p][q]);\n        const sign = theta >= 0 ? 1 : -1;\n        const tVal = sign / (Math.abs(theta) + Math.sqrt(theta * theta + 1));\n        const c = 1 / Math.sqrt(tVal * tVal + 1);\n        const s = tVal * c;\n        const app = A[p][p];\n        const aqq = A[q][q];\n        const apq = A[p][q];\n        A[p][p] = c * c * app - 2 * s * c * apq + s * s * aqq;\n        A[q][q] = s * s * app + 2 * s * c * apq + c * c * aqq;\n        A[p][q] = 0;\n        A[q][p] = 0;\n        for (let i = 0; i < n; i++) {\n          if (i !== p && i !== q) {\n            const aip = A[i][p];\n            const aiq = A[i][q];\n            A[i][p] = c * aip - s * aiq;\n            A[p][i] = A[i][p];\n            A[i][q] = s * aip + c * aiq;\n            A[q][i] = A[i][q];\n          }\n        }\n        for (let i = 0; i < n; i++) {\n          const vip = V[i][p];\n          const viq = V[i][q];\n          V[i][p] = c * vip - s * viq;\n          V[i][q] = s * vip + c * viq;\n        }\n      }\n    }\n  }\n  const values = Array.from({ length: n }, (_, i) => A[i][i]);\n  return { values, vectors: V };\n}\n\nconst { values: eigenvalues, vectors: eigenvectors } = jacobiEigen(corr, nFeatures);\nconst order = eigenvalues.map((_, i) => i).sort((a, b) => eigenvalues[b] - eigenvalues[a]);\nconst eigen1 = eigenvalues[order[0]];\nconst eigen2 = eigenvalues[order[1]];\nconst pc1Vec = eigenvectors.map((row) => row[order[0]]);\nconst pc2Vec = eigenvectors.map((row) => row[order[1]]);\nconst totalVariance = eigenvalues.reduce((s, v) => s + v, 0);\nconst pc1Pct = (eigen1 / totalVariance) * 100;\nconst pc2Pct = (eigen2 / totalVariance) * 100;\n\n// --- Scores (observations projected onto PC1/PC2) --------------------------\nconst scores = standardized.map((row) => [\n  row.reduce((s, v, k) => s + v * pc1Vec[k], 0),\n  row.reduce((s, v, k) => s + v * pc2Vec[k], 0),\n]);\n\n// --- Loadings (correlation between variable and component) -----------------\nconst loadings = featureNames.map((name, k) => ({\n  name,\n  x: pc1Vec[k] * Math.sqrt(eigen1),\n  y: pc2Vec[k] * Math.sqrt(eigen2),\n}));\n\n// --- Scale loadings so arrows read alongside the score cloud ---------------\nconst maxScoreAbs = Math.max(...scores.flat().map(Math.abs));\nconst maxLoadingMag = Math.max(...loadings.map((l) => Math.hypot(l.x, l.y)));\nconst loadingScale = (maxScoreAbs * 0.85) / maxLoadingMag;\nconst scaledLoadings = loadings.map((l) => ({\n  name: l.name,\n  x: l.x * loadingScale,\n  y: l.y * loadingScale,\n}));\nconst axisMax = Math.max(maxScoreAbs, loadingScale) * 1.2;\n\n// --- Dominant loading + its most-aligned cultivar (interpretive callout) ---\nconst loadingMagnitudes = loadings.map((l) => Math.hypot(l.x, l.y));\nconst dominantIdx = loadingMagnitudes.indexOf(Math.max(...loadingMagnitudes));\nconst dominantUnit = {\n  x: loadings[dominantIdx].x / loadingMagnitudes[dominantIdx],\n  y: loadings[dominantIdx].y / loadingMagnitudes[dominantIdx],\n};\nconst groupCentroids = groups.map((g, gi) => {\n  const groupScores = scores.filter((_, i) => groupIndex[i] === gi);\n  const n = groupScores.length;\n  return {\n    name: g.name,\n    meanX: groupScores.reduce((s, r) => s + r[0], 0) / n,\n    meanY: groupScores.reduce((s, r) => s + r[1], 0) / n,\n  };\n});\nconst alignedGroup = groupCentroids.reduce((best, c) => {\n  const proj = c.meanX * dominantUnit.x + c.meanY * dominantUnit.y;\n  return proj > best.proj ? { name: c.name, proj } : best;\n}, { name: null, proj: -Infinity });\n\n// --- Group score series ------------------------------------------------------\nconst seriesData = groups.map((g, gi) => ({\n  name: g.name,\n  type: \"scatter\",\n  color: t.palette[gi],\n  data: scores.filter((_, i) => groupIndex[i] === gi),\n  marker: { radius: 6, symbol: \"circle\" },\n}));\n\n// --- Chart -------------------------------------------------------------------\nHighcharts.chart(\"container\", {\n  chart: {\n    type: \"scatter\",\n    backgroundColor: \"transparent\",\n    animation: false,\n    style: { fontFamily: \"inherit\" },\n    events: {\n      load: function () {\n        const chart = this;\n        const xAxis = chart.xAxis[0];\n        const yAxis = chart.yAxis[0];\n        const renderer = chart.renderer;\n\n        // Reference unit circle (correlation biplot scaling)\n        const circlePoints = 72;\n        const circlePath = [];\n        for (let i = 0; i <= circlePoints; i++) {\n          const angle = (i / circlePoints) * 2 * Math.PI;\n          const px = xAxis.toPixels(loadingScale * Math.cos(angle), false);\n          const py = yAxis.toPixels(loadingScale * Math.sin(angle), false);\n          circlePath.push(i === 0 ? \"M\" : \"L\", px, py);\n        }\n        circlePath.push(\"Z\");\n        renderer\n          .path(circlePath)\n          .attr({ stroke: t.inkSoft, \"stroke-width\": 1, \"stroke-dasharray\": \"4,4\", fill: \"none\", opacity: 0.5 })\n          .add();\n\n        // Loading vectors — arrows drawn from the origin\n        const originX = xAxis.toPixels(0, false);\n        const originY = yAxis.toPixels(0, false);\n\n        // Precompute label anchors, then resolve vertical collisions per side\n        // (right-side labels otherwise crowd together at mobile width).\n        const labelAnchors = scaledLoadings.map((l) => ({\n          x: xAxis.toPixels(l.x * 1.12, false),\n          y: yAxis.toPixels(l.y * 1.12, false),\n          align: l.x >= 0 ? \"left\" : \"right\",\n        }));\n        const minGap = 16;\n        [\"left\", \"right\"].forEach((side) => {\n          const group = labelAnchors\n            .map((a, idx) => ({ ...a, idx }))\n            .filter((a) => a.align === side)\n            .sort((a, b) => a.y - b.y);\n          for (let i = 1; i < group.length; i++) {\n            if (group[i].y - group[i - 1].y < minGap) {\n              group[i].y = group[i - 1].y + minGap;\n            }\n            labelAnchors[group[i].idx].y = group[i].y;\n          }\n        });\n\n        scaledLoadings.forEach((l, idx) => {\n          const isDominant = idx === dominantIdx;\n          const strokeColor = isDominant ? t.amber : t.ink;\n          const strokeWidth = isDominant ? 3 : 2;\n\n          const tipX = xAxis.toPixels(l.x, false);\n          const tipY = yAxis.toPixels(l.y, false);\n          renderer\n            .path([\"M\", originX, originY, \"L\", tipX, tipY])\n            .attr({ stroke: strokeColor, \"stroke-width\": strokeWidth })\n            .add();\n\n          const angle = Math.atan2(tipY - originY, tipX - originX);\n          const headLen = isDominant ? 12 : 10;\n          const headAngle = 0.45;\n          const h1x = tipX - headLen * Math.cos(angle - headAngle);\n          const h1y = tipY - headLen * Math.sin(angle - headAngle);\n          const h2x = tipX - headLen * Math.cos(angle + headAngle);\n          const h2y = tipY - headLen * Math.sin(angle + headAngle);\n          renderer\n            .path([\"M\", h1x, h1y, \"L\", tipX, tipY, \"L\", h2x, h2y])\n            .attr({ stroke: strokeColor, \"stroke-width\": strokeWidth, fill: \"none\" })\n            .add();\n\n          const anchor = labelAnchors[idx];\n          renderer\n            .text(l.name, anchor.x, anchor.y)\n            .attr({ align: anchor.align })\n            .css({\n              color: isDominant ? t.amber : t.ink,\n              fontSize: \"13px\",\n              fontWeight: isDominant ? \"700\" : \"600\",\n            })\n            .add();\n        });\n      },\n    },\n  },\n  credits: { enabled: false },\n  colors: t.palette,\n  title: {\n    text: \"biplot-pca · javascript · highcharts · anyplot.ai\",\n    style: { color: t.ink, fontSize: \"22px\", fontWeight: \"600\" },\n  },\n  subtitle: {\n    text: \"Points: standardized wine-cultivar scores · Arrows: variable loadings (scaled) · Dashed: unit circle\",\n    style: { color: t.inkSoft, fontSize: \"14px\" },\n  },\n  caption: {\n    text:\n      \"Strongest signal: \" +\n      featureNames[dominantIdx] +\n      \" loads most heavily, aligning most closely with \" +\n      alignedGroup.name,\n    align: \"left\",\n    style: { color: t.amber, fontSize: \"13px\", fontStyle: \"italic\" },\n  },\n  xAxis: {\n    title: {\n      text: \"PC1 (\" + pc1Pct.toFixed(1) + \"%)\",\n      style: { color: t.inkSoft, fontSize: \"16px\" },\n    },\n    min: -axisMax,\n    max: axisMax,\n    lineWidth: 0,\n    tickLength: 0,\n    gridLineColor: t.grid,\n    gridLineWidth: 1,\n    gridLineDashStyle: \"Dot\",\n    labels: { style: { color: t.inkSoft, fontSize: \"14px\" } },\n    plotLines: [{ value: 0, color: t.inkSoft, width: 1.5, zIndex: 1 }],\n  },\n  yAxis: {\n    title: {\n      text: \"PC2 (\" + pc2Pct.toFixed(1) + \"%)\",\n      style: { color: t.inkSoft, fontSize: \"16px\" },\n    },\n    min: -axisMax,\n    max: axisMax,\n    lineWidth: 0,\n    tickLength: 0,\n    gridLineColor: t.grid,\n    gridLineWidth: 1,\n    gridLineDashStyle: \"Dot\",\n    labels: { style: { color: t.inkSoft, fontSize: \"14px\" } },\n    plotLines: [{ value: 0, color: t.inkSoft, width: 1.5, zIndex: 1 }],\n  },\n  legend: {\n    enabled: true,\n    verticalAlign: \"bottom\",\n    itemStyle: { color: t.inkSoft, fontSize: \"14px\" },\n    itemHoverStyle: { color: t.ink },\n  },\n  tooltip: { enabled: false },\n  plotOptions: {\n    scatter: {\n      animation: false,\n      states: { hover: { enabled: false } },\n    },\n    series: { animation: false },\n  },\n  series: seriesData,\n});\n"}