Orfeo Toolbox  4.0
itkTDistribution.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  *=========================================================================*/
18 #include "itkTDistribution.h"
20 #include "vnl/vnl_erf.h"
21 
22 extern "C" double dbetai_(double *x, double *pin, double *qin);
23 
24 extern "C" double dgamma_(double *x);
25 
26 namespace itk
27 {
28 namespace Statistics
29 {
32 {
33  m_Parameters = ParametersType(1);
34  m_Parameters[0] = 1.0;
35 }
36 
37 void
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 
67 {
68  if ( m_Parameters.GetSize() != 1 )
69  {
70  itkGenericExceptionMacro(
71  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
72  << m_Parameters.size()
73  << " parameters.");
74  }
75  return static_cast< SizeValueType >( m_Parameters[0] );
76 }
77 
78 double
80 ::PDF(double x, SizeValueType degreesOfFreedom)
81 {
82  double dof = static_cast< double >( degreesOfFreedom );
83  double dofplusoneon2 = 0.5 * ( dof + 1.0 );
84  double dofon2 = 0.5 * dof;
85  double pdf;
86 
87  pdf = ( dgamma_(&dofplusoneon2) / dgamma_(&dofon2) )
88  / ( vcl_sqrt(dof * vnl_math::pi) * vcl_pow(1.0 + ( ( x * x ) / dof ), dofplusoneon2) );
89 
90  return pdf;
91 }
92 
93 double
95 ::PDF(double x, const ParametersType & p)
96 {
97  if ( p.GetSize() != 1 )
98  {
99  itkGenericExceptionMacro(
100  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
101  << p.size()
102  << " parameters.");
103  }
104  return TDistribution::PDF( x, static_cast< SizeValueType >( p[0] ) );
105 }
106 
107 double
109 ::CDF(double x, SizeValueType degreesOfFreedom)
110 {
111  double bx;
112  double pin, qin;
113  double dof;
114 
115  // Based on Abramowitz and Stegun 26.7.1, which gives the probability
116  // that the absolute value of a random variable with a Student-t
117  // distribution is less than or equal to a specified t.
118  //
119  // P[|x| <= t] = 1 - Ix(v/2, 1/2)
120  //
121  // where Ix is the incomplete beta function and v is the number of
122  // degrees of freedom in the Student-t distribution and x is
123  // v / (v + t^2).
124  //
125  // To calculate the cdf of the Student-t we need to convert
126  // this formula. For an x >= 0,
127  //
128  // P[|x| <= t] = \int_{-t}^{t} p(x) dx
129  // = 2 \int_0^t p(x) dx
130  //
131  // The cdf of the Student-t is
132  //
133  // P[x <= t] = \int_{-\inf}^t p(x) dx
134  // = 0.5 + \int_0^t p(x) dx (for x >= 0)
135  // = 0.5 + 0.5 * P[|x| < t] (from above)
136  // = 0.5 + 0.5 * (1 - Ix(v/2. 1/2))
137  // = 1 - 0.5 * Ix(v/2, 1/2)
138  //
139  dof = static_cast< double >( degreesOfFreedom );
140  bx = dof / ( dof + ( x * x ) );
141  pin = dof / 2.0;
142  qin = 0.5;
143 
144  if ( x >= 0.0 )
145  {
146  return 1.0 - 0.5 * dbetai_(&bx, &pin, &qin);
147  }
148  else
149  {
150  return 0.5 * dbetai_(&bx, &pin, &qin);
151  }
152 }
153 
154 double
156 ::CDF(double x, const ParametersType & p)
157 {
158  if ( p.GetSize() != 1 )
159  {
160  itkGenericExceptionMacro(
161  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
162  << p.size()
163  << " parameters.");
164  }
165  return TDistribution::CDF( x, static_cast< SizeValueType >( p[0] ) );
166 }
167 
168 double
170 ::InverseCDF(double p, SizeValueType degreesOfFreedom)
171 {
172  if ( p <= 0.0 )
173  {
175  }
176  else if ( p >= 1.0 )
177  {
179  }
180 
181  double x;
182  double dof, dof2, dof3, dof4;
183  double gaussX, gaussX3, gaussX5, gaussX7, gaussX9;
184 
185  // Based on Abramowitz and Stegun 26.7.5
186  dof = static_cast< double >( degreesOfFreedom );
187  dof2 = dof * dof;
188  dof3 = dof * dof2;
189  dof4 = dof * dof3;
190 
192  gaussX3 = vcl_pow(gaussX, 3.0);
193  gaussX5 = vcl_pow(gaussX, 5.0);
194  gaussX7 = vcl_pow(gaussX, 7.0);
195  gaussX9 = vcl_pow(gaussX, 9.0);
196 
197  x = gaussX
198  + ( gaussX3 + gaussX ) / ( 4.0 * dof )
199  + ( 5.0 * gaussX5 + 16.0 * gaussX3 + 3 * gaussX ) / ( 96.0 * dof2 )
200  + ( 3.0 * gaussX7 + 19.0 * gaussX5 + 17.0 * gaussX3 - 15.0 * gaussX ) / ( 384.0 * dof3 )
201  + ( 79.0 * gaussX9
202  + 776.0 * gaussX7
203  + 1482.0 * gaussX5
204  - 1920.0 * gaussX3
205  - 945.0 * gaussX ) / ( 92160.0 * dof4 );
206 
207  // The polynomial approximation above is only accurate for large degrees
208  // of freedom. We'll improve the approximation by a few Newton
209  // iterations.
210  //
211  // 0 iterations, error = 1 at 1 degree of freedom
212  // 3 iterations, error = 10^-10 at 1 degree of freedom
213  // 100 iterations, erorr = 10^-12 at 1 degree of freedom
214  //
215  // 0 iterations, error = 10^-2 at 11 degrees of freedom
216  // 3 iterations, error = 10^-11 at 11 degrees of freedom
217  // 100 iterations, erorr = 10^-12 at 11 degrees of freedom
218  //
219  //
220  // We are trying to find the zero of
221  //
222  // f(x) = p - tcdf(x) = 0;
223  //
224  // So,
225  //
226  // x(n+1) = x(n) - f(x(n)) / f'(x(n))
227  // = x(n) + (p - tcdf(x)) / tpdf(x)
228  //
229  // Note that f'(x) = - tpdf(x)
230  //
231  double delta;
232  for ( unsigned int newt = 0; newt < 3; ++newt )
233  {
234  delta = ( p - TDistribution::CDF(x, degreesOfFreedom) )
235  / TDistribution::PDF(x, degreesOfFreedom);
236  x += delta;
237  }
238 
239  return x;
240 }
241 
242 double
244 ::InverseCDF(double p, const ParametersType & params)
245 {
246  if ( params.GetSize() != 1 )
247  {
248  itkGenericExceptionMacro(
249  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
250  << params.size()
251  << " parameters.");
252  }
253  return TDistribution::InverseCDF( p, static_cast< SizeValueType >( params[0] ) );
254 }
255 
256 double
258 ::EvaluatePDF(double x) const
259 {
260  if ( m_Parameters.GetSize() != 1 )
261  {
262  itkGenericExceptionMacro(
263  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
264  << m_Parameters.size()
265  << " parameters.");
266  }
267  return TDistribution::PDF( x, static_cast< SizeValueType >( m_Parameters[0] ) );
268 }
269 
270 double
272 ::EvaluatePDF(double x, const ParametersType & p) const
273 {
274  if ( p.GetSize() != 1 )
275  {
276  itkGenericExceptionMacro(
277  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
278  << p.size()
279  << " parameters.");
280  }
281  return TDistribution::PDF( x, static_cast< SizeValueType >( p[0] ) );
282 }
283 
284 double
286 ::EvaluatePDF(double x, SizeValueType degreesOfFreedom) const
287 {
288  return TDistribution::PDF(x, degreesOfFreedom);
289 }
290 
291 double
293 ::EvaluateCDF(double x) const
294 {
295  if ( m_Parameters.GetSize() != 1 )
296  {
297  itkGenericExceptionMacro(
298  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
299  << m_Parameters.size()
300  << " parameters.");
301  }
302  return TDistribution::CDF( x, static_cast< SizeValueType >( m_Parameters[0] ) );
303 }
304 
305 double
307 ::EvaluateCDF(double x, const ParametersType & p) const
308 {
309  if ( p.GetSize() != 1 )
310  {
311  itkGenericExceptionMacro(
312  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
313  << p.size()
314  << " parameters.");
315  }
316  return TDistribution::CDF( x, static_cast< SizeValueType >( p[0] ) );
317 }
318 
319 double
321 ::EvaluateCDF(double x, SizeValueType degreesOfFreedom) const
322 {
323  return TDistribution::CDF(x, degreesOfFreedom);
324 }
325 
326 double
328 ::EvaluateInverseCDF(double p) const
329 {
330  if ( m_Parameters.GetSize() != 1 )
331  {
332  itkGenericExceptionMacro(
333  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
334  << m_Parameters.size()
335  << " parameters.");
336  }
337  return TDistribution::InverseCDF( p, static_cast< SizeValueType >( m_Parameters[0] ) );
338 }
339 
340 double
342 ::EvaluateInverseCDF(double p, const ParametersType & params) const
343 {
344  if ( params.GetSize() != 1 )
345  {
346  itkGenericExceptionMacro(
347  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
348  << params.size()
349  << " parameters.");
350  }
351  return TDistribution::InverseCDF( p, static_cast< SizeValueType >( params[0] ) );
352 }
353 
354 double
356 ::EvaluateInverseCDF(double p, SizeValueType degreesOfFreedom) const
357 {
358  return TDistribution::InverseCDF(p, degreesOfFreedom);
359 }
360 
361 bool
364 {
365  if ( m_Parameters.GetSize() == 1 )
366  {
367  if ( m_Parameters[0] > 2 )
368  {
369  return true;
370  }
371  }
372  else
373  {
374  itkGenericExceptionMacro(
375  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
376  << m_Parameters.size()
377  << " parameters.");
378  }
379 
380  return false;
381 }
382 
383 double
385 ::GetMean() const
386 {
387  return 0.0;
388 }
389 
390 double
393 {
394  if ( m_Parameters.GetSize() == 1 )
395  {
396  if ( m_Parameters[0] > 2 )
397  {
398  double dof = static_cast< double >( m_Parameters[0] );
399 
400  return dof / ( dof - 2.0 );
401  }
402  else
403  {
405  }
406  }
407  else
408  {
409  itkGenericExceptionMacro(
410  "Invalid number of parameters to describe distribution. Expected 1 parameter, but got "
411  << m_Parameters.size()
412  << " parameters.");
413  }
414 
416 }
417 
418 void
420 ::PrintSelf(std::ostream & os, Indent indent) const
421 {
422  Superclass::PrintSelf(os, indent);
423 
424  if ( m_Parameters.GetSize() > 0 )
425  {
426  os << indent << "Degrees of freedom: "
427  << static_cast< SizeValueType >( m_Parameters[0] ) << std::endl;
428  }
429  else
430  {
431  os << indent << "Degrees of freedom: (unknown)"
432  << std::endl;
433  }
434 }
435 } // end of namespace Statistics
436 } // end namespace itk

Generated at Sat Mar 8 2014 15:37:00 for Orfeo Toolbox with doxygen 1.8.3.1