Orfeo Toolbox  3.16
itkChiSquareDistribution.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkChiSquareDistribution.cxx,v $
5  Language: C++
6  Date: $Date: 2009-04-06 11:15:09 $
7  Version: $Revision: 1.3 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 
18 #include "itkChiSquareDistribution.h"
19 #include "itkGaussianDistribution.h"
20 #include "itkNumericTraits.h"
21 #include "vnl/vnl_math.h"
22 #include "vnl/vnl_erf.h"
23 
24 extern "C" double dgami_(double *a, double *x);
25 extern "C" double dgamma_(double *x);
26 
27 namespace itk {
28 namespace Statistics {
29 
30 ChiSquareDistribution
31 ::ChiSquareDistribution()
32 {
33  m_Parameters = ParametersType(1);
34  m_Parameters[0] = 1.0;
35 }
36 
37 void
38 ChiSquareDistribution
39 ::SetDegreesOfFreedom(long dof)
40 {
41  bool modified = false;
42 
43  if (m_Parameters.GetSize() > 0)
44  {
45  if (m_Parameters[0] != static_cast<double>(dof) )
46  {
47  modified = true;
48  }
49  }
50 
51  if (m_Parameters.GetSize() != 1)
52  {
53  m_Parameters = ParametersType(1);
54  }
55 
56  m_Parameters[0] = static_cast<double>(dof);
57 
58  if (modified)
59  {
60  this->Modified();
61  }
62 }
63 
64 long
65 ChiSquareDistribution
66 ::GetDegreesOfFreedom() const
67 {
68  if (m_Parameters.GetSize() == 1)
69  {
70  return static_cast<long>(m_Parameters[0]);
71  }
72  else
73  {
74  InvalidArgumentError excp(__FILE__, __LINE__);
75  ::itk::OStringStream message;
76  message << "itk::ERROR: " << this->GetNameOfClass()
77  << "(" << this << "): "
78  << "Invalid number of parameters to describe distribution.";
79  excp.SetDescription(message.str());
80  excp.SetLocation(ITK_LOCATION);
81  throw excp;
82  }
83 
84  return 1;
85 }
86 
87 double
88 ChiSquareDistribution
89 ::PDF(double x, long degreesOfFreedom)
90 {
91  double dof = static_cast<double>(degreesOfFreedom);
92  double dofon2 = 0.5*dof;
93  double pdf = 0.0;
94 
95  if (x >= 0.0)
96  {
97  pdf = vcl_exp(-0.5*x) * vcl_pow(x, dofon2 - 1.0)
98  / (vcl_pow(2.0, dofon2) * dgamma_(&dofon2));
99  }
100 
101  return pdf;
102 }
103 
104 
105 double
106 ChiSquareDistribution
107 ::PDF(double x, const ParametersType& p)
108 {
109  if (p.GetSize() == 1)
110  {
111  return ChiSquareDistribution::PDF(x, static_cast<long>(p[0]));
112  }
113  else
114  {
115  InvalidArgumentError excp(__FILE__, __LINE__);
116  ::itk::OStringStream message;
117  message << "itk::ERROR: "
118  << "ChiSquareDistribution: "
119  << "Invalid number of parameters to describe distribution.";
120  excp.SetDescription(message.str());
121  excp.SetLocation(ITK_LOCATION);
122  throw excp;
123  }
124 
125  return 0.0;
126 }
127 
128 double
129 ChiSquareDistribution
130 ::CDF(double x, long degreesOfFreedom)
131 {
132  // Based on Abramowitz and Stegun 26.4.19 which relates the
133  // cumulative of the chi-square to incomplete (and complete) gamma
134  // function.
135  if (x < 0)
136  {
137  return 0.0;
138  }
139 
140  double dofon2 = 0.5*degreesOfFreedom;
141  double xon2 = 0.5*x;
142 
143  return dgami_(&dofon2, &xon2) / dgamma_(&dofon2);
144 }
145 
146 
147 double
148 ChiSquareDistribution
149 ::CDF(double x, const ParametersType& p)
150 {
151  if (p.GetSize() == 1)
152  {
153  return ChiSquareDistribution::CDF(x, (long) p[0]);
154  }
155  else
156  {
157  InvalidArgumentError excp(__FILE__, __LINE__);
158  ::itk::OStringStream message;
159  message << "itk::ERROR: "
160  << "ChiSquareDistribution: "
161  << "Invalid number of parameters to describe distribution.";
162  excp.SetDescription(message.str());
163  excp.SetLocation(ITK_LOCATION);
164  throw excp;
165  }
166 
167  return 0.0;
168 }
169 
170 double
171 ChiSquareDistribution
172 ::InverseCDF(double p, long degreesOfFreedom)
173 {
174  if (p <= 0.0)
175  {
176  return itk::NumericTraits<double>::Zero;
177  }
178  else if (p >= 1.0)
179  {
180  return itk::NumericTraits<double>::max();
181  }
182 
183  double x;
184  double dof;
185  double nx;
186 
187  // Based on Abramowitz and Stegun 26.4.17
188  dof = static_cast<double>(degreesOfFreedom);
189  nx = GaussianDistribution::InverseCDF(p);
190 
191  double f = 2.0 / (9.0*dof);
192  x = dof*vcl_pow(1.0 - f + nx*vcl_sqrt(f), 3.0);
193 
194 
195  // The approximation above is only accurate for large degrees of
196  // freedom. We'll improve the approximation by a few Newton iterations.
197  //
198  // 0 iterations, error = 10^-1 at 1 degree of freedom
199  // 3 iterations, error = 10^-11 at 1 degree of freedom
200  // 10 iterations, erorr = 10^-13 at 1 degree of freedom
201  // 20 iterations, erorr = 10^-13 at 1 degree of freedom
202  //
203  // 0 iterations, error = 10^-1 at 11 degrees of freedom
204  // 3 iterations, error = 10^-8 at 11 degrees of freedom
205  // 10 iterations, erorr = 10^-13 at 11 degrees of freedom
206  // 20 iterations, erorr = 10^-13 at 11 degrees of freedom
207  //
208  // 0 iterations, error = 10^-1 at 100 degrees of freedom
209  // 3 iterations, error = 10^-9 at 100 degrees of freedom
210  // 10 iterations, erorr = 10^-10 at 100 degrees of freedom
211  // 20 iterations, erorr = 10^-9 at 100 degrees of freedom
212 
213 
214  // We are trying to find the zero of
215  //
216  // f(x) = p - chisquarecdf(x) = 0;
217  //
218  // So,
219  //
220  // x(n+1) = x(n) - f(x(n)) / f'(x(n))
221  // = x(n) + (p - chisquarecdf(x)) / chisquarepdf(x)
222  //
223  // Note that f'(x) = - chisquarepdf(x)
224  //
225  double delta;
226  for (unsigned int newt = 0; newt < 10; ++newt)
227  {
228  delta = (p - ChiSquareDistribution::CDF(x, degreesOfFreedom))
229  / ChiSquareDistribution::PDF(x, degreesOfFreedom);
230  x += delta;
231  }
232 
233  return x;
234 }
235 
236 
237 double
238 ChiSquareDistribution
239 ::InverseCDF(double p, const ParametersType& params)
240 {
241  if (params.GetSize() == 1)
242  {
243  return ChiSquareDistribution::InverseCDF(p, static_cast<long>(params[0]));
244  }
245  else
246  {
247  InvalidArgumentError excp(__FILE__, __LINE__);
248  ::itk::OStringStream message;
249  message << "itk::ERROR: "
250  << "ChiSquareDistribution: "
251  << "Invalid number of parameters to describe distribution.";
252  excp.SetDescription(message.str());
253  excp.SetLocation(ITK_LOCATION);
254  throw excp;
255  }
256 
257  return 0.0;
258 }
259 
260 double
261 ChiSquareDistribution
262 ::EvaluatePDF(double x) const
263 {
264  if (m_Parameters.GetSize() == 1)
265  {
266  return ChiSquareDistribution::PDF(x, static_cast<long>(m_Parameters[0]));
267  }
268  else
269  {
270  InvalidArgumentError excp(__FILE__, __LINE__);
271  ::itk::OStringStream message;
272  message << "itk::ERROR: " << this->GetNameOfClass()
273  << "(" << this << "): "
274  << "Invalid number of parameters to describe distribution.";
275  excp.SetDescription(message.str());
276  excp.SetLocation(ITK_LOCATION);
277  throw excp;
278  }
279  return 0.0;
280 }
281 
282 double
283 ChiSquareDistribution
284 ::EvaluatePDF(double x, const ParametersType& p) const
285 {
286  if (p.GetSize() == 1)
287  {
288  return ChiSquareDistribution::PDF(x, static_cast<long>(p[0]));
289  }
290  else
291  {
292  InvalidArgumentError excp(__FILE__, __LINE__);
293  ::itk::OStringStream message;
294  message << "itk::ERROR: " << this->GetNameOfClass()
295  << "(" << this << "): "
296  << "Invalid number of parameters to describe distribution.";
297  excp.SetDescription(message.str());
298  excp.SetLocation(ITK_LOCATION);
299  throw excp;
300  }
301  return 0.0;
302 }
303 
304 double
305 ChiSquareDistribution
306 ::EvaluatePDF(double x, long degreesOfFreedom) const
307 {
308  return ChiSquareDistribution::PDF(x, degreesOfFreedom);
309 }
310 
311 
312 double
313 ChiSquareDistribution
314 ::EvaluateCDF(double x) const
315 {
316  if (m_Parameters.GetSize() == 1)
317  {
318  return ChiSquareDistribution::CDF(x, static_cast<long>(m_Parameters[0]));
319  }
320  else
321  {
322  InvalidArgumentError excp(__FILE__, __LINE__);
323  ::itk::OStringStream message;
324  message << "itk::ERROR: " << this->GetNameOfClass()
325  << "(" << this << "): "
326  << "Invalid number of parameters to describe distribution.";
327  excp.SetDescription(message.str());
328  excp.SetLocation(ITK_LOCATION);
329  throw excp;
330  }
331  return 0.0;
332 }
333 
334 double
335 ChiSquareDistribution
336 ::EvaluateCDF(double x, const ParametersType& p) const
337 {
338  if (p.GetSize() == 1)
339  {
340  return ChiSquareDistribution::CDF(x, static_cast<long>(p[0]));
341  }
342  else
343  {
344  InvalidArgumentError excp(__FILE__, __LINE__);
345  ::itk::OStringStream message;
346  message << "itk::ERROR: " << this->GetNameOfClass()
347  << "(" << this << "): "
348  << "Invalid number of parameters to describe distribution.";
349  excp.SetDescription(message.str());
350  excp.SetLocation(ITK_LOCATION);
351  throw excp;
352  }
353  return 0.0;
354 }
355 
356 double
357 ChiSquareDistribution
358 ::EvaluateCDF(double x, long degreesOfFreedom) const
359 {
360  return ChiSquareDistribution::CDF(x, degreesOfFreedom);
361 }
362 
363 
364 double
365 ChiSquareDistribution
366 ::EvaluateInverseCDF(double p) const
367 {
368  if (m_Parameters.GetSize() == 1)
369  {
370  return ChiSquareDistribution::InverseCDF(p,
371  static_cast<long>(m_Parameters[0]));
372  }
373  else
374  {
375  InvalidArgumentError excp(__FILE__, __LINE__);
376  ::itk::OStringStream message;
377  message << "itk::ERROR: " << this->GetNameOfClass()
378  << "(" << this << "): "
379  << "Invalid number of parameters to describe distribution.";
380  excp.SetDescription(message.str());
381  excp.SetLocation(ITK_LOCATION);
382  throw excp;
383  }
384  return 0.0;
385 }
386 
387 double
388 ChiSquareDistribution
389 ::EvaluateInverseCDF(double p, const ParametersType& params) const
390 {
391  if (params.GetSize() == 1)
392  {
393  return ChiSquareDistribution::InverseCDF(p, static_cast<long>(params[0]));
394  }
395  else
396  {
397  InvalidArgumentError excp(__FILE__, __LINE__);
398  ::itk::OStringStream message;
399  message << "itk::ERROR: " << this->GetNameOfClass()
400  << "(" << this << "): "
401  << "Invalid number of parameters to describe distribution.";
402  excp.SetDescription(message.str());
403  excp.SetLocation(ITK_LOCATION);
404  throw excp;
405  }
406  return 0.0;
407 }
408 
409 double
410 ChiSquareDistribution
411 ::EvaluateInverseCDF(double p, long degreesOfFreedom) const
412 {
413  return ChiSquareDistribution::InverseCDF(p, degreesOfFreedom);
414 }
415 
416 
417 double
418 ChiSquareDistribution
419 ::GetMean() const
420 {
421  if (m_Parameters.GetSize() == 1)
422  {
423  return m_Parameters[0];
424  }
425  else
426  {
427  InvalidArgumentError excp(__FILE__, __LINE__);
428  ::itk::OStringStream message;
429  message << "itk::ERROR: " << this->GetNameOfClass()
430  << "(" << this << "): "
431  << "Invalid number of parameters to describe distribution.";
432  excp.SetDescription(message.str());
433  excp.SetLocation(ITK_LOCATION);
434  throw excp;
435  }
436 
437  return NumericTraits<double>::quiet_NaN();
438 }
439 
440 double
441 ChiSquareDistribution
442 ::GetVariance() const
443 {
444  if (m_Parameters.GetSize() == 1)
445  {
446  return 2.0*m_Parameters[0];
447  }
448  else
449  {
450  InvalidArgumentError excp(__FILE__, __LINE__);
451  ::itk::OStringStream message;
452  message << "itk::ERROR: " << this->GetNameOfClass()
453  << "(" << this << "): "
454  << "Invalid number of parameters to describe distribution.";
455  excp.SetDescription(message.str());
456  excp.SetLocation(ITK_LOCATION);
457  throw excp;
458  }
459 
460  return NumericTraits<double>::quiet_NaN();
461 }
462 
463 void
464 ChiSquareDistribution
465 ::PrintSelf(std::ostream& os, Indent indent) const
466 {
467  Superclass::PrintSelf(os,indent);
468 
469  if (m_Parameters.GetSize() > 0)
470  {
471  os << indent << "Degrees of freedom: "
472  << static_cast<long>(m_Parameters[0]) << std::endl;
473  }
474  else
475  {
476  os << indent << "Degrees of freedom: (unknown)"
477  << std::endl;
478  }
479 }
480 
481 } // end of namespace Statistics
482 } // end namespace itk
Generated at Sat Feb 2 2013 23:32:10 for Orfeo Toolbox with doxygen 1.8.1.1