/**
 * 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.");
  }

  return samples.reduce((acc, observed, i) => {
    const exp = expected[i];
    if (exp === 0) {
      throw new Error("Expected frequency cannot be zero.");
    }
    return acc + Math.pow(observed - exp, 2) / exp;
  }, 0);
}