1 /*
2 *========================================================================
3 * $Id: Vtest.c 223 2006-08-17 06:19:38Z rgb $
4 *
5 * See copyright in copyright.h and the accompanying file COPYING
6 *========================================================================
7 */
8
9 /*
10 *========================================================================
11 * Vtest.c contains a set of routines for evaluating a p-value, the
12 * probability that a given test result was obtained IF the underlying
13 * random number generator was a "good" one (the null hypothesis), given a
14 * vector x of observed results and a vector y of expected results. It
15 * determines the p-value using Pearson's chisq, which does not require
16 * the independent input of the expected sigma for each "bin" (vector
17 * position).
18 *========================================================================
19 */
20 #include <dieharder/libdieharder.h>
21
Vtest_create(Vtest * vtest,unsigned int nvec)22 void Vtest_create(Vtest *vtest, unsigned int nvec)
23 {
24
25 int i;
26 MYDEBUG(D_VTEST){
27 printf("#==================================================================\n");
28 printf("# Vtest_create(): Creating test struct for %u nvec.\n",nvec);
29 }
30 vtest->x = (double *) malloc(sizeof(double)*nvec); /* sample results */
31 vtest->y = (double *) malloc(sizeof(double)*nvec); /* expected sample results */
32 /* zero or set everything */
33 for(i=0;i<nvec;i++){
34 vtest->x[i] = 0.0;
35 vtest->y[i] = 0.0;
36 }
37 vtest->nvec = nvec;
38 vtest->ndof = 0; /* The user must enter this, or it will try to compute it */
39 vtest->chisq = 0.0;
40 vtest->pvalue = 0.0;
41 vtest->cutoff = 5;
42 MYDEBUG(D_VTEST){
43 printf("# Vtest_create(): Done.\n");
44 }
45
46
47 }
48
Vtest_destroy(Vtest * vtest)49 void Vtest_destroy(Vtest *vtest)
50 {
51
52 free(vtest->x);
53 free(vtest->y);
54
55 }
56
Vtest_eval(Vtest * vtest)57 void Vtest_eval(Vtest *vtest)
58 {
59
60 uint i,ndof,itail;
61 double delchisq,chisq;
62 double x_tot,y_tot;
63
64
65 /*
66 * This routine evaluates chi-squared, where:
67 * vtest->x is the trial vector
68 * vtest->y is the exact vector
69 * vtest->sigma is the vector of expected error for y
70 * (for the exact/true distribution)
71 * vtest->nvec is the vector length(s).
72 * vtest->ndof is the number of degrees of freedom (default nvec-1)
73 * vtest->cutoff is the minimum expected count for a cell to be
74 * included in the chisq sum (it should be at least 5, in general,
75 * probably higher in some cases).
76 *
77 * x, y, sigma, nvec all must be filled in my the calling routine.
78 * Be sure to override the default value of ndof if it is known to
79 * the caller.
80 *
81 * Note well that chisq is KNOWN to do poorly -- sometimes very
82 * poorly -- if ndof=1 (two mutually exclusive and exhaustive parameters,
83 * e.g. a normal approximation to the binomial) or if y (the expected
84 * value) for any given cell is less than a cutoff usually set to around
85 * 5. This test will therefore routinely bundle all cells with expected
86 * returns less than the user-defined cutoff AUTOMATICALLY into a single
87 * cell (itail) and use the total number of cells EXCLUSIVE of this
88 * "garbage" cell as the number of degrees of freedom unless ndof is
89 * overridden.
90 */
91 /* verbose=1; */
92 MYDEBUG(D_VTEST){
93 printf("Evaluating chisq and pvalue for %d points\n",vtest->nvec);
94 printf("Using a cutoff of %f\n",vtest->cutoff);
95 }
96
97 chisq = 0.0;
98 x_tot = 0.0;
99 y_tot = 0.0;
100 ndof = 0;
101 itail = -1;
102 MYDEBUG(D_VTEST){
103 printf("# %7s %3s %3s %10s %10s %10s %9s\n",
104 "bit/bin","DoF","X","Y","sigma","del-chisq","chisq");
105 printf("#==================================================================\n");
106 }
107 /*
108 * If vtest->ndof is nonzero, we use it to compute chisq. If not, we try
109 * to estimate it based on a vtest->cutoff that can be set by the caller.
110 * If vtest->ndof is set, the cutoff should probably not be.
111 */
112 for (i=0;i<vtest->nvec;i++) {
113 if(vtest->y[i] >= vtest->cutoff) {
114 x_tot += vtest->x[i];
115 y_tot += vtest->y[i];
116 delchisq = (vtest->x[i] - vtest->y[i])*(vtest->x[i] - vtest->y[i])/vtest->y[i];
117 /* Alternative way of evaluating chisq for binomial only.
118 delchisq = (vtest->x[i] - vtest->y[i])/vtest->sigma[i];
119 delchisq *= delchisq;
120 */
121 chisq += delchisq;
122 MYDEBUG(D_VTEST){
123 printf("# %5u\t%3u\t%12.4f\t%12.4f\t%8.4f\t%10.4f\n",
124 i,vtest->ndof,vtest->x[i],vtest->y[i],delchisq,chisq);
125 }
126 /* increment only if the data is substantial */
127 if(vtest->ndof == 0) ndof++;
128 } else {
129 if(itail == -1){
130 itail = i; /* Do nothing; just remember the index */
131 MYDEBUG(D_VTEST){
132 printf(" Saving itail = %u because vtest->x[i] = %f <= %f\n",itail,vtest->x[i],vtest->cutoff);
133 }
134 } else {
135 /*
136 * Accumulate all the tail expectation here.
137 */
138 vtest->y[itail] += vtest->y[i];
139 vtest->x[itail] += vtest->x[i];
140 }
141 }
142 }
143 /*
144 * At the end, ALL the counts that are statistically weak should sum into
145 * a statistically significant tail count, but the tail count still has
146 * to make the cutoff! Sometimes it won't! Note that the toplevel
147 * conditional guards against itail still being -1 because Vtest did nothing
148 * in its last pass through the code above.
149 */
150 if(itail != -1){
151 if(vtest->y[itail] >= vtest->cutoff){
152 delchisq = (vtest->x[itail] - vtest->y[itail])*
153 (vtest->x[itail] - vtest->y[itail])/vtest->y[itail];
154 chisq += delchisq;
155 /* increment only if the data is substantial */
156 if(vtest->ndof == 0) ndof++;
157 MYDEBUG(D_VTEST){
158 printf("# %5u\t%3u\t%12.4f\t%12.4f\t%8.4f\t%10.4f\n",
159 itail,vtest->ndof,vtest->x[itail],vtest->y[itail],delchisq,chisq);
160 }
161 }
162 }
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179 /*
180 * Interestingly, one simply cannot make the tail "work" as a
181 * contribution to the ndof. The number of degrees of freedom is one
182 * less than the number that make the cutoff, although it does seem
183 * useful to add the last chunk from the tail before doing the
184 * computation of the chisq p-value.
185 */
186 if(vtest->ndof == 0){
187 vtest->ndof = ndof-1;
188 /*
189 * David Bauer: TODO BUG?? The returned ndof is correct, but the
190 * wrong value is used.
191 *
192 * RGB comment: Really, if ndof = 0, the test fails, does it not?
193 * How can you fit a curve with no degrees of freedom. Perhaps
194 * an error and exit, or some other signal of failure? (This should
195 * almost never happen...)
196 */
197 }
198
199 MYDEBUG(D_VTEST){
200 printf("Total: %10.4f %10.4f\n",x_tot,y_tot);
201 printf("#==================================================================\n");
202 printf("Evaluated chisq = %f for %u degrees of freedom\n",chisq,vtest->ndof);
203 }
204 vtest->chisq = chisq;
205
206 /*
207 * Now evaluate the corresponding pvalue. The only real question
208 * is what is the correct number of degrees of freedom. I'd argue we
209 * did use a constraint when we set expected = binomial*nsamp, so we'll
210 * go for ndof (count of nvec tallied) - 1.
211 */
212 vtest->pvalue = gsl_sf_gamma_inc_Q((double)(vtest->ndof)/2.0,chisq/2.0);
213 /* printf("Evaluted pvalue = %6.4f in Vtest_eval() with %u ndof.\n",vtest->pvalue,vtest->ndof); */
214 MYDEBUG(D_VTEST){
215 printf("Evaluted pvalue = %6.4f in Vtest_eval().\n",vtest->pvalue);
216 }
217 /* verbose=0; */
218
219 }
220
221