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