Main MRPT website > C++ reference for MRPT 1.5.7
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-2017, 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 <mrpt/utils/utils_defs.h>
13 #include <mrpt/math/math_frwds.h>
15 
16 #include <mrpt/math/ops_matrices.h>
17 #include <mrpt/math/ops_vectors.h>
18 
19 /*---------------------------------------------------------------
20  Namespace
21  ---------------------------------------------------------------*/
22 namespace mrpt
23 {
24  namespace math
25  {
26  /** \addtogroup stats_grp Statistics functions, probability distributions
27  * \ingroup mrpt_base_grp
28  * @{ */
29 
30  /** Evaluates the univariate normal (Gaussian) distribution at a given point "x".
31  */
32  double BASE_IMPEXP normalPDF(double x, double mu, double std);
33 
34  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
35  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
36  * \param mu A vector or column or row matrix with the Gaussian mean.
37  * \param cov_inv The inverse covariance (information) matrix of the Gaussian.
38  * \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.
39  */
40  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
41  inline typename MATRIXLIKE::Scalar
43  const VECTORLIKE1 & x,
44  const VECTORLIKE2 & mu,
45  const MATRIXLIKE & cov_inv,
46  const bool scaled_pdf = false )
47  {
49  typedef typename MATRIXLIKE::Scalar T;
50  ASSERTDEB_(cov_inv.isSquare())
51  ASSERTDEB_(size_t(cov_inv.getColCount())==size_t(x.size()) && size_t(cov_inv.getColCount())==size_t(mu.size()))
52  T ret = ::exp( static_cast<T>(-0.5) * mrpt::math::multiply_HCHt_scalar((x-mu).eval(), cov_inv ) );
53  return scaled_pdf ? ret : ret * ::sqrt(cov_inv.det() / ::pow(static_cast<T>(M_2PI),static_cast<T>( size(cov_inv,1) )) );
54  MRPT_END
55  }
56 
57  /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
58  * \param x A vector or column or row matrix with the point at which to evaluate the pdf.
59  * \param mu A vector or column or row matrix with the Gaussian mean.
60  * \param cov The covariance matrix of the Gaussian.
61  * \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.
62  */
63  template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
64  inline typename MATRIXLIKE::Scalar
66  const VECTORLIKE1 & x,
67  const VECTORLIKE2 & mu,
68  const MATRIXLIKE & cov,
69  const bool scaled_pdf = false )
70  {
71  return normalPDFInf(x,mu,cov.inverse(),scaled_pdf);
72  }
73 
74  /** Evaluates the multivariate normal (Gaussian) distribution at a given point given its distance vector "d" from the Gaussian mean.
75  */
76  template <typename VECTORLIKE,typename MATRIXLIKE>
77  typename MATRIXLIKE::Scalar
78  normalPDF(const VECTORLIKE &d,const MATRIXLIKE &cov)
79  {
81  ASSERTDEB_(cov.isSquare())
82  ASSERTDEB_(size_t(cov.getColCount())==size_t(d.size()))
83  return std::exp( static_cast<typename MATRIXLIKE::Scalar>(-0.5)*mrpt::math::multiply_HCHt_scalar(d,cov.inverse()))
84  / (::pow(
85  static_cast<typename MATRIXLIKE::Scalar>(M_2PI),
86  static_cast<typename MATRIXLIKE::Scalar>(0.5*cov.getColCount()))
87  * ::sqrt(cov.det()));
88  MRPT_END
89  }
90 
91  /** Kullback-Leibler divergence (KLD) between two independent multivariate Gaussians.
92  *
93  * \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$
94  */
95  template <typename VECTORLIKE1,typename MATRIXLIKE1,typename VECTORLIKE2,typename MATRIXLIKE2>
96  double KLD_Gaussians(
97  const VECTORLIKE1 &mu0, const MATRIXLIKE1 &cov0,
98  const VECTORLIKE2 &mu1, const MATRIXLIKE2 &cov1)
99  {
100  MRPT_START
101  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() )
102  const size_t N = mu0.size();
103  MATRIXLIKE2 cov1_inv;
104  cov1.inv(cov1_inv);
105  const VECTORLIKE1 mu_difs = mu0-mu1;
106  return 0.5*( log(cov1.det()/cov0.det()) + (cov1_inv*cov0).trace() + multiply_HCHt_scalar(mu_difs,cov1_inv) - N );
107  MRPT_END
108  }
109 
110 
111  /** The complementary error function of a Normal distribution
112  */
113  double BASE_IMPEXP erfc(const double x);
114 
115  /** The error function of a Normal distribution
116  */
117  double BASE_IMPEXP erf(const double x);
118 
119  /** Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
120  * The employed approximation is that from Peter J. Acklam (pjacklam@online.no),
121  * freely available in http://home.online.no/~pjacklam.
122  */
123  double BASE_IMPEXP normalQuantile(double p);
124 
125  /** Evaluates the Gaussian cumulative density function.
126  * The employed approximation is that from W. J. Cody
127  * freely available in http://www.netlib.org/specfun/erf
128  * \note Equivalent to MATLAB normcdf(x,mu,s) with p=(x-mu)/s
129  */
130  double BASE_IMPEXP normalCDF(double p);
131 
132  /** The "quantile" of the Chi-Square distribution, for dimension "dim" and probability 0<P<1 (the inverse of chi2CDF)
133  * An aproximation from the Wilson-Hilferty transformation is used.
134  * \note Equivalent to MATLAB chi2inv(), but note that this is just an approximation, which becomes very poor for small values of "P".
135  */
136  double BASE_IMPEXP chi2inv(double P, unsigned int dim=1);
137 
138  /*! Cumulative non-central chi square distribution (approximate).
139 
140  Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom,
141  and noncentrality parameter \a noncentrality at the given argument
142  \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
143  It uses the approximate transform into a normal distribution due to Wilson and Hilferty
144  (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32).
145  The algorithm's running time is independent of the inputs. The accuracy is only
146  about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
147 
148  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
149  * \sa noncentralChi2PDF_CDF
150  */
151  double BASE_IMPEXP noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg);
152 
153  /*! Cumulative chi square distribution.
154 
155  Computes the cumulative density of a chi square distribution with \a degreesOfFreedom
156  and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
157  a random number drawn from the distribution is below \a arg
158  by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
159 
160  \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
161  */
162  double BASE_IMPEXP chi2CDF(unsigned int degreesOfFreedom, double arg);
163 
164  /*! Chi square distribution PDF.
165  * Computes the density of a chi square distribution with \a degreesOfFreedom
166  * and tolerance \a accuracy at the given argument \a arg
167  * by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
168  * \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
169  *
170  * \note Equivalent to MATLAB's chi2pdf(arg,degreesOfFreedom)
171  */
172  double BASE_IMPEXP chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7);
173 
174  /** Returns the 'exact' PDF (first) and CDF (second) of a Non-central chi-squared probability distribution, using an iterative method.
175  * \note Equivalent to MATLAB's ncx2cdf(arg,degreesOfFreedom,noncentrality)
176  */
177  std::pair<double, double> BASE_IMPEXP noncentralChi2PDF_CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double eps = 1e-7);
178 
179  /** 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.
180  * The container can be any MRPT container (CArray, matrices, vectors).
181  * \param confidenceInterval A number in the range (0,1) such as the confidence interval will be [100*confidenceInterval, 100*(1-confidenceInterval)].
182  */
183  template <typename CONTAINER>
185  const CONTAINER &data,
187  typename mrpt::math::ContainerType<CONTAINER>::element_t &out_lower_conf_interval,
188  typename mrpt::math::ContainerType<CONTAINER>::element_t &out_upper_conf_interval,
189  const double confidenceInterval = 0.1,
190  const size_t histogramNumBins = 1000 )
191  {
192  MRPT_START
193  ASSERT_(data.size()!=0) // don't use .empty() here to allow using matrices
194  ASSERT_(confidenceInterval>0 && confidenceInterval<1)
195 
196  out_mean = mean(data);
198  minimum_maximum(data,x_min,x_max);
199 
200  const typename mrpt::math::ContainerType<CONTAINER>::element_t binWidth = (x_max-x_min)/histogramNumBins;
201 
202  const std::vector<double> H = mrpt::math::histogram(data,x_min,x_max,histogramNumBins);
203  std::vector<double> Hc;
204  cumsum(H,Hc); // CDF
205  Hc*=1.0/ mrpt::math::maximum(Hc);
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 chi2CDF(unsigned int degreesOfFreedom, double arg)
Definition: math.cpp:2098
size_t size(const MATRIXLIKE &m, const int dim)
Definition: bits.h:43
double BASE_IMPEXP normalCDF(double p)
Evaluates the Gaussian cumulative density function.
Definition: math.cpp:1108
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
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
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:96
#define M_2PI
Definition: mrpt_macros.h:380
double BASE_IMPEXP normalQuantile(double p)
Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
Definition: math.cpp:1046
double BASE_IMPEXP erfc(const double x)
The complementary error function of a Normal distribution.
Definition: math.cpp:917
#define MRPT_END
double BASE_IMPEXP erf(const double x)
The error function of a Normal distribution.
Definition: math.cpp:955
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
const double eps
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:42
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...
Definition: math.cpp:967
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...
Definition: math.cpp:2120
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...
GLenum GLint x
Definition: glext.h:3516
double BASE_IMPEXP normalPDF(double x, double mu, double std)
Evaluates the univariate normal (Gaussian) distribution at a given point "x".
Definition: math.cpp:908
GLsizei GLsizei GLenum GLenum const GLvoid * data
Definition: glext.h:3520
double Scalar
Definition: KmUtils.h:41
GLfloat GLfloat p
Definition: glext.h:5587
double BASE_IMPEXP chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Definition: math.cpp:2189
double BASE_IMPEXP distance(const TPoint2D &p1, const TPoint2D &p2)
Gets the distance between two points in a 2D space.
Definition: geometry.cpp:1504
double BASE_IMPEXP noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg)
Definition: math.cpp:2194



Page generated by Doxygen 1.8.14 for MRPT 1.5.7 Git: 5902e14cc Wed Apr 24 15:04:01 2019 +0200 at lun oct 28 01:39:17 CET 2019