00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00039
00040 namespace mrpt
00041 {
00042 namespace math
00043 {
00044 using namespace mrpt::utils;
00045
00046
00047
00048
00049
00050
00051 double BASE_IMPEXP normalPDF(double x, double mu, double std);
00052
00053
00054
00055
00056
00057
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
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
00094
00095
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
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
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
00129
00130
00131
00132 double BASE_IMPEXP normalQuantile(double p);
00133
00134
00135
00136
00137
00138 double BASE_IMPEXP normalCDF(double p);
00139
00140
00141
00142
00143 double BASE_IMPEXP chi2inv(double P, unsigned int dim=1);
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
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
00167
00168
00169
00170
00171
00172
00173
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
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
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
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;
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 }
00268
00269
00270
00271
00272
00273
00274
00275
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
00283
00284
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)
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
00304
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);
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 }
00325
00326 }
00327
00328
00329 #endif