Fix test for uniformity of distribution for generation of integers

This lead me into a rabbit hole to understand how one should assert
the shape of the distribution of some empirical data w.r.t. to the
"true" expected distribution which uses the Pearson coefficient
test. Unfortunately for very large number of samples deviations sum up
quickly to a point where the computed p-value is basically 0, thus
rejecting the null hypothesis.

While this is theoretically true, there are a lot of ways to introduce
biases that will lead to rejection, one of them being the use of
modulus to select numbers from the generator's output!

1 files changed, 154 insertions, 2 deletions

diff --git a/pbt/ts/src/gen.test.ts b/pbt/ts/src/gen.test.ts
index eaba36f..75c41c8 100644
--- a/pbt/ts/src/gen.test.ts
+++ b/pbt/ts/src/gen.test.ts
@@ -1,7 +1,7 @@
 import Prando from 'prando';
-import {none} from 'fp-ts/lib/Option.js';
+import { none } from 'fp-ts/lib/Option.js';
 import * as G from './gen.js';
-import {test, expect} from 'vitest';
+import { test, expect } from 'vitest';
 
 test('can generate and shrink bytes', async () => {
   const prng = new Prando(Math.random() * 1000);
@@ -17,3 +17,155 @@ test('can generate and shrink bytes', async () => {
     expect(replayed, `Replaying at ${i}`).toEqual(bytes[i - 1]);
   }
 });
+
+test('can generate an arbitrary integer', () => {
+  const prng = new Prando(Math.random() * 1000);
+  const gen = G.genInt();
+  const value = gen.generate();
+  expect(typeof value).toBe('number');
+});
+
+test('can generate N arbitrary integers', () => {
+  const prng = new Prando(Math.random() * 1000);
+  const gen = G.genInt();
+  const value = gen.generateN(5);
+  expect(value.length).toBe(5);
+});
+
+test('can generate arbitrary integer within bounds', () => {
+  const seed = Math.random() * 1000;
+  const gen = G.genIntWithin(1000, 200000, seed);
+  for (let i = 0; i < 100; i++) {
+    const value = gen.generate();
+    expect(value).toBeGreaterThanOrEqual(1000);
+    expect(value).toBeLessThanOrEqual(200000);
+  }
+});
+
+test('generated integers are uniformly distributed', () => {
+  const seed = Math.random() * 1000;
+  const lb = 0;
+  const ub = 65536;
+  const gen = G.genIntWithin(lb, ub, seed);
+  const range = ub - lb;
+  const buckets = 128;
+  const numSamples = 100000;
+  // indexed by lower bound
+  const counts: Record<number, number> = {};
+  for (let i = 1; i <= buckets; i++) {
+    counts[Math.ceil(range / buckets) * i + lb] = 0;
+  }
+  const values = [];
+  for (let i = 0; i < numSamples; i++) {
+    const value = gen.generate();
+    values.push(value);
+    const bucket = Object.keys(counts).find((b) => value < parseInt(b));
+    if (!bucket) {
+      console.log('ignoring bucket for value', value);
+    } else {
+      const b = parseInt(bucket);
+      counts[b] = (counts[b] || 0) + 1;
+    }
+  }
+
+  // compute χ² statistic
+  const expectedCount = numSamples / buckets;
+  const chiSquared = Object.values(counts).reduce((acc, count) => {
+    const diff = count - expectedCount;
+    return acc + (diff * diff) / expectedCount;
+  });
+  // compute p-value and Cramér's V.
+  // In our case the p-value will
+  // systematically be negligible and therefore will reject the null
+  // hypothesis, eg. that the frequency distribution of numbers we
+  // generate is truly uniform. However, Cramér's V is a better
+  // measure of effect size and will be small if the distribution is
+  // uniform.
+  const cramerv = computeCramersV(chiSquared, numSamples, buckets);
+  const pValue = computeChiSquaredPValue(chiSquared, buckets);
+  console.log(`Chi-squared: ${chiSquared}, p-value: ${pValue}, Cramer V: ${cramerv}`);
+  expect(cramerv).toBeLessThan(0.01);
+});
+
+
+/**
+ * Computes the right-tailed p-value for a Chi-Squared distribution
+ * using the Wilson-Hilferty approximation.
+ *
+ * @param chiSquareValue The calculated empirical chi-square statistic (e.g., 1190.372)
+ * @param numBuckets The total number of bins/categories in your test (e.g., 100)
+ * @returns The p-value between 0.0 and 1.0
+ */
+function computeChiSquaredPValue(chiSquareValue: number, numBuckets: number): number {
+  // Degrees of Freedom = Number of buckets - 1
+  const df = numBuckets - 1;
+
+  if (df <= 0) {
+    throw new Error("Number of buckets must be greater than 1.");
+  }
+  if (chiSquareValue < 0) {
+    return 1.0;
+  }
+
+  // 1. Compute the Wilson-Hilferty Normal Approximation Z-score
+  const term1 = Math.pow(chiSquareValue / df, 1 / 3);
+  const term2 = 1 - 2 / (9 * df);
+  const denominator = Math.sqrt(2 / (9 * df));
+
+  const zScore = (term1 - term2) / denominator;
+
+  // 2. Compute right-tailed probability from Z-score using standard Normal CDF
+  return standardNormalRightTail(zScore);
+}
+
+/**
+ * Standard Normal Right-Tail Probability Q(z) = P(Z > z)
+ * Uses a highly precise numerical approximation for the error function (erf).
+ */
+function standardNormalRightTail(z: number): number {
+  // Standard Normal Cumulative Distribution Function (CDF)
+  const cdf = 0.5 * (1 + errorFunction(z / Math.sqrt(2)));
+  // Right tail is 1 - CDF
+  return 1 - cdf;
+}
+
+/**
+ * Abramowitz and Stegun approximation for the Error Function (erf),
+ * accurate to 1.5 x 10^-7
+ *
+ * See https://en.wikipedia.org/wiki/Error_function#Bounds_and_numerical_approximations
+ */
+function errorFunction(x: number): number {
+  const sign = x < 0 ? -1 : 1;
+  const absX = Math.abs(x);
+
+  // Constants for the numerical approximation
+  const p = 0.3275911;
+  const a1 = 0.254829592;
+  const a2 = -0.284496736;
+  const a3 = 1.421413741;
+  const a4 = -1.453152027;
+  const a5 = 1.061405429;
+
+  const t = 1.0 / (1.0 + p * absX);
+  const y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-absX * absX);
+
+  return sign * y;
+}
+
+/**
+ * Calculates Cramér's V to determine the practical effect size (strength of bias).
+ * See https://en.wikipedia.org/wiki/Cramér%27s_V
+ *
+ * @param chiSquare The empirical chi-square value (e.g., 1190.372)
+ * @param totalSamples The total number of data points (e.g., 100000)
+ * @param numBuckets The total number of bins/buckets (e.g., 100)
+ * @returns A value between 0.0 (no bias) and 1.0 (extreme bias)
+ */
+function computeCramersV(chiSquare: number, totalSamples: number, numBuckets: number): number {
+  const df = numBuckets - 1;
+  if (df <= 0 || totalSamples <= 0) return 0;
+
+  // Core Formula: sqrt( chi2 / (N * df) )
+  return Math.sqrt(chiSquare / (totalSamples * df));
+}