← Back to branch summary

~ubuntu-branches/ubuntu/utopic/paml/utopic-proposed

« back to all changes in this revision

Viewing changes to src/chi2.c

  • Committer: Bazaar Package Importer
  • Author(s): Pjotr Prins
  • Date: 2010-09-11 23:01:37 UTC
  • Revision ID: james.westby@ubuntu.com-20100911230137-jjf5d0blx5p0m9ba
Tags: upstream-4.4c
Import upstream version 4.4c

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/chi2.c
 
1
/* chi2.c 
 
2
          cc -o chi2 chi2.c -lm
 
3
 
 
4
   This program calculates the chi-square significance values for given 
 
5
   degrees of freedom and the tail probability (type I error rate) for 
 
6
   given observed chi-square statistic and degree of freedom.
 
7
 
 
8
      Ziheng Yang,  October 1993.
 
9
*/
 
10
 
 
11
#include <stdio.h>
 
12
#include <stdlib.h>
 
13
#include <math.h>
 
14
 
 
15
double PointChi2 (double prob, double v);
 
16
 
 
17
#define PointGamma(prob,alpha,beta) PointChi2(prob,2.0*(alpha))/(2.0*(beta))
 
18
#define CDFGamma(x,alpha,beta) IncompleteGamma((beta)*(x),alpha,LnGammaFunction(alpha))
 
19
#define CDFChi2(x,v) CDFGamma(x,(v)/2.0,0.5)
 
20
 
 
21
double PointNormal (double prob);
 
22
double CDFNormal (double x);
 
23
double LnGammaFunction (double alpha);
 
24
double IncompleteGamma (double x, double alpha, double ln_gamma_alpha);
 
25
 
 
26
int main(int argc, char*argv[])
 
27
{
 
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};
 
30
 
 
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 {
 
42
      printf ("\n\nChi-square critical values\n");
 
43
      for (i=0; i<ndf; i++) {
 
44
         if (i%15==0) {
 
45
            printf ("\n\t\t\t\tSignificance level\n");
 
46
            printf ("\n DF ");  
 
47
            for (j=0; j<nprob; j++) printf ("%9.4f", 1-prob[j]);
 
48
            printf ("\n");  
 
49
         }
 
50
         printf ("\n%3d ", i+1);
 
51
         for (j=0; j<nprob; j++)
 
52
            printf ("%9.4f", PointChi2(prob[j],(double)(i+1)));
 
53
         if (i%5==4) printf ("\n");
 
54
         if (i%15==14) {
 
55
            printf ("\nENTER for more, (q+ENTER) for quit... ");
 
56
            if (getchar()=='q') break;
 
57
         }
 
58
      }
 
59
   }
 
60
   printf ("\n");
 
61
   return (0);
 
62
}
 
63
 
 
64
double PointNormal (double prob)
 
65
{
 
66
/* returns z so that Prob{x<z}=prob where x ~ N(0,1) and (1e-12)<prob<1-(1e-12)
 
67
   returns (-9999) if in error
 
68
   Odeh RE & Evans JO (1974) The percentage points of the normal distribution.
 
69
   Applied Statistics 22: 96-97 (AS70)
 
70
*/
 
71
   double a0=-.322232431088, a1=-1, a2=-.342242088547, a3=-.0204231210245;
 
72
   double a4=-.453642210148e-4, b0=.0993484626060, b1=.588581570495;
 
73
   double b2=.531103462366, b3=.103537752850, b4=.0038560700634;
 
74
   double y, z=0, p=prob, p1;
 
75
 
 
76
   p1 = (p<0.5 ? p : 1-p);
 
77
   if (p1<1e-20) return (-9999);
 
78
 
 
79
   y = sqrt (log(1/(p1*p1)));   
 
80
   z = y + ((((y*a4+a3)*y+a2)*y+a1)*y+a0) / ((((y*b4+b3)*y+b2)*y+b1)*y+b0);
 
81
   return (p<0.5 ? -z : z);
 
82
}
 
83
 
 
84
double CDFNormal (double x)
 
85
{
 
86
/*  Hill ID  (1973)  The normal integral.  Applied Statistics, 22:424-427.
 
87
    Algorithm AS 66.
 
88
    adapted by Z. Yang, March 1994.  Hill's routine is quite bad, and I 
 
89
    haven't consulted 
 
90
      Adams AG  (1969)  Algorithm 39.  Areas under the normal curve.
 
91
      Computer J. 12: 197-198.
 
92
*/
 
93
    int invers=0;
 
94
    double p, limit=10, t=1.28, y=x*x/2;
 
95
 
 
96
    if (x<0) {  invers=1;  x*=-1; }
 
97
    if (x>limit)  return (invers?0:1);
 
98
    if (x<1.28)  
 
99
       p = .5 - x * (    .398942280444 - .399903438504 * y
 
100
                   /(y + 5.75885480458 - 29.8213557808
 
101
                   /(y + 2.62433121679 + 48.6959930692
 
102
                   /(y + 5.92885724438))));
 
103
    else 
 
104
       p = 0.398942280385 * exp(-y) /
 
105
           (x - 3.8052e-8 + 1.00000615302 /
 
106
           (x + 3.98064794e-4 + 1.98615381364 /
 
107
           (x - 0.151679116635 + 5.29330324926 /
 
108
           (x + 4.8385912808 - 15.1508972451 /
 
109
           (x + 0.742380924027 + 30.789933034 /
 
110
           (x + 3.99019417011))))));
 
111
 
 
112
    return  invers ? p : 1-p;
 
113
}
 
114
 
 
115
double LnGammaFunction (double alpha)
 
116
{
 
117
/* returns ln(gamma(alpha)) for alpha>0, accurate to 10 decimal places.  
 
118
   Stirling's formula is used for the central polynomial part of the procedure.
 
119
   Pike MC & Hill ID (1966) Algorithm 291: Logarithm of the gamma function.
 
120
   Communications of the Association for Computing Machinery, 9:684
 
121
*/
 
122
   double x=alpha, f=0, z;
 
123
 
 
124
   if (x<7) {
 
125
       f=1;  z=x-1;
 
126
       while (++z<7)  f*=z;
 
127
       x=z;   f=-log(f);
 
128
   }
 
129
   z = 1/(x*x);
 
130
   return  f + (x-0.5)*log(x) - x + .918938533204673 
 
131
          + (((-.000595238095238*z+.000793650793651)*z-.002777777777778)*z
 
132
               +.083333333333333)/x;  
 
133
}
 
134
 
 
135
double IncompleteGamma (double x, double alpha, double ln_gamma_alpha)
 
136
{
 
137
/* returns the incomplete gamma ratio I(x,alpha) where x is the upper 
 
138
           limit of the integration and alpha is the shape parameter.
 
139
   returns (-1) if in error
 
140
   ln_gamma_alpha = ln(Gamma(alpha)), is almost redundant.
 
141
   (1) series expansion     if (alpha>x || x<=1)
 
142
   (2) continued fraction   otherwise
 
143
   RATNEST FORTRAN by
 
144
   Bhattacharjee GP (1970) The incomplete gamma integral.  Applied Statistics,
 
145
   19: 285-287 (AS32)
 
146
*/
 
147
   int i;
 
148
   double p=alpha, g=ln_gamma_alpha;
 
149
   double accurate=1e-8, overflow=1e30;
 
150
   double factor, gin=0, rn=0, a=0,b=0,an=0,dif=0, term=0, pn[6];
 
151
 
 
152
   if (x==0) return (0);
 
153
   if (x<0 || p<=0) return (-1);
 
154
 
 
155
   factor=exp(p*log(x)-x-g);   
 
156
   if (x>1 && x>=p) goto l30;
 
157
   /* (1) series expansion */
 
158
   gin=1;  term=1;  rn=p;
 
159
 l20:
 
160
   rn++;
 
161
   term*=x/rn;   gin+=term;
 
162
 
 
163
   if (term > accurate) goto l20;
 
164
   gin*=factor/p;
 
165
   goto l50;
 
166
 l30:
 
167
   /* (2) continued fraction */
 
168
   a=1-p;   b=a+x+1;  term=0;
 
169
   pn[0]=1;  pn[1]=x;  pn[2]=x+1;  pn[3]=x*b;
 
170
   gin=pn[2]/pn[3];
 
171
 l32:
 
172
   a++;  b+=2;  term++;   an=a*term;
 
173
   for (i=0; i<2; i++) pn[i+4]=b*pn[i+2]-an*pn[i];
 
174
   if (pn[5] == 0) goto l35;
 
175
   rn=pn[4]/pn[5];   dif=fabs(gin-rn);
 
176
   if (dif>accurate) goto l34;
 
177
   if (dif<=accurate*rn) goto l42;
 
178
 l34:
 
179
   gin=rn;
 
180
 l35:
 
181
   for (i=0; i<4; i++) pn[i]=pn[i+2];
 
182
   if (fabs(pn[4]) < overflow) goto l32;
 
183
   for (i=0; i<4; i++) pn[i]/=overflow;
 
184
   goto l32;
 
185
 l42:
 
186
   gin=1-factor*gin;
 
187
 
 
188
 l50:
 
189
   return (gin);
 
190
}
 
191
 
 
192
 
 
193
double PointChi2 (double prob, double v)
 
194
{
 
195
/* returns z so that Prob{x<z}=prob where x is Chi2 distributed with df=v
 
196
   returns -1 if in error.   0.000002<prob<0.999998
 
197
   RATNEST FORTRAN by
 
198
       Best DJ & Roberts DE (1975) The percentage points of the 
 
199
       Chi2 distribution.  Applied Statistics 24: 385-388.  (AS91)
 
200
   Converted into C by Ziheng Yang, Oct. 1993.
 
201
*/
 
202
   double e=.5e-6, aa=.6931471805, p=prob, g;
 
203
   double xx, c, ch, a=0,q=0,p1=0,p2=0,t=0,x=0,b=0,s1,s2,s3,s4,s5,s6;
 
204
 
 
205
   if (p<.000002 || p>.999998 || v<=0) return (-1);
 
206
 
 
207
   g = LnGammaFunction (v/2);
 
208
   xx=v/2;   c=xx-1;
 
209
   if (v >= -1.24*log(p)) goto l1;
 
210
 
 
211
   ch=pow((p*xx*exp(g+xx*aa)), 1/xx);
 
212
   if (ch-e<0) return (ch);
 
213
   goto l4;
 
214
l1:
 
215
   if (v>.32) goto l3;
 
216
   ch=0.4;   a=log(1-p);
 
217
l2:
 
218
   q=ch;  p1=1+ch*(4.67+ch);  p2=ch*(6.73+ch*(6.66+ch));
 
219
   t=-0.5+(4.67+2*ch)/p1 - (6.73+ch*(13.32+3*ch))/p2;
 
220
   ch-=(1-exp(a+g+.5*ch+c*aa)*p2/p1)/t;
 
221
   if (fabs(q/ch-1)-.01 <= 0) goto l4;
 
222
   else                       goto l2;
 
223
  
 
224
l3: 
 
225
   x=PointNormal (p);
 
226
   p1=0.222222/v;   ch=v*pow((x*sqrt(p1)+1-p1), 3.0);
 
227
   if (ch>2.2*v+6)  ch=-2*(log(1-p)-c*log(.5*ch)+g);
 
228
l4:
 
229
   q=ch;   p1=.5*ch;
 
230
   if ((t=IncompleteGamma (p1, xx, g))<0) {
 
231
      printf ("\nerr IncompleteGamma");
 
232
      return (-1);
 
233
   }
 
234
   p2=p-t;
 
235
   t=p2*exp(xx*aa+g+p1-c*log(ch));   
 
236
   b=t/ch;  a=0.5*t-b*c;
 
237
 
 
238
   s1=(210+a*(140+a*(105+a*(84+a*(70+60*a))))) / 420;
 
239
   s2=(420+a*(735+a*(966+a*(1141+1278*a))))/2520;
 
240
   s3=(210+a*(462+a*(707+932*a)))/2520;
 
241
   s4=(252+a*(672+1182*a)+c*(294+a*(889+1740*a)))/5040;
 
242
   s5=(84+264*a+c*(175+606*a))/2520;
 
243
   s6=(120+c*(346+127*c))/5040;
 
244
   ch+=t*(1+0.5*t*s1-b*c*(s1-b*(s2-b*(s3-b*(s4-b*(s5-b*s6))))));
 
245
   if (fabs(q/ch-1) > e) goto l4;
 
246
 
 
247
   return (ch);
 
248
}