Orfeo Toolbox  4.0
itkChiSquareDistribution.cxx
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
20 
21 extern "C" double dgami_(double *a, double *x);
22 
23 extern "C" double dgamma_(double *x);
24 
25 namespace itk
26 {
27 namespace Statistics
28 {
31 {
32  m_Parameters = ParametersType(1);
33  m_Parameters[0] = 1.0;
34 }
35 
36 void
39 {
40  bool modified = false;
41 
42  if ( m_Parameters.GetSize() > 0 )
43  {
44  if ( m_Parameters[0] != static_cast< double >( dof ) )
45  {
46  modified = true;
47  }
48  }
49 
50  if ( m_Parameters.GetSize() != 1 )
51  {
52  m_Parameters = ParametersType(1);
53  }
54 
55  m_Parameters[0] = static_cast< double >( dof );
56 
57  if ( modified )
58  {
59  this->Modified();
60  }
61 }
62 
66 {
67  if ( m_Parameters.GetSize() != 1 )
68  {
69  itkGenericExceptionMacro(
70  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
71  << m_Parameters.size()
72  << " parameters.");
73  }
74  return static_cast< SizeValueType >( m_Parameters[0] );
75 }
76 
77 double
79 ::PDF(double x, SizeValueType degreesOfFreedom)
80 {
81  double dof = static_cast< double >( degreesOfFreedom );
82  double dofon2 = 0.5 * dof;
83  double pdf = 0.0;
84 
85  if ( x >= 0.0 )
86  {
87  pdf = vcl_exp(-0.5 * x) * vcl_pow(x, dofon2 - 1.0)
88  / ( vcl_pow(2.0, dofon2) * dgamma_(&dofon2) );
89  }
90 
91  return pdf;
92 }
93 
94 double
96 ::PDF(double x, const ParametersType & p)
97 {
98  if ( p.GetSize() != 1 )
99  {
100  itkGenericExceptionMacro(
101  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
102  << p.size()
103  << " parameters.");
104  }
105  return ChiSquareDistribution::PDF( x, static_cast< SizeValueType >( p[0] ) );
106 }
107 
108 double
110 ::CDF(double x, SizeValueType degreesOfFreedom)
111 {
112  // Based on Abramowitz and Stegun 26.4.19 which relates the
113  // cumulative of the chi-square to incomplete (and complete) gamma
114  // function.
115  if ( x < 0 )
116  {
117  return 0.0;
118  }
119 
120  double dofon2 = 0.5 * degreesOfFreedom;
121  double xon2 = 0.5 * x;
122 
123  return dgami_(&dofon2, &xon2) / dgamma_(&dofon2);
124 }
125 
126 double
128 ::CDF(double x, const ParametersType & p)
129 {
130  if ( p.GetSize() != 1 )
131  {
132  itkGenericExceptionMacro(
133  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
134  << p.size()
135  << " parameters.");
136  }
137  return ChiSquareDistribution::CDF(x, (SizeValueType)p[0]);
138 }
139 
140 double
142 ::InverseCDF(double p, SizeValueType degreesOfFreedom)
143 {
144  if ( p <= 0.0 )
145  {
147  }
148  else if ( p >= 1.0 )
149  {
151  }
152 
153  double x;
154  double dof;
155  double nx;
156 
157  // Based on Abramowitz and Stegun 26.4.17
158  dof = static_cast< double >( degreesOfFreedom );
160 
161  double f = 2.0 / ( 9.0 * dof );
162  x = dof * vcl_pow(1.0 - f + nx * vcl_sqrt(f), 3.0);
163 
164  // The approximation above is only accurate for large degrees of
165  // freedom. We'll improve the approximation by a few Newton iterations.
166  //
167  // 0 iterations, error = 10^-1 at 1 degree of freedom
168  // 3 iterations, error = 10^-11 at 1 degree of freedom
169  // 10 iterations, erorr = 10^-13 at 1 degree of freedom
170  // 20 iterations, erorr = 10^-13 at 1 degree of freedom
171  //
172  // 0 iterations, error = 10^-1 at 11 degrees of freedom
173  // 3 iterations, error = 10^-8 at 11 degrees of freedom
174  // 10 iterations, erorr = 10^-13 at 11 degrees of freedom
175  // 20 iterations, erorr = 10^-13 at 11 degrees of freedom
176  //
177  // 0 iterations, error = 10^-1 at 100 degrees of freedom
178  // 3 iterations, error = 10^-9 at 100 degrees of freedom
179  // 10 iterations, erorr = 10^-10 at 100 degrees of freedom
180  // 20 iterations, erorr = 10^-9 at 100 degrees of freedom
181 
182  // We are trying to find the zero of
183  //
184  // f(x) = p - chisquarecdf(x) = 0;
185  //
186  // So,
187  //
188  // x(n+1) = x(n) - f(x(n)) / f'(x(n))
189  // = x(n) + (p - chisquarecdf(x)) / chisquarepdf(x)
190  //
191  // Note that f'(x) = - chisquarepdf(x)
192  //
193  double delta;
194  for ( unsigned int newt = 0; newt < 10; ++newt )
195  {
196  delta = ( p - ChiSquareDistribution::CDF(x, degreesOfFreedom) )
197  / ChiSquareDistribution::PDF(x, degreesOfFreedom);
198  x += delta;
199  }
200 
201  return x;
202 }
203 
204 double
206 ::InverseCDF(double p, const ParametersType & params)
207 {
208  if ( params.GetSize() != 1 )
209  {
210  itkGenericExceptionMacro(
211  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
212  << params.size()
213  << " parameters.");
214  }
215  return ChiSquareDistribution::InverseCDF( p, static_cast< SizeValueType >( params[0] ) );
216 }
217 
218 double
220 ::EvaluatePDF(double x) const
221 {
222  if ( m_Parameters.GetSize() != 1 )
223  {
224  itkGenericExceptionMacro(
225  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
226  << m_Parameters.size()
227  << " parameters.");
228  }
229  return ChiSquareDistribution::PDF( x, static_cast< SizeValueType >( m_Parameters[0] ) );
230 }
231 
232 double
234 ::EvaluatePDF(double x, const ParametersType & p) const
235 {
236  if ( p.GetSize() != 1 )
237  {
238  itkGenericExceptionMacro(
239  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
240  << p.size()
241  << " parameters.");
242  }
243  return ChiSquareDistribution::PDF( x, static_cast< SizeValueType >( p[0] ) );
244 }
245 
246 double
248 ::EvaluatePDF(double x, SizeValueType degreesOfFreedom) const
249 {
250  return ChiSquareDistribution::PDF(x, degreesOfFreedom);
251 }
252 
253 double
255 ::EvaluateCDF(double x) const
256 {
257  if ( m_Parameters.GetSize() != 1 )
258  {
259  itkGenericExceptionMacro(
260  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
261  << m_Parameters.size()
262  << " parameters.");
263  }
264  return ChiSquareDistribution::CDF( x, static_cast< SizeValueType >( m_Parameters[0] ) );
265 }
266 
267 double
269 ::EvaluateCDF(double x, const ParametersType & p) const
270 {
271  if ( p.GetSize() != 1 )
272  {
273  itkGenericExceptionMacro(
274  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
275  << p.size()
276  << " parameters.");
277  }
278  return ChiSquareDistribution::CDF( x, static_cast< SizeValueType >( p[0] ) );
279 }
280 
281 double
283 ::EvaluateCDF(double x, SizeValueType degreesOfFreedom) const
284 {
285  return ChiSquareDistribution::CDF(x, degreesOfFreedom);
286 }
287 
288 double
290 ::EvaluateInverseCDF(double p) const
291 {
292  if ( m_Parameters.GetSize() != 1 )
293  {
294  itkGenericExceptionMacro(
295  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
296  << m_Parameters.size()
297  << " parameters.");
298  }
299  return ChiSquareDistribution::InverseCDF( p, static_cast< SizeValueType >( m_Parameters[0] ) );
300 }
301 
302 double
304 ::EvaluateInverseCDF(double p, const ParametersType & params) const
305 {
306  if ( params.GetSize() != 1 )
307  {
308  itkGenericExceptionMacro(
309  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
310  << params.size()
311  << " parameters.");
312  }
313  return ChiSquareDistribution::InverseCDF( p, static_cast< SizeValueType >( params[0] ) );
314 }
315 
316 double
318 ::EvaluateInverseCDF(double p, SizeValueType degreesOfFreedom) const
319 {
320  return ChiSquareDistribution::InverseCDF(p, degreesOfFreedom);
321 }
322 
323 double
325 ::GetMean() const
326 {
327  if ( m_Parameters.GetSize() != 1 )
328  {
329  itkGenericExceptionMacro(
330  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
331  << m_Parameters.size()
332  << " parameters.");
333  }
334  return m_Parameters[0];
335 }
336 
337 double
340 {
341  if ( m_Parameters.GetSize() != 1 )
342  {
343  itkGenericExceptionMacro(
344  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
345  << m_Parameters.size()
346  << " parameters.");
347  }
348  return 2.0 * m_Parameters[0];
349 }
350 
351 void
353 ::PrintSelf(std::ostream & os, Indent indent) const
354 {
355  Superclass::PrintSelf(os, indent);
356 
357  if ( m_Parameters.GetSize() > 0 )
358  {
359  os << indent << "Degrees of freedom: "
360  << static_cast< SizeValueType >( m_Parameters[0] ) << std::endl;
361  }
362  else
363  {
364  os << indent << "Degrees of freedom: (unknown)"
365  << std::endl;
366  }
367 }
368 } // end of namespace Statistics
369 } // end namespace itk

Generated at Sat Mar 8 2014 14:27:11 for Orfeo Toolbox with doxygen 1.8.3.1