using System; using System.Collections.Generic; using System.Linq; using System.Text; namespace Science.Statistics.BasicStatistics { public class OneSampleChiSquareTest { private double[] of; private double[] ef; public OneSampleChiSquareTest(FrequencyTable table, int degreesOfFreedom) { int cc = table.ObservedFrequency.Count; of = new double[cc]; ef = new double[cc]; df = degreesOfFreedom; for (int i = 0; i < cc; i++) { of[i] = table.ObservedFrequency[i]; ef[i] = table.ExpectedFrequency[i]; } FindpValue(); } private void FindpValue() { int nbins = of.Length; Gamma gam = new Gamma(); int j; double temp; chsq = 0.0; for (j = 0; j < nbins; j++) { if (ef[j] <= 0.0) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Bad expected frequency."); } temp = of[j] - ef[j]; chsq += temp * temp / ef[j]; } pv = gam.gammq(0.5 * (double)df, 0.5 * chsq) * 100.0; } private double chsq; public double ChiSquare { get { return chsq; } } private double pv; public double PValue { get { return pv; } } private int df; public int DegreesOfFreedom { get { return df; } } } class Gauleg18 { protected int ngau = 18; protected double[] y = {0.0021695375159141994, 0.011413521097787704,0.027972308950302116,0.051727015600492421, 0.082502225484340941, 0.12007019910960293,0.16415283300752470, 0.21442376986779355, 0.27051082840644336, 0.33199876341447887, 0.39843234186401943, 0.46931971407375483, 0.54413605556657973, 0.62232745288031077, 0.70331500465597174, 0.78649910768313447, 0.87126389619061517, 0.95698180152629142}; protected double[] w = {0.0055657196642445571, 0.012915947284065419,0.020181515297735382,0.027298621498568734, 0.034213810770299537,0.040875750923643261,0.047235083490265582, 0.053244713977759692,0.058860144245324798,0.064039797355015485, 0.068745323835736408,0.072941885005653087,0.076598410645870640, 0.079687828912071670,0.082187266704339706,0.084078218979661945, 0.085346685739338721,0.085983275670394821}; } class Gamma : Gauleg18 { private int ASWITCH = 100; private double EPS = 2.22044604925031E-10; private double FPMIN = 2.2250738585072E-60; private double gln; public double gammp(double a, double x) { if (x < 0.0 || a <= 0.0) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("bad args in gammp."); } if (x == 0.0) return 0.0; else if ((int)a >= ASWITCH) return gammpapprox(a, x, 1); else if (x < a + 1.0) return gser(a, x); else return 1.0 - gcf(a, x); } public double gammq(double a, double x) { if (x < 0.0 || a <= 0.0) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("bad args in gammq."); } if (x == 0.0) return 1.0; else if ((int)a >= ASWITCH) return gammpapprox(a, x, 0); else if (x < a + 1.0) return 1.0 - gser(a, x); else return gcf(a, x); } public double gser(double a, double x) { Gammln obj = new Gammln(); double sum, del, ap; gln = obj.gammln(a); ap = a; del = sum = 1.0 / a; for (; ; ) { ++ap; del *= x / ap; sum += del; if (Math.Abs(del) < Math.Abs(sum) * EPS) { return sum * Math.Exp(-x + a * Math.Log(x) - gln); } } } public double gcf(double a, double x) { Gammln obj = new Gammln(); int i; double an, b, c, d, del, h; gln = obj.gammln(a); b = x + 1.0 - a; c = 1.0 / FPMIN; d = 1.0 / b; h = d; for (i = 1; ; i++) { an = -i * (i - a); b += 2.0; d = an * d + b; if (Math.Abs(d) < FPMIN) d = FPMIN; c = b + an / c; if (Math.Abs(c) < FPMIN) c = FPMIN; d = 1.0 / d; del = d * c; h *= del; if (Math.Abs(del - 1.0) <= EPS) break; } return Math.Exp(-x + a * Math.Log(x) - gln) * h; } public double gammpapprox(double a, double x, int psig) { Gammln obj = new Gammln(); int j; double xu, t, sum, ans; double a1 = a - 1.0, lna1 = Math.Log(a1), sqrta1 = Math.Sqrt(a1); gln = obj.gammln(a); if (x > a1) xu = Math.Max(a1 + 11.5 * sqrta1, x + 6.0 * sqrta1); else xu = Math.Max(0.0, Math.Min(a1 - 7.5 * sqrta1, x - 5.0 * sqrta1)); sum = 0; for (j = 0; j < ngau; j++) { t = x + (xu - x) * y[j]; sum += w[j] * Math.Exp(-(t - a1) + a1 * (Math.Log(t) - lna1)); } ans = sum * (xu - x) * Math.Exp(a1 * (lna1 - 1.0) - gln); return (psig != 0 ? (ans > 0.0 ? 1.0 - ans : -ans) : (ans >= 0.0 ? ans : 1.0 + ans)); } public double invgammp(double p, double a) { Gammln obj = new Gammln(); int j; double x, err, t, u, pp, lna1 = 0.0, afac = 0.0, a1 = a - 1; double EPS = 1.0e-8; gln = obj.gammln(a); if (a <= 0.0) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("a must be pos in invgammap."); } if (p >= 1.0) return Math.Max(100.0, a + 100.0 * Math.Sqrt(a)); if (p <= 0.0) return 0.0; if (a > 1.0) { lna1 = Math.Log(a1); afac = Math.Exp(a1 * (lna1 - 1.0) - gln); pp = (p < 0.5) ? p : 1.0 - p; t = Math.Sqrt(-2.0 * Math.Log(pp)); x = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t; if (p < 0.5) x = -x; x = Math.Max(1.0e-3, a * Math.Pow(1.0 - 1.0 / (9.0 * a) - x / (3.0 * Math.Sqrt(a)), 3.0)); } else { t = 1.0 - a * (0.253 + a * 0.12); if (p < t) x = Math.Pow(p / t, 1.0 / a); else x = 1.0 - Math.Log(1.0 - (p - t) / (1.0 - t)); } for (j = 0; j < 12; j++) { if (x <= 0.0) return 0.0; err = gammp(a, x) - p; if (a > 1.0) t = afac * Math.Exp(-(x - a1) + a1 * (Math.Log(x) - lna1)); else t = Math.Exp(-x + a1 * Math.Log(x) - gln); u = err / t; x -= (t = u / (1.0 - 0.5 * Math.Min(1.0, u * ((a - 1.0) / x - 1)))); if (x <= 0.0) x = 0.5 * (x + t); if (Math.Abs(t) < EPS * x) break; } return x; } } class Gammln { public Gammln() { } public double gammln(double xx) { int j; double x, tmp, y, ser; double[] cof ={57.1562356658629235,-59.5979603554754912, 14.1360979747417471,-0.491913816097620199,0.339946499848118887e-4, 0.465236289270485756e-4,-0.983744753048795646e-4,0.158088703224912494e-3, -0.210264441724104883e-3,0.217439618115212643e-3,-0.164318106536763890e-3, 0.844182239838527433e-4,-0.261908384015814087e-4,0.368991826595316234e-5}; if (xx <= 0) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("bad arg in gammln."); } y = x = xx; tmp = x + 5.24218750000000000; tmp = (x + 0.5) * Math.Log(tmp) - tmp; ser = 0.999999999999997092; for (j = 0; j < 14; j++) ser += cof[j] / ++y; return tmp + Math.Log(2.5066282746310005 * ser / x); } } }