aGrUM  0.13.2
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< Size >& var_modalities, double confidence_proba) :
56  __modalities(var_modalities), __confidence_proba(confidence_proba) {
57  // for debugging purposes
58  GUM_CONSTRUCTOR(Chi2);
59  }
60 
61  // destructor
63  // for debugging purposes
64  GUM_DESTRUCTOR(Chi2);
65  }
66 
67  // computes the probability of normal z value (used by the cache)
68  double Chi2::__probaZValue(double z) {
69  double y, x, w;
70 
71  if (z == 0.0)
72  x = 0.0;
73  else {
74  y = 0.5 * std::fabs(z);
75 
76  if (y >= (GUM_Z_MAX * 0.5))
77  x = 1.0;
78  else if (y < 1.0) {
79  w = y * y;
80  x = ((((((((0.000124818987 * w - 0.001075204047) * w + 0.005198775019) * w
81  - 0.019198292004)
82  * w
83  + 0.059054035642)
84  * w
85  - 0.151968751364)
86  * w
87  + 0.319152932694)
88  * w
89  - 0.531923007300)
90  * w
91  + 0.797884560593)
92  * y * 2.0;
93  } else {
94  y -= 2.0;
95  x =
96  (((((((((((((-0.000045255659 * y + 0.000152529290) * y - 0.000019538132)
97  * y
98  - 0.000676904986)
99  * y
100  + 0.001390604284)
101  * y
102  - 0.000794620820)
103  * y
104  - 0.002034254874)
105  * y
106  + 0.006549791214)
107  * y
108  - 0.010557625006)
109  * y
110  + 0.011630447319)
111  * y
112  - 0.009279453341)
113  * y
114  + 0.005353579108)
115  * y
116  - 0.002141268741)
117  * y
118  + 0.000535310849)
119  * y
120  + 0.999936657524;
121  }
122  }
123 
124  return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
125  }
126 
127  // computes the probability of chi2 value (used by the cache)
128  double Chi2::__probaChi2(double x, Size df) {
129  double a, y = 0, s;
130  double e, c, z;
131  int even; /* true if df is an even number */
132 
133  if ((x <= 0.0) || (df < 1)) return (1.0);
134 
135  a = 0.5 * x;
136 
137  even = (2 * (df / 2)) == df;
138 
139  if (df > 1) y = gum__ex(-a);
140 
141  s = (even ? y : (2.0 * __probaZValue(-std::sqrt(x))));
142 
143  if (df > 2) {
144  x = 0.5 * (df - 1.0);
145  z = (even ? 1.0 : 0.5);
146 
147  if (a > GUM_BIGX) {
148  e = (even ? 0.0 : GUM_LOG_SQRT_PI);
149  c = std::log(a);
150 
151  while (z <= x) {
152  e = std::log(z) + e;
153  s += gum__ex(c * z - a - e);
154  z += 1.0;
155  }
156 
157  return (s);
158  } else {
159  e = (even ? 1.0 : (GUM_I_SQRT_PI / std::sqrt(a)));
160  c = 0.0;
161 
162  while (z <= x) {
163  e = e * (a / z);
164  c = c + e;
165  z += 1.0;
166  }
167 
168  return (c * y + s);
169  }
170  } else
171  return (s);
172  }
173 
174  // computes the critical value of a given chi2 test (used by the cache)
175  double Chi2::__criticalValue(double proba, Size df) {
176  double minchisq = 0.0;
177  double maxchisq = GUM_CHI_MAX;
178  double chisqval;
179 
180  if (proba <= 0.0)
181  return (maxchisq);
182  else if (proba >= 1.0)
183  return (0.0);
184 
185  chisqval = df / std::sqrt(proba); /* fair first value */
186 
187  while (maxchisq - minchisq > GUM_CHI_EPSILON) {
188  if (__probaChi2(chisqval, df) < proba)
189  maxchisq = chisqval;
190  else
191  minchisq = chisqval;
192 
193  chisqval = (maxchisq + minchisq) * 0.5;
194  }
195 
196  return (chisqval);
197  }
198 
199 } /* namespace gum */
unsigned long Size
In aGrUM, hashed values are unsigned long int.
Definition: types.h:50
static double __criticalValue(double proba, Size df)
Computes the critical value of a given chi2 test (used by the cache).
Definition: chi2.cpp:175
gum is the global namespace for all aGrUM entities
Definition: agrum.h:25
~Chi2()
Class destructor.
Definition: chi2.cpp:62
static double __probaChi2(double x, Size df)
Computes the probability of chi2 value.
Definition: chi2.cpp:128
Represent the chi2 distribution.
Definition: chi2.h:56
The class that represents the chi2 distribution.
The class that represents the chi2 distribution.
static double __probaZValue(double z)
Computes the probability of normal z value.
Definition: chi2.cpp:68
Chi2(const std::vector< Idx > &var_modalities, double confidence_proba=GUM_LEARNING_CONFIDENCE_PROBA)
Default constructor.
Definition: chi2.cpp:55