{"spec_id":"biplot-pca","library":"makie","language":"julia","code":"# anyplot.ai\n# biplot-pca: PCA Biplot with Scores and Loading Vectors\n# Library: makie 0.21.9 | Julia 1.11.9\n# Quality: 88/100 | Created: 2026-09-01\n\nusing CairoMakie\nusing Colors\nusing RDatasets\nusing Statistics\nusing LinearAlgebra\n\n# --- Theme tokens ------------------------------------------------------------\nconst THEME       = get(ENV, \"ANYPLOT_THEME\", \"light\")\nconst PAGE_BG     = THEME == \"light\" ? colorant\"#FAF8F1\" : colorant\"#1A1A17\"\nconst INK         = THEME == \"light\" ? colorant\"#1A1A17\" : colorant\"#F0EFE8\"\nconst INK_SOFT    = THEME == \"light\" ? colorant\"#4A4A44\" : colorant\"#B8B7B0\"\nconst IMPRINT_PALETTE = [\n    colorant\"#009E73\", colorant\"#C475FD\", colorant\"#4467A3\", colorant\"#BD8233\",\n    colorant\"#AE3030\", colorant\"#2ABCCD\", colorant\"#954477\", colorant\"#99B314\",\n]\n\n# --- Data ---------------------------------------------------------------------\niris = RDatasets.dataset(\"datasets\", \"iris\")\nfeature_names = [\"Sepal Length\", \"Sepal Width\", \"Petal Length\", \"Petal Width\"]\nfeatures = Matrix(iris[:, [:SepalLength, :SepalWidth, :PetalLength, :PetalWidth]])\nspecies = String.(iris.Species)\ngroups = unique(species)\n\n# Standardize so PCA decomposes the correlation matrix, not the covariance matrix\nfeature_means = mean(features; dims = 1)\nfeature_stds = std(features; dims = 1)\nstandardized = (features .- feature_means) ./ feature_stds\n\n# Eigen-decompose the correlation matrix, sorted by descending eigenvalue\ncorrelation = Symmetric(cor(standardized))\neigen_result = eigen(correlation)\norder = sortperm(eigen_result.values; rev = true)\neigenvalues = eigen_result.values[order]\neigenvectors = eigen_result.vectors[:, order]\n\nvariance_explained = eigenvalues ./ sum(eigenvalues) .* 100\nscores = standardized * eigenvectors[:, 1:2]\nloadings = eigenvectors[:, 1:2] .* sqrt.(eigenvalues[1:2])'\n\n# Scale correlation loadings so the arrows sit alongside the score cloud\nscore_radius = maximum(sqrt.(scores[:, 1] .^ 2 .+ scores[:, 2] .^ 2))\nloading_radius = maximum(sqrt.(loadings[:, 1] .^ 2 .+ loadings[:, 2] .^ 2))\narrow_scale = 0.85 * score_radius / loading_radius\narrow_xy = loadings .* arrow_scale\n\n# --- Plot -----------------------------------------------------------------\nfig = Figure(\n    resolution      = (1200, 1200),\n    fontsize        = 14,\n    backgroundcolor = PAGE_BG,\n    figure_padding  = (10, 10, 10, 10),\n)\n\npc1_label = \"PC1 ($(round(variance_explained[1]; digits = 1))%)\"\npc2_label = \"PC2 ($(round(variance_explained[2]; digits = 1))%)\"\n\nax = Axis(\n    fig[1, 1];\n    title              = \"biplot-pca · julia · makie · anyplot.ai\",\n    titlesize          = 20,\n    titlecolor         = INK,\n    xlabel             = pc1_label,\n    ylabel             = pc2_label,\n    xlabelcolor        = INK,\n    ylabelcolor        = INK,\n    xlabelsize         = 14,\n    ylabelsize         = 14,\n    xticklabelcolor    = INK_SOFT,\n    yticklabelcolor    = INK_SOFT,\n    xticklabelsize     = 12,\n    yticklabelsize     = 12,\n    backgroundcolor    = PAGE_BG,\n    topspinevisible    = false,\n    rightspinevisible  = false,\n    leftspinecolor     = INK_SOFT,\n    bottomspinecolor   = INK_SOFT,\n    xgridcolor         = RGBAf(INK.r, INK.g, INK.b, 0.15),\n    ygridcolor         = RGBAf(INK.r, INK.g, INK.b, 0.15),\n    xautolimitmargin   = (0.03, 0.03),\n    yautolimitmargin   = (0.03, 0.03),\n    aspect             = DataAspect(),\n)\n\n# Unit circle reference for the correlation-scaled loadings\ncircle_theta = range(0, 2π; length = 200)\nlines!(ax, cos.(circle_theta) .* arrow_scale, sin.(circle_theta) .* arrow_scale;\n    color = INK_SOFT, linestyle = :dash, linewidth = 1.5)\n\nhlines!(ax, [0]; color = RGBAf(INK.r, INK.g, INK.b, 0.3), linewidth = 1)\nvlines!(ax, [0]; color = RGBAf(INK.r, INK.g, INK.b, 0.3), linewidth = 1)\n\n# Subtle per-group 1.5σ density ellipses give the species separation a\n# deliberate focal point instead of relying on point color alone\nellipse_theta = range(0, 2π; length = 100)\nfor (i, group) in enumerate(groups)\n    mask = species .== group\n    group_scores = scores[mask, :]\n    center = vec(mean(group_scores; dims = 1))\n    group_evals, group_evecs = eigen(Symmetric(cov(group_scores)))\n    group_order = sortperm(group_evals; rev = true)\n    radii = 1.5 .* sqrt.(max.(group_evals[group_order], 0))\n    circle_pts = radii .* permutedims(hcat(cos.(ellipse_theta), sin.(ellipse_theta)))\n    ellipse_pts = group_evecs[:, group_order] * circle_pts\n    poly!(ax, Point2f.(ellipse_pts[1, :] .+ center[1], ellipse_pts[2, :] .+ center[2]);\n        color = (IMPRINT_PALETTE[i], 0.10), strokecolor = (IMPRINT_PALETTE[i], 0.35), strokewidth = 1)\nend\n\nfor (i, group) in enumerate(groups)\n    mask = species .== group\n    scatter!(ax, scores[mask, 1], scores[mask, 2];\n        color = IMPRINT_PALETTE[i], markersize = 12, strokewidth = 0, label = group)\nend\n\n# Heavier line weight for the dominant loadings makes the strongest drivers\n# of variance the clear focal point rather than four uniform arrows\nloading_magnitude = sqrt.(arrow_xy[:, 1] .^ 2 .+ arrow_xy[:, 2] .^ 2)\narrow_linewidth = 1.8 .+ 1.8 .* (loading_magnitude ./ maximum(loading_magnitude))\n\narrows!(ax, zeros(length(feature_names)), zeros(length(feature_names)),\n    arrow_xy[:, 1], arrow_xy[:, 2];\n    color = INK, linewidth = arrow_linewidth, arrowsize = 18)\n\n# Nudge labels apart vertically when their arrow tips sit too close to read\nlabel_offset = zeros(length(feature_names), 2)\nfor i in 1:length(feature_names), j in (i + 1):length(feature_names)\n    tip_distance = hypot(arrow_xy[i, 1] - arrow_xy[j, 1], arrow_xy[i, 2] - arrow_xy[j, 2])\n    if tip_distance < 0.4\n        label_offset[i, 2] += 0.35\n        label_offset[j, 2] -= 0.35\n    end\nend\n\n# Anchor each label on whichever horizontal side has more clearance from the\n# score cloud, so the text itself never runs into a nearby point cluster\n# (a fixed left-anchor would run \"Sepal Width\" straight into the setosa\n# points sitting just beyond its arrow tip)\ntext_halfwidth = [0.025 * length(name) for name in feature_names]\nh_align = Vector{Symbol}(undef, length(feature_names))\nfor i in 1:length(feature_names)\n    tip = arrow_xy[i, :] .+ label_offset[i, :]\n    right_center = tip .+ [text_halfwidth[i], 0]\n    left_center = tip .- [text_halfwidth[i], 0]\n    dist_right = minimum(hypot.(scores[:, 1] .- right_center[1], scores[:, 2] .- right_center[2]))\n    dist_left = minimum(hypot.(scores[:, 1] .- left_center[1], scores[:, 2] .- left_center[2]))\n    h_align[i] = dist_right >= dist_left ? :left : :right\nend\n\nfor (i, name) in enumerate(feature_names)\n    text!(ax, arrow_xy[i, 1] + label_offset[i, 1], arrow_xy[i, 2] + label_offset[i, 2]; text = name,\n        color = INK, fontsize = 13, align = (h_align[i], :bottom))\nend\n\nLegend(fig[1, 2], ax, \"Species\";\n    framevisible = false, labelcolor = INK, titlecolor = INK)\ncolsize!(fig.layout, 2, Relative(0.15))\ncolgap!(fig.layout, 1, 8)\n\n# --- Save -------------------------------------------------------------------\nsave(\"plot-$(THEME).png\", fig; px_per_unit = 2)\n"}