Chi-Square Test

using System;
using System.Diagnostics;
using System.Linq;
using System.Numerics;

public class Program
{
	public static void Main()
	{
		double[] dset1 =
		{
			199809.0,
			200665.0,
			199607.0,
			200270.0,
			199649.0
		};
		double[] dset2 =
		{
			522573.0,
			244456.0,
			139979.0,
			71531.0,
			21461.0
		};
		double[][] dsets =
		{
			dset1,
			dset2
		};
		int[] dslens =
		{
			5,
			5
		};
		for (var k = 0; k < 2; k++)
		{
			Console.Write("Dataset: [ ");
			int l;
			for (l = 0; l < dslens[k]; l++)
				Console.Write("{0}, ", dsets[k][l]);

			Console.Write("]\n");
			var dist = Chi2UniformDistance(dsets[k]);
			var dof = dslens[k] - 1;
			Console.Write("dof: {0}  distance: {1}", dof, dist);
			var prob = Chi2Probability(dof, dist);
			Console.Write(" probability: {0}", prob);
			Console.Write(" uniform? {0}\n", ChiIsUniform(dsets[k], 0.05) ? "Yes" : "No");
		}
		Console.ReadLine();
	}

	private static double _aa1;

	static double Simpson3_8(Func<double, double> f, double a, double b, int n)
	{
		var h = (b - a) / n;
		var sum = f(a) + f(b);
		for (var j = 3 * n - 1; j > 0; j--)
		{
			sum += (j % 3 == 0 ? 3.0 : 2.0) * f(a + h / 3.0 * j);
		}
		return h * sum / 8.0;
	}


	public static double GammaInc_Q(double a, double x)
	{
		Func<double, double> f0 = t => Math.Pow(t, _aa1) * Math.Exp(-t);
		const double h = 1.5e-2;

		var y = _aa1 = a - 1;
		while ((f0(y) * (x - y) > 2.0e-8) && (y < x))
			y += .4;
		if (y > x)
			y = x;

		return 1.0 - Simpson3_8(f0, 0, y, (int)(y / h)) / Gamma(a);
	}

	#region Gamma with Lanczos approximation

	static int g = 7;
	static double[] p = {0.99999999999980993, 676.5203681218851, -1259.1392167224028,
	     771.32342877765313, -176.61502916214059, 12.507343278686905,
	     -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7};

	static double Gamma(double z)
	{
		// Reflection formula
		//if (z < 0.5)
		//{
		//	return Math.PI / (Math.Sin(Math.PI * z) * Gamma(1 - z));
		//}
		//z -= 1;
		//var x = p[0];
		//for (var i = 1; i < g + 2; i++)
		//{
		//	x += p[i] / (z + i);
		//}
		//var t = z + g + 0.5;
		//return Math.Sqrt(2 * Math.PI) * (Math.Pow(t, z + 0.5)) * Math.Exp(-t) * x;
		return Math.Sqrt(2 * Math.PI / z) * Math.Pow((z / Math.E), z);
	}

	#endregion

	public static double Chi2UniformDistance(double[] ds)
	{
		var expected = ds.Sum() / ds.Length;
		return ds.Sum(t => Math.Pow(t - expected, 2)) / expected;
	}

	public static double Chi2Probability(int dof, double distance)
	{
		return 1.0 - GammaInc_Q(0.5 * dof, 0.5 * distance);
	}

	static bool ChiIsUniform(double[] dset, double significance)
	{
		return Chi2Probability(dset.Length - 1, Chi2UniformDistance(dset)) > significance;
	}

}