1 files changed, 100 insertions, 0 deletions
diff --git a/pbt/ts/src/stats.ts b/pbt/ts/src/stats.ts
new file mode 100644
index 0000000..8c50856
--- /dev/null
+++ b/pbt/ts/src/stats.ts
@@ -0,0 +1,100 @@
+/**
+ * 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
+ */
+export 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)
+ */
+export 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));
+}
+
+
+export function chiSquare (samples : number [], expected : number []) : number {
+  if (samples.length !== expected.length) {
+    throw new Error("Samples and expected arrays must have the same length.");
+  }
+
+  let chiSquareValue = 0;
+  for (let i = 0; i < samples.length; i++) {
+    const observed = samples[i];
+    const exp = expected[i];
+    if (exp === 0) {
+      throw new Error("Expected frequency cannot be zero.");
+    }
+    chiSquareValue += Math.pow(observed - exp, 2) / exp;
+  }
+
+  return chiSquareValue;
+}