Tools/CIC/utl/Math.cc

Go to the documentation of this file.

#include <utl/Math.h>
#include <utl/Accumulator.h>
#include <numeric>

using namespace std;
using namespace utl::Accumulator;


namespace utl {

inline
TRandom3*
GetRandom3()
{
  gRandom->SetSeed(0);
  return (TRandom3*)(gRandom);
}


TRandom3* const gRandom3 = GetRandom3();


// Kolmogorov-Smirnov measure of uniformity
/*double
UniformKolmogorovSmirnov(const vector<double>& cdf)
{
  // cdf is sorted list of discrete points where jump occurs
  // the values are assumed to be normalized into [0, 1]
  const double n = cdf.size();
  Max<double> maxm(0);
  Max<double> maxp(0);
  unsigned int i = 0;
  double c = 0;
  for (const auto x : cdf) {
    maxm(x - c);
    c = (++i) / n;
    maxp(c - x);
  }
  return sqrt(n) * (maxp.GetMax() + maxm.GetMax());
}*/


double
UniformKolmogorovSmirnov(const vector<pair<double, double>>& wcdf,
                         double wtot)
{
  // wcdf a list of (x,w) points, sorted in x, assuming x in [0,1].
  // weights w can be arbitrary but positive
  if (wcdf.empty())
    return 0;
  if (!wtot) {
    for (const auto& xw : wcdf)
      wtot += xw.second;
    if (!wtot)
      return 0;
  }
  Max<double> maxm(0);
  Max<double> maxp(0);
  double wsum = 0;
  for (const auto& xw : wcdf) {
    maxm(xw.first - wsum/wtot);
    wsum += xw.second;
    maxp(wsum/wtot - xw.first);
  }
  return sqrt(wtot) * (maxp.GetMax() + maxm.GetMax());
}


double
UniformAndersonDarling(const vector<pair<double, double>>& wcdf,
                       double wtot)
{
  // wcdf a list of (x,w) points, sorted in x, assuming x in [0,1].
  // weights w can be arbitrary but positive
  if (wcdf.empty())
    return 0;
  if (!wtot) {
    for (const auto& xw : wcdf)
      wtot += xw.second;
    if (!wtot)
      return 0;
  }
  double wsum = 0;
  double a = wcdf.front().first;
  double a1 = 1 - a;
  double anderson = -log(a1);
  for (unsigned int i = 0, n = wcdf.size() - 1; i < n; ++i) {
    const double b = wcdf[i+1].first;
    const double b1 = 1 - b;
    wsum += wcdf[i].second;
    const double d = wsum/wtot;
    anderson += Sqr(1 - d)*log(a1/b1) + Sqr(d)*log(b/a);
    a = b;
    a1 = b1;
  }
  anderson -= log(a) + 1;
  return wtot * anderson;
}

}