MRPT  1.9.9
math.cpp
Go to the documentation of this file.
1 /* +------------------------------------------------------------------------+
2  | Mobile Robot Programming Toolkit (MRPT) |
3  | https://www.mrpt.org/ |
4  | |
5  | Copyright (c) 2005-2019, Individual contributors, see AUTHORS file |
6  | See: https://www.mrpt.org/Authors - All rights reserved. |
7  | Released under BSD License. See: https://www.mrpt.org/License |
8  +------------------------------------------------------------------------+ */
9 
10 #include "math-precomp.h" // Precompiled headers
11 
12 #include <mrpt/math/CMatrixD.h>
13 #include <mrpt/math/data_utils.h>
15 //#include <mrpt/math/eigen_extensions.h>
16 #include <mrpt/math/ops_matrices.h>
17 #include <mrpt/math/utils.h>
18 #include <mrpt/math/wrap2pi.h>
20 #include <Eigen/Dense>
21 #include <algorithm>
22 #include <cmath> // erf(), ...
23 #include <iostream>
24 #include <mrpt/math/interp_fit.hpp>
25 
26 using namespace mrpt;
27 using namespace mrpt::math;
28 using namespace std;
30 
31 /*---------------------------------------------------------------
32  normalPDF
33  ---------------------------------------------------------------*/
34 double math::normalPDF(double x, double mu, double std)
35 {
36  return ::exp(-0.5 * square((x - mu) / std)) /
37  (std * 2.506628274631000502415765284811);
38 }
39 
40 /*---------------------------------------------------------------
41  chi2inv
42  ---------------------------------------------------------------*/
43 double math::chi2inv(double P, unsigned int dim)
44 {
45  ASSERT_(P >= 0 && P < 1);
46  if (P == 0)
47  return 0;
48  else
49  return dim * pow(1.0 - 2.0 / (9 * dim) +
50  sqrt(2.0 / (9 * dim)) * normalQuantile(P),
51  3);
52 }
53 
54 /*---------------------------------------------------------------
55  factorial64
56  ---------------------------------------------------------------*/
58 {
59  uint64_t ret = 1;
60 
61  for (unsigned int i = 2; i <= n; i++) ret *= i;
62 
63  return ret;
64 }
65 
66 /*---------------------------------------------------------------
67  factorial
68  ---------------------------------------------------------------*/
69 double math::factorial(unsigned int n)
70 {
71  double retLog = 0;
72 
73  for (unsigned int i = 2; i <= n; i++) retLog += ::log((double)n);
74 
75  return ::exp(retLog);
76 }
77 
78 /*---------------------------------------------------------------
79  normalQuantile
80  ---------------------------------------------------------------*/
81 double math::normalQuantile(double p)
82 {
83  double q, t, u;
84 
85  static const double a[6] = {-3.969683028665376e+01, 2.209460984245205e+02,
86  -2.759285104469687e+02, 1.383577518672690e+02,
87  -3.066479806614716e+01, 2.506628277459239e+00};
88  static const double b[5] = {-5.447609879822406e+01, 1.615858368580409e+02,
89  -1.556989798598866e+02, 6.680131188771972e+01,
90  -1.328068155288572e+01};
91  static const double c[6] = {-7.784894002430293e-03, -3.223964580411365e-01,
92  -2.400758277161838e+00, -2.549732539343734e+00,
93  4.374664141464968e+00, 2.938163982698783e+00};
94  static const double d[4] = {7.784695709041462e-03, 3.224671290700398e-01,
95  2.445134137142996e+00, 3.754408661907416e+00};
96 
97  ASSERT_(!std::isnan(p));
99 
100  q = min(p, 1 - p);
101 
102  if (q > 0.02425)
103  {
104  /* Rational approximation for central region. */
105  u = q - 0.5;
106  t = u * u;
107  u = u *
108  (((((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4]) * t +
109  a[5]) /
110  (((((b[0] * t + b[1]) * t + b[2]) * t + b[3]) * t + b[4]) * t + 1);
111  }
112  else
113  {
114  /* Rational approximation for tail region. */
115  t = sqrt(-2 * ::log(q));
116  u = (((((c[0] * t + c[1]) * t + c[2]) * t + c[3]) * t + c[4]) * t +
117  c[5]) /
118  ((((d[0] * t + d[1]) * t + d[2]) * t + d[3]) * t + 1);
119  }
120 
121  /* The relative error of the approximation has absolute value less
122  than 1.15e-9. One iteration of Halley's rational method (third
123  order) gives full machine precision... */
124  t = normalCDF(u) - q; /* error */
125  t = t * 2.506628274631000502415765284811 *
126  ::exp(u * u / 2); /* f(u)/df(u) */
127  u = u - t / (1 + u * t / 2); /* Halley's method */
128 
129  return (p > 0.5 ? -u : u);
130 }
131 
132 /*---------------------------------------------------------------
133  normalCDF
134  ---------------------------------------------------------------*/
135 double math::normalCDF(double u)
136 {
137  static const double a[5] = {1.161110663653770e-002, 3.951404679838207e-001,
138  2.846603853776254e+001, 1.887426188426510e+002,
139  3.209377589138469e+003};
140  static const double b[5] = {1.767766952966369e-001, 8.344316438579620e+000,
141  1.725514762600375e+002, 1.813893686502485e+003,
142  8.044716608901563e+003};
143  static const double c[9] = {
144  2.15311535474403846e-8, 5.64188496988670089e-1, 8.88314979438837594e00,
145  6.61191906371416295e01, 2.98635138197400131e02, 8.81952221241769090e02,
146  1.71204761263407058e03, 2.05107837782607147e03, 1.23033935479799725E03};
147  static const double d[9] = {
148  1.00000000000000000e00, 1.57449261107098347e01, 1.17693950891312499e02,
149  5.37181101862009858e02, 1.62138957456669019e03, 3.29079923573345963e03,
150  4.36261909014324716e03, 3.43936767414372164e03, 1.23033935480374942e03};
151  static const double p[6] = {1.63153871373020978e-2, 3.05326634961232344e-1,
152  3.60344899949804439e-1, 1.25781726111229246e-1,
153  1.60837851487422766e-2, 6.58749161529837803e-4};
154  static const double q[6] = {1.00000000000000000e00, 2.56852019228982242e00,
155  1.87295284992346047e00, 5.27905102951428412e-1,
156  6.05183413124413191e-2, 2.33520497626869185e-3};
157  double y, z;
158 
159  ASSERT_(!std::isnan(u));
160  ASSERT_(std::isfinite(u));
161 
162  y = fabs(u);
163  if (y <= 0.46875 * 1.4142135623730950488016887242097)
164  {
165  /* evaluate erf() for |u| <= sqrt(2)*0.46875 */
166  z = y * y;
167  y = u * ((((a[0] * z + a[1]) * z + a[2]) * z + a[3]) * z + a[4]) /
168  ((((b[0] * z + b[1]) * z + b[2]) * z + b[3]) * z + b[4]);
169  return 0.5 + y;
170  }
171 
172  z = ::exp(-y * y / 2) / 2;
173  if (y <= 4.0)
174  {
175  /* evaluate erfc() for sqrt(2)*0.46875 <= |u| <= sqrt(2)*4.0 */
176  y = y / 1.4142135623730950488016887242097;
177  y = ((((((((c[0] * y + c[1]) * y + c[2]) * y + c[3]) * y + c[4]) * y +
178  c[5]) *
179  y +
180  c[6]) *
181  y +
182  c[7]) *
183  y +
184  c[8])
185 
186  / ((((((((d[0] * y + d[1]) * y + d[2]) * y + d[3]) * y + d[4]) * y +
187  d[5]) *
188  y +
189  d[6]) *
190  y +
191  d[7]) *
192  y +
193  d[8]);
194 
195  y = z * y;
196  }
197  else
198  {
199  /* evaluate erfc() for |u| > sqrt(2)*4.0 */
200  z = z * 1.4142135623730950488016887242097 / y;
201  y = 2 / (y * y);
202  y = y *
203  (((((p[0] * y + p[1]) * y + p[2]) * y + p[3]) * y + p[4]) * y +
204  p[5]) /
205  (((((q[0] * y + q[1]) * y + q[2]) * y + q[3]) * y + q[4]) * y +
206  q[5]);
207  y = z * (0.564189583547756286948 - y);
208  }
209  return (u < 0.0 ? y : 1 - y);
210 }
211 
212 // Loads a vector from a text file:
213 bool math::loadVector(std::istream& f, ::std::vector<int>& d)
214 {
215  MRPT_START
216 
217  std::string str;
218  if (!std::getline(f, str)) return false;
219 
220  const char* s = str.c_str();
221 
222  char *nextTok, *context;
223  const char* delim = " \t";
224 
225  d.clear();
226  nextTok = mrpt::system::strtok((char*)s, delim, &context);
227  while (nextTok != nullptr)
228  {
229  d.push_back(atoi(nextTok));
230  nextTok = mrpt::system::strtok(nullptr, delim, &context);
231  };
232 
233  return true;
234  MRPT_END
235 }
236 
237 bool math::loadVector(std::istream& f, ::std::vector<double>& d)
238 {
239  MRPT_START
240 
241  std::string str;
242  if (!std::getline(f, str)) return false;
243 
244  const char* s = str.c_str();
245 
246  char *nextTok, *context;
247  const char* delim = " \t";
248 
249  d.clear();
250  nextTok = mrpt::system::strtok((char*)s, delim, &context);
251  while (nextTok != nullptr)
252  {
253  d.push_back(atof(nextTok));
254  nextTok = mrpt::system::strtok(nullptr, delim, &context);
255  };
256 
257  return true;
258  MRPT_END
259 }
260 
261 // See declaration for the documentation
263  const CVectorDouble& logWeights, const CVectorDouble& logLikelihoods)
264 {
265  MRPT_START
266 
267  // Explained in:
268  // https://www.mrpt.org/Averaging_Log-Likelihood_Values:Numerical_Stability
269  ASSERT_(logWeights.size() == logLikelihoods.size());
270 
271  if (!logWeights.size())
272  THROW_EXCEPTION("ERROR: logWeights vector is empty!");
273 
275  double lw_max = math::maximum(logWeights);
276  double ll_max = math::maximum(logLikelihoods);
277  double SUM1 = 0, SUM2 = 0, tmpVal;
278 
279  for (itLW = logWeights.begin(), itLL = logLikelihoods.begin();
280  itLW != logWeights.end(); itLW++, itLL++)
281  {
282  tmpVal = *itLW - lw_max;
283  SUM1 += std::exp(tmpVal);
284  SUM2 += std::exp(tmpVal + *itLL - ll_max);
285  }
286 
287  double res = -std::log(SUM1) + std::log(SUM2) + ll_max;
289  return res;
290  MRPT_END
291 }
292 
293 // Unweighted version:
294 double math::averageLogLikelihood(const CVectorDouble& logLikelihoods)
295 {
296  MRPT_START
297 
298  // Explained in:
299  // https://www.mrpt.org/Averaging_Log-Likelihood_Values:Numerical_Stability
300  size_t N = logLikelihoods.size();
301  if (!N) THROW_EXCEPTION("ERROR: logLikelihoods vector is empty!");
302 
303  double ll_max = math::maximum(logLikelihoods);
304  double SUM1 = 0;
305 
306  for (size_t i = 0; i < N; i++) SUM1 += exp(logLikelihoods[i] - ll_max);
307 
308  double res = log(SUM1) - log(static_cast<double>(N)) + ll_max;
309 
311  return res;
312  MRPT_END
313 }
314 
315 // Wrapped angles average:
316 double math::averageWrap2Pi(const CVectorDouble& angles)
317 {
318  if (angles.empty()) return 0;
319 
320  int W_phi_R = 0, W_phi_L = 0;
321  double phi_R = 0, phi_L = 0;
322 
323  // First: XY
324  // -----------------------------------
325  for (CVectorDouble::Index i = 0; i < angles.size(); i++)
326  {
327  double phi = angles[i];
328  if (abs(phi) > 1.5707963267948966192313216916398)
329  {
330  // LEFT HALF: 0,2pi
331  if (phi < 0) phi = (M_2PI + phi);
332 
333  phi_L += phi;
334  W_phi_L++;
335  }
336  else
337  {
338  // RIGHT HALF: -pi,pi
339  phi_R += phi;
340  W_phi_R++;
341  }
342  }
343 
344  // Next: PHI
345  // -----------------------------------
346  // The mean value from each side:
347  if (W_phi_L) phi_L /= static_cast<double>(W_phi_L); // [0,2pi]
348  if (W_phi_R) phi_R /= static_cast<double>(W_phi_R); // [-pi,pi]
349 
350  // Left side to [-pi,pi] again:
351  if (phi_L > M_PI) phi_L -= M_2PI;
352 
353  // The total mean:
354  return ((phi_L * W_phi_L + phi_R * W_phi_R) / (W_phi_L + W_phi_R));
355 }
356 
358  const CMatrixFloat& cov, const CVectorFloat& mean, const float& stdCount,
359  const string& style, const size_t& nEllipsePoints)
360 {
361  MRPT_START
362  CMatrixD cov2(cov);
363  CVectorDouble mean2(2);
364  mean2[0] = mean[0];
365  mean2[1] = mean[1];
366 
368  cov2, mean2, stdCount, style, nEllipsePoints);
369 
370  MRPT_END
371 }
372 
374  const CMatrixDouble& cov, const CVectorDouble& mean, const float& stdCount,
375  const string& style, const size_t& nEllipsePoints)
376 {
377  MRPT_START
378 
379  ASSERT_(cov.cols() == cov.rows() && cov.cols() == 2);
380  ASSERT_(cov(0, 1) == cov(1, 0));
381  ASSERT_(!((cov(0, 0) == 0) ^ (cov(1, 1) == 0))); // Both or none 0
382  ASSERT_(mean.size() == 2);
383 
384  std::vector<float> X, Y, COS, SIN;
385  std::vector<float>::iterator x, y, Cos, Sin;
386  double ang;
387  string str;
388 
389  X.resize(nEllipsePoints);
390  Y.resize(nEllipsePoints);
391  COS.resize(nEllipsePoints);
392  SIN.resize(nEllipsePoints);
393 
394  // Fill the angles:
395  for (Cos = COS.begin(), Sin = SIN.begin(), ang = 0; Cos != COS.end();
396  ++Cos, ++Sin, ang += (M_2PI / (nEllipsePoints - 1)))
397  {
398  *Cos = (float)cos(ang);
399  *Sin = (float)sin(ang);
400  }
401 
403  const auto eigVec = eigensolver.eigenvectors();
404  auto eigVal = eigensolver.eigenvalues();
405 
406  eigVal = eigVal.array().sqrt().matrix();
407  Eigen::MatrixXd M = eigVal * eigVec.transpose();
408 
409  // Compute the points of the ellipsoid:
410  // ----------------------------------------------
411  for (x = X.begin(), y = Y.begin(), Cos = COS.begin(), Sin = SIN.begin();
412  x != X.end(); ++x, ++y, ++Cos, ++Sin)
413  {
414  *x = (float)(mean[0] + stdCount * (*Cos * M(0, 0) + *Sin * M(1, 0)));
415  *y = (float)(mean[1] + stdCount * (*Cos * M(0, 1) + *Sin * M(1, 1)));
416  }
417 
418  // Save the code to plot the ellipsoid:
419  // ----------------------------------------------
420  str += string("plot([ ");
421  for (x = X.begin(); x != X.end(); ++x)
422  {
423  str += format("%.4f", *x);
424  if (x != (X.end() - 1)) str += format(",");
425  }
426  str += string("],[ ");
427  for (y = Y.begin(); y != Y.end(); ++y)
428  {
429  str += format("%.4f", *y);
430  if (y != (Y.end() - 1)) str += format(",");
431  }
432 
433  str += format("],'%s');\n", style.c_str());
434 
435  return str;
436  MRPT_END_WITH_CLEAN_UP(std::cerr << "The matrix that led to error was: "
437  << std::endl
438  << cov.asEigen() << std::endl;)
439 }
440 
442  const double x, const double x0, const double y0, const double x1,
443  const double y1, bool wrap2pi)
444 {
445  MRPT_START
446  if (x0 == x1)
447  THROW_EXCEPTION_FMT("ERROR: Both x0 and x1 are equal (=%f)", x0);
448 
449  const double Ax = x1 - x0;
450  const double Ay = y1 - y0;
451 
452  double r = y0 + Ay * (x - x0) / Ax;
453  if (!wrap2pi)
454  return r;
455  else
456  return mrpt::math::wrapToPi(r);
457 
458  MRPT_END
459 }
460 
461 /*---------------------------------------------------------------
462  median filter of a vector
463  ---------------------------------------------------------------*/
464 // template<typename VECTOR>
466  const std::vector<double>& inV, std::vector<double>& outV,
467  const int& _winSize, const int& numberOfSigmas)
468 {
469  MRPT_UNUSED_PARAM(numberOfSigmas);
470  ASSERT_((int)inV.size() >= _winSize);
471  ASSERT_(_winSize >= 2); // The minimum window size is 3 elements
472  size_t winSize = _winSize;
473 
474  if (!(winSize % 2)) // We use an odd number of elements for the window size
475  winSize++;
476 
477  size_t sz = inV.size();
478  outV.resize(sz);
479 
480  std::vector<double> aux(winSize);
481  size_t mpoint = winSize / 2;
482  for (size_t k = 0; k < sz; ++k)
483  {
484  aux.clear();
485 
486  size_t idx_to_start =
487  std::max(size_t(0), k - mpoint); // Dealing with the boundaries
488  size_t n_elements =
489  std::min(std::min(winSize, sz + mpoint - k), k + mpoint + 1);
490 
491  aux.resize(n_elements);
492  for (size_t m = idx_to_start, n = 0; m < idx_to_start + n_elements;
493  ++m, ++n)
494  aux[n] = inV[m];
495 
496  std::sort(aux.begin(), aux.end());
497 
498  size_t auxSz = aux.size();
499  size_t auxMPoint = auxSz / 2;
500  outV[k] = (auxSz % 2) ? (aux[auxMPoint])
501  : (0.5 * (aux[auxMPoint - 1] +
502  aux[auxMPoint])); // If the window is
503  // even, take the
504  // mean value of the
505  // middle points
506  } // end-for
507 } // end medianFilter
508 
509 double mrpt::math::chi2CDF(unsigned int degreesOfFreedom, double arg)
510 {
511  return noncentralChi2CDF(degreesOfFreedom, 0.0, arg);
512 }
513 
514 template <class T>
516  T arg, T& lans, T& dans, T& pans, unsigned int& j)
517 {
518  double tol = -50.0;
519  if (lans < tol)
520  {
521  lans = lans + std::log(arg / j);
522  dans = std::exp(lans);
523  }
524  else
525  {
526  dans = dans * arg / j;
527  }
528  pans = pans - dans;
529  j += 2;
530 }
531 
532 std::pair<double, double> mrpt::math::noncentralChi2PDF_CDF(
533  unsigned int degreesOfFreedom, double noncentrality, double arg, double eps)
534 {
535  ASSERTMSG_(
536  noncentrality >= 0.0 && arg >= 0.0 && eps > 0.0,
537  "noncentralChi2PDF_CDF(): parameters must be positive.");
538 
539  if (arg == 0.0 && degreesOfFreedom > 0) return std::make_pair(0.0, 0.0);
540 
541  // Determine initial values
542  double b1 = 0.5 * noncentrality, ao = std::exp(-b1), eps2 = eps / ao,
543  lnrtpi2 = 0.22579135264473, probability, density, lans, dans, pans,
544  sum, am, hold;
545  unsigned int maxit = 500, i, m;
546  if (degreesOfFreedom % 2)
547  {
548  i = 1;
549  lans = -0.5 * (arg + std::log(arg)) - lnrtpi2;
550  dans = std::exp(lans);
551  pans = std::erf(std::sqrt(arg / 2.0));
552  }
553  else
554  {
555  i = 2;
556  lans = -0.5 * arg;
557  dans = std::exp(lans);
558  pans = 1.0 - dans;
559  }
560 
561  // Evaluate first term
562  if (degreesOfFreedom == 0)
563  {
564  m = 1;
565  degreesOfFreedom = 2;
566  am = b1;
567  sum = 1.0 / ao - 1.0 - am;
568  density = am * dans;
569  probability = 1.0 + am * pans;
570  }
571  else
572  {
573  m = 0;
574  degreesOfFreedom = degreesOfFreedom - 1;
575  am = 1.0;
576  sum = 1.0 / ao - 1.0;
577  while (i < degreesOfFreedom)
578  noncentralChi2OneIteration(arg, lans, dans, pans, i);
579  degreesOfFreedom = degreesOfFreedom + 1;
580  density = dans;
581  probability = pans;
582  }
583  // Evaluate successive terms of the expansion
584  for (++m; m < maxit; ++m)
585  {
586  am = b1 * am / m;
587  noncentralChi2OneIteration(arg, lans, dans, pans, degreesOfFreedom);
588  sum = sum - am;
589  density = density + am * dans;
590  hold = am * pans;
591  probability = probability + hold;
592  if ((pans * sum < eps2) && (hold < eps2)) break; // converged
593  }
594  if (m == maxit) THROW_EXCEPTION("noncentralChi2PDF_CDF(): no convergence.");
595  return std::make_pair(
596  0.5 * ao * density, std::min(1.0, std::max(0.0, ao * probability)));
597 }
598 
600  unsigned int degreesOfFreedom, double arg, double accuracy)
601 {
603  degreesOfFreedom, 0.0, arg, accuracy)
604  .first;
605 }
606 
608  unsigned int degreesOfFreedom, double noncentrality, double arg)
609 {
610  const double a = degreesOfFreedom + noncentrality;
611  const double b = (a + noncentrality) / square(a);
612  const double t =
613  (std::pow((double)arg / a, 1.0 / 3.0) - (1.0 - 2.0 / 9.0 * b)) /
614  std::sqrt(2.0 / 9.0 * b);
615  return 0.5 * (1.0 + std::erf(t / std::sqrt(2.0)));
616 }
617 
619 {
620  v.resize(s.ReadAs<uint32_t>());
621  if (v.size() > 0) s.ReadBufferFixEndianness(&v[0], v.size());
622  return s;
623 }
625 {
626  v.resize(s.ReadAs<uint32_t>());
627  if (v.size() > 0) s.ReadBufferFixEndianness(&v[0], v.size());
628  return s;
629 }
631 {
632  s.WriteAs<uint32_t>(v.size());
633  if (v.size() > 0) s.WriteBufferFixEndianness(&v[0], v.size());
634  return s;
635 }
637 {
638  s.WriteAs<uint32_t>(v.size());
639  if (v.size() > 0) s.WriteBufferFixEndianness(&v[0], v.size());
640  return s;
641 }
double averageWrap2Pi(const CVectorDouble &angles)
Computes the average of a sequence of angles in radians taking into account the correct wrapping in t...
Definition: math.cpp:316
std::string MATLAB_plotCovariance2D(const CMatrixFloat &cov22, const CVectorFloat &mean, const float &stdCount, const std::string &style=std::string("b"), const size_t &nEllipsePoints=30)
Generates a string with the MATLAB commands required to plot an confidence interval (ellipse) for a 2...
Definition: math.cpp:357
#define MRPT_START
Definition: exceptions.h:241
GLdouble GLdouble t
Definition: glext.h:3695
GLdouble GLdouble z
Definition: glext.h:3879
This class is a "CSerializable" wrapper for "CMatrixDynamic<double>".
Definition: CMatrixD.h:23
#define min(a, b)
#define M_2PI
Definition: common.h:58
#define THROW_EXCEPTION(msg)
Definition: exceptions.h:67
double normalCDF(double p)
Evaluates the Gaussian cumulative density function.
Definition: math.cpp:135
GLdouble GLdouble GLdouble GLdouble q
Definition: glext.h:3727
GLenum GLsizei n
Definition: glext.h:5136
size_type size() const
Get a 2-vector with [NROWS NCOLS] (as in MATLAB command size(x))
This file implements miscelaneous matrix and matrix/vector operations, and internal functions in mrpt...
STL namespace.
void WriteAs(const TYPE_FROM_ACTUAL &value)
Definition: CArchive.h:157
#define MRPT_END_WITH_CLEAN_UP(stuff)
Definition: exceptions.h:247
GLdouble s
Definition: glext.h:3682
T square(const T x)
Inline function for the square of a number.
double factorial(unsigned int n)
Computes the factorial of an integer number and returns it as a double value (internally it uses loga...
Definition: math.cpp:69
double normalQuantile(double p)
Evaluates the Gaussian distribution quantile for the probability value p=[0,1].
Definition: math.cpp:81
#define ASSERT_(f)
Defines an assertion mechanism.
Definition: exceptions.h:120
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:532
#define MRPT_CHECK_NORMAL_NUMBER(v)
Throws an exception if the number is NaN, IND, or +/-INF, or return the same number otherwise...
Definition: exceptions.h:125
double averageLogLikelihood(const CVectorDouble &logLikelihoods)
A numerically-stable method to compute average likelihood values with strongly different ranges (unwe...
Definition: math.cpp:294
const GLubyte * c
Definition: glext.h:6406
bool loadVector(std::istream &f, std::vector< int > &d)
Loads one row of a text file as a numerical std::vector.
double chi2PDF(unsigned int degreesOfFreedom, double arg, double accuracy=1e-7)
Definition: math.cpp:599
char * strtok(char *str, const char *strDelimit, char **context) noexcept
An OS-independent method for tokenizing a string.
CONTAINER::Scalar sum(const CONTAINER &v)
Computes the sum of all the elements.
uint64_t factorial64(unsigned int n)
Computes the factorial of an integer number and returns it as a 64-bit integer number.
Definition: math.cpp:57
CONTAINER::Scalar maximum(const CONTAINER &v)
GLubyte GLubyte b
Definition: glext.h:6372
const double eps
#define ASSERTMSG_(f, __ERROR_MSG)
Defines an assertion mechanism.
Definition: exceptions.h:108
typename vec_t::const_iterator const_iterator
CMatrixDouble 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
GLsizei const GLchar ** string
Definition: glext.h:4116
T wrapToPi(T a)
Modifies the given angle to translate it into the ]-pi,pi] range.
Definition: wrap2pi.h:50
size_type rows() const
Number of rows in the matrix.
size_type cols() const
Number of columns in the matrix.
mrpt::serialization::CArchive & operator>>(mrpt::serialization::CArchive &in, CMatrixD::Ptr &pObj)
double chi2CDF(unsigned int degreesOfFreedom, double arg)
Definition: math.cpp:509
void medianFilter(const std::vector< double > &inV, std::vector< double > &outV, const int &winSize, const int &numberOfSigmas=2)
Definition: math.cpp:465
void noncentralChi2OneIteration(T arg, T &lans, T &dans, T &pans, unsigned int &j)
Definition: math.cpp:515
unsigned __int64 uint64_t
Definition: rptypes.h:53
const GLdouble * v
Definition: glext.h:3684
This is the global namespace for all Mobile Robot Programming Toolkit (MRPT) libraries.
Virtual base class for "archives": classes abstracting I/O streams.
Definition: CArchive.h:53
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:43
GLdouble GLdouble GLdouble r
Definition: glext.h:3711
mrpt::serialization::CArchive & operator<<(mrpt::serialization::CArchive &s, const CVectorFloat &a)
Definition: math.cpp:630
std::string format(const char *fmt,...) MRPT_printf_format_check(1
A std::string version of C sprintf.
Definition: format.cpp:16
#define MRPT_END
Definition: exceptions.h:245
double mean(const CONTAINER &v)
Computes the mean value of a vector.
GLenum GLint GLint y
Definition: glext.h:3542
double noncentralChi2CDF(unsigned int degreesOfFreedom, double noncentrality, double arg)
Definition: math.cpp:607
double interpolate2points(const double x, const double x0, const double y0, const double x1, const double y1, bool wrap2pi=false)
Linear interpolation/extrapolation: evaluates at "x" the line (x0,y0)-(x1,y1).
Definition: math.cpp:441
GLuint res
Definition: glext.h:7385
GLenum GLint x
Definition: glext.h:3542
#define THROW_EXCEPTION_FMT(_FORMAT_STRING,...)
Definition: exceptions.h:69
unsigned __int32 uint32_t
Definition: rptypes.h:50
size_t ReadBufferFixEndianness(T *ptr, size_t ElementCount)
Reads a sequence of elemental datatypes, taking care of reordering their bytes from the MRPT stream s...
Definition: CArchive.h:89
double normalPDF(double x, double mu, double std)
Evaluates the univariate normal (Gaussian) distribution at a given point "x".
Definition: math.cpp:34
EIGEN_MAP asEigen()
Get as an Eigen-compatible Eigen::Map object.
This template class provides the basic functionality for a general 2D any-size, resizable container o...
GLubyte GLubyte GLubyte a
Definition: glext.h:6372
GLfloat GLfloat p
Definition: glext.h:6398
#define MRPT_UNUSED_PARAM(a)
Determines whether this is an X86 or AMD64 platform.
Definition: common.h:186



Page generated by Doxygen 1.8.14 for MRPT 1.9.9 Git: ce1a28c9f Fri Aug 23 08:02:09 2019 +0200 at vie ago 23 08:10:11 CEST 2019