Main MRPT website > C++ reference
MRPT logo

distributions.h

Go to the documentation of this file.
00001 /* +---------------------------------------------------------------------------+
00002    |          The Mobile Robot Programming Toolkit (MRPT) C++ library          |
00003    |                                                                           |
00004    |                   http://mrpt.sourceforge.net/                            |
00005    |                                                                           |
00006    |   Copyright (C) 2005-2011  University of Malaga                           |
00007    |                                                                           |
00008    |    This software was written by the Machine Perception and Intelligent    |
00009    |      Robotics Lab, University of Malaga (Spain).                          |
00010    |    Contact: Jose-Luis Blanco  <jlblanco@ctima.uma.es>                     |
00011    |                                                                           |
00012    |  This file is part of the MRPT project.                                   |
00013    |                                                                           |
00014    |     MRPT is free software: you can redistribute it and/or modify          |
00015    |     it under the terms of the GNU General Public License as published by  |
00016    |     the Free Software Foundation, either version 3 of the License, or     |
00017    |     (at your option) any later version.                                   |
00018    |                                                                           |
00019    |   MRPT is distributed in the hope that it will be useful,                 |
00020    |     but WITHOUT ANY WARRANTY; without even the implied warranty of        |
00021    |     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         |
00022    |     GNU General Public License for more details.                          |
00023    |                                                                           |
00024    |     You should have received a copy of the GNU General Public License     |
00025    |     along with MRPT.  If not, see <http://www.gnu.org/licenses/>.         |
00026    |                                                                           |
00027    +---------------------------------------------------------------------------+ */
00028 #ifndef  mrpt_math_distributions_H
00029 #define  mrpt_math_distributions_H
00030 
00031 #include <mrpt/utils/utils_defs.h>
00032 #include <mrpt/math/math_frwds.h>
00033 #include <mrpt/math/CMatrixTemplateNumeric.h>
00034 
00035 #include <mrpt/math/ops_matrices.h>
00036 
00037 /*---------------------------------------------------------------
00038                 Namespace
00039   ---------------------------------------------------------------*/
00040 namespace mrpt
00041 {
00042         namespace math
00043         {
00044                 using namespace mrpt::utils;
00045 
00046                 /** @name Statistics functions
00047                 @{ */
00048 
00049                 /** Evaluates the univariate normal (Gaussian) distribution at a given point "x".
00050                   */
00051                 double BASE_IMPEXP  normalPDF(double x, double mu, double std);
00052 
00053                 /** Evaluates the multivariate normal (Gaussian) distribution at a given point "x".
00054                   *  \param  x   A vector or column or row matrix with the point at which to evaluate the pdf.
00055                   *  \param  mu  A vector or column or row matrix with the Gaussian mean.
00056                   *  \param  cov  The covariance matrix of the Gaussian.
00057                   *  \param  scaled_pdf If set to true, the PDF will be scaled to be in the range [0,1].
00058                   */
00059                 template <class VECTORLIKE1,class VECTORLIKE2,class MATRIXLIKE>
00060                 inline typename MATRIXLIKE::value_type
00061                         normalPDF(
00062                                 const VECTORLIKE1  & x,
00063                                 const VECTORLIKE2  & mu,
00064                                 const MATRIXLIKE   & cov,
00065                                 const bool scaled_pdf = false )
00066                 {
00067                         MRPT_START
00068                         typedef typename MATRIXLIKE::value_type T;
00069                         ASSERTDEB_(cov.isSquare())
00070                         ASSERTDEB_(size_t(cov.getColCount())==size_t(x.size()) && size_t(cov.getColCount())==size_t(mu.size()))
00071                         T ret = ::exp( static_cast<T>(-0.5) * mrpt::math::multiply_HCHt_scalar((x-mu), cov.inverse() ) );
00072                         return scaled_pdf ? ret : ret / (::pow(static_cast<T>(M_2PI),static_cast<T>( size(cov,1) )) * ::sqrt(cov.det()));
00073                         MRPT_END
00074                 }
00075 
00076                 /** Evaluates the multivariate normal (Gaussian) distribution at a given point given its distance vector "d" from the Gaussian mean.
00077                   */
00078                 template <typename VECTORLIKE,typename MATRIXLIKE>
00079                 typename MATRIXLIKE::value_type
00080                 normalPDF(const VECTORLIKE &d,const MATRIXLIKE &cov)
00081                 {
00082                         MRPT_START
00083                         ASSERTDEB_(cov.isSquare())
00084                         ASSERTDEB_(size_t(cov.getColCount())==size_t(d.size()))
00085                         return std::exp( static_cast<typename MATRIXLIKE::value_type>(-0.5)*mrpt::math::multiply_HCHt_scalar(d,cov.inverse()))
00086                         / (::pow(
00087                                         static_cast<typename MATRIXLIKE::value_type>(M_2PI),
00088                                         static_cast<typename MATRIXLIKE::value_type>(cov.getColCount()))
00089                                 * ::sqrt(cov.det()));
00090                         MRPT_END
00091                 }
00092 
00093                 /** Kullback-Leibler divergence (KLD) between two independent multivariate Gaussians.
00094                   *
00095                   * \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$
00096                   */
00097                 template <typename VECTORLIKE1,typename MATRIXLIKE1,typename VECTORLIKE2,typename MATRIXLIKE2>
00098                 double KLD_Gaussians(
00099                         const VECTORLIKE1 &mu0, const MATRIXLIKE1 &cov0,
00100                         const VECTORLIKE2 &mu1, const MATRIXLIKE2 &cov1)
00101                 {
00102                         MRPT_START
00103                         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() )
00104                         const size_t N = mu0.size();
00105                         MATRIXLIKE2 cov1_inv;
00106                         cov1.inv(cov1_inv);
00107                         const VECTORLIKE1 mu_difs = mu0-mu1;
00108                         return 0.5*( log(cov1.det()/cov0.det()) + (cov1_inv*cov0).trace() + multiply_HCHt_scalar(mu_difs,cov1_inv) - N );
00109                         MRPT_END
00110                 }
00111 
00112 
00113                 /** The complementary error function of a Normal distribution
00114                   */
00115 #ifdef HAVE_ERF
00116                 inline double erfc(double x) { return ::erfc(x); }
00117 #else
00118                 double BASE_IMPEXP  erfc(double x);
00119 #endif
00120 
00121                 /** The error function of a Normal distribution
00122                   */
00123 #ifdef HAVE_ERF
00124                 inline double erf(double x) { return ::erf(x); }
00125 #else
00126                 double BASE_IMPEXP   erf(double x);
00127 #endif
00128                 /** Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
00129                   *  The employed approximation is that from Peter J. Acklam (pjacklam@online.no),
00130                   *  freely available in http://home.online.no/~pjacklam.
00131                   */
00132                 double BASE_IMPEXP  normalQuantile(double p);
00133 
00134                 /** Evaluates the Gaussian cumulative density function.
00135                   *  The employed approximation is that from W. J. Cody
00136                   *  freely available in http://www.netlib.org/specfun/erf
00137                   */
00138                 double  BASE_IMPEXP normalCDF(double p);
00139 
00140                 /** The "quantile" of the Chi-Square distribution, for dimension "dim" and probability 0<P<1 (the inverse of chi2CDF)
00141                   * An aproximation from the Wilson-Hilferty transformation is used.
00142                   */
00143                 double  BASE_IMPEXP chi2inv(double P, unsigned int dim=1);
00144 
00145                 /*! Cumulative non-central chi square distribution (approximate).
00146 
00147                         Computes approximate values of the cumulative density of a chi square distribution with \a degreesOfFreedom,
00148                         and noncentrality parameter \a noncentrality at the given argument
00149                         \a arg, i.e. the probability that a random number drawn from the distribution is below \a arg
00150                         It uses the approximate transform into a normal distribution due to Wilson and Hilferty
00151                         (see Abramovitz, Stegun: "Handbook of Mathematical Functions", formula 26.3.32).
00152                         The algorithm's running time is independent of the inputs. The accuracy is only
00153                         about 0.1 for few degrees of freedom, but reaches about 0.001 above dof = 5.
00154 
00155                         \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
00156                 */
00157                 template <class T>
00158                 double noncentralChi2CDF(unsigned int degreesOfFreedom, T noncentrality, T arg)
00159                 {
00160                         const double a = degreesOfFreedom + noncentrality;
00161                         const double b = (a + noncentrality) / square(a);
00162                         const double t = (std::pow((double)arg / a, 1.0/3.0) - (1.0 - 2.0 / 9.0 * b)) / std::sqrt(2.0 / 9.0 * b);
00163                         return 0.5*(1.0 + mrpt::math::erf(t/std::sqrt(2.0)));
00164                 }
00165 
00166                 /*! Cumulative chi square distribution.
00167 
00168                         Computes the cumulative density of a chi square distribution with \a degreesOfFreedom
00169                         and tolerance \a accuracy at the given argument \a arg, i.e. the probability that
00170                         a random number drawn from the distribution is below \a arg
00171                         by calling <tt>noncentralChi2CDF(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
00172 
00173                         \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
00174                 */
00175                 inline double chi2CDF(unsigned int degreesOfFreedom, double arg)
00176                 {
00177                         return noncentralChi2CDF(degreesOfFreedom, 0.0, arg);
00178                 }
00179 
00180                 namespace detail
00181                 {
00182                         template <class T>
00183                         void noncentralChi2OneIteration(T arg, T & lans, T & dans, T & pans, unsigned int & j)
00184                         {
00185                                 double tol = -50.0;
00186                                 if(lans < tol)
00187                                 {
00188                                         lans = lans + std::log(arg / j);
00189                                         dans = std::exp(lans);
00190                                 }
00191                                 else
00192                                 {
00193                                         dans = dans * arg / j;
00194                                 }
00195                                 pans = pans - dans;
00196                                 j += 2;
00197                         }
00198 
00199                         template <class T>
00200                         std::pair<double, double> noncentralChi2CDF_exact(unsigned int degreesOfFreedom, T noncentrality, T arg, T eps)
00201                         {
00202                                 ASSERTMSG_(noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,"noncentralChi2P(): parameters must be positive.");
00203                                 if (arg == 0.0 && degreesOfFreedom > 0)
00204                                         return std::make_pair(0.0, 0.0);
00205 
00206                                 // Determine initial values
00207                                 double b1 = 0.5 * noncentrality,
00208                                            ao = std::exp(-b1),
00209                                            eps2 = eps / ao,
00210                                            lnrtpi2 = 0.22579135264473,
00211                                            probability, density, lans, dans, pans, sum, am, hold;
00212                                 unsigned int maxit = 500,
00213                                         i, m;
00214                                 if(degreesOfFreedom % 2)
00215                                 {
00216                                         i = 1;
00217                                         lans = -0.5 * (arg + std::log(arg)) - lnrtpi2;
00218                                         dans = std::exp(lans);
00219                                         pans = erf(std::sqrt(arg/2.0));
00220                                 }
00221                                 else
00222                                 {
00223                                         i = 2;
00224                                         lans = -0.5 * arg;
00225                                         dans = std::exp(lans);
00226                                         pans = 1.0 - dans;
00227                                 }
00228 
00229                                 // Evaluate first term
00230                                 if(degreesOfFreedom == 0)
00231                                 {
00232                                         m = 1;
00233                                         degreesOfFreedom = 2;
00234                                         am = b1;
00235                                         sum = 1.0 / ao - 1.0 - am;
00236                                         density = am * dans;
00237                                         probability = 1.0 + am * pans;
00238                                 }
00239                                 else
00240                                 {
00241                                         m = 0;
00242                                         degreesOfFreedom = degreesOfFreedom - 1;
00243                                         am = 1.0;
00244                                         sum = 1.0 / ao - 1.0;
00245                                         while(i < degreesOfFreedom)
00246                                                 detail::noncentralChi2OneIteration(arg, lans, dans, pans, i);
00247                                         degreesOfFreedom = degreesOfFreedom + 1;
00248                                         density = dans;
00249                                         probability = pans;
00250                                 }
00251                                 // Evaluate successive terms of the expansion
00252                                 for(++m; m<maxit; ++m)
00253                                 {
00254                                         am = b1 * am / m;
00255                                         detail::noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
00256                                         sum = sum - am;
00257                                         density = density + am * dans;
00258                                         hold = am * pans;
00259                                         probability = probability + hold;
00260                                         if((pans * sum < eps2) && (hold < eps2))
00261                                                 break; // converged
00262                                 }
00263                                 if(m == maxit)
00264                                         THROW_EXCEPTION("noncentralChi2P(): no convergence.");
00265                                 return std::make_pair(0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
00266                         }
00267                 } // namespace detail
00268 
00269                 /*! Chi square distribution.
00270 
00271                         Computes the density of a chi square distribution with \a degreesOfFreedom
00272                         and tolerance \a accuracy at the given argument \a arg
00273                         by calling <tt>noncentralChi2(degreesOfFreedom, 0.0, arg, accuracy)</tt>.
00274 
00275                         \note Function code from the Vigra project (http://hci.iwr.uni-heidelberg.de/vigra/); code under "MIT X11 License", GNU GPL-compatible.
00276                 */
00277                 inline double chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy = 1e-7)
00278                 {
00279                         return detail::noncentralChi2CDF_exact(degreesOfFreedom, 0.0, arg, accuracy).first;
00280                 }
00281 
00282                 /** 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.
00283                   *  The container can be any MRPT container (CArray, matrices, vectors).
00284                   * \param confidenceInterval A number in the range (0,1) such as the confidence interval will be [100*confidenceInterval, 100*(1-confidenceInterval)].
00285                   */
00286                 template <typename CONTAINER>
00287                 void condidenceIntervals(
00288                         const CONTAINER &data,
00289                         typename CONTAINER::value_type &out_mean,
00290                         typename CONTAINER::value_type &out_lower_conf_interval,
00291                         typename CONTAINER::value_type &out_upper_conf_interval,
00292                         const double confidenceInterval = 0.1,
00293                         const size_t histogramNumBins = 1000 )
00294                 {
00295                         MRPT_START
00296                         ASSERT_(data.size()!=0)  // don't use .empty() here to allow using matrices
00297                         ASSERT_(confidenceInterval>0 && confidenceInterval<1)
00298 
00299                         out_mean = mean(data);
00300                         typename CONTAINER::value_type x_min,x_max;
00301                         minimum_maximum(data,x_min,x_max);
00302 
00303                         //std::vector<typename CONTAINER::value_type> xs;
00304                         //linspace(x_min,x_max,histogramNumBins, xs);
00305                         const typename CONTAINER::value_type binWidth = (x_max-x_min)/histogramNumBins;
00306 
00307                         const vector_double H = mrpt::math::histogram(data,x_min,x_max,histogramNumBins);
00308                         vector_double Hc;
00309                         cumsum(H,Hc); // CDF
00310                         Hc*=1.0/Hc.maximum();
00311 
00312                         vector_double::iterator it_low  = std::lower_bound(Hc.begin(),Hc.end(),confidenceInterval);   ASSERT_(it_low!=Hc.end())
00313                         vector_double::iterator it_high = std::upper_bound(Hc.begin(),Hc.end(),1-confidenceInterval); ASSERT_(it_high!=Hc.end())
00314                         const size_t idx_low = std::distance(Hc.begin(),it_low);
00315                         const size_t idx_high = std::distance(Hc.begin(),it_high);
00316                         out_lower_conf_interval = x_min + idx_low * binWidth;
00317                         out_upper_conf_interval = x_min + idx_high * binWidth;
00318 
00319                         MRPT_END
00320                 }
00321 
00322                 /** @} */
00323 
00324         } // End of MATH namespace
00325 
00326 } // End of namespace
00327 
00328 
00329 #endif



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN:exported at Tue Jan 25 21:56:31 UTC 2011