Main MRPT website > C++ reference for MRPT 1.3.2
distributions.h
Go to the documentation of this file.
1 /* +---------------------------------------------------------------------------+
2  | Mobile Robot Programming Toolkit (MRPT) |
3  | http://www.mrpt.org/ |
4  | |
5  | Copyright (c) 2005-2015, Individual contributors, see AUTHORS file |
6  | See: http://www.mrpt.org/Authors - All rights reserved. |
7  | Released under BSD License. See details in http://www.mrpt.org/License |
8  +---------------------------------------------------------------------------+ */
9 #ifndef mrpt_math_distributions_H
10 #define mrpt_math_distributions_H
11 
12 #include <algorithm>
13 #include <mrpt/utils/utils_defs.h>
14 #include <mrpt/math/math_frwds.h>
16 
17 #include <mrpt/math/ops_matrices.h>
18 /*---------------------------------------------------------------
19  Namespace
20  ---------------------------------------------------------------*/
21 namespace mrpt
22 {
23  namespace math
24  {
25  /** \addtogroup stats_grp Statistics functions, probability distributions
26  * \ingroup mrpt_base_grp
27  * @{ */
28 
29  /** Evaluates the univariate normal (Gaussian) distribution at a given point "x".
30  */
31  double BASE_IMPEXP normalPDF(double x, double mu, double std);
32 
33  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
34  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
35  * \param mu A vector or column or row matrix with the Gaussian mean.
36  * \param cov_inv The inverse covariance (information) matrix of the Gaussian.
37  * \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.
38  */
39  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
40  inline typename MATRIXLIKE::Scalar
42  const VECTORLIKE1 & x,
43  const VECTORLIKE2 & mu,
44  const MATRIXLIKE & cov_inv,
45  const bool scaled_pdf = false )
46  {
48  typedef typename MATRIXLIKE::Scalar T;
49  ASSERTDEB_(cov_inv.isSquare())
50  ASSERTDEB_(size_t(cov_inv.getColCount())==size_t(x.size()) && size_t(cov_inv.getColCount())==size_t(mu.size()))
51  T ret = ::exp( static_cast<T>(-0.5) * mrpt::math::multiply_HCHt_scalar((x-mu).eval(), cov_inv ) );
52  return scaled_pdf ? ret : ret * ::sqrt(cov_inv.det() / ::pow(static_cast<T>(M_2PI),static_cast<T>( size(cov_inv,1) )) );
53  MRPT_END
54  }
55 
56  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
57  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
58  * \param mu A vector or column or row matrix with the Gaussian mean.
59  * \param cov The covariance matrix of the Gaussian.
60  * \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.
61  */
62  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
63  inline typename MATRIXLIKE::Scalar
65  const VECTORLIKE1 & x,
66  const VECTORLIKE2 & mu,
67  const MATRIXLIKE & cov,
68  const bool scaled_pdf = false )
69  {
70  return normalPDFInf(x,mu,cov.inverse(),scaled_pdf);
71  }
72 
73  /** Evaluates the multivariate normal (Gaussian) distribution at a given point given its distance vector "d" from the Gaussian mean.
74  */
75  template <typename VECTORLIKE,typename MATRIXLIKE>
76  typename MATRIXLIKE::Scalar
77  normalPDF(const VECTORLIKE &d,const MATRIXLIKE &cov)
78  {
80  ASSERTDEB_(cov.isSquare())
81  ASSERTDEB_(size_t(cov.getColCount())==size_t(d.size()))
82  return std::exp( static_cast<typename MATRIXLIKE::Scalar>(-0.5)*mrpt::math::multiply_HCHt_scalar(d,cov.inverse()))
83  / (::pow(
84  static_cast<typename MATRIXLIKE::Scalar>(M_2PI),
85  static_cast<typename MATRIXLIKE::Scalar>(0.5*cov.getColCount()))
86  * ::sqrt(cov.det()));
87  MRPT_END
88  }
89 
90  /** Kullback-Leibler divergence (KLD) between two independent multivariate Gaussians.
91  *
92  * \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$
93  */
94  template <typename VECTORLIKE1,typename MATRIXLIKE1,typename VECTORLIKE2,typename MATRIXLIKE2>
95  double KLD_Gaussians(
96  const VECTORLIKE1 &mu0, const MATRIXLIKE1 &cov0,
97  const VECTORLIKE2 &mu1, const MATRIXLIKE2 &cov1)
98  {
100  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() )
101  const size_t N = mu0.size();
102  MATRIXLIKE2 cov1_inv;
103  cov1.inv(cov1_inv);
104  const VECTORLIKE1 mu_difs = mu0-mu1;
105  return 0.5*( log(cov1.det()/cov0.det()) + (cov1_inv*cov0).trace() + multiply_HCHt_scalar(mu_difs,cov1_inv) - N );
106  MRPT_END
107  }
108 
109 
110  /** The complementary error function of a Normal distribution
111  */
112  double BASE_IMPEXP erfc(const double x);
113 
114  /** The error function of a Normal distribution
115  */
116  double BASE_IMPEXP erf(const double x);
117 
118  /** Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
119  * The employed approximation is that from Peter J. Acklam (pjacklam@online.no),
120  * freely available in http://home.online.no/~pjacklam.
121  */
122  double BASE_IMPEXP normalQuantile(double p);
123 
124  /** Evaluates the Gaussian cumulative density function.
125  * The employed approximation is that from W. J. Cody
126  * freely available in http://www.netlib.org/specfun/erf
127  * \note Equivalent to MATLAB normcdf(x,mu,s) with p=(x-mu)/s
128  */
129  double BASE_IMPEXP normalCDF(double p);
130 
131  /** The "quantile" of the Chi-Square distribution, for dimension "dim" and probability 0<P<1 (the inverse of chi2CDF)
132  * An aproximation from the Wilson-Hilferty transformation is used.
133  * \note Equivalent to MATLAB chi2inv(), but note that this is just an approximation, which becomes very poor for small values of "P".
134  */
135  double BASE_IMPEXP chi2inv(double P, unsigned int dim=1);
136 
137  /*! Cumulative non-central chi square distribution (approximate).
138 
139  Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom,
140  and noncentrality parameter \a noncentrality at the given argument
141  \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
142  It uses the approximate transform into a normal distribution due to Wilson and Hilferty
143  (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32).
144  The algorithm's running time is independent of the inputs. The accuracy is only
145  about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
146 
147  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
148  * \sa noncentralChi2PDF_CDF
149  */
150  double BASE_IMPEXP noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg);
151 
152  /*! Cumulative chi square distribution.
153 
154  Computes the cumulative density of a chi square distribution with \a degreesOfFreedom
155  and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
156  a random number drawn from the distribution is below \a arg
157  by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
158 
159  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
160  */
161  double BASE_IMPEXP chi2CDF(unsigned int degreesOfFreedom, double arg);
162 
163  /*! Chi square distribution PDF.
164  * Computes the density of a chi square distribution with \a degreesOfFreedom
165  * and tolerance \a accuracy at the given argument \a arg
166  * by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
167  * \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
168  *
169  * \note Equivalent to MATLAB's chi2pdf(arg,degreesOfFreedom)
170  */
171  double BASE_IMPEXP chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7);
172 
173  /** Returns the 'exact' PDF (first) and CDF (second) of a Non-central chi-squared probability distribution, using an iterative method.
174  * \note Equivalent to MATLAB's ncx2cdf(arg,degreesOfFreedom,noncentrality)
175  */
176  std::pair<double, double> BASE_IMPEXP noncentralChi2PDF_CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double eps = 1e-7);
177 
178  /** 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.
179  * The container can be any MRPT container (CArray, matrices, vectors).
180  * \param confidenceInterval A number in the range (0,1) such as the confidence interval will be [100*confidenceInterval, 100*(1-confidenceInterval)].
181  */
182  template <typename CONTAINER>
184  const CONTAINER &data,
186  typename mrpt::math::ContainerType<CONTAINER>::element_t &out_lower_conf_interval,
187  typename mrpt::math::ContainerType<CONTAINER>::element_t &out_upper_conf_interval,
188  const double confidenceInterval = 0.1,
189  const size_t histogramNumBins = 1000 )
190  {
191  MRPT_START
192  ASSERT_(data.size()!=0) // don't use .empty() here to allow using matrices
193  ASSERT_(confidenceInterval>0 && confidenceInterval<1)
194 
195  out_mean = mean(data);
197  minimum_maximum(data,x_min,x_max);
198 
199  const typename mrpt::math::ContainerType<CONTAINER>::element_t binWidth = (x_max-x_min)/histogramNumBins;
200 
201  const std::vector<double> H = mrpt::math::histogram(data,x_min,x_max,histogramNumBins);
202  std::vector<double> Hc;
203  cumsum(H,Hc); // CDF
204  double factor = 1.0/mrpt::math::maximum(Hc);
205  std::transform(Hc.begin(), Hc.end(), Hc.begin(), std::bind1st(std::multiplies<double>(), factor));
206 
207  std::vector<double>::iterator it_low = std::lower_bound(Hc.begin(),Hc.end(),confidenceInterval); ASSERT_(it_low!=Hc.end())
208  std::vector<double>::iterator it_high = std::upper_bound(Hc.begin(),Hc.end(),1-confidenceInterval); ASSERT_(it_high!=Hc.end())
209  const size_t idx_low = std::distance(Hc.begin(),it_low);
210  const size_t idx_high = std::distance(Hc.begin(),it_high);
211  out_lower_conf_interval = x_min + idx_low * binWidth;
212  out_upper_conf_interval = x_min + idx_high * binWidth;
213 
214  MRPT_END
215  }
216 
217  /** @} */
218 
219  } // End of MATH namespace
220 
221 } // End of namespace
222 
223 
224 #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].
double BASE_IMPEXP normalCDF(double p)
Evaluates the Gaussian cumulative density function.
std::vector< double > histogram(const CONTAINER &v, double limit_min, double limit_max, size_t number_bins, bool do_normalization=false, std::vector< double > *out_bin_centers=NULL)
Computes the normalized or normal histogram of a sequence of numbers given the number of bins and the...
Scalar * iterator
Definition: eigen_plugins.h:23
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, and internal functions in mrpt...
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:62
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. ...
Definition: distributions.h:95
#define M_2PI
Definition: mrpt_macros.h:363
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...
double BASE_IMPEXP chi2CDF(unsigned int degreesOfFreedom, double arg)
#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:135
CONTAINER::Scalar maximum(const CONTAINER &v)
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:41
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:38
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)
CONTAINER::value_type element_t
Definition: math_frwds.h:85
void confidenceIntervals(const CONTAINER &data, typename mrpt::math::ContainerType< CONTAINER >::element_t &out_mean, typename mrpt::math::ContainerType< CONTAINER >::element_t &out_lower_conf_interval, typename mrpt::math::ContainerType< CONTAINER >::element_t &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 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.12 for MRPT 1.3.2 SVN: at Thu Nov 10 13:46:27 UTC 2016