Determine if two distributions are the same

Question

I develop GNU Parallel. For each release I would like to test whether the version I am about to release has a memory leak.

I can measure the memory usage. Unfortunately the usage varies slightly given the same input. The memory usage is: base + leak. The base varies slightly and there is an expected leak of 1 bit per job. What I want to detect is if there is a bigger leak than the expected leak.

For one release the base is between, say, 13050 and 13150 kB and for another release it can be between 15000 and 15500 kB. They are more or less normally distributed.

It is very easy to see if there is a memory leak: Just compare runs of 1000 jobs to 1000000 jobs and if there is a leak of 1 byte/job then the usage will increase by 1000000 bytes which is bigger than the range of the base memory usage.

But running 1000000 jobs will take around 10000 seconds (1 job takes around 0.01 sec to run). So is there a faster way in which I can determine if there is likely a bigger than expected leak?

I am thinking of running 10 runs with 1000 jobs and 10 runs with 2000 jobs and comparing the distributions of memory usage of these. If the distributions are significantly different then there is a leak. How do I do that?

Intuitively I can see that you need fewer runs and fewer jobs if the leak is big (say, 500 bytes per job) than if the leak is small (say, 2 bytes per job), and thus it should be possible to stop much earlier. How do I figure out how many runs and how many jobs I should run to get a 99% certainty?

(I am a novice at statistics, R and Python, and advanced in Perl. So a solution will have to show real code, not just refer to some statistical methods; also code that gives a "yes/no" is preferred over "eyeballing it").

((Ideally I want 2 pieces of code: One that I can tell "You have 3 minutes to run. After that I want to know how certain you are about a memleak of 1 byte, 10 bytes, 100 bytes, and 1000 bytes", and the second version that show me the certainties continously while it is running more runs to gather more datapoints.))

Answer 1

Kolmogorov–Smirnov statistic may help you in this case.
Following is an implementation which uses Kolmogorov-Smirnov statistic and the function returns the probability of similarity.

#include <math.h>
#define EPS1 0.001
#define EPS2 1.0e-8

float kstest(float alam) {
    int j;
    float a2, fac = 2.0, sum = 0.0, term, termbf = 0.0;
    a2 = -2.0 * alam * alam;
    for (j = 1; j <= 100; j++) {
    term = fac * exp(a2 * j * j);
    sum += term;
    if (fabs(term) <= EPS1 * termbf || fabs(term) <= EPS2 * sum)
        return sum;
    fac = -fac;
    termbf = fabs(term);
    }
    return 1.0;
}

void checkSameDist(float data1[], unsigned long n1, float data2[],
          unsigned long n2, float *d, float *prob) {
    float           kstest(float alam);
    void            sort(unsigned long n, float arr[]);
    unsigned long   j1 = 1,
    j2 = 1;
    float d1, d2, dt, en1, en2, en, fn1 = 0.0, fn2 = 0.0;
    sort(n1, data1);
    sort(n2, data2);
    en1 = n1;
    en2 = n2;
    *d = 0.0;
    while (j1 <= n1 && j2 <= n2) {
        if ((d1 = data1[j1]) <= (d2 = data2[j2]))
            fn1 = j1++ / en1;
        if (d2 <= d1)
            fn2 = j2++ / en2;
        if ((dt = fabs(fn2 - fn1)) > *d)
            *d = dt;
    }
    en = sqrt(en1 * en2 / (en1 + en2));
    *prob = kstest((en + 0.12 + 0.11 / en) * (*d));
}

Also check, following function checks if a particular distribution is normal, you could modify it a little bit (this would give you more intuition about the statistic and how can you implement it from scratch
(https://walteis.wordpress.com/2012/04/26/a-kolmogorov-smirnov-implementation/)

public bool IsNormal
{
    get
    {
        // This method uses the Kolmogorov-Smirnov test to determine a normal distribution.
        // The level of significance (alpha) used is .05, and the critical values used are from Table 1 of: 
        // The Kolmogorov-Smirnov Test for Goodness of Fit
        // Frank J. Massey, Jr.
        // Journal of the American Statistical Association
        // Vol. 46, No. 253 (Mar., 1951) (pp. 68-78)

        if (DataSet.Count == 0)
            return false;

        List<double> vals = DataSet.Values.ToList();
        Accumulator acc = new Accumulator(vals.ToArray());
        double dmax = double.MinValue;
        double cv = 0;

        MathNet.Numerics.Distributions.NormalDistribution test = new MathNet.Numerics.Distributions.NormalDistribution(acc.Mean, acc.Sigma);

        // the 0 entry is to force the list to be a base 1 index table.
        List<double> cvTable = new List<double>() { 0, .975, .842, .708, .624, .565,
                                                .521, .486, .457, .432, .410,
                                                .391, .375, .361, .349, .338,
                                                .328, .318, .309, .301, .294};

        test.EstimateDistributionParameters(DataSet.Values.ToArray());
        vals.Sort();

        for (int i = 0; i < vals.Count; i++)
        {
            double dr = Math.Abs(((i + 1) / (double)vals.Count) - test.CumulativeDistribution(vals[i]));
            double dl = Math.Abs(test.CumulativeDistribution(vals[i]) - (i / (double)vals.Count));
            dmax = Math.Max(dmax, Math.Max(dl, dr));
        }

        // get critical value and compare to d(N)
        if (vals.Count <= 10)
            cv = cvTable[vals.Count];
        else if (vals.Count > 10)
            cv = 1.36 / Math.Sqrt(vals.Count);

        return (dmax < cv);
    }
}

Best of Luck