Bootstrapper.cc

Go to the documentation of this file.

#include "Bootstrapper.h"

#include <utl/TabulatedFunctionErrors.h>
#include <utl/MultiTabulatedFunctionErrors.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerException.h>

#undef NDEBUG
//#define NDEBUG
#include <assert.h>
#include <cmath>
#include <string>
#include <sstream>
#include <iostream>
#include <vector>
#include <algorithm>

using namespace FdLightCollectionEfficiencyKG;
using namespace utl;
using namespace std;


/*************************************************************************/
Bootstrapper::Bootstrapper(const std::list<utl::TabulatedFunctionErrors>& data, const int verbosity) :
  fData(data),
  fVerbosity(verbosity)
{
  Bootstrap();
}

/*************************************************************************/
void
Bootstrapper::Bootstrap()
{
  const unsigned int nSamples = fData.size();
  assert(nSamples);

  // for each bin i, calculate the uncertainty on the mean

  // one vector of efficiencies for each bin i
  vector<vector<double> > efficiencies;
  vector<vector<double> > efficiencyUncertainties;
  const unsigned int nBins = fData.begin()->GetNPoints();
  efficiencies.resize(nBins);
  efficiencyUncertainties.resize(nBins);

  for (unsigned int i = 0; i < nBins; ++i) {
    efficiencies[i].reserve(nSamples);
    efficiencyUncertainties[i].reserve(nSamples);
  }

  vector<double> xvalues;
  xvalues.resize(nBins, 0.);

  // iterate over the efficiency functions
  for (std::list<utl::TabulatedFunctionErrors>::const_iterator iter = fData.begin();
       iter != fData.end(); ++iter)
  {
    const utl::TabulatedFunctionErrors& lceff = *iter;
    assert(lceff.GetNPoints() > 0);

    for (unsigned int i = 0; i < lceff.GetNPoints(); ++i) {
      const double eff = lceff.GetY(i);
      const double effUnc = lceff.GetYErr(i);
      assert(isfinite(eff));
      efficiencies[i].push_back(eff);
      efficiencyUncertainties[i].push_back(effUnc);
      xvalues[i] = lceff.GetX(i); // FIXME what a waste
    }
  } // end foreach efficiency function

  fMean.Clear();
  for (unsigned int i = 0; i < efficiencies.size(); ++i) {
    if (efficiencies[i].size() != 0) {
      //pair<double, double> meanWUncert = MeanWithUncertainty(efficiencies[i], efficiencyUncertainties[i], eSampleVariance);
      //pair<double, double> meanWUncert = MeanWithUncertainty(efficiencies[i], efficiencyUncertainties[i], eMedianAbsoluteDeviation);
      pair<double, double> meanWUncert = MeanWithUncertainty(efficiencies[i], efficiencyUncertainties[i], eSampleVarianceWeighted);

      assert(isfinite(meanWUncert.first));
      assert(isfinite(meanWUncert.second));
      //cout << meanWUncert.first << " " << meanWUncert.second << endl;
      fMean.PushBack( xvalues[i], 0., meanWUncert.first, meanWUncert.second );
    }
    else
      fMean.PushBack( xvalues[i], 0., 0., 1.e3 );
  }

  ostringstream msg;
  msg << "Bootstrapped using " << nSamples << " efficiencies";
  INFO(msg);
}


/*************************************************************************/
bool
Bootstrapper::MaxRelUncertaintyBelowThreshold(const double relUncertaintyThreshold, const double absValueThreshold)
{
  const utl::TabulatedFunctionErrors& lcEff = fMean;

  for (unsigned int i = 0; i < lcEff.GetNPoints(); ++i) {
    const double val = lcEff.GetY(i);
    const double unc = lcEff.GetYErr(i);
    assert(isfinite(val));
    assert(isfinite(unc));
    //cout << val << " +/- " << unc << " " << unc/val*100. << endl;
    if (val >= absValueThreshold && unc/val > relUncertaintyThreshold) {
      if (fVerbosity >= 0) {
        ostringstream msg;
        msg << "Uncertainty still too large at"
            << " index " << i << ", value " << unc/val*100. << "%";
        INFO(msg);
      }
      return false;
    }
  }

  return true;
}


/*************************************************************************/
std::pair<double, double>
Bootstrapper::MeanWithUncertainty(const std::vector<double>& samples,
                                  const std::vector<double>& uncertainties,
                                  const VarianceEstimator est)
{
  std::pair<double, double> retval; // mean, unc. of mean
  if (est == eSampleVarianceWeighted) {
    const double mean = WeightedMean(samples, uncertainties);
    retval = std::make_pair( mean, sqrt( WeightedSampleVariance(samples, uncertainties, mean) / (double)samples.size() ) );
  }
  else {
    const double mean = Mean(samples);
    assert(isfinite(mean));
    const double kMAD = 1.4826;

    switch (est) {
    case eSampleVariance:
      retval = std::make_pair( mean, sqrt( SampleVariance(samples, mean) / (double)samples.size() ) );
      break;
    case eMedianAbsoluteDeviation:
      retval = std::make_pair( mean, kMAD*MedianAbsoluteDeviation(samples, mean) / sqrt((double)samples.size()) );
      break;
    default:
      throw utl::DoesNotComputeException("Bad variance estimator");
    };
  }

  assert(isfinite(retval.first));
  assert(isfinite(retval.second));
  return retval;
}

/*************************************************************************/
double
Bootstrapper::Mean(const std::vector<double>& samples)
{
  double sum = 0;
  const unsigned int n = samples.size();
  assert(n > 0);
  for (unsigned int i = 0; i < n; ++i) {
    assert(isfinite(samples[i]));
    sum += samples[i];
  }
  return sum / (double)n;
}

/*************************************************************************/
double
Bootstrapper::WeightedMean(const std::vector<double>& samples, const std::vector<double>& uncertainties)
{
  const unsigned int n = samples.size();
  assert(n > 0);
  assert(uncertainties.size() == n);

  double sum = 0, sumW = 0;
  for (unsigned int i = 0; i < n; ++i) {
    assert(isfinite(samples[i]));
    assert(isfinite(uncertainties[i]));
    const double unc = uncertainties[i];
    const double w = 1./(unc*unc);
    sum  += samples[i]*w;
    sumW += w;
  }
  return sum / sumW;
}

/*************************************************************************/
double
Bootstrapper::SampleVariance(const std::vector<double>& samples, const double mean)
{
  double variance = 0.;
  const unsigned int n = samples.size();
  assert(n > 1);
  for (unsigned int i = 0; i < n; i++) {
    assert(isfinite(samples[i]));
    variance += pow(samples[i]-mean, 2.);
  }
  return variance / ((double)n - 1.);
}

/*************************************************************************/
double
Bootstrapper::WeightedSampleVariance(const std::vector<double>& samples,
                                     const std::vector<double>& uncertainties,
                                     const double mean)
{
  const unsigned int n = samples.size();
  assert(n > 1);
  assert(uncertainties.size() == n);

  double sumW = 0., sumW2 = 0., wvariance = 0.;
  for (unsigned int i = 0; i < n; ++i) {
    const double unc = uncertainties[i];
    assert(isfinite(samples[i]));
    assert(isfinite(unc));
    assert(unc > 0.);
    const double w = 1./(unc*unc);
    sumW      += w;
    sumW2     += w*w;
    wvariance += w*pow(samples[i]-mean, 2.);
  }
  return (sumW / (sumW*sumW - sumW2)) * wvariance;
}

/*************************************************************************/
double
Bootstrapper::Median(const std::vector<double>& sortedSamples)
{
  const double medianElem = ((double)sortedSamples.size()+1.)/2.;
  const unsigned int idx  = (unsigned int)medianElem;
  double median;
  if ( (medianElem-(double)idx > 0.1) )
    median = 0.5 * (sortedSamples[idx] + sortedSamples[idx+1]);
  else
    median = sortedSamples[idx];

  assert(isfinite(median));
  return median;
}

/*************************************************************************/
double
Bootstrapper::MedianAbsoluteDeviation(const std::vector<double>& samples, const double /*mean*/)
{
  // MAD=median_i(|X_i-median_j(X_j)|)
  vector<double> sortedSamples = samples; // Note: Consider the non-const case as an optimization?
  const unsigned int nSamples = sortedSamples.size();

  std::sort(sortedSamples.begin(), sortedSamples.end());
  const double sampleMedian = Median(sortedSamples);
  assert(isfinite(sampleMedian));

  for (unsigned int i = 0; i < nSamples; ++i) {
    assert(isfinite(sortedSamples[i]));
    sortedSamples[i] = fabs(sortedSamples[i] - sampleMedian);
  }

  std::sort(sortedSamples.begin(), sortedSamples.end());

  const double mad = Median(sortedSamples);
  assert(isfinite(mad));
  return mad;
}

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: