OnePassDescriptiveStats.hpp

#pragma once
/*
   Copyright 2015 Peter Gaultney

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
*/
#include <limits>

// This class is not thread safe.
class OnePassDescriptiveStats
{
public:
	void addValue(double value);
	double getVariance() const;
	double getStddev() const;

	OnePassDescriptiveStats aggregateWithSet(const OnePassDescriptiveStats& B) const;

	unsigned long long count = 0;
	double mean = 0.0;

	double min = std::numeric_limits<double>::max();
	double max = std::numeric_limits<double>::min();
		
	double M2 = 0.0;

public:
	static OnePassDescriptiveStats aggregateSets(
		const OnePassDescriptiveStats& A,
		const OnePassDescriptiveStats& B);
};

/*
 * Implementation sourced from http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
 * which is itself (mostly) due to Donald Knuth, according to Wikipedia.
 */
#include <cmath>

inline void OnePassDescriptiveStats::addValue(double value)
{
	if (value > this->max) {
		this->max = value;
	}
	if (value < this->min) {
		this->min = value;
	}

	++ this->count;
	double delta = value - this->mean;
	this->mean = this->mean + delta / this->count;
	this->M2 = this->M2 + delta * (value - this->mean);
}

/**
 * This is the sample, not population, variance.
 */
inline double OnePassDescriptiveStats::getVariance() const
{
	if (this->count < 2) {
		return 0.0;
	} else {
		return this->M2 / (this->count - 1);
	}
}

inline double OnePassDescriptiveStats::getStddev() const
{
	return sqrt(this->getVariance());
}

inline OnePassDescriptiveStats OnePassDescriptiveStats::aggregateWithSet(
	const OnePassDescriptiveStats& B) const
{
	return OnePassDescriptiveStats::aggregateSets(*this, B);
}

inline OnePassDescriptiveStats OnePassDescriptiveStats::aggregateSets(
	const OnePassDescriptiveStats& A,
	const OnePassDescriptiveStats& B)
{
	OnePassDescriptiveStats combined;
	// this algorithm due to Chan et al., also from Wikipedia page above
	double delta = B.mean - A.mean;
	combined.count = A.count + B.count;
	if (combined.count > 0) {
		combined.mean = (A.count * A.mean + B.count * B.mean) / combined.count;
		combined.M2 = A.M2 + B.M2 + delta * delta * ((double)(A.count * B.count) / combined.count);
	} else { // if we're combining two empty sets, we don't want floating point exceptions.
		combined.mean = 0;
		combined.M2 = 0;
	}

	// mins, maxes
	if (B.max > A.max) {
		combined.max = B.max;
	} else {
		combined.max = A.max;
	}
	if (B.min < A.min) {
		combined.min = B.min;
	} else {
		combined.min = A.min;
	}

	return combined;
}