Main MRPT website > C++ reference for MRPT 1.4.0
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-2016, 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
15
18
19/*---------------------------------------------------------------
20 Namespace
21 ---------------------------------------------------------------*/
22namespace 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) )) );
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()));
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>
97 const VECTORLIKE1 &mu0, const MATRIXLIKE1 &cov0,
98 const VECTORLIKE2 &mu1, const MATRIXLIKE2 &cov1)
99 {
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 );
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 {
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
215 }
216
217 /** @} */
218
219 } // End of MATH namespace
220
221} // End of namespace
222
223
224#endif
Scalar * iterator
Definition: eigen_plugins.h:23
double BASE_IMPEXP distance(const TPoint2D &p1, const TPoint2D &p2)
Gets the distance between two points in a 2D space.
double BASE_IMPEXP normalCDF(double p)
Evaluates the Gaussian cumulative density function.
double BASE_IMPEXP chi2CDF(unsigned int degreesOfFreedom, double arg)
std::pair< double, double > BASE_IMPEXP noncentralChi2PDF_CDF(unsigned int degreesOfFreedom, double noncentrality, double arg, double eps=1e-7)
Returns the 'exact' PDF (first) and CDF (second) of a Non-central chi-squared probability distributio...
double BASE_IMPEXP normalPDF(double x, double mu, double std)
Evaluates the univariate normal (Gaussian) distribution at a given point "x".
double BASE_IMPEXP erf(const double x)
The error function of a Normal distribution.
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 chi2inv(double P, unsigned int dim=1)
The "quantile" of the Chi-Square distribution, for dimension "dim" and probability 0<P<1 (the inverse...
double BASE_IMPEXP noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg)
double BASE_IMPEXP erfc(const double x)
The complementary error function of a Normal distribution.
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
double BASE_IMPEXP chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
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
double BASE_IMPEXP normalQuantile(double p)
Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
#define MRPT_START
Definition: mrpt_macros.h:349
#define ASSERT_(f)
Definition: mrpt_macros.h:261
#define M_2PI
Definition: mrpt_macros.h:363
#define ASSERTDEB_(f)
Defines an assertion mechanism - only when compiled in debug.
Definition: mrpt_macros.h:283
#define MRPT_END
Definition: mrpt_macros.h:353
double mean(const CONTAINER &v)
Computes the mean value of a vector.
size_t size(const MATRIXLIKE &m, int dim)
Definition: bits.h:38
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...
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
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
void minimum_maximum(const std::vector< T > &V, T &curMin, T &curMax)
Return the maximum and minimum values of a std::vector.
void cumsum(const CONTAINER1 &in_data, CONTAINER2 &out_cumsum)
CONTAINER::Scalar maximum(const CONTAINER &v)
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
STL namespace.
This file implements miscelaneous matrix and matrix/vector operations, and internal functions in mrpt...
CONTAINER::value_type element_t
Definition: math_frwds.h:85



Page generated by Doxygen 1.9.2 for MRPT 1.4.0 SVN: at Mon Sep 20 00:47:55 UTC 2021