Main MRPT website > C++ reference for MRPT 1.9.9
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-2018, Individual contributors, see AUTHORS file |
8  +------------------------------------------------------------------------+ */
9 #pragma once
10
11 #include <mrpt/math/math_frwds.h>
13
14 #include <mrpt/math/ops_matrices.h>
15 #include <mrpt/math/ops_vectors.h>
16
17 namespace mrpt
18 {
19 namespace math
20 {
21 /** \addtogroup stats_grp Statistics functions, probability distributions
22  * \ingroup mrpt_math_grp
23  * @{ */
24
25 /** Evaluates the univariate normal (Gaussian) distribution at a given point
26  * "x".
27  */
28 double normalPDF(double x, double mu, double std);
29
30 /** Evaluates the multivariate normal (Gaussian) distribution at a given point
31  * "x".
32  * \param x A vector or column or row matrix with the point at which to
33  * evaluate the pdf.
34  * \param mu A vector or column or row matrix with the Gaussian mean.
35  * \param cov_inv The inverse covariance (information) matrix of the
36  * Gaussian.
37  * \param scaled_pdf If set to true, the PDF will be scaled to be in the
38  * range [0,1], in contrast to its integral from [-inf,+inf] being 1.
39  */
40 template <class VECTORLIKE1, class VECTORLIKE2, class MATRIXLIKE>
42  const VECTORLIKE1& x, const VECTORLIKE2& mu, const MATRIXLIKE& cov_inv,
43  const bool scaled_pdf = false)
44 {
46  using T = typename MATRIXLIKE::Scalar;
47  ASSERTDEB_(cov_inv.isSquare());
48  ASSERTDEB_(
49  size_t(cov_inv.cols()) == size_t(x.size()) &&
50  size_t(cov_inv.cols()) == size_t(mu.size()));
51  T ret = ::exp(
52  static_cast<T>(-0.5) *
53  mrpt::math::multiply_HCHt_scalar((x - mu).eval(), cov_inv));
54  return scaled_pdf
55  ? ret
56  : ret * ::sqrt(
57  cov_inv.det() / ::pow(
58  static_cast<T>(M_2PI),
59  static_cast<T>(cov_inv.rows())));
60  MRPT_END
61 }
62
63 /** Evaluates the multivariate normal (Gaussian) distribution at a given point
64  * "x".
65  * \param x A vector or column or row matrix with the point at which to
66  * evaluate the pdf.
67  * \param mu A vector or column or row matrix with the Gaussian mean.
68  * \param cov The covariance matrix of the Gaussian.
69  * \param scaled_pdf If set to true, the PDF will be scaled to be in the
70  * range [0,1], in contrast to its integral from [-inf,+inf] being 1.
71  */
72 template <class VECTORLIKE1, class VECTORLIKE2, class MATRIXLIKE>
73 inline typename MATRIXLIKE::Scalar normalPDF(
74  const VECTORLIKE1& x, const VECTORLIKE2& mu, const MATRIXLIKE& cov,
75  const bool scaled_pdf = false)
76 {
77  return normalPDFInf(x, mu, cov.inverse(), scaled_pdf);
78 }
79
80 /** Evaluates the multivariate normal (Gaussian) distribution at a given point
81  * given its distance vector "d" from the Gaussian mean.
82  */
83 template <typename VECTORLIKE, typename MATRIXLIKE>
85  const VECTORLIKE& d, const MATRIXLIKE& cov)
86 {
88  ASSERTDEB_(cov.isSquare());
89  ASSERTDEB_(size_t(cov.cols()) == size_t(d.size()));
90  return std::exp(
91  static_cast<typename MATRIXLIKE::Scalar>(-0.5) *
92  mrpt::math::multiply_HCHt_scalar(d, cov.inverse())) /
93  (::pow(
94  static_cast<typename MATRIXLIKE::Scalar>(M_2PI),
95  static_cast<typename MATRIXLIKE::Scalar>(0.5 * cov.cols())) *
96  ::sqrt(cov.det()));
97  MRPT_END
98 }
99
100 /** Kullback-Leibler divergence (KLD) between two independent multivariate
101  * Gaussians.
102  *
103  * \f$D_\mathrm{KL}(\mathcal{N}_0 \| \mathcal{N}_1) = { 1 \over 2 } ( \log_e ( 104 * { \det \Sigma_1 \over \det \Sigma_0 } ) + \mathrm{tr} ( \Sigma_1^{-1} 105 * \Sigma_0 ) + ( \mu_1 - \mu_0 )^\top \Sigma_1^{-1} ( \mu_1 - \mu_0 ) - N ) \f$
106  */
107 template <typename VECTORLIKE1, typename MATRIXLIKE1, typename VECTORLIKE2,
108  typename MATRIXLIKE2>
110  const VECTORLIKE1& mu0, const MATRIXLIKE1& cov0, const VECTORLIKE2& mu1,
111  const MATRIXLIKE2& cov1)
112 {
113  MRPT_START
114  ASSERT_(
115  size_t(mu0.size()) == size_t(mu1.size()) &&
116  size_t(mu0.size()) == size_t(cov0.rows()) &&
117  size_t(mu0.size()) == size_t(cov1.cols()) && cov0.isSquare() &&
118  cov1.isSquare());
119  const size_t N = mu0.size();
120  MATRIXLIKE2 cov1_inv;
121  cov1.inv(cov1_inv);
122  const VECTORLIKE1 mu_difs = mu0 - mu1;
123  return 0.5 * (log(cov1.det() / cov0.det()) + (cov1_inv * cov0).trace() +
124  multiply_HCHt_scalar(mu_difs, cov1_inv) - N);
125  MRPT_END
126 }
127
128 /** Evaluates the Gaussian distribution quantile for the probability value
129  * p=[0,1].
130  * The employed approximation is that from Peter J. Acklam
131  * (pjacklam@online.no),
132  * freely available in http://home.online.no/~pjacklam.
133  */
134 double normalQuantile(double p);
135
136 /** Evaluates the Gaussian cumulative density function.
137  * The employed approximation is that from W. J. Cody
138  * freely available in http://www.netlib.org/specfun/erf
139  * \note Equivalent to MATLAB normcdf(x,mu,s) with p=(x-mu)/s
140  */
141 double normalCDF(double p);
142
143 /** The "quantile" of the Chi-Square distribution, for dimension "dim" and
144  * probability 0<P<1 (the inverse of chi2CDF)
145  * An aproximation from the Wilson-Hilferty transformation is used.
146  * \note Equivalent to MATLAB chi2inv(), but note that this is just an
147  * approximation, which becomes very poor for small values of "P".
148  */
149 double chi2inv(double P, unsigned int dim = 1);
150
151 /*! Cumulative non-central chi square distribution (approximate).
152
153  Computes approximate values of the cumulative density of a chi square
154 distribution with \a degreesOfFreedom,
155  and noncentrality parameter \a noncentrality at the given argument
156  \a arg, i.e. the probability that a random number drawn from the
157 distribution is below \a arg
158  It uses the approximate transform into a normal distribution due to Wilson
159 and Hilferty
160  (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula
161 26.3.32).
162  The algorithm's running time is independent of the inputs. The accuracy is
163 only
164  about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
165
166  \note Function code from the Vigra project
167 (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU
168 GPL-compatible.
169 * \sa noncentralChi2PDF_CDF
170 */
171 double noncentralChi2CDF(
172  unsigned int degreesOfFreedom, double noncentrality, double arg);
173
174 /*! Cumulative chi square distribution.
175
176  Computes the cumulative density of a chi square distribution with \a
177  degreesOfFreedom
178  and tolerance \a accuracy at the given argument \a arg, i.e. the probability
179  that
180  a random number drawn from the distribution is below \a arg
181  by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
182
183  \note Function code from the Vigra project
184  (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU
185  GPL-compatible.
186 */
187 double chi2CDF(unsigned int degreesOfFreedom, double arg);
188
189 /*! Chi square distribution PDF.
190  * Computes the density of a chi square distribution with \a degreesOfFreedom
191  * and tolerance \a accuracy at the given argument \a arg
192  * by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
193  * \note Function code from the Vigra project
194  *(http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU
195  *GPL-compatible.
196  *
197  * \note Equivalent to MATLAB's chi2pdf(arg,degreesOfFreedom)
198  */
199 double chi2PDF(
200  unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7);
201
202 /** Returns the 'exact' PDF (first) and CDF (second) of a Non-central
203  * chi-squared probability distribution, using an iterative method.
204  * \note Equivalent to MATLAB's ncx2cdf(arg,degreesOfFreedom,noncentrality)
205  */
206 std::pair<double, double> noncentralChi2PDF_CDF(
207  unsigned int degreesOfFreedom, double noncentrality, double arg,
208  double eps = 1e-7);
209
210 /** Return the mean and the 10%-90% confidence points (or with any other
211  * confidence value) of a set of samples by building the cummulative CDF of all
212  * the elements of the container.
213  * The container can be any MRPT container (CArray, matrices, vectors).
214  * \param confidenceInterval A number in the range (0,1) such as the confidence
215  * interval will be [100*confidenceInterval, 100*(1-confidenceInterval)].
216  */
217 template <typename CONTAINER>
219  const CONTAINER& data,
222  out_lower_conf_interval,
224  out_upper_conf_interval,
225  const double confidenceInterval = 0.1, const size_t histogramNumBins = 1000)
226 {
227  MRPT_START
228  ASSERT_(
229  data.size() != 0); // don't use .empty() here to allow using matrices
230  ASSERT_(confidenceInterval > 0 && confidenceInterval < 1);
231
232  out_mean = mean(data);
234  minimum_maximum(data, x_min, x_max);
235
236  const typename mrpt::math::ContainerType<CONTAINER>::element_t binWidth =
237  (x_max - x_min) / histogramNumBins;
238
239  const std::vector<double> H =
240  mrpt::math::histogram(data, x_min, x_max, histogramNumBins);
241  std::vector<double> Hc;
242  cumsum(H, Hc); // CDF
243  Hc *= 1.0 / mrpt::math::maximum(Hc);
244
246  std::lower_bound(Hc.begin(), Hc.end(), confidenceInterval);
247  ASSERT_(it_low != Hc.end());
249  std::upper_bound(Hc.begin(), Hc.end(), 1 - confidenceInterval);
250  ASSERT_(it_high != Hc.end());
251  const size_t idx_low = std::distance(Hc.begin(), it_low);
252  const size_t idx_high = std::distance(Hc.begin(), it_high);
253  out_lower_conf_interval = x_min + idx_low * binWidth;
254  out_upper_conf_interval = x_min + idx_high * binWidth;
255
256  MRPT_END
257 }
258
259 /** @} */
260
261 } // End of MATH namespace
262 } // End of namespace
Scalar * iterator
Definition: eigen_plugins.h:26
double Scalar
Definition: KmUtils.h:44
#define MRPT_START
Definition: exceptions.h:262
#define M_2PI
Definition: common.h:58
double normalCDF(double p)
Evaluates the Gaussian cumulative density function.
Definition: math.cpp:133
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:69
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. ...
double normalQuantile(double p)
Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
Definition: math.cpp:80
#define ASSERT_(f)
Defines an assertion mechanism.
Definition: exceptions.h:113
std::pair< double, double > 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:527
double chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Definition: math.cpp:594
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:41
void cumsum(const CONTAINER1 &in_data, CONTAINER2 &out_cumsum)
double chi2CDF(unsigned int degreesOfFreedom, double arg)
Definition: math.cpp:504
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
double 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:42
#define ASSERTDEB_(f)
Defines an assertion mechanism - only when compiled in debug.
Definition: exceptions.h:205
typename CONTAINER::value_type element_t
Definition: math_frwds.h:92
#define MRPT_END
Definition: exceptions.h:266
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=nullptr)
Computes the normalized or normal histogram of a sequence of numbers given the number of bins and the...
double mean(const CONTAINER &v)
Computes the mean value of a vector.
double noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg)
Definition: math.cpp:602
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:148
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:3538
double normalPDF(double x, double mu, double std)
Evaluates the univariate normal (Gaussian) distribution at a given point "x".
Definition: math.cpp:33
GLsizei GLsizei GLenum GLenum const GLvoid * data
Definition: glext.h:3546
GLfloat GLfloat p
Definition: glext.h:6305
double distance(const TPoint2D &p1, const TPoint2D &p2)
Gets the distance between two points in a 2D space.
Definition: geometry.cpp:1891

 Page generated by Doxygen 1.8.14 for MRPT 1.9.9 Git: ad3a9d8ae Tue May 1 23:10:22 2018 -0700 at lun oct 28 00:14:14 CET 2019