Orfeo Toolbox  4.0
itkGaussianDistribution.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  *=========================================================================*/
19 #include "vnl/vnl_erf.h"
20 
21 namespace itk
22 {
23 namespace Statistics
24 {
27 {
28  m_Parameters = ParametersType(2);
29  m_Parameters[0] = 0.0;
30  m_Parameters[1] = 1.0;
31 }
32 
33 double
35 ::GetMean() const
36 {
37  if ( m_Parameters.GetSize() != 2 )
38  {
39  itkGenericExceptionMacro(
40  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
41  << m_Parameters.size()
42  << " parameters.");
43  }
44  return m_Parameters[0];
45 }
46 
47 void
49 ::SetMean(double mean)
50 {
51  bool modified = false;
52 
53  if ( m_Parameters.GetSize() > 0 )
54  {
55  if ( m_Parameters[0] != mean )
56  {
57  modified = true;
58  }
59  }
60 
61  if ( m_Parameters.GetSize() != 2 )
62  {
63  bool cache = false;
64  double t = 1.0;
65 
66  // cache current variance if there is one
67  if ( m_Parameters.GetSize() > 1 )
68  {
69  t = m_Parameters[1];
70  cache = true;
71  }
72 
73  // create a parameters array of the right length
74  m_Parameters = ParametersType(2);
75 
76  // reapply variance if there was one, otherwise use a default value
77  if ( cache )
78  {
79  m_Parameters[1] = t;
80  }
81  else
82  {
83  m_Parameters[1] = 1.0;
84  }
85 
86  modified = true;
87  }
88 
89  m_Parameters[0] = mean;
90 
91  if ( modified )
92  {
93  this->Modified();
94  }
95 }
96 
97 double
100 {
101  if ( m_Parameters.GetSize() != 2 )
102  {
103  itkGenericExceptionMacro(
104  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
105  << m_Parameters.size()
106  << " parameters.");
107  }
108  return m_Parameters[1];
109 }
110 
111 void
113 ::SetVariance(double variance)
114 {
115  bool modified = false;
116 
117  if ( m_Parameters.GetSize() > 1 )
118  {
119  if ( m_Parameters[1] != variance )
120  {
121  modified = true;
122  }
123  }
124 
125  if ( m_Parameters.GetSize() != 2 )
126  {
127  bool cache = false;
128  double t = 0.0;
129 
130  // cache current mean if there is one
131  if ( m_Parameters.GetSize() > 0 )
132  {
133  t = m_Parameters[0];
134  cache = true;
135  }
136 
137  // create a parameters array of the right length
138  m_Parameters = ParametersType(2);
139 
140  // reapply mean if there was one, otherwise use a default value
141  if ( cache )
142  {
143  m_Parameters[0] = t;
144  }
145  else
146  {
147  m_Parameters[0] = 0.0;
148  }
149 
150  modified = true;
151  }
152 
153  m_Parameters[1] = variance;
154 
155  if ( modified )
156  {
157  this->Modified();
158  }
159 }
160 
161 double
163 ::PDF(double x)
164 {
165  return vnl_math::one_over_sqrt2pi * vcl_exp(-0.5 * x * x);
166 }
167 
168 double
170 ::PDF(double x, double mean, double variance)
171 {
172  double xminusmean = x - mean;
173 
174  return ( vnl_math::one_over_sqrt2pi / vcl_sqrt(variance) )
175  * vcl_exp(-0.5 * xminusmean * xminusmean / variance);
176 }
177 
178 double
180 ::PDF(double x, const ParametersType & p)
181 {
182  // verify the parameter vector length
183  if ( p.GetSize() == 2 )
184  {
185  return PDF(x, p[0], p[1]);
186  }
187  else
188  {
189  itkGenericExceptionMacro(
190  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
191  << p.size()
192  << " parameters.");
193  }
194 }
195 
196 double
198 ::CDF(double x)
199 {
200  return 0.5 * ( vnl_erf(vnl_math::sqrt1_2 * x) + 1.0 );
201 }
202 
203 double
205 ::CDF(double x, double mean, double variance)
206 {
207  // convert to zero mean unit variance
208  double u = ( x - mean ) / vcl_sqrt(variance);
209 
210  return 0.5 * ( vnl_erf(vnl_math::sqrt1_2 * u) + 1.0 );
211 }
212 
213 double
215 ::CDF(double x, const ParametersType & p)
216 {
217  if ( p.GetSize() != 2 )
218  {
219  itkGenericExceptionMacro(
220  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
221  << p.size()
222  << " parameters.");
223  }
224  return GaussianDistribution::CDF(x, p[0], p[1]);
225 }
226 
227 double
229 ::InverseCDF(double p)
230 {
231  double dp, dx, dt, ddq, dq;
232  int newt;
233 
234  dp = ( p <= 0.5 ) ? ( p ) : ( 1.0 - p ); /* make between 0 and 0.5 */
235 
236  // if original value is invalid, return +infinity or -infinity
237  // changed from original code to reflect the fact that the
238  // the inverse of P(x) not Q(x) is desired.
239  //
240  // Original line: used for inverse of Q(x)
241  // if( dp <= 0.0 ){ dx = 13.0; return ( (p <= 0.5) ? (dx) : (-dx) ); }
242 
243  // replaced with this if construct for the inverse of P(x)
244  if ( p <= 0.0 )
245  {
247  }
248  else if ( p >= 1.0 )
249  {
251  }
252 
255  dt = vcl_sqrt( -2.0 * vcl_log(dp) );
256  dx = dt
257  - ( ( .010328e+0 * dt + .802853e+0 ) * dt + 2.515517e+0 )
258  / ( ( ( .001308e+0 * dt + .189269e+0 ) * dt + 1.432788e+0 ) * dt + 1.e+0 );
259 
262  for ( newt = 0; newt < 3; newt++ )
263  {
264  dq = 0.5e+0 * vnl_erfc(dx * vnl_math::sqrt1_2) - dp;
265  ddq = vcl_exp(-0.5e+0 * dx * dx) / 2.506628274631000e+0;
266  dx = dx + dq / ddq;
267  }
268 
269  // original line when computing the inverse of Q(x)
270  // return ( (p <= 0.5) ? (dx) : (-dx) ); /* return with correct sign */
271  //
272  // Note that P(-x) = Q(x), so whatever x was calculated for Q(x) = p,
273  // we simply need to return the negative of x to get P(xp) = p.
274  //
275  dx = ( ( p <= 0.5 ) ? ( -dx ) : ( dx ) ); // return with correct sign
276 
277  return dx;
278 }
279 
280 double
282 ::InverseCDF(double p, double mean, double variance)
283 {
284  double x = GaussianDistribution::InverseCDF(p);
285 
287  {
288  return x;
289  }
290  else if ( x == itk::NumericTraits< double >::max() )
291  {
292  return x;
293  }
294  else
295  {
296  // apply the mean and variance to provide the value for the
297  // prescribed Gaussian
298  x = x * vcl_sqrt(variance) + mean;
299  return x;
300  }
301 }
302 
303 double
305 ::InverseCDF(double p, const ParametersType & params)
306 {
307  if ( params.GetSize() != 2 )
308  {
309  itkGenericExceptionMacro(
310  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
311  << params.size()
312  << " parameters.");
313  }
314  return GaussianDistribution::InverseCDF(p, params[0], params[1]);
315 }
316 
317 double
319 ::EvaluatePDF(double x) const
320 {
321  if ( m_Parameters.GetSize() == 2 )
322  {
323  if ( m_Parameters[0] == 0.0 && m_Parameters[1] == 1.0 )
324  {
325  return GaussianDistribution::PDF(x);
326  }
327 
328  return GaussianDistribution::PDF(x, m_Parameters[0], m_Parameters[1]);
329  }
330  else
331  {
332  itkGenericExceptionMacro(
333  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
334  << m_Parameters.size()
335  << " parameters.");
336  }
337 }
338 
339 double
341 ::EvaluatePDF(double x, const ParametersType & p) const
342 {
343  if ( p.GetSize() == 2 )
344  {
345  if ( p[0] == 0.0 && p[1] == 1.0 )
346  {
347  return GaussianDistribution::PDF(x);
348  }
349 
350  return GaussianDistribution::PDF(x, p[0], p[1]);
351  }
352  else
353  {
354  itkGenericExceptionMacro(
355  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
356  << p.size()
357  << " parameters.");
358  }
359 }
360 
361 double
363 ::EvaluatePDF(double x, double mean, double variance) const
364 {
365  if ( mean == 0.0 && variance == 1.0 )
366  {
367  return GaussianDistribution::PDF(x);
368  }
369 
370  return GaussianDistribution::PDF(x, mean, variance);
371 }
372 
373 double
375 ::EvaluateCDF(double x) const
376 {
377  if ( m_Parameters.GetSize() == 2 )
378  {
379  if ( m_Parameters[0] == 0.0 && m_Parameters[1] == 1.0 )
380  {
381  return GaussianDistribution::CDF(x);
382  }
383 
384  return GaussianDistribution::CDF(x, m_Parameters[0], m_Parameters[1]);
385  }
386  else
387  {
388  itkGenericExceptionMacro(
389  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
390  << m_Parameters.size()
391  << " parameters.");
392  }
393 }
394 
395 double
397 ::EvaluateCDF(double x, const ParametersType & p) const
398 {
399  if ( p.GetSize() == 2 )
400  {
401  if ( p[0] == 0.0 && p[1] == 1.0 )
402  {
403  return GaussianDistribution::CDF(x);
404  }
405 
406  return GaussianDistribution::CDF(x, p[0], p[1]);
407  }
408  else
409  {
410  itkGenericExceptionMacro(
411  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
412  << p.size()
413  << " parameters.");
414  }
415 }
416 
417 double
419 ::EvaluateCDF(double x, double mean, double variance) const
420 {
421  if ( mean == 0.0 && variance == 1.0 )
422  {
423  return GaussianDistribution::CDF(x);
424  }
425 
426  return GaussianDistribution::CDF(x, mean, variance);
427 }
428 
429 double
431 ::EvaluateInverseCDF(double p) const
432 {
433  if ( m_Parameters.GetSize() == 2 )
434  {
435  if ( m_Parameters[0] == 0.0 && m_Parameters[1] == 1.0 )
436  {
438  }
439 
440  return GaussianDistribution::InverseCDF(p, m_Parameters[0], m_Parameters[1]);
441  }
442  else
443  {
444  itkGenericExceptionMacro(
445  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
446  << m_Parameters.size()
447  << " parameters.");
448  }
449 }
450 
451 double
453 ::EvaluateInverseCDF(double p, const ParametersType & params) const
454 {
455  if ( params.GetSize() == 2 )
456  {
457  if ( params[0] == 0.0 && params[1] == 1.0 )
458  {
460  }
461 
462  return GaussianDistribution::InverseCDF(p, params[0], params[1]);
463  }
464  else
465  {
466  itkGenericExceptionMacro(
467  "Invalid number of parameters to describe distribution. Expected 2 parameters, but got "
468  << params.size()
469  << " parameters.");
470  }
471 }
472 
473 double
475 ::EvaluateInverseCDF(double p, double mean, double variance) const
476 {
477  if ( mean == 0.0 && variance == 1.0 )
478  {
480  }
481 
482  return GaussianDistribution::InverseCDF(p, mean, variance);
483 }
484 
485 void
487 ::PrintSelf(std::ostream & os, Indent indent) const
488 {
489  Superclass::PrintSelf(os, indent);
490 
491  if ( m_Parameters.GetSize() > 0 )
492  {
493  os << indent << "Mean: " << m_Parameters[0] << std::endl;
494  }
495  else
496  {
497  os << indent << "Mean: (unknown)" << std::endl;
498  }
499 
500  if ( m_Parameters.GetSize() > 1 )
501  {
502  os << indent << "Variance: " << m_Parameters[1] << std::endl;
503  }
504  else
505  {
506  os << indent << "Variance: (unknown)" << std::endl;
507  }
508 }
509 } // end of namespace Statistics
510 } // end namespace itk

Generated at Sat Mar 8 2014 14:39:58 for Orfeo Toolbox with doxygen 1.8.3.1