Stats.cpp

#include "Stats.h"

//-----------------------------------------------------------------------------

// If you want to compute these two statistics, uncomment the code and link with
// the GSL library.

/*
extern "C"
{
	double gsl_sf_gamma_inc_P(const double a, const double x);
	double gsl_sf_gamma_inc_Q(const double a, const double x);
};

// P-val for a set of binomial distributions

void pval_binomial ( int * buckets, int len, int n, double p, double & sdev, double & pval )
{
	double c = 0;

	double u = n*p;
	double s = sqrt(n*p*(1-p));

	for(int i = 0; i < len; i++)
	{
		double x = buckets[i];

		double n = (x-u)/s;

		c += n*n;
	}

	sdev = sqrt(c / len);

	pval = gsl_sf_gamma_inc_P( len/2, c/2 );
}

// P-val for a histogram - K keys distributed between N buckets
// Note the (len-1) due to the degree-of-freedom reduction

void pval_pearson ( int * buckets, int len, int keys, double & sdev, double & pval )
{
	double c = 0;

	double n = keys;
	double p = 1.0 / double(len);

	double u = n*p;
	double s = sqrt(n*p*(1-p));

	for(int i = 0; i < len; i++)
	{
		double x = buckets[i];

		double n = (x-u)/s;

		c += n*n;
	}

	sdev = sqrt(c / len);

	pval = gsl_sf_gamma_inc_P( (len-1)/2, c/2 );
}
*/

//----------------------------------------------------------------------------

double erf2 ( double x )
{
    const double a1 =  0.254829592;
    const double a2 = -0.284496736;
    const double a3 =  1.421413741;
    const double a4 = -1.453152027;
    const double a5 =  1.061405429;
    const double p  =  0.3275911;

    double sign = 1;
    if(x < 0) sign = -1;

    x = abs(x);

    double t = 1.0/(1.0 + p*x);
    double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

    return sign*y;
}

double normal_cdf ( double u, double s2, double x )
{
	x = (x - u) / sqrt(2*s2);

	double c = (1 + erf2(x)) / 2;

	return c;
}

double binom_cdf ( double n, double p, double k )
{
	double u = n*p;
	double s2 = n*p*(1-p);

	return normal_cdf(u,s2,k);
}

// return the probability that a random variable from distribution A is greater than a random variable from distribution B

double comparenorms ( double uA, double sA, double uB, double sB )
{
	double c = 1.0 - normal_cdf(uA-uB,sA*sA+sB*sB,0);

	return c;
}

// convert beta distribution to normal distribution approximation

void beta2norm ( double a, double b, double & u, double & s )
{
	u = a / (a+b);

	double t1 = a*b;
	double t2 = a+b;
	double t3 = t2*t2*(t2+1);

	s = sqrt( t1 / t3 );
}

#pragma warning(disable : 4189)

double comparecoins ( double hA, double tA, double hB, double tB )
{
	double uA,sA,uB,sB;

	beta2norm(hA+1,tA+1,uA,sA);
	beta2norm(hB+1,tB+1,uB,sB);

	// this is not the right way to handle the discontinuity at 0.5, but i don't want to deal with truncated normal distributions...

	if(uA < 0.5) uA = 1.0 - uA;
	if(uB < 0.5) uB = 1.0 - uB;

	return 1.0 - comparenorms(uA,sA,uB,sB);
}

// Binomial distribution using the normal approximation

double binom2 ( double n, double p, double k )
{
	double u = n*p;
	double s2 = n*p*(1-p);

	double a = k-u;

	const double pi = 3.14159265358979323846264338327950288419716939937510;

	a = a*a / (-2.0*s2);
	a = exp(a) / sqrt(s2*2.0*pi);

	return a;
}

double RandWork ( double bucketcount, double keycount )
{
	double avgload = keycount / bucketcount; 

	double total = 0;

	if(avgload <= 16)
	{
		// if the load is low enough we can compute the expected work directly

		double p = pow((bucketcount-1)/bucketcount,keycount);

		double work = 0;

		for(double i = 0; i < 50; i++)
		{
			work  += i;
			total += work * p;

			p *= (keycount-i) / ( (i+1) * (bucketcount-1) );
		}
	}
	else
	{
		// otherwise precision errors screw up the calculation, and so we fall back
		// to the normal approxmation to the binomial distribution

		double min = avgload / 5.0;
		double max = avgload * 5.0;

		for(double i = min; i <= max; i++)
		{
			double p = binom2(keycount,1.0 / bucketcount,i);

			total += double((i*i+i) / 2) * p;
		}
	}

	return total / avgload;
}

// Normalized standard deviation.

double nsdev ( int * buckets, int len, int keys )
{
	double n = len;
	double k = keys;
	double p = 1.0/n;

	double u = k*p;
	double s = sqrt(k*p*(1-p));

	double c = 0;

	for(int i = 0; i < len; i++)
	{
		double d = buckets[i];

		d = (d-u)/s;

		c += d*d;
	}

	double nsd = sqrt(c / n);

	return nsd;
}


double chooseK ( int n, int k )
{
    if(k > (n - k)) k = n - k;

	double c = 1;

	for(int i = 0; i < k; i++)
	{
		double t = double(n-i) / double(i+1);

		c *= t;
	}

    return c;
}

double chooseUpToK ( int n, int k )
{
	double c = 0;

	for(int i = 1; i <= k; i++)
	{
		c += chooseK(n,i);
	}

	return c;
}

//-----------------------------------------------------------------------------

uint32_t bitrev ( uint32_t v )
{
	v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
	v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
	v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
	v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
	v = ( v >> 16             ) | ( v               << 16);

	return v;
}

//-----------------------------------------------------------------------------

// Distribution "score"
// TODO - big writeup of what this score means

// Basically, we're computing a constant that says "The test distribution is as
// uniform, RMS-wise, as a random distribution restricted to (1-X)*100 percent of
// the bins. This makes for a nice uniform way to rate a distribution that isn't
// dependent on the number of bins or the number of keys

// (as long as # keys > # bins * 3 or so, otherwise random fluctuations show up
// as distribution weaknesses)

double calcScore ( std::vector<int> const & bins, int keys )
{
	double n = (int)bins.size();
	double k = keys;

	// compute rms value

	double r = 0;

	for(size_t i = 0; i < bins.size(); i++)
	{
		double b = bins[i];

		r += b*b;
	}

	r = sqrt(r / n);

	// compute fill factor

	double f = (k*k - 1) / (n*r*r - k);

	// rescale to (0,1) with 0 = good, 1 = bad

	return 1 - (f / n);
}


//----------------------------------------------------------------------------

void plot ( double n )
{
	double n2 = n * 1;

	if(n2 < 0) n2 = 0;

	n2 *= 100;

	if(n2 > 64) n2 = 64;

	int n3 = (int)floor(n2 + 0.5);

	if(n3 == 0)
		printf(".");
	else
	{
		char x = '0' + char(n3);

		if(x > '9') x = 'X';

		printf("%c",x);
	}
}

//-----------------------------------------------------------------------------