Main MRPT website > C++ reference
MRPT logo
distributions.h
Go to the documentation of this file.
1 /* +---------------------------------------------------------------------------+
2  | The Mobile Robot Programming Toolkit (MRPT) |
3  | |
4  | http://www.mrpt.org/ |
5  | |
6  | Copyright (c) 2005-2013, Individual contributors, see AUTHORS file |
7  | Copyright (c) 2005-2013, MAPIR group, University of Malaga |
8  | Copyright (c) 2012-2013, University of Almeria |
9  | All rights reserved. |
10  | |
11  | Redistribution and use in source and binary forms, with or without |
12  | modification, are permitted provided that the following conditions are |
13  | met: |
14  | * Redistributions of source code must retain the above copyright |
15  | notice, this list of conditions and the following disclaimer. |
16  | * Redistributions in binary form must reproduce the above copyright |
17  | notice, this list of conditions and the following disclaimer in the |
18  | documentation and/or other materials provided with the distribution. |
19  | * Neither the name of the copyright holders nor the |
20  | names of its contributors may be used to endorse or promote products |
21  | derived from this software without specific prior written permission.|
22  | |
23  | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
24  | 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED |
25  | TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR|
26  | PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE |
27  | FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL|
28  | DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR|
29  | SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
30  | HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, |
31  | STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN |
32  | ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
33  | POSSIBILITY OF SUCH DAMAGE. |
34  +---------------------------------------------------------------------------+ */
35 #ifndef mrpt_math_distributions_H
36 #define mrpt_math_distributions_H
37 
38 #include <mrpt/utils/utils_defs.h>
39 #include <mrpt/math/math_frwds.h>
41 
42 #include <mrpt/math/ops_matrices.h>
43 
44 /*---------------------------------------------------------------
45  Namespace
46  ---------------------------------------------------------------*/
47 namespace mrpt
48 {
49  namespace math
50  {
51  using namespace mrpt::utils;
52 
53  /** \addtogroup stats_grp Statistics functions, probability distributions
54  * \ingroup mrpt_base_grp
55  * @{ */
56 
57  /** Evaluates the univariate normal (Gaussian) distribution at a given point "x".
58  */
59  double BASE_IMPEXP normalPDF(double x, double mu, double std);
60 
61  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
62  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
63  * \param mu A vector or column or row matrix with the Gaussian mean.
64  * \param cov_inv The inverse covariance (information) matrix of the Gaussian.
65  * \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.
66  */
67  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
68  inline typename MATRIXLIKE::Scalar
70  const VECTORLIKE1 & x,
71  const VECTORLIKE2 & mu,
72  const MATRIXLIKE & cov_inv,
73  const bool scaled_pdf = false )
74  {
76  typedef typename MATRIXLIKE::Scalar T;
77  ASSERTDEB_(cov_inv.isSquare())
78  ASSERTDEB_(size_t(cov_inv.getColCount())==size_t(x.size()) && size_t(cov_inv.getColCount())==size_t(mu.size()))
79  T ret = ::exp( static_cast<T>(-0.5) * mrpt::math::multiply_HCHt_scalar((x-mu).eval(), cov_inv ) );
80  return scaled_pdf ? ret : ret * ::sqrt(cov_inv.det() / ::pow(static_cast<T>(M_2PI),static_cast<T>( size(cov_inv,1) )) );
81  MRPT_END
82  }
83 
84  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
85  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
86  * \param mu A vector or column or row matrix with the Gaussian mean.
87  * \param cov The covariance matrix of the Gaussian.
88  * \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.
89  */
90  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
91  inline typename MATRIXLIKE::Scalar
93  const VECTORLIKE1 & x,
94  const VECTORLIKE2 & mu,
95  const MATRIXLIKE & cov,
96  const bool scaled_pdf = false )
97  {
98  return normalPDFInf(x,mu,cov.inverse(),scaled_pdf);
99  }
100 
101  /** Evaluates the multivariate normal (Gaussian) distribution at a given point given its distance vector "d" from the Gaussian mean.
102  */
103  template <typename VECTORLIKE,typename MATRIXLIKE>
104  typename MATRIXLIKE::Scalar
105  normalPDF(const VECTORLIKE &d,const MATRIXLIKE &cov)
106  {
107  MRPT_START
108  ASSERTDEB_(cov.isSquare())
109  ASSERTDEB_(size_t(cov.getColCount())==size_t(d.size()))
110  return std::exp( static_cast<typename MATRIXLIKE::Scalar>(-0.5)*mrpt::math::multiply_HCHt_scalar(d,cov.inverse()))
111  / (::pow(
112  static_cast<typename MATRIXLIKE::Scalar>(M_2PI),
113  static_cast<typename MATRIXLIKE::Scalar>(0.5*cov.getColCount()))
114  * ::sqrt(cov.det()));
115  MRPT_END
116  }
117 
118  /** Kullback-Leibler divergence (KLD) between two independent multivariate Gaussians.
119  *
120  * \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$
121  */
122  template <typename VECTORLIKE1,typename MATRIXLIKE1,typename VECTORLIKE2,typename MATRIXLIKE2>
124  const VECTORLIKE1 &mu0, const MATRIXLIKE1 &cov0,
125  const VECTORLIKE2 &mu1, const MATRIXLIKE2 &cov1)
126  {
127  MRPT_START
128  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() )
129  const size_t N = mu0.size();
130  MATRIXLIKE2 cov1_inv;
131  cov1.inv(cov1_inv);
132  const VECTORLIKE1 mu_difs = mu0-mu1;
133  return 0.5*( log(cov1.det()/cov0.det()) + (cov1_inv*cov0).trace() + multiply_HCHt_scalar(mu_difs,cov1_inv) - N );
134  MRPT_END
135  }
136 
137 
138  /** The complementary error function of a Normal distribution
139  */
140  double BASE_IMPEXP erfc(const double x);
141 
142  /** The error function of a Normal distribution
143  */
144  double BASE_IMPEXP erf(const double x);
145 
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  * \note Equivalent to MATLAB normcdf(x,mu,s) with p=(x-mu)/s
156  */
157  double BASE_IMPEXP normalCDF(double p);
158 
159  /** The "quantile" of the Chi-Square distribution, for dimension "dim" and probability 0<P<1 (the inverse of chi2CDF)
160  * An aproximation from the Wilson-Hilferty transformation is used.
161  * \note Equivalent to MATLAB chi2inv(), but note that this is just an approximation, which becomes very poor for small values of "P".
162  */
163  double BASE_IMPEXP chi2inv(double P, unsigned int dim=1);
164 
165  /*! Cumulative non-central chi square distribution (approximate).
166 
167  Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom,
168  and noncentrality parameter \a noncentrality at the given argument
169  \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
170  It uses the approximate transform into a normal distribution due to Wilson and Hilferty
171  (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32).
172  The algorithm's running time is independent of the inputs. The accuracy is only
173  about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
174 
175  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
176  * \sa noncentralChi2PDF_CDF
177  */
178  double BASE_IMPEXP noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg);
179 
180  /*! Cumulative chi square distribution.
181 
182  Computes the cumulative density of a chi square distribution with \a degreesOfFreedom
183  and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
184  a random number drawn from the distribution is below \a arg
185  by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
186 
187  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
188  */
189  double BASE_IMPEXP chi2CDF(unsigned int degreesOfFreedom, double arg);
190 
191  /*! Chi square distribution PDF.
192  * Computes the density of a chi square distribution with \a degreesOfFreedom
193  * and tolerance \a accuracy at the given argument \a arg
194  * by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
195  * \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
196  *
197  * \note Equivalent to MATLAB's chi2pdf(arg,degreesOfFreedom)
198  */
199  double BASE_IMPEXP chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7);
200 
201  /** Returns the 'exact' PDF (first) and CDF (second) of a Non-central chi-squared probability distribution, using an iterative method.
202  * \note Equivalent to MATLAB's ncx2cdf(arg,degreesOfFreedom,noncentrality)
203  */
204  std::pair<double, double> BASE_IMPEXP noncentralChi2PDF_CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double eps = 1e-7);
205 
206  /** 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.
207  * The container can be any MRPT container (CArray, matrices, vectors).
208  * \param confidenceInterval A number in the range (0,1) such as the confidence interval will be [100*confidenceInterval, 100*(1-confidenceInterval)].
209  */
210  template <typename CONTAINER>
212  const CONTAINER &data,
213  typename CONTAINER::value_type &out_mean,
214  typename CONTAINER::value_type &out_lower_conf_interval,
215  typename CONTAINER::value_type &out_upper_conf_interval,
216  const double confidenceInterval = 0.1,
217  const size_t histogramNumBins = 1000 )
218  {
219  MRPT_START
220  ASSERT_(data.size()!=0) // don't use .empty() here to allow using matrices
221  ASSERT_(confidenceInterval>0 && confidenceInterval<1)
222 
223  out_mean = mean(data);
224  typename CONTAINER::value_type x_min,x_max;
225  minimum_maximum(data,x_min,x_max);
226 
227  const typename CONTAINER::value_type binWidth = (x_max-x_min)/histogramNumBins;
228 
229  const vector_double H = mrpt::math::histogram(data,x_min,x_max,histogramNumBins);
230  vector_double Hc;
231  cumsum(H,Hc); // CDF
232  Hc*=1.0/Hc.maximum();
233 
234  vector_double::iterator it_low = std::lower_bound(Hc.begin(),Hc.end(),confidenceInterval); ASSERT_(it_low!=Hc.end())
235  vector_double::iterator it_high = std::upper_bound(Hc.begin(),Hc.end(),1-confidenceInterval); ASSERT_(it_high!=Hc.end())
236  const size_t idx_low = std::distance(Hc.begin(),it_low);
237  const size_t idx_high = std::distance(Hc.begin(),it_high);
238  out_lower_conf_interval = x_min + idx_low * binWidth;
239  out_upper_conf_interval = x_min + idx_high * binWidth;
240 
241  MRPT_END
242  }
243 
244  /** @} */
245 
246  } // End of MATH namespace
247 
248 } // End of namespace
249 
250 
251 #endif
double BASE_IMPEXP erfc(const double x)
The complementary error function of a Normal distribution.
double BASE_IMPEXP normalQuantile(double p)
Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
Classes for serialization, sockets, ini-file manipulation, streams, list of properties-values, timewatch, extensions to STL.
Definition: zip.h:42
double BASE_IMPEXP normalCDF(double p)
Evaluates the Gaussian cumulative density function.
Scalar * iterator
Definition: eigen_plugins.h:49
double BASE_IMPEXP normalPDF(double x, double mu, double std)
Evaluates the univariate normal (Gaussian) distribution at a given point "x".
This file implements miscelaneous matrix and matrix/vector operations, plus internal functions in mrp...
MAT_C::Scalar multiply_HCHt_scalar(const VECTOR_H &H, const MAT_C &C)
r (a scalar) = H * C * H^t (with a vector H and a symmetric matrix C)
Definition: ops_matrices.h:167
double BASE_IMPEXP noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg)
STL namespace.
double KLD_Gaussians(const VECTORLIKE1 &mu0, const MATRIXLIKE1 &cov0, const VECTORLIKE2 &mu1, const MATRIXLIKE2 &cov1)
Kullback-Leibler divergence (KLD) between two independent multivariate Gaussians. ...
#define M_2PI
Definition: mrpt_macros.h:412
std::pair< double, double > BASE_IMPEXP noncentralChi2PDF_CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double eps=1e-7)
Returns the &#39;exact&#39; PDF (first) and CDF (second) of a Non-central chi-squared probability distributio...
vector_double histogram(const CONTAINER &v, double limit_min, double limit_max, size_t number_bins, bool do_normalization=false, vector_double *out_bin_centers=NULL)
Computes the normalized or normal histogram of a sequence of numbers given the number of bins and the...
double BASE_IMPEXP chi2CDF(unsigned int degreesOfFreedom, double arg)
The base class of MRPT vectors, actually, Eigen column matrices of dynamic size with specialized cons...
#define MRPT_END
Eigen::Matrix< typename MATRIX::Scalar, MATRIX::ColsAtCompileTime, MATRIX::ColsAtCompileTime > cov(const MATRIX &v)
Computes the covariance matrix from a list of samples in an NxM matrix, where each row is a sample...
Definition: ops_matrices.h:240
void minimum_maximum(const std::vector< T > &V, T &curMin, T &curMax)
Return the maximum and minimum values of a std::vector.
MATRIXLIKE::Scalar normalPDFInf(const VECTORLIKE1 &x, const VECTORLIKE2 &mu, const MATRIXLIKE &cov_inv, const bool scaled_pdf=false)
Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
Definition: distributions.h:69
void cumsum(const CONTAINER1 &in_data, CONTAINER2 &out_cumsum)
#define MRPT_START
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
#define ASSERTDEB_(f)
Defines an assertion mechanism - only when compiled in debug.
double BASE_IMPEXP chi2inv(double P, unsigned int dim=1)
The "quantile" of the Chi-Square distribution, for dimension "dim" and probability 0<P<1 (the inverse...
size_t size(const MATRIXLIKE &m, int dim)
Definition: bits.h:72
void confidenceIntervals(const CONTAINER &data, typename CONTAINER::value_type &out_mean, typename CONTAINER::value_type &out_lower_conf_interval, typename CONTAINER::value_type &out_upper_conf_interval, const double confidenceInterval=0.1, const size_t histogramNumBins=1000)
Return the mean and the 10%-90% confidence points (or with any other confidence value) of a set of sa...
double BASE_IMPEXP chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
double mean(const CONTAINER &v)
Computes the mean value of a vector.
#define ASSERT_(f)
double BASE_IMPEXP erf(const double x)
The error function of a Normal distribution.
double BASE_IMPEXP distance(const TPoint2D &p1, const TPoint2D &p2)
Gets the distance between two points in a 2D space.



Page generated by Doxygen 1.8.14 for MRPT 1.0.2 SVN: at lun oct 28 00:52:41 CET 2019 Hosted on:
SourceForge.net Logo