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