← Back to branch summary

~ubuntu-branches/ubuntu/raring/paml/raring

« back to all changes in this revision

Viewing changes to src/chi2.c

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2011-09-30 21:48:20 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20110930214820-34le7vftkcq97n9i
Tags: 4.4e-1
http://bugs.debian.org/643843
* Team upload.
* New upstream release.
* Upload to unstable (Closes: #643843).
* Corrected debian/watch for downloading sources with uscan.
* Removed links to inexsistant Git repositories in debian/README.source.
* Corrected VCS URLs (debian/control).
* Conforms to Policy 3.9.2 (debian/control, no other changes needed).
* Use Debhelper 8 (debian/control, debian/compat).
* Converted debian/copyright to machine readable format version 1.0.
* Compress binary packages with xz.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/chi2.c
1
1
/* chi2.c 
2
2
          cc -o chi2 chi2.c -lm
 
3
          cl -O2 chi2.c
3
4
 
4
5
   This program calculates the chi-square significance values for given 
5
6
   degrees of freedom and the tail probability (type I error rate) for 
12
13
#include <stdlib.h>
13
14
#include <math.h>
14
15
 
15
 
double PointChi2 (double prob, double v);
 
16
double QuantileChi2 (double prob, double v);
16
17
 
17
 
#define PointGamma(prob,alpha,beta) PointChi2(prob,2.0*(alpha))/(2.0*(beta))
 
18
#define QuantileGamma(prob,alpha,beta) QuantileChi2(prob,2.0*(alpha))/(2.0*(beta))
18
19
#define CDFGamma(x,alpha,beta) IncompleteGamma((beta)*(x),alpha,LnGammaFunction(alpha))
19
20
#define CDFChi2(x,v) CDFGamma(x,(v)/2.0,0.5)
20
21
 
21
 
double PointNormal (double prob);
 
22
double QuantileNormal (double prob);
22
23
double CDFNormal (double x);
23
24
double LnGammaFunction (double alpha);
24
25
double IncompleteGamma (double x, double alpha, double ln_gamma_alpha);
25
26
 
26
27
int main(int argc, char*argv[])
27
28
{
28
 
   int i,j, n=20, ndf=200, nprob=8, option=0; /* 0:table, 1:prob */
29
 
   double df,chi2, d=1.0/n, prob[]={.005, .025, .1, .5, .90, .95, .99, .999};
 
29
   int i,j, n=20, ndf=200, nprob=8, option=0;
 
30
   double df, chi2, d=1.0/n, prob[]={.005, .025, .1, .5, .90, .95, .99, .999};
30
31
 
31
 
   option=(argc>1);
32
 
   if (option) {
33
 
      for (; ; ) {
34
 
         printf ("\nd.f. & Chi^2 value (Ctrl-c to break)? ");
35
 
         scanf ("%lf%lf", &df, &chi2);
36
 
         if(df<1 || chi2<0) break;
37
 
         prob[0]=1-CDFChi2(chi2,df);
38
 
         printf ("\nprob = %.9f = %.3e\n", prob[0],prob[0]);
39
 
      }
40
 
   }
41
 
   else {
 
32
   if (argc!=2 && argc!=3) {
42
33
      printf ("\n\nChi-square critical values\n");
43
34
      for (i=0; i<ndf; i++) {
44
35
         if (i%15==0) {
49
40
         }
50
41
         printf ("\n%3d ", i+1);
51
42
         for (j=0; j<nprob; j++)
52
 
            printf ("%9.4f", PointChi2(prob[j],(double)(i+1)));
 
43
            printf ("%9.4f", QuantileChi2(prob[j],(double)(i+1)));
53
44
         if (i%5==4) printf ("\n");
54
45
         if (i%15==14) {
55
46
            printf ("\nENTER for more, (q+ENTER) for quit... ");
57
48
         }
58
49
      }
59
50
   }
 
51
   else if(argc==2) {
 
52
      for (; ; ) {
 
53
         printf ("\nd.f. & Chi^2 value (Ctrl-c to break)? ");
 
54
         scanf ("%lf%lf", &df, &chi2);
 
55
         if(df<1 || chi2<0) break;
 
56
         prob[0] = 1-CDFChi2(chi2,df);
 
57
         printf ("\ndf = %2.0f  prob = %.9f = %.3e\n", df, prob[0], prob[0]);
 
58
      }
 
59
   }
 
60
   else if(argc==3) {
 
61
      df = atoi(argv[1]);
 
62
      chi2 = atof(argv[2]);
 
63
      if(df<1 || chi2<0) {
 
64
         printf("df = %d  ch2 = %.4f invalid");
 
65
         exit(-1);
 
66
      }
 
67
      prob[0] = 1 - CDFChi2(chi2, df);
 
68
      printf ("\ndf = %2.0f  prob = %.9f = %.3e\n", df, prob[0], prob[0]);
 
69
   }
60
70
   printf ("\n");
61
71
   return (0);
62
72
}
63
73
 
64
 
double PointNormal (double prob)
 
74
double QuantileNormal (double prob)
65
75
{
66
76
/* returns z so that Prob{x<z}=prob where x ~ N(0,1) and (1e-12)<prob<1-(1e-12)
67
77
   returns (-9999) if in error
190
200
}
191
201
 
192
202
 
193
 
double PointChi2 (double prob, double v)
 
203
double QuantileChi2 (double prob, double v)
194
204
{
195
205
/* returns z so that Prob{x<z}=prob where x is Chi2 distributed with df=v
196
206
   returns -1 if in error.   0.000002<prob<0.999998
197
207
   RATNEST FORTRAN by
198
 
       Best DJ & Roberts DE (1975) The percentage points of the 
 
208
       Best DJ & Roberts DE (1975) The percentage Quantiles of the 
199
209
       Chi2 distribution.  Applied Statistics 24: 385-388.  (AS91)
200
210
   Converted into C by Ziheng Yang, Oct. 1993.
201
211
*/
222
232
   else                       goto l2;
223
233
  
224
234
l3: 
225
 
   x=PointNormal (p);
 
235
   x = QuantileNormal (p);
226
236
   p1=0.222222/v;   ch=v*pow((x*sqrt(p1)+1-p1), 3.0);
227
237
   if (ch>2.2*v+6)  ch=-2*(log(1-p)-c*log(.5*ch)+g);
228
238
l4: