Simultaneous estimation of Several Percentiles

percentile.cc

#include <algorithm>
#include <functional>
#include <iomanip>
#include <iostream>

#include "percentile.h"

using namespace std;
using namespace stat;

const float_t Percentile::DEFAULT_PERCENTILE[] = {0.05F, 0.5F, 0.683F, 0.75F, 0.85F, 0.954F, 0.99F};

Percentile::Percentile() : name("Percentile"), quantiles(DEFAULT_PERCENTILE), m(7) {

  const int N = 2 * m + 3;

  q = new float_t[N + 1];
  n = new int[N + 1];
  f = new float_t[N + 1];
  d = new float_t[N + 1];

  e = new float_t[m];

  clear();
}

Percentile::Percentile(const char name[]) : name(name), quantiles(DEFAULT_PERCENTILE), m(7) {

  const int N = 2 * m + 3;

  q = new float_t[N + 1];
  n = new int[N + 1];
  f = new float_t[N + 1];
  d = new float_t[N + 1];

  e = new float_t[m];

  clear();
}

Percentile::Percentile(const char name[], const float_t quantiles[], const int nQuantiles) : name(name), m(nQuantiles) {

  this->quantiles = quantiles;

  const int N = 2 * m + 3;

  q = new float_t[N + 1];
  n = new int[N + 1];
  f = new float_t[N + 1];
  d = new float_t[N + 1];

  e = new float_t[m];

  clear();
}

Percentile::~Percentile() {
  delete[] q;
  delete[] n;
  delete[] f;
  delete[] d;
  delete[] e;
}

void Percentile::clear() {
  for (int i = 1; i <= 2 * m + 3; i++) {
    n[i] = i + 1;
  }

  f[1] = 0.0;
  f[2 * m + 3] = 1.0;

  for (int i = 1; i <= m; i++) {
    f[2 * i + 1] = quantiles[i - 1];
  }

  for (int i = 1; i <= m + 1; i++) {
    f[2 * i] = (f[2 * i - 1] + f[2 * i + 1]) / 2.0;
  }

  for (int i = 1; i <= 2 * m + 3; i++) {
    d[i] = 1.0 + 2 * (m + 1) * f[i];
  }
  isInitializing = true;
  ni = 1;
}

int Percentile::getNQuantiles() const { return m; }

void Percentile::add(const float_t x) {
  if (isInitializing) {
    q[ni++] = x;

    const int N = 2 * m + 3;

    if (ni == N + 1) {
      sort(q, q + N, greater<float_t>());
      isInitializing = false;
    }
  } else {
    addMeasurement(x);
  }
}

const float_t* Percentile::getQuantiles() const { return quantiles; }

bool Percentile::isReady() const { return !isInitializing; }

int Percentile::getNSamples() const {
  if (!isInitializing)
    return n[2 * m + 3] - 1;
  else {
    return ni - 1;
  }
}

const float_t* Percentile::getEstimates() const {
  if (!isInitializing) {
    for (int i = 1; i <= m; i++) {
      e[i - 1] = q[2 * i + 1];
    }
    return e;
  } else {
    return NULL;
  }
}

float_t Percentile::getMin() const { return q[1]; }

float_t Percentile::getMax() const { return q[2 * m + 3]; }

void Percentile::addMeasurement(const float_t x) {
  int k = 1;

  if (x < q[1]) {
    k = 1;
    q[1] = x;
  } else if (x >= q[2 * m + 3]) {
    k = 2 * m + 2;
    q[2 * m + 3] = x;
  } else {
    for (int i = 1; i <= 2 * m + 2; i++) {
      if ((q[i] <= x) && (x < q[i + 1])) {
        k = i;
        break;
      }
    }
  }

  for (int i = k + 1; i <= 2 * m + 3; i++) {
    n[i] = n[i] + 1;
  }

  for (int i = 1; i <= 2 * m + 3; i++) {
    d[i] = d[i] + f[i];
  }

  for (int i = 2; i <= 2 * m + 2; i++) {
    const float_t dval = d[i] - n[i];
    const float_t dp = n[i + 1] - n[i];
    const float_t dm = n[i - 1] - n[i];
    const float_t qp = (q[i + 1] - q[i]) / dp;
    const float_t qm = (q[i - 1] - q[i]) / dm;
    if ((dval >= 1.0) && (dp > 1.0)) {

      const float_t qt = q[i] + ((1.0 - dm) * qp + (dp - 1.0) * qm) / (dp - dm);

      if ((q[i - 1] < qt) && (qt < q[i + 1])) {
        q[i] = qt;
      } else {
        q[i] = q[i] + qp;
      }
      n[i] = n[i] + 1;
    } else if ((dval <= -1) && dm < -1) {
      const float_t qt = q[i] - ((1.0 + dp) * qm - (dm + 1.0) * qp) / (dp - dm);
      if ((q[i - 1] < qt) && (qt < q[i + 1])) {
        q[i] = qt;
      } else {
        q[i] = q[i] - qm;
      }
      n[i] = n[i] - 1;
    }
  }
}

ostream& Percentile::print(ostream& os) const {
  if (isReady()) {
    const float_t* q = getQuantiles();
    const float_t* e = getEstimates();
    const int SCREENWIDTH = 59;
    const int n = getNQuantiles();

    os << name << ", min(" << getMin() << "), max(" << getMax() << "), samples(" << getNSamples() << ")" << endl;

    const float_t max = e[n - 1];
    for (int i = 0; i < n; i++) {
      os << fixed << setprecision(3) << q[i] << ": ";

      const int len = (int)(e[i] * SCREENWIDTH / max);
      for (int j = 0; j < len; j++) {
        os << '#';
      }
      os << " ";
      os << fixed << setprecision(3) << e[i] << endl;
    }
  }
  return os;
}

percentile.h

#pragma once

#include <iostream>

namespace stat {

using float_t = double;

/**
 * Implementation of "Simulatenous Estimation of Several Persentiles," by Kimmo E. E. Raatikainen
 *
 * This is very useful for profiling the performance of timing characteristics
 *
 * Created by jcairns on 5/28/14.
 */
class Percentile {
public:
  const static float_t DEFAULT_PERCENTILE[];

private:
  const char* name;

  const float_t* quantiles;

  const int m;

  float_t* q; // heights

  int* n; // actual positions

  float_t* f; // increments of desired positions

  float_t* d; // desired positions

  float_t* e; // estimates

  bool isInitializing;

  int ni; // which x is initialized so far

  void addMeasurement(const float_t x);

public:
  Percentile();
  Percentile(const char* name);
  Percentile(const char* name, const float_t quantiles[], const int nQuantiles);

  ~Percentile();

  /**
   * clear existing samples
   */
  void clear();

  /**
   * Add a measurement to estimate
   *
   * @param x - the value of the measurement
   */
  void add(const float_t x);

  /**
   * @return float_t[] - percentiles requested at initialization
   */
  const float_t* getQuantiles() const;

  /**
   * return the number of quantiles
   */
  int getNQuantiles() const;

  /**
   * @return boolean - true if sufficient samples have been seen to form an estimate
   */
  bool isReady() const;

  /**
   * @return int - the number of samples in the estimate
   */
  int getNSamples() const;

  /**
   * get the estimates based on the last sample
   *
   * @return float_t[]
   *
   * @throws InsufficientSamplesException - if no estimate is currently available due to insufficient data
   */
  const float_t* getEstimates() const;

  /**
   * @return float_t - the minimum sample seen in the distribution
   */
  float_t getMin() const;

  /**
   * @return float_t - the maximum sample seen in the distribution
   */
  float_t getMax() const;

  /**
   * print a nice histogram of percentiles
   *
   * @param out - output stream
   * @param name - data set name
   * @param p - percentile
   *
   */
  std::ostream& print(std::ostream& os) const;
};

inline std::ostream& operator<<(std::ostream& os, const Percentile& p) { return p.print(os); }

}; // namespace stat