1 package org.lcsim.math.chisq;
2 /** A utility class to return the probability that a value of chi-squared <b> x </b> measured
3 * for a system with <b> n </b>degrees of freedom would be exceeded by chance. Based on routines
4 * in the book <em> "Numerical Recipes: The Art of Scientific Computing"</em>.
5 * @author Norman Graf
6 * @version $Id: ChisqProb.java,v 1.1.1.1 2010/11/30 21:32:00 jeremy Exp $
7 *
8 */
9
10 public class ChisqProb
11 {
12
13 /**Returns the incomplete gamma function P(a, x).
14 *
15 * <br clear="all" /><table border="0" width="100%"><tr><td>
16 * <table align="center"><tr><td nowrap="nowrap" align="center">
17 * <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
18 * > </td><td nowrap="nowrap" align="center">
19 * <font face="symbol">g</font
20 * >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
21 * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
22 * <font face="symbol">�</font
23 * > </td><td nowrap="nowrap" align="center">
24 * 1<hr noshade="noshade" /><font face="symbol">G</font
25 * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
26 * </td><td nowrap="nowrap" align="center">
27 * <font size="-1"><i>x</i></font><!--sup
28 * --><br /><font face="symbol">�<br />�<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
29 * <i>e</i><sup> <font face="symbol">-</font
30 * > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
31 * > 1</sup> <i>dt</i> </td></tr></table>
32 * </td></tr></table>
33 *
34 * This is the probability that the observed chi-square for a correct
35 * model should be less than a value of chi-square.
36 * @param a the number of degrees of freedom
37 * @param x the value of chi-squared
38 * @return The incomplete gamma function P(a, x).
39 */
40 public static double gammp(double a, double x)
41 {
42 if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine gammp");
43 if (x < (a+1.0))
44 { //Use the series representation.
45 return gser(a/2.,x/2.);
46 }
47 else
48 { //Use the continued fraction representation and take its complement.
49 return 1.0-gcf(a/2.,x/2.);
50 }
51 }
52
53 /**Returns the incomplete gamma function Q(a, x)= 1 - P(a, x).
54 *
55 * <br clear="all" /><table border="0" width="100%"><tr><td>
56 * <table align="center"><tr><td nowrap="nowrap" align="center">
57 * <i>Q</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
58 * > 1 <font face="symbol">-</font
59 * > <i>P</i>(<i>a</i>,<i>x</i>) <font face="symbol">�</font
60 * > </td><td nowrap="nowrap" align="center">
61 * <font face="symbol">G</font
62 * >(<i>a</i>,<i>x</i>)<hr noshade="noshade" /><font face="symbol">G</font
63 * >(<i>a</i>)<br /></td><td nowrap="nowrap" align="center">
64 * </td><td nowrap="nowrap" align="center">
65 * <font size="-1"><font face="symbol">�</font
66 * ></font><!--sup
67 * --><br /><font face="symbol">�<br />�<br /></font><font size="-1"><i>x</i></font> <br /></td><td nowrap="nowrap" align="center">
68 * <i>e</i><sup> <font face="symbol">-</font
69 * > <i>t</i></sup> <i>t</i><sup><i>a</i> <font face="symbol">-</font
70 * > 1</sup> <i>dt</i> </td></tr></table>
71 * </td></tr></table>
72 *
73 * This is the probability that the observed chi-square
74 * will exceed the value of chi-square by chance even for a correct model.
75 * @param a the number of degrees of freedom
76 * @param x the value of chi-squared
77 * @return The incomplete gamma function Q(a, x)= 1 - P(a, x).
78 */
79 public static double gammq(double a, double x)
80 {
81
82 if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine gammq");
83 if (x < (a+1.0))
84 { //Use the series representation and take its complement.
85 return 1.0-gser(a/2.,x/2.);
86 }
87 else
88 { //Use the continued fraction representation.
89 return gcf(a/2.,x/2.);
90 }
91 }
92
93 /**Returns the incomplete gamma function P(a, x) evaluated by its series representation.
94 *
95 * <br clear="all" /><table border="0" width="100%"><tr><td>
96 * <table align="center"><tr><td nowrap="nowrap" align="center">
97 * <font face="symbol">g</font
98 * >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
99 * > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td nowrap="nowrap" align="center">
100 * <font size="-1"><font face="symbol">�</font
101 * ></font><!--sup
102 * --><br /><font size="+3"><font face="symbol">�<br />
103 * </font></font><font size="-1"><i>n</i> = 0</font> <br /></td><td nowrap="nowrap" align="center">
104 * </td><td nowrap="nowrap" align="center">
105 * <font face="symbol">G</font
106 * >(<i>a</i>)
107 * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><font face="symbol">G</font
108 * >(<i>a</i> + 1 + <i>n</i>)<br /></td><td nowrap="nowrap" align="center">
109 * <i>x</i><sup><i>n</i></sup></td></tr></table>
110 * </td></tr></table>
111 *
112 * @param a the number of degrees of freedom/2
113 * @param x the value of chi-squared/2
114 * @return The incomplete gamma P(a, x) evaluated by its series representation.
115 */
116 public static double gser(double a, double x)
117 {
118 int ITMAX=100;
119 double EPS=3.0e-7;
120
121 int n;
122 double sum,del,ap;
123 double gln=gammln(a);
124 if (x <= 0.0)
125 {
126 if (x < 0.0) System.out.println("x less than 0 in routine gser");
127 return 0.0;
128 }
129 else
130 {
131 ap=a;
132 del=sum=1.0/a;
133 for (n=1;n<=ITMAX;n++)
134 {
135 ++ap;
136 del *= x/ap;
137 sum += del;
138 if (Math.abs(del) < Math.abs(sum)*EPS)
139 {
140 return sum*Math.exp(-x+a*Math.log(x)-(gln));
141 }
142 }
143 System.out.println("a too large, ITMAX too small in routine gser");
144 return 0.;
145 }
146 }
147
148 /**Returns the incomplete gamma function Q(a, x) evaluated by its continued fraction representation.
149 *
150 * <br clear="all" /><table border="0" width="100%"><tr><td>
151 * <table align="center"><tr><td nowrap="nowrap" align="center">
152 * <font face="symbol">G</font
153 * >(<i>a</i>,<i>x</i>) = <i>e</i><sup> <font face="symbol">-</font
154 * > <i>x</i></sup> <i>x</i><sup><i>a</i></sup> </td><td align="left" class="cl"><font face="symbol">
155 * �<br />�
156 * </font> </td><td nowrap="nowrap" align="center">
157 * 1
158 * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
159 * </td><td nowrap="nowrap" align="center">
160 * 1 <font face="symbol">-</font
161 * > <i>a</i>
162 * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>1 + <br /></td><td nowrap="nowrap" align="center">
163 * </td><td nowrap="nowrap" align="center">
164 * 1
165 * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
166 * </td><td nowrap="nowrap" align="center">
167 * 2 <font face="symbol">-</font
168 * > <i>a</i>
169 * <div class="hrcomp"><hr noshade="noshade" size="1"/></div>1 + <br /></td><td nowrap="nowrap" align="center">
170 * </td><td nowrap="nowrap" align="center">
171 * 2
172 * <div class="hrcomp"><hr noshade="noshade" size="1"/></div><i>x</i> + <br /></td><td nowrap="nowrap" align="center">
173 * <font face="symbol">�</font
174 * > </td><td align="left" class="cl"><font face="symbol">
175 * �<br />�
176 * </font></td><td nowrap="nowrap" align="center">
177 * (<i>x</i> > 0)</td></tr></table>
178 * </td></tr></table>
179 *
180 * @param a the number of degrees of freedom/2
181 * @param x the value of chi-squared/2
182 * @return The incomplete gamma Q(a, x) evaluated by its continued fraction representation.
183 */
184 public static double gcf(double a, double x)
185 {
186 int ITMAX=100; //Maximum allowed number of iteration
187 double EPS=3.0e-7; //Relative accuracy.
188 double FPMIN=1.0e-30; //Number near the smallest representable doubleing-point number.
189 {
190 int i;
191 double an,b,c,d,del,h;
192 double gln=gammln(a);
193 b=x+1.0-a; //Set up for evaluating continued fraction by modified Lentz's method with b0 = 0.
194 c=1.0/FPMIN;
195 d=1.0/b;
196 h=d;
197 for (i=1;i<=ITMAX;i++)
198 { //Iterate to convergence.
199 an = -i*(i-a);
200 b += 2.0;
201 d=an*d+b;
202 if (Math.abs(d) < FPMIN) d=FPMIN;
203 c=b+an/c;
204 if (Math.abs(c) < FPMIN) c=FPMIN;
205 d=1.0/d;
206 del=d*c;
207 h *= del;
208 if (Math.abs(del-1.0) < EPS) break;
209 }
210 if (i > ITMAX) System.out.println("a too large, ITMAX too small in gcf");
211 return Math.exp(-x+a*Math.log(x)-(gln))*h; //Put factors in front.
212 }
213 }
214
215 /**Returns the logarithm of the Gamma function of x, for x>0.
216 *
217 * <br clear="all" /><table border="0" width="100%"><tr><td>
218 * <table align="center"><tr><td nowrap="nowrap" align="center">
219 * <font face="symbol">G</font
220 * >(<i>z</i>) = </td><td nowrap="nowrap" align="center">
221 * <font size="-1"><font face="symbol">�</font
222 * ></font><!--sup
223 * --><br /><font face="symbol">�<br />�<br /></font><font size="-1">0</font> <br /></td><td nowrap="nowrap" align="center">
224 * <i>t</i><sup><i>z</i> <font face="symbol">-</font
225 * > 1</sup> <i>e</i><sup> <font face="symbol">-</font
226 * > <i>t</i></sup> <i>dt</i> </td></tr></table>
227 * </td></tr></table>
228 *
229 *@param x
230 *@return The value ln(Gamma(x))
231 */
232 public static double gammln(double x)
233 {
234 double y,tmp,ser;
235 double[] cof=
236 {76.18009172947146,-86.50532032941677,
237 24.01409824083091,-1.231739572450155,
238 0.1208650973866179e-2,-0.5395239384953e-5};
239 int j;
240 y=x;
241 tmp=x+5.5;
242 tmp -= (x+0.5)*Math.log(tmp);
243 ser=1.000000000190015;
244 for (j=0;j<=5;j++) ser += cof[j]/++y;
245 return -tmp+Math.log(2.5066282746310005*ser/x);
246 }
247 }