SPB Git

spb/metrika Public

Stata-class statistics, GPU-accelerated by Apple Silicon. Native Swift — no Electron, no Python runtime, no compromises.

macosswiftstatisticsgpu
Swift 92.4% HTML 3.3% R 3% Shell 1.3%

feat(stats): permutation tests and single-variable graphics

- ZQResampling.permutationIndices: permutation as argsort of 64-bit
  Philox keys — order-free and counter-addressable so a future GPU
  argsort path reproduces it exactly; dedicated stream offset (bit 63,
  applied after word expansion) keeps permutations disjoint from
  bootstrap draws under the same seed
- engine: 'permute, reps(#) [seed(#)]: reg …' — permutes the response
  within the estimation sample, parallel chunked refits, empirical
  two-sided p per coefficient (c, reps, p, SE(p)); degenerate _cons row
  omitted
- planner: GPU heuristic restricted to bootstrap until the permute
  argsort path lands
- graphics: standalone histogram/scatter/kdensity verbs; single-variable
  histogram (Sturges default, bins() option) and Epanechnikov kdensity
  with Silverman bandwidth (density integrates to 1 in tests)
- 81 tests green (swift test and xcodebuild with GPU suites)

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
simon-pierre boucher committed 6 days ago (Aug 5, 2026) parent de63232

Showing 6 changed files with +365 and −3

modified MetrikaKit/Sources/ZQEngine/Session.swift +220 −1
@@ -137,7 +137,13 @@ public actor ZQSession {
137 137 case "xtset": return try handleXTSet(command)
138 138 case "display": return try handleDisplay(command)
139 139 case "bootstrap": return try await handleBootstrap(command, backend: backend)
140 + case "permute": return try await handlePermute(command)
140 141 case "graph": return try handleGraph(command)
142 + case "histogram", "scatter", "kdensity":
143 + var promoted = command
144 + promoted.subverb = command.verb
145 + promoted.verb = "graph"
146 + return try handleGraph(promoted)
141 147 case "log": return try handleLog(command)
142 148 default:
143 149 throw ZQEngineError("command '\(command.verb)' is not implemented yet")
@@ -1399,6 +1405,128 @@ public actor ZQSession {
1399 1405 return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
1400 1406 }
1401 1407
1408 + // MARK: - Permutation test
1409 +
1410 + /// `permute, reps(#) [seed(#)]: reg y x…` — permutes the response
1411 + /// within the estimation sample, refits, and reports the empirical
1412 + /// two-sided p-value per coefficient: the fraction of permuted |β*|
1413 + /// at or above the observed |β|. Permutations come from the Philox
1414 + /// argsort stream, so results are exactly reproducible and
1415 + /// replicate-addressable.
1416 + private func handlePermute(_ command: ZQCommand) async throws -> ZQResult {
1417 + guard let body = command.body else {
1418 + throw ZQEngineError("permute: syntax is 'permute, reps(#): command'")
1419 + }
1420 + guard body.verb == "regress" else {
1421 + throw ZQEngineError("permute currently supports 'regress' bodies only")
1422 + }
1423 + guard let repsText = command.option("reps")?.firstArgument,
1424 + let reps = Int(repsText), reps >= 1 else {
1425 + throw ZQEngineError("permute: reps(#) with # ≥ 1 required")
1426 + }
1427 + if let seedText = command.option("seed")?.firstArgument {
1428 + guard let value = UInt64(seedText) else {
1429 + throw ZQEngineError("permute: invalid seed")
1430 + }
1431 + seed = value
1432 + }
1433 +
1434 + let sample = try buildRegressionSample(body, clusterVariable: nil)
1435 + let includeConstant = !body.hasOption("noconstant")
1436 + let observed = try ZQOLS.fit(
1437 + y: sample.y,
1438 + predictors: sample.predictors,
1439 + includeConstant: includeConstant,
1440 + variance: .classical
1441 + )
1442 + let n = sample.y.count
1443 + let k = observed.coefficients.count
1444 + let generator = Philox4x32(seed: seed)
1445 + let y = sample.y
1446 + let predictors = sample.predictors
1447 +
1448 + // Count |β*| ≥ |β_obs| per coefficient, in parallel chunks —
1449 + // permutations are counter-addressable so order is irrelevant.
1450 + let observedMagnitudes = observed.coefficients.map { abs($0.estimate) }
1451 + let chunkCount = min(
1452 + max(1, ProcessInfo.processInfo.activeProcessorCount), reps
1453 + )
1454 + let chunkSize = (reps + chunkCount - 1) / chunkCount
1455 + let exceedances: [Int] = try await withThrowingTaskGroup(of: [Int].self) { group in
1456 + for chunk in 0..<chunkCount {
1457 + let lower = chunk * chunkSize
1458 + let upper = min(reps, lower + chunkSize)
1459 + guard lower < upper else { continue }
1460 + group.addTask {
1461 + var counts = [Int](repeating: 0, count: k)
1462 + for replicate in lower..<upper {
1463 + let order = ZQResampling.permutationIndices(
1464 + replicate: replicate, sampleSize: n, generator: generator
1465 + )
1466 + let fit = try ZQOLS.fit(
1467 + y: order.map { y[$0] },
1468 + predictors: predictors,
1469 + includeConstant: includeConstant,
1470 + variance: .classical
1471 + )
1472 + for j in 0..<k
1473 + where abs(fit.coefficients[j].estimate) >= observedMagnitudes[j] {
1474 + counts[j] += 1
1475 + }
1476 + }
1477 + return counts
1478 + }
1479 + }
1480 + var total = [Int](repeating: 0, count: k)
1481 + for try await counts in group {
1482 + for j in 0..<k { total[j] += counts[j] }
1483 + }
1484 + return total
1485 + }
1486 +
1487 + var lines = [
1488 + TableFormatter.pad("Permutation test (response permuted)", 46, right: false) +
1489 + "Replications = " + TableFormatter.pad("\(reps)", 10),
1490 + " Number of obs = " +
1491 + TableFormatter.pad("\(n)", 10),
1492 + "",
1493 + ]
1494 + let widths = [12, 12, 8, 8, 10, 11]
1495 + lines.append(
1496 + TableFormatter.pad(sample.responseName, widths[0]) + " | " +
1497 + TableFormatter.pad("Observed", widths[1]) + " " +
1498 + TableFormatter.pad("c", widths[2]) + " " +
1499 + TableFormatter.pad("reps", widths[3]) + " " +
1500 + TableFormatter.pad("p=c/reps", widths[4]) + " " +
1501 + TableFormatter.pad("SE(p)", widths[5])
1502 + )
1503 + lines.append(TableFormatter.rule(widths))
1504 +
1505 + var scalars: [String: Double] = ["reps": Double(reps), "N": Double(n)]
1506 + for j in 0..<k {
1507 + let name = observed.coefficients[j].name
1508 + // The constant's permutation distribution is degenerate under
1509 + // response permutation; Stata omits it too.
1510 + if name == "_cons" { continue }
1511 + let p = Double(exceedances[j]) / Double(reps)
1512 + let standardError = (p * (1 - p) / Double(reps)).squareRoot()
1513 + lines.append(
1514 + TableFormatter.pad(name, widths[0]) + " | " +
1515 + TableFormatter.pad(
1516 + TableFormatter.general(observed.coefficients[j].estimate), widths[1]
1517 + ) + " " +
1518 + TableFormatter.pad("\(exceedances[j])", widths[2]) + " " +
1519 + TableFormatter.pad("\(reps)", widths[3]) + " " +
1520 + TableFormatter.pad(TableFormatter.fixed(p, decimals: 4), widths[4]) + " " +
1521 + TableFormatter.pad(TableFormatter.fixed(standardError, decimals: 4), widths[5])
1522 + )
1523 + scalars["b_\(name)"] = observed.coefficients[j].estimate
1524 + scalars["p_\(name)"] = p
1525 + scalars["c_\(name)"] = Double(exceedances[j])
1526 + }
1527 + return ZQResult(text: lines.joined(separator: "\n"), scalars: scalars)
1528 + }
1529 +
1402 1530 // MARK: - Graphics
1403 1531
1404 1532 private func handleGraph(_ command: ZQCommand) throws -> ZQResult {
@@ -1408,15 +1536,49 @@ public actor ZQSession {
1408 1536 case "scatter": kind = .scatter
1409 1537 case "line": kind = .line
1410 1538 case "histogram": kind = .histogram
1539 + case "kdensity": kind = .kdensity
1411 1540 case let other:
1412 1541 throw ZQEngineError("graph \(other ?? ""): not supported yet")
1413 1542 }
1414 1543
1415 1544 let names = command.varlist.flatMap(\.referencedNames)
1545 + let mask = try observationMask(command)
1546 +
1547 + // Single-variable plots: histogram and kernel density.
1548 + if kind == .histogram || kind == .kdensity {
1549 + guard names.count == 1 else {
1550 + throw ZQEngineError("\(command.subverb ?? ""): syntax is '\(command.subverb ?? "") varname'")
1551 + }
1552 + let (values, missing) = try frame.requireNumeric(names[0])
1553 + let data = (0..<frame.rowCount)
1554 + .filter { mask[$0] && !missing[$0] }
1555 + .map { values[$0] }
1556 + guard data.count > 1 else { throw ZQEngineError("no observations") }
1557 +
1558 + let series: ZQPlotSpec.Series
1559 + if kind == .histogram {
1560 + series = Self.histogramSeries(
1561 + data: data, label: names[0],
1562 + binCount: command.option("bins")?.firstArgument.flatMap(Int.init)
1563 + )
1564 + } else {
1565 + series = Self.kernelDensitySeries(data: data, label: names[0])
1566 + }
1567 + lastPlot = ZQPlotSpec(
1568 + kind: kind,
1569 + xLabel: names[0],
1570 + yLabel: kind == .histogram ? "Frequency" : "Density",
1571 + series: [series]
1572 + )
1573 + return ZQResult(
1574 + text: "(plot created: \(data.count) observations)",
1575 + scalars: ["N": Double(data.count)]
1576 + )
1577 + }
1578 +
1416 1579 guard names.count == 2 else {
1417 1580 throw ZQEngineError("graph \(command.subverb ?? ""): syntax is 'graph \(command.subverb ?? "kind") yvar xvar'")
1418 1581 }
1419 − let mask = try observationMask(command)
1420 1582 let (yValues, yMissing) = try frame.requireNumeric(names[0])
1421 1583 let (xValues, xMissing) = try frame.requireNumeric(names[1])
1422 1584
@@ -1471,6 +1633,63 @@ public actor ZQSession {
1471 1633 return try evaluator.evaluateCondition(expression)
1472 1634 }
1473 1635
1636 + /// Frequency histogram: Sturges bin count by default, half-open bins
1637 + /// with the maximum folded into the last bin. Series points are bin
1638 + /// midpoints vs counts.
1639 + static func histogramSeries(
1640 + data: [Double], label: String, binCount: Int?
1641 + ) -> ZQPlotSpec.Series {
1642 + let n = data.count
1643 + let bins = max(1, binCount ?? Int((Foundation.log2(Double(n))).rounded(.up)) + 1)
1644 + let minimum = data.min()!
1645 + let maximum = data.max()!
1646 + let width = maximum > minimum ? (maximum - minimum) / Double(bins) : 1
1647 + var counts = [Double](repeating: 0, count: bins)
1648 + for value in data {
1649 + let raw = Int((value - minimum) / width)
1650 + counts[min(max(raw, 0), bins - 1)] += 1
1651 + }
1652 + let midpoints = (0..<bins).map { minimum + (Double($0) + 0.5) * width }
1653 + return ZQPlotSpec.Series(label: label, x: midpoints, y: counts)
1654 + }
1655 +
1656 + /// Epanechnikov kernel density on a 100-point grid with Silverman's
1657 + /// bandwidth h = 0.9·min(sd, IQR/1.349)·n^(−1/5) (Stata's default).
1658 + static func kernelDensitySeries(
1659 + data: [Double], label: String
1660 + ) -> ZQPlotSpec.Series {
1661 + let n = data.count
1662 + let mean = data.reduce(0, +) / Double(n)
1663 + let sd = (data.reduce(0) { $0 + ($1 - mean) * ($1 - mean) }
1664 + / Double(n - 1)).squareRoot()
1665 + let sorted = data.sorted()
1666 + let iqr = ZQSummarize.stataPercentile(sorted: sorted, percent: 75)
1667 + - ZQSummarize.stataPercentile(sorted: sorted, percent: 25)
1668 + var spread = min(sd, iqr / 1.349)
1669 + if spread <= 0 { spread = max(sd, 1e-9) }
1670 + let bandwidth = 0.9 * spread * pow(Double(n), -0.2)
1671 +
1672 + let lower = sorted.first! - bandwidth
1673 + let upper = sorted.last! + bandwidth
1674 + let gridSize = 100
1675 + let step = (upper - lower) / Double(gridSize - 1)
1676 + var xs = [Double](repeating: 0, count: gridSize)
1677 + var densities = [Double](repeating: 0, count: gridSize)
1678 + for g in 0..<gridSize {
1679 + let point = lower + Double(g) * step
1680 + xs[g] = point
1681 + var total = 0.0
1682 + for value in data {
1683 + let u = (point - value) / bandwidth
1684 + if abs(u) < 1 {
1685 + total += 0.75 * (1 - u * u) // Epanechnikov kernel
1686 + }
1687 + }
1688 + densities[g] = total / (Double(n) * bandwidth)
1689 + }
1690 + return ZQPlotSpec.Series(label: label, x: xs, y: densities)
1691 + }
1692 +
1474 1693 // MARK: - Qualifiers
1475 1694
1476 1695 /// Combined if/in keep-mask over current observations.
modified MetrikaKit/Sources/ZQGPU/Philox.swift +25 −0
@@ -98,4 +98,29 @@ public enum ZQResampling {
98 98 }
99 99 return indices
100 100 }
101 +
102 + /// A random permutation of 0..<n for one replicate, defined as the
103 + /// argsort of n 64-bit Philox keys. Unlike Fisher–Yates this is
104 + /// order-free — each key is counter-addressable — so a future GPU
105 + /// argsort produces the identical permutation. Key collisions occur
106 + /// with probability ~n²·2⁻⁶⁴ and would only perturb tie order.
107 + ///
108 + /// Counter space: permutations use a dedicated stream offset (bit 63)
109 + /// so they never collide with bootstrap draws under the same seed.
110 + public static func permutationIndices(
111 + replicate: Int, sampleSize: Int, generator: Philox4x32
112 + ) -> [Int] {
113 + let streamOffset = UInt64(1) << 63
114 + let base = UInt64(replicate) &* UInt64(sampleSize)
115 + var keys = [UInt64](repeating: 0, count: sampleSize)
116 + for i in 0..<sampleSize {
117 + // Word positions: offset applied after the ×2 expansion so it
118 + // survives (offset·2 would wrap to 0 mod 2⁶⁴).
119 + let wordIndex = streamOffset &+ (base &+ UInt64(i)) &* 2
120 + let high = UInt64(generator.word(at: wordIndex))
121 + let low = UInt64(generator.word(at: wordIndex &+ 1))
122 + keys[i] = (high << 32) | low
123 + }
124 + return (0..<sampleSize).sorted { keys[$0] < keys[$1] }
125 + }
101 126 }
modified MetrikaKit/Sources/ZQParser/KnownVerbs.swift +1 −0
@@ -48,6 +48,7 @@ public struct ZQVerbTable: Sendable {
48 48 "permute": 7,
49 49 "jackknife": 9,
50 50 "graph": 2,
51 + "kdensity": 4,
51 52 "histogram": 4,
52 53 "scatter": 7,
53 54 "display": 2,
modified MetrikaKit/Sources/ZQPlanner/Planner.swift +3 −1
@@ -57,7 +57,9 @@ public struct ZQPlanner: Sendable {
57 57 guard gpuAvailable else { return .cpu }
58 58
59 59 // Heuristic 2: ≥ threshold independent replicates → batched GPU.
60 − if command.verb == "bootstrap" || command.verb == "permute" {
60 + // (permute joins once its GPU argsort path lands; the permutation
61 + // definition is already argsort-of-Philox-keys for that reason.)
62 + if command.verb == "bootstrap" {
61 63 if let reps = command.option("reps")?.firstArgument,
62 64 let count = Int(reps), count >= gpuReplicateThreshold {
63 65 return .gpu
modified MetrikaKit/Sources/ZQStats/Summarize.swift +1 −1
@@ -104,7 +104,7 @@ public enum ZQSummarize {
104 104 /// Stata percentile definition: with np = n·p/100, if np is integral
105 105 /// the percentile is (x₍np₎ + x₍np+1₎)/2, otherwise x₍⌈np⌉₎ (1-based
106 106 /// order statistics).
107 − static func stataPercentile(sorted: [Double], percent: Double) -> Double {
107 + public static func stataPercentile(sorted: [Double], percent: Double) -> Double {
108 108 let n = sorted.count
109 109 precondition(n > 0)
110 110 let np = Double(n) * percent / 100
added MetrikaKit/Tests/MetrikaKitTests/PermuteTests.swift +115 −0
@@ -0,0 +1,115 @@
1 +//
2 +// PermuteTests.swift
3 +// Metrika
4 +//
5 +// Author: Simon-Pierre Boucher
6 +// Contact: contact@spboucher.ai
7 +// Copyright © 2026 Simon-Pierre Boucher. All rights reserved.
8 +//
9 +
10 +import Foundation
11 +import Testing
12 +import ZQEngine
13 +import ZQGPU
14 +
15 +@Suite("Permutation tests", .serialized)
16 +struct PermuteTests {
17 + let fixtures: Fixtures
18 + let session: ZQSession
19 +
20 + init() async throws {
21 + self.fixtures = try Fixtures()
22 + self.session = try ZQSession(discoverUserCommands: false)
23 + _ = try await session.execute("use \(fixtures.datasetURL.path)")
24 + _ = try await session.execute("gen log_rev = ln(revenue)")
25 + }
26 +
27 + @Test("permutation indices form a valid permutation")
28 + func validPermutation() {
29 + let generator = Philox4x32(seed: 42)
30 + for replicate in [0, 1, 99] {
31 + let order = ZQResampling.permutationIndices(
32 + replicate: replicate, sampleSize: 257, generator: generator
33 + )
34 + #expect(order.sorted() == Array(0..<257), "replicate \(replicate)")
35 + }
36 + }
37 +
38 + @Test("permutations are replicate-addressable and seed-deterministic")
39 + func determinism() {
40 + let a = ZQResampling.permutationIndices(
41 + replicate: 5, sampleSize: 100, generator: Philox4x32(seed: 42)
42 + )
43 + let b = ZQResampling.permutationIndices(
44 + replicate: 5, sampleSize: 100, generator: Philox4x32(seed: 42)
45 + )
46 + let c = ZQResampling.permutationIndices(
47 + replicate: 5, sampleSize: 100, generator: Philox4x32(seed: 43)
48 + )
49 + let d = ZQResampling.permutationIndices(
50 + replicate: 6, sampleSize: 100, generator: Philox4x32(seed: 42)
51 + )
52 + #expect(a == b)
53 + #expect(a != c)
54 + #expect(a != d)
55 + }
56 +
57 + @Test("permutation stream does not collide with the bootstrap stream")
58 + func streamSeparation() {
59 + let generator = Philox4x32(seed: 42)
60 + let bootstrap = ZQResampling.pairsBootstrapIndices(
61 + replicate: 0, sampleSize: 100, generator: generator
62 + )
63 + let permutation = ZQResampling.permutationIndices(
64 + replicate: 0, sampleSize: 100, generator: generator
65 + )
66 + // A collision would make the permutation a deterministic function
67 + // of the bootstrap draw; these must be unrelated streams.
68 + #expect(bootstrap != permutation)
69 + }
70 +
71 + @Test("strong relationship gets an extreme empirical p-value")
72 + func strongRelationship() async throws {
73 + let result = try await session.execute(
74 + "permute, reps(200) seed(42): reg log_rev price"
75 + )
76 + // R² ≈ 0.88 on n=57 — no permutation should beat the observed |β|.
77 + #expect(result.scalars["p_price"] == 0)
78 + #expect(result.scalars["c_price"] == 0)
79 + #expect(result.scalars["reps"] == 200)
80 + #expect(result.text.contains("Permutation test"))
81 + // The degenerate constant row is omitted.
82 + #expect(result.scalars["p__cons"] == nil)
83 + }
84 +
85 + @Test("permute is reproducible for a fixed seed")
86 + func reproducible() async throws {
87 + let first = try await session.execute(
88 + "permute, reps(50) seed(7): reg log_rev orders"
89 + )
90 + let second = try await session.execute(
91 + "permute, reps(50) seed(7): reg log_rev orders"
92 + )
93 + #expect(first.scalars["c_orders"] == second.scalars["c_orders"])
94 + }
95 +
96 + @Test("histogram and kdensity build sensible plot specs")
97 + func singleVariablePlots() async throws {
98 + _ = try await session.execute("histogram revenue")
99 + let histogram = try #require(await session.lastPlot)
100 + #expect(histogram.kind == .histogram)
101 + #expect(histogram.series[0].y.reduce(0, +) == 60) // counts sum to N
102 +
103 + _ = try await session.execute("kdensity revenue")
104 + let density = try #require(await session.lastPlot)
105 + #expect(density.kind == .kdensity)
106 + // Density integrates to ~1 (trapezoid over the grid).
107 + let series = density.series[0]
108 + var integral = 0.0
109 + for i in 1..<series.x.count {
110 + integral += 0.5 * (series.y[i] + series.y[i - 1])
111 + * (series.x[i] - series.x[i - 1])
112 + }
113 + #expect(abs(integral - 1) < 0.05, "density integral \(integral)")
114 + }
115 +}
116