Main MRPT website > C++ reference
MRPT logo
distributions.h
Go to the documentation of this file.
1 /* +---------------------------------------------------------------------------+
2  | The Mobile Robot Programming Toolkit (MRPT) C++ library |
3  | |
4  | http://www.mrpt.org/ |
5  | |
6  | Copyright (C) 2005-2012 University of Malaga |
7  | |
8  | This software was written by the Machine Perception and Intelligent |
9  | Robotics Lab, University of Malaga (Spain). |
10  | Contact: Jose-Luis Blanco <jlblanco@ctima.uma.es> |
11  | |
12  | This file is part of the MRPT project. |
13  | |
14  | MRPT is free software: you can redistribute it and/or modify |
15  | it under the terms of the GNU General Public License as published by |
16  | the Free Software Foundation, either version 3 of the License, or |
17  | (at your option) any later version. |
18  | |
19  | MRPT is distributed in the hope that it will be useful, |
20  | but WITHOUT ANY WARRANTY; without even the implied warranty of |
21  | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
22  | GNU General Public License for more details. |
23  | |
24  | You should have received a copy of the GNU General Public License |
25  | along with MRPT. If not, see <http://www.gnu.org/licenses/>. |
26  | |
27  +---------------------------------------------------------------------------+ */
28 #ifndef mrpt_math_distributions_H
29 #define mrpt_math_distributions_H
30 
31 #include <mrpt/utils/utils_defs.h>
32 #include <mrpt/math/math_frwds.h>
33 #include <mrpt/math/CMatrixTemplateNumeric.h>
34 
35 #include <mrpt/math/ops_matrices.h>
36 
37 /*---------------------------------------------------------------
38  Namespace
39  ---------------------------------------------------------------*/
40 namespace mrpt
41 {
42  namespace math
43  {
44  using namespace mrpt::utils;
45 
46  /** \addtogroup stats_grp Statistics functions, probability distributions
47  * \ingroup mrpt_base_grp
48  * @{ */
49 
50  /** Evaluates the univariate normal (Gaussian) distribution at a given point "x".
51  */
52  double BASE_IMPEXP normalPDF(double x, double mu, double std);
53 
54  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
55  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
56  * \param mu A vector or column or row matrix with the Gaussian mean.
57  * \param cov_inv The inverse covariance (information) matrix of the Gaussian.
58  * \param scaled_pdf If set to true, the PDF will be scaled to be in the range [0,1], in contrast to its integral from [-inf,+inf] being 1.
59  */
60  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
61  inline typename MATRIXLIKE::value_type
62  normalPDFInf(
63  const VECTORLIKE1 & x,
64  const VECTORLIKE2 & mu,
65  const MATRIXLIKE & cov_inv,
66  const bool scaled_pdf = false )
67  {
68  MRPT_START
69  typedef typename MATRIXLIKE::value_type T;
70  ASSERTDEB_(cov_inv.isSquare())
71  ASSERTDEB_(size_t(cov_inv.getColCount())==size_t(x.size()) && size_t(cov_inv.getColCount())==size_t(mu.size()))
72  T ret = ::exp( static_cast<T>(-0.5) * mrpt::math::multiply_HCHt_scalar((x-mu), cov_inv ) );
73  return scaled_pdf ? ret : ret * ::sqrt(cov_inv.det()) / ::pow(static_cast<T>(M_2PI),static_cast<T>( size(cov_inv,1) ));
74  MRPT_END
75  }
76 
77  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
78  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
79  * \param mu A vector or column or row matrix with the Gaussian mean.
80  * \param cov The covariance matrix of the Gaussian.
81  * \param scaled_pdf If set to true, the PDF will be scaled to be in the range [0,1], in contrast to its integral from [-inf,+inf] being 1.
82  */
83  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
84  inline typename MATRIXLIKE::value_type
85  normalPDF(
86  const VECTORLIKE1 & x,
87  const VECTORLIKE2 & mu,
88  const MATRIXLIKE & cov,
89  const bool scaled_pdf = false )
90  {
91  return normalPDFInf(x,mu,cov.inverse(),scaled_pdf);
92  }
93 
94  /** Evaluates the multivariate normal (Gaussian) distribution at a given point given its distance vector "d" from the Gaussian mean.
95  */
96  template <typename VECTORLIKE,typename MATRIXLIKE>
97  typename MATRIXLIKE::value_type
98  normalPDF(const VECTORLIKE &d,const MATRIXLIKE &cov)
99  {
100  MRPT_START
101  ASSERTDEB_(cov.isSquare())
102  ASSERTDEB_(size_t(cov.getColCount())==size_t(d.size()))
103  return std::exp( static_cast<typename MATRIXLIKE::value_type>(-0.5)*mrpt::math::multiply_HCHt_scalar(d,cov.inverse()))
104  / (::pow(
105  static_cast<typename MATRIXLIKE::value_type>(M_2PI),
106  static_cast<typename MATRIXLIKE::value_type>(cov.getColCount()))
107  * ::sqrt(cov.det()));
108  MRPT_END
109  }
110 
111  /** Kullback-Leibler divergence (KLD) between two independent multivariate Gaussians.
112  *
113  * \f$ D_\mathrm{KL}(\mathcal{N}_0 \| \mathcal{N}_1) = { 1 \over 2 } ( \log_e ( { \det \Sigma_1 \over \det \Sigma_0 } ) + \mathrm{tr} ( \Sigma_1^{-1} \Sigma_0 ) + ( \mu_1 - \mu_0 )^\top \Sigma_1^{-1} ( \mu_1 - \mu_0 ) - N ) \f$
114  */
115  template <typename VECTORLIKE1,typename MATRIXLIKE1,typename VECTORLIKE2,typename MATRIXLIKE2>
116  double KLD_Gaussians(
117  const VECTORLIKE1 &mu0, const MATRIXLIKE1 &cov0,
118  const VECTORLIKE2 &mu1, const MATRIXLIKE2 &cov1)
119  {
120  MRPT_START
121  ASSERT_(size_t(mu0.size())==size_t(mu1.size()) && size_t(mu0.size())==size_t(size(cov0,1)) && size_t(mu0.size())==size_t(size(cov1,1)) && cov0.isSquare() && cov1.isSquare() )
122  const size_t N = mu0.size();
123  MATRIXLIKE2 cov1_inv;
124  cov1.inv(cov1_inv);
125  const VECTORLIKE1 mu_difs = mu0-mu1;
126  return 0.5*( log(cov1.det()/cov0.det()) + (cov1_inv*cov0).trace() + multiply_HCHt_scalar(mu_difs,cov1_inv) - N );
127  MRPT_END
128  }
129 
130 
131  /** The complementary error function of a Normal distribution
132  */
133 #ifdef HAVE_ERF
134  inline double erfc(double x) { return ::erfc(x); }
135 #else
136  double BASE_IMPEXP erfc(double x);
137 #endif
138 
139  /** The error function of a Normal distribution
140  */
141 #ifdef HAVE_ERF
142  inline double erf(double x) { return ::erf(x); }
143 #else
144  double BASE_IMPEXP erf(double x);
145 #endif
146  /** Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
147  * The employed approximation is that from Peter J. Acklam (pjacklam@online.no),
148  * freely available in http://home.online.no/~pjacklam.
149  */
150  double BASE_IMPEXP normalQuantile(double p);
151 
152  /** Evaluates the Gaussian cumulative density function.
153  * The employed approximation is that from W. J. Cody
154  * freely available in http://www.netlib.org/specfun/erf
155  */
156  double BASE_IMPEXP normalCDF(double p);
157 
158  /** The "quantile" of the Chi-Square distribution, for dimension "dim" and probability 0<P<1 (the inverse of chi2CDF)
159  * An aproximation from the Wilson-Hilferty transformation is used.
160  */
161  double BASE_IMPEXP chi2inv(double P, unsigned int dim=1);
162 
163  /*! Cumulative non-central chi square distribution (approximate).
164 
165  Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom,
166  and noncentrality parameter \a noncentrality at the given argument
167  \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
168  It uses the approximate transform into a normal distribution due to Wilson and Hilferty
169  (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32).
170  The algorithm's running time is independent of the inputs. The accuracy is only
171  about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
172 
173  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
174  */
175  template <class T>
176  double noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg)
177  {
178  const double a = degreesOfFreedom + noncentrality;
179  const double b = (a + noncentrality) / square(a);
180  const double t = (std::pow((double)arg / a, 1.0/3.0) - (1.0 - 2.0 / 9.0 * b)) / std::sqrt(2.0 / 9.0 * b);
181  return 0.5*(1.0 + mrpt::math::erf(t/std::sqrt(2.0)));
182  }
183 
184  /*! Cumulative chi square distribution.
185 
186  Computes the cumulative density of a chi square distribution with \a degreesOfFreedom
187  and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
188  a random number drawn from the distribution is below \a arg
189  by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
190 
191  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
192  */
193  inline double chi2CDF(unsigned int degreesOfFreedom, double arg)
194  {
195  return noncentralChi2CDF(degreesOfFreedom, 0.0, arg);
196  }
197 
198  namespace detail
199  {
200  template <class T>
201  void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j)
202  {
203  double tol = -50.0;
204  if(lans < tol)
205  {
206  lans = lans + std::log(arg / j);
207  dans = std::exp(lans);
208  }
209  else
210  {
211  dans = dans * arg / j;
212  }
213  pans = pans - dans;
214  j += 2;
215  }
216 
217  template <class T>
218  std::pair<double, double> noncentralChi2CDF_exact(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
219  {
220  ASSERTMSG_(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,"noncentralChi2P(): parameters must be positive.");
221  if (arg == 0.0 && degreesOfFreedom > 0)
222  return std::make_pair(0.0, 0.0);
223 
224  // Determine initial values
225  double b1 = 0.5 * noncentrality,
226  ao = std::exp(-b1),
227  eps2 = eps / ao,
228  lnrtpi2 = 0.22579135264473,
229  probability, density, lans, dans, pans, sum, am, hold;
230  unsigned int maxit = 500,
231  i, m;
232  if(degreesOfFreedom % 2)
233  {
234  i = 1;
235  lans = -0.5 * (arg + std::log(arg)) - lnrtpi2;
236  dans = std::exp(lans);
237  pans = erf(std::sqrt(arg/2.0));
238  }
239  else
240  {
241  i = 2;
242  lans = -0.5 * arg;
243  dans = std::exp(lans);
244  pans = 1.0 - dans;
245  }
246 
247  // Evaluate first term
248  if(degreesOfFreedom == 0)
249  {
250  m = 1;
251  degreesOfFreedom = 2;
252  am = b1;
253  sum = 1.0 / ao - 1.0 - am;
254  density = am * dans;
255  probability = 1.0 + am * pans;
256  }
257  else
258  {
259  m = 0;
260  degreesOfFreedom = degreesOfFreedom - 1;
261  am = 1.0;
262  sum = 1.0 / ao - 1.0;
263  while(i < degreesOfFreedom)
264  detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
265  degreesOfFreedom = degreesOfFreedom + 1;
266  density = dans;
267  probability = pans;
268  }
269  // Evaluate successive terms of the expansion
270  for(++m; m<maxit; ++m)
271  {
272  am = b1 * am / m;
273  detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
274  sum = sum - am;
275  density = density + am * dans;
276  hold = am * pans;
277  probability = probability + hold;
278  if((pans * sum < eps2) && (hold < eps2))
279  break; // converged
280  }
281  if(m == maxit)
282  THROW_EXCEPTION("noncentralChi2P(): no convergence.");
283  return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
284  }
285  } // namespace detail
286 
287  /*! Chi square distribution.
288 
289  Computes the density of a chi square distribution with \a degreesOfFreedom
290  and tolerance \a accuracy at the given argument \a arg
291  by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
292 
293  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
294  */
295  inline double chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
296  {
297  return detail::noncentralChi2CDF_exact(degreesOfFreedom, 0.0, arg, accuracy).first;
298  }
299 
300  /** Return the mean and the 10%-90% confidence points (or with any other confidence value) of a set of samples by building the cummulative CDF of all the elements of the container.
301  * The container can be any MRPT container (CArray, matrices, vectors).
302  * \param confidenceInterval A number in the range (0,1) such as the confidence interval will be [100*confidenceInterval, 100*(1-confidenceInterval)].
303  */
304  template <typename CONTAINER>
305  void confidenceIntervals(
306  const CONTAINER &data,
307  typename CONTAINER::value_type &out_mean,
308  typename CONTAINER::value_type &out_lower_conf_interval,
309  typename CONTAINER::value_type &out_upper_conf_interval,
310  const double confidenceInterval = 0.1,
311  const size_t histogramNumBins = 1000 )
312  {
313  MRPT_START
314  ASSERT_(data.size()!=0) // don't use .empty() here to allow using matrices
315  ASSERT_(confidenceInterval>0 && confidenceInterval<1)
316 
317  out_mean = mean(data);
318  typename CONTAINER::value_type x_min,x_max;
319  minimum_maximum(data,x_min,x_max);
320 
321  //std::vector<typename CONTAINER::value_type> xs;
322  //linspace(x_min,x_max,histogramNumBins, xs);
323  const typename CONTAINER::value_type binWidth = (x_max-x_min)/histogramNumBins;
324 
325  const vector_double H = mrpt::math::histogram(data,x_min,x_max,histogramNumBins);
326  vector_double Hc;
327  cumsum(H,Hc); // CDF
328  Hc*=1.0/Hc.maximum();
329 
330  vector_double::iterator it_low = std::lower_bound(Hc.begin(),Hc.end(),confidenceInterval); ASSERT_(it_low!=Hc.end())
331  vector_double::iterator it_high = std::upper_bound(Hc.begin(),Hc.end(),1-confidenceInterval); ASSERT_(it_high!=Hc.end())
332  const size_t idx_low = std::distance(Hc.begin(),it_low);
333  const size_t idx_high = std::distance(Hc.begin(),it_high);
334  out_lower_conf_interval = x_min + idx_low * binWidth;
335  out_upper_conf_interval = x_min + idx_high * binWidth;
336 
337  MRPT_END
338  }
339 
340  /** @} */
341 
342  } // End of MATH namespace
343 
344 } // End of namespace
345 
346 
347 #endif


Page generated by Doxygen 1.8.3 for MRPT 0.9.6 SVN: at Fri Feb 15 22:05:02 EST 2013