StRoot  1
probChiSquared.cc
1 /***************************************************************************
2  *
3  * \$Id: probChiSquared.cc,v 1.3 2003/09/02 17:59:35 perev Exp \$
4  *
5  * Author: Thomas Ullrich, Apr 2000
6  ***************************************************************************
7  *
8  * Description:
9  *
10  * Computes the probability that a random variable having a
11  * chi2-distribution with N >= 1 degrees of freedom assumes
12  * a value which is larger than a given value chi2 >= 0.
13  *
14  * The algorithm was taken from CERNLIB prob(G100).
15  *
16  ***************************************************************************
17  *
18  * \$Log: probChiSquared.cc,v \$
19  * Revision 1.3 2003/09/02 17:59:35 perev
20  * gcc 3.2 updates + WarnOff
21  *
22  * Revision 1.2 2000/09/25 19:30:19 ullrich
23  * Fixed code to avoid signed/unsigned comparison warning.
24  *
25  * Revision 1.1 2000/04/06 22:23:32 ullrich
26  * Initial Revision
27  *
28  **************************************************************************/
29 #include <math.h>
30
31 double probChiSquared(double x, unsigned int n)
32 {
33  if(x < 0) return 0;
34
35  const unsigned int nmax = 300;
36  const double r1 = 1;
37  const double hf = r1/2;
38  const double th = r1/3;
39  const double f1 = 2*r1/9;
40  const double c1 = 1.128379167095513;
41  const double chipdf = 100;
42  const double xmax = 174.673;
43  const double xmax2 = 2*xmax;
44  const double xlim = 24;
45  const double eps = 1e-30;
46
47  int m, i;
48  double h, w, s, t, fi, e;
49  double u = hf*x;
50
51  if (x == 0 || n/20 > x)
52  h=1;
53  else if (n == 1) {
54  w = ::sqrt(u);
55  h = w < xlim ? erfc(w) : 0;
56  }
57  else if (n > nmax) {
58  s = r1/n;
59  t = f1*s;
60  w = (::pow(x*s,th)-(1-t))/::sqrt(2*t);
61  if (w < -xlim)
62  h=1;
63  else if (w < xlim)
64  h=hf*erfc(w);
65  else
66  h=0;
67  }
68  else {
69  m=n/2;
70  if(u < xmax2 && (x/n) <= chipdf ) {
71  s=exp(-hf*u);
72  t=s;
73  e=s;
74  if (2*m == (int)n) {
75  fi=0;
76  for (i=1; i<m; i++) {
77  fi += 1;
78  t=u*t/fi;
79  s=s+t;
80  }
81  h=s*e;
82  }
83  else {
84  fi=1;
85  for (i=1; i<m; i++) {
86  fi+=2;
87  t=t*x/fi;
88  s=s+t;
89  }
90  w=::sqrt(u);
91  h = w < xlim ? c1*w*s*e+erfc(w) : 0;
92  }
93  }
94  else
95  h=0;
96  }
97  return h > eps ? h : 0;
98 }