Orfeo Toolbox  3.16
itkTDistribution.cxx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkTDistribution.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 "itkTDistribution.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 dbetai_(double *x, double *pin, double *qin);
25 extern "C" double dgamma_(double *x);
26 
27 namespace itk {
28 namespace Statistics {
29 
30 TDistribution
31 ::TDistribution()
32 {
33  m_Parameters = ParametersType(1);
34  m_Parameters[0] = 1.0;
35 }
36 
37 void
38 TDistribution
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 TDistribution
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 TDistribution
89 ::PDF(double x, long degreesOfFreedom)
90 {
91  double dof = static_cast<double>(degreesOfFreedom);
92  double dofplusoneon2 = 0.5*(dof+1.0);
93  double dofon2 = 0.5*dof;
94  double pdf;
95 
96  pdf = (dgamma_(&dofplusoneon2) / dgamma_(&dofon2))
97  / (vcl_sqrt(dof*vnl_math::pi) * vcl_pow(1.0 + ((x*x)/dof), dofplusoneon2));
98 
99  return pdf;
100 }
101 
102 double
103 TDistribution
104 ::PDF(double x, const ParametersType& p)
105 {
106  if (p.GetSize() == 1)
107  {
108  return TDistribution::PDF(x, static_cast<long>(p[0]));
109  }
110  else
111  {
112  InvalidArgumentError excp(__FILE__, __LINE__);
113  ::itk::OStringStream message;
114  message << "itk::ERROR: "
115  << "TDistribution: "
116  << "Invalid number of parameters to describe distribution.";
117  excp.SetDescription(message.str());
118  excp.SetLocation(ITK_LOCATION);
119  throw excp;
120  }
121 
122  return 0.0;
123 }
124 
125 double
126 TDistribution
127 ::CDF(double x, long degreesOfFreedom)
128 {
129  double bx;
130  double pin, qin;
131  double dof;
132 
133 
134  // Based on Abramowitz and Stegun 26.7.1, which gives the probability
135  // that the absolute value of a random variable with a Student-t
136  // distribution is less than or equal to a specified t.
137  //
138  // P[|x| <= t] = 1 - Ix(v/2, 1/2)
139  //
140  // where Ix is the incomplete beta function and v is the number of
141  // degrees of freedom in the Student-t distribution and x is
142  // v / (v + t^2).
143  //
144  // To calculate the cdf of the Student-t we need to convert
145  // this formula. For an x >= 0,
146  //
147  // P[|x| <= t] = \int_{-t}^{t} p(x) dx
148  // = 2 \int_0^t p(x) dx
149  //
150  // The cdf of the Student-t is
151  //
152  // P[x <= t] = \int_{-\inf}^t p(x) dx
153  // = 0.5 + \int_0^t p(x) dx (for x >= 0)
154  // = 0.5 + 0.5 * P[|x| < t] (from above)
155  // = 0.5 + 0.5 * (1 - Ix(v/2. 1/2))
156  // = 1 - 0.5 * Ix(v/2, 1/2)
157  //
158  dof = static_cast<double>(degreesOfFreedom);
159  bx = dof / (dof + (x*x));
160  pin = dof / 2.0;
161  qin = 0.5;
162 
163  if (x >= 0.0)
164  {
165  return 1.0 - 0.5 * dbetai_(&bx, &pin, &qin);
166  }
167  else
168  {
169  return 0.5 * dbetai_(&bx, &pin, &qin);
170  }
171 }
172 
173 double
174 TDistribution
175 ::CDF(double x, const ParametersType& p)
176 {
177  if (p.GetSize() == 1)
178  {
179  return TDistribution::CDF(x, static_cast<long>(p[0]));
180  }
181  else
182  {
183  InvalidArgumentError excp(__FILE__, __LINE__);
184  ::itk::OStringStream message;
185  message << "itk::ERROR: "
186  << "TDistribution: "
187  << "Invalid number of parameters to describe distribution.";
188  excp.SetDescription(message.str());
189  excp.SetLocation(ITK_LOCATION);
190  throw excp;
191  }
192 
193  return 0.0;
194 }
195 
196 
197 double
198 TDistribution
199 ::InverseCDF(double p, long degreesOfFreedom)
200 {
201  if (p <= 0.0)
202  {
203  return itk::NumericTraits<double>::NonpositiveMin();
204  }
205  else if (p >= 1.0)
206  {
207  return itk::NumericTraits<double>::max();
208  }
209 
210  double x;
211  double dof, dof2, dof3, dof4;
212  double gaussX, gaussX3, gaussX5, gaussX7, gaussX9;
213 
214  // Based on Abramowitz and Stegun 26.7.5
215  dof = static_cast<double>(degreesOfFreedom);
216  dof2 = dof*dof;
217  dof3 = dof*dof2;
218  dof4 = dof*dof3;
219 
220  gaussX = GaussianDistribution::InverseCDF(p);
221  gaussX3 = vcl_pow(gaussX, 3.0);
222  gaussX5 = vcl_pow(gaussX, 5.0);
223  gaussX7 = vcl_pow(gaussX, 7.0);
224  gaussX9 = vcl_pow(gaussX, 9.0);
225 
226  x = gaussX
227  + (gaussX3 + gaussX) / (4.0 * dof)
228  + (5.0*gaussX5 + 16.0*gaussX3 + 3*gaussX) / (96.0 * dof2)
229  + (3.0*gaussX7 + 19.0*gaussX5 + 17.0*gaussX3 - 15.0*gaussX) / (384.0*dof3)
230  + (79.0*gaussX9
231  + 776.0*gaussX7
232  + 1482.0*gaussX5
233  - 1920.0*gaussX3
234  - 945.0*gaussX) / (92160.0 * dof4);
235 
236  // The polynomial approximation above is only accurate for large degrees
237  // of freedom. We'll improve the approximation by a few Newton
238  // iterations.
239  //
240  // 0 iterations, error = 1 at 1 degree of freedom
241  // 3 iterations, error = 10^-10 at 1 degree of freedom
242  // 100 iterations, erorr = 10^-12 at 1 degree of freedom
243  //
244  // 0 iterations, error = 10^-2 at 11 degrees of freedom
245  // 3 iterations, error = 10^-11 at 11 degrees of freedom
246  // 100 iterations, erorr = 10^-12 at 11 degrees of freedom
247  //
248  //
249  // We are trying to find the zero of
250  //
251  // f(x) = p - tcdf(x) = 0;
252  //
253  // So,
254  //
255  // x(n+1) = x(n) - f(x(n)) / f'(x(n))
256  // = x(n) + (p - tcdf(x)) / tpdf(x)
257  //
258  // Note that f'(x) = - tpdf(x)
259  //
260  double delta;
261  for (unsigned int newt = 0; newt < 3; ++newt)
262  {
263  delta = (p - TDistribution::CDF(x, degreesOfFreedom))
264  / TDistribution::PDF(x, degreesOfFreedom);
265  x += delta;
266  }
267 
268 
269  return x;
270 }
271 
272 double
273 TDistribution
274 ::InverseCDF(double p, const ParametersType& params)
275 {
276  if (params.GetSize() == 1)
277  {
278  return TDistribution::InverseCDF(p, static_cast<long>(params[0]));
279  }
280  else
281  {
282  InvalidArgumentError excp(__FILE__, __LINE__);
283  ::itk::OStringStream message;
284  message << "itk::ERROR: "
285  << "TDistribution: "
286  << "Invalid number of parameters to describe distribution.";
287  excp.SetDescription(message.str());
288  excp.SetLocation(ITK_LOCATION);
289  throw excp;
290  }
291 
292  return 0.0;
293 }
294 
295 
296 double
297 TDistribution
298 ::EvaluatePDF(double x) const
299 {
300  if (m_Parameters.GetSize() == 1)
301  {
302  return TDistribution::PDF(x, static_cast<long>(m_Parameters[0]));
303  }
304  else
305  {
306  InvalidArgumentError excp(__FILE__, __LINE__);
307  ::itk::OStringStream message;
308  message << "itk::ERROR: " << this->GetNameOfClass()
309  << "(" << this << "): "
310  << "Invalid number of parameters to describe distribution.";
311  excp.SetDescription(message.str());
312  excp.SetLocation(ITK_LOCATION);
313  throw excp;
314  }
315  return 0.0;
316 }
317 
318 double
319 TDistribution
320 ::EvaluatePDF(double x, const ParametersType& p) const
321 {
322  if (p.GetSize() == 1)
323  {
324  return TDistribution::PDF(x, static_cast<long>(p[0]));
325  }
326  else
327  {
328  InvalidArgumentError excp(__FILE__, __LINE__);
329  ::itk::OStringStream message;
330  message << "itk::ERROR: " << this->GetNameOfClass()
331  << "(" << this << "): "
332  << "Invalid number of parameters to describe distribution.";
333  excp.SetDescription(message.str());
334  excp.SetLocation(ITK_LOCATION);
335  throw excp;
336  }
337  return 0.0;
338 }
339 
340 double
341 TDistribution
342 ::EvaluatePDF(double x, long degreesOfFreedom) const
343 {
344  return TDistribution::PDF(x, degreesOfFreedom);
345 }
346 
347 
348 double
349 TDistribution
350 ::EvaluateCDF(double x) const
351 {
352  if (m_Parameters.GetSize() == 1)
353  {
354  return TDistribution::CDF(x, static_cast<long>(m_Parameters[0]));
355  }
356  else
357  {
358  InvalidArgumentError excp(__FILE__, __LINE__);
359  ::itk::OStringStream message;
360  message << "itk::ERROR: " << this->GetNameOfClass()
361  << "(" << this << "): "
362  << "Invalid number of parameters to describe distribution.";
363  excp.SetDescription(message.str());
364  excp.SetLocation(ITK_LOCATION);
365  throw excp;
366  }
367  return 0.0;
368 }
369 
370 double
371 TDistribution
372 ::EvaluateCDF(double x, const ParametersType& p) const
373 {
374  if (p.GetSize() == 1)
375  {
376  return TDistribution::CDF(x, static_cast<long>(p[0]));
377  }
378  else
379  {
380  InvalidArgumentError excp(__FILE__, __LINE__);
381  ::itk::OStringStream message;
382  message << "itk::ERROR: " << this->GetNameOfClass()
383  << "(" << this << "): "
384  << "Invalid number of parameters to describe distribution.";
385  excp.SetDescription(message.str());
386  excp.SetLocation(ITK_LOCATION);
387  throw excp;
388  }
389  return 0.0;
390 }
391 
392 double
393 TDistribution
394 ::EvaluateCDF(double x, long degreesOfFreedom) const
395 {
396  return TDistribution::CDF(x, degreesOfFreedom);
397 }
398 
399 
400 double
401 TDistribution
402 ::EvaluateInverseCDF(double p) const
403 {
404  if (m_Parameters.GetSize() == 1)
405  {
406  return TDistribution::InverseCDF(p, static_cast<long>(m_Parameters[0]));
407  }
408  else
409  {
410  InvalidArgumentError excp(__FILE__, __LINE__);
411  ::itk::OStringStream message;
412  message << "itk::ERROR: " << this->GetNameOfClass()
413  << "(" << this << "): "
414  << "Invalid number of parameters to describe distribution.";
415  excp.SetDescription(message.str());
416  excp.SetLocation(ITK_LOCATION);
417  throw excp;
418  }
419  return 0.0;
420 }
421 
422 double
423 TDistribution
424 ::EvaluateInverseCDF(double p, const ParametersType& params) const
425 {
426  if (params.GetSize() == 1)
427  {
428  return TDistribution::InverseCDF(p, static_cast<long>(params[0]));
429  }
430  else
431  {
432  InvalidArgumentError excp(__FILE__, __LINE__);
433  ::itk::OStringStream message;
434  message << "itk::ERROR: " << this->GetNameOfClass()
435  << "(" << this << "): "
436  << "Invalid number of parameters to describe distribution.";
437  excp.SetDescription(message.str());
438  excp.SetLocation(ITK_LOCATION);
439  throw excp;
440  }
441  return 0.0;
442 }
443 
444 double
445 TDistribution
446 ::EvaluateInverseCDF(double p, long degreesOfFreedom) const
447 {
448  return TDistribution::InverseCDF(p, degreesOfFreedom);
449 }
450 
451 bool
452 TDistribution
453 ::HasVariance() const
454 {
455  if (m_Parameters.GetSize() == 1)
456  {
457  if (m_Parameters[0] > 2)
458  {
459  return true;
460  }
461  }
462  else
463  {
464  InvalidArgumentError excp(__FILE__, __LINE__);
465  ::itk::OStringStream message;
466  message << "itk::ERROR: " << this->GetNameOfClass()
467  << "(" << this << "): "
468  << "Invalid number of parameters to describe distribution.";
469  excp.SetDescription(message.str());
470  excp.SetLocation(ITK_LOCATION);
471  throw excp;
472  }
473 
474  return false;
475 }
476 
477 double
478 TDistribution
479 ::GetMean() const
480 {
481  return 0.0;
482 }
483 
484 double
485 TDistribution
486 ::GetVariance() const
487 {
488  if (m_Parameters.GetSize() == 1)
489  {
490  if (m_Parameters[0] > 2)
491  {
492  double dof = static_cast<double>(m_Parameters[0]);
493 
494  return dof / (dof - 2.0);
495  }
496  else
497  {
498  return NumericTraits<double>::quiet_NaN();
499  }
500  }
501  else
502  {
503  InvalidArgumentError excp(__FILE__, __LINE__);
504  ::itk::OStringStream message;
505  message << "itk::ERROR: " << this->GetNameOfClass()
506  << "(" << this << "): "
507  << "Invalid number of parameters to describe distribution.";
508  excp.SetDescription(message.str());
509  excp.SetLocation(ITK_LOCATION);
510  throw excp;
511  }
512 
513 
514  return NumericTraits<double>::quiet_NaN();
515 }
516 
517 void
518 TDistribution
519 ::PrintSelf(std::ostream& os, Indent indent) const
520 {
521  Superclass::PrintSelf(os,indent);
522 
523  if (m_Parameters.GetSize() > 0)
524  {
525  os << indent << "Degrees of freedom: "
526  << static_cast<long>(m_Parameters[0]) << std::endl;
527  }
528  else
529  {
530  os << indent << "Degrees of freedom: (unknown)"
531  << std::endl;
532  }
533 }
534 
535 } // end of namespace Statistics
536 } // end namespace itk
Generated at Sun Feb 3 2013 00:09:29 for Orfeo Toolbox with doxygen 1.8.1.1