aGrUM  0.14.1
chi2.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * Copyright (C) 2005 by Christophe GONZALES and Pierre-Henri WUILLEMIN *
3  * {prenom.nom}_at_lip6.fr *
4  * *
5  * This program is free software; you can redistribute it and/or modify *
6  * it under the terms of the GNU General Public License as published by *
7  * the Free Software Foundation; either version 2 of the License, or *
8  * (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public License *
16  * along with this program; if not, write to the *
17  * Free Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ***************************************************************************/
30 #include <agrum/core/math/chi2.h>
31 
32 #ifndef DOXYGEN_SHOULD_SKIP_THIS
33 
34 // constants used by Gary Perlman for his code for computing chi2 critical
35 // values
36 # define GUM_Z_MAX 6.0 // maximum meaningful z value
37 # define GUM_CHI_EPSILON 0.000001 // accuracy of critchi approximation
38 # define GUM_CHI_MAX 99999.0 // maximum chi square value
39 # define GUM_LOG_SQRT_PI \
40  0.5723649429247000870717135 // std::log (std::sqrt (pi))
41 # define GUM_I_SQRT_PI 0.5641895835477562869480795 // 1 / std::sqrt (pi)
42 # define GUM_BIGX 20.0 // max value to represent exp (x)
43 # define gum__ex(x) (((x) < -GUM_BIGX) ? 0.0 : std::exp(x))
44 
45 #endif /* DOXYGEN_SHOULD_SKIP_THIS */
46 
47 // include the inlined functions if necessary
48 #ifdef GUM_NO_INLINE
50 #endif /* GUM_NO_INLINE */
51 
52 namespace gum {
53 
54  // default constructor
55  Chi2::Chi2(const std::vector< std::size_t >& var_modalities,
56  double confidence_proba) :
57  __modalities(var_modalities),
58  __confidence_proba(confidence_proba) {
59  // for debugging purposes
60  GUM_CONSTRUCTOR(Chi2);
61  }
62 
63  // destructor
65  // for debugging purposes
66  GUM_DESTRUCTOR(Chi2);
67  }
68 
69  // computes the probability of normal z value (used by the cache)
70  double Chi2::__probaZValue(double z) {
71  double y, x, w;
72 
73  if (z == 0.0)
74  x = 0.0;
75  else {
76  y = 0.5 * std::fabs(z);
77 
78  if (y >= (GUM_Z_MAX * 0.5))
79  x = 1.0;
80  else if (y < 1.0) {
81  w = y * y;
82  x = ((((((((0.000124818987 * w - 0.001075204047) * w + 0.005198775019) * w
83  - 0.019198292004)
84  * w
85  + 0.059054035642)
86  * w
87  - 0.151968751364)
88  * w
89  + 0.319152932694)
90  * w
91  - 0.531923007300)
92  * w
93  + 0.797884560593)
94  * y * 2.0;
95  } else {
96  y -= 2.0;
97  x =
98  (((((((((((((-0.000045255659 * y + 0.000152529290) * y - 0.000019538132)
99  * y
100  - 0.000676904986)
101  * y
102  + 0.001390604284)
103  * y
104  - 0.000794620820)
105  * y
106  - 0.002034254874)
107  * y
108  + 0.006549791214)
109  * y
110  - 0.010557625006)
111  * y
112  + 0.011630447319)
113  * y
114  - 0.009279453341)
115  * y
116  + 0.005353579108)
117  * y
118  - 0.002141268741)
119  * y
120  + 0.000535310849)
121  * y
122  + 0.999936657524;
123  }
124  }
125 
126  return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
127  }
128 
129  // computes the probability of chi2 value (used by the cache)
130  double Chi2::probaChi2(double x, Size df) {
131  double a, y = 0, s;
132  double e, c, z;
133  int even; /* true if df is an even number */
134 
135  if ((x <= 0.0) || (df < 1)) return (1.0);
136 
137  a = 0.5 * x;
138 
139  even = (2 * (df / 2)) == df;
140 
141  if (df > 1) y = gum__ex(-a);
142 
143  s = (even ? y : (2.0 * __probaZValue(-std::sqrt(x))));
144 
145  if (df > 2) {
146  x = 0.5 * (df - 1.0);
147  z = (even ? 1.0 : 0.5);
148 
149  if (a > GUM_BIGX) {
150  e = (even ? 0.0 : GUM_LOG_SQRT_PI);
151  c = std::log(a);
152 
153  while (z <= x) {
154  e = std::log(z) + e;
155  s += gum__ex(c * z - a - e);
156  z += 1.0;
157  }
158 
159  return (s);
160  } else {
161  e = (even ? 1.0 : (GUM_I_SQRT_PI / std::sqrt(a)));
162  c = 0.0;
163 
164  while (z <= x) {
165  e = e * (a / z);
166  c = c + e;
167  z += 1.0;
168  }
169 
170  return (c * y + s);
171  }
172  } else
173  return (s);
174  }
175 
176  // computes the critical value of a given chi2 test (used by the cache)
177  double Chi2::__criticalValue(double proba, Size df) {
178  double minchisq = 0.0;
179  double maxchisq = GUM_CHI_MAX;
180  double chisqval;
181 
182  if (proba <= 0.0)
183  return (maxchisq);
184  else if (proba >= 1.0)
185  return (0.0);
186 
187  chisqval = df / std::sqrt(proba); /* fair first value */
188 
189  while (maxchisq - minchisq > GUM_CHI_EPSILON) {
190  if (probaChi2(chisqval, df) < proba)
191  maxchisq = chisqval;
192  else
193  minchisq = chisqval;
194 
195  chisqval = (maxchisq + minchisq) * 0.5;
196  }
197 
198  return (chisqval);
199  }
200 
201 } /* namespace gum */
static double __criticalValue(double proba, Size df)
Computes the critical value of a given chi2 test (used by the cache).
Definition: chi2.cpp:177
gum is the global namespace for all aGrUM entities
Definition: agrum.h:25
~Chi2()
Class destructor.
Definition: chi2.cpp:64
Represent the chi2 distribution.
Definition: chi2.h:55
static double probaChi2(double x, Size df)
Computes the probability of chi2 value.
Definition: chi2.cpp:130
The class that represents the chi2 distribution.
The class that represents the chi2 distribution.
Chi2(const std::vector< std::size_t > &var_modalities, double confidence_proba=GUM_LEARNING_CONFIDENCE_PROBA)
Default constructor.
Definition: chi2.cpp:55
std::size_t Size
In aGrUM, hashed values are unsigned long int.
Definition: types.h:45
static double __probaZValue(double z)
Computes the probability of normal z value.
Definition: chi2.cpp:70