I am trying to learn Bayesian methods, and to that end, I picked up an application of interest to me to develop the concepts in practice.


Suppose I wrote an initial version of a performance-sensitive piece of software, and want to optimize its execution time. I may have a baseline version and an "improved" version (or at least, I suspect it may be an improvement -- I need to measure).

I'm looking to quantify how likely it is that this new version is actually be an improvement (as opposed to being equivalent or possibly even worse than the baseline), as well as how much -- is it 20% faster? 100% faster? 10% slower? Also I'd like to give credible intervals rather than only point estimates of the speedup.

To that end, I time a number of runs of the two version of the software, trying to keep all other factors the same (input data, hardware, OS, etc.) I also try to kill every running app and service, and even turn off networking, to make sure that, to the extent possible by modern feature-heavy code, these apps have the CPU all to themselves. I also disable Turbo Boost on my CPU to prevent CPU clock rate changes over time and temperature, and run my fans at maximum to minimize the change of CPU thermal throttling (and in practice my computer's thermal solution is good enough that I've never seen this happen). I've tried to restrict the portion of the code being measured to the computational part only, so no I/O to add variability.

Despite my best efforts, this is not an embedded system with a single-core processor running on bare metal, so there is some variability, possibly due to OS processes that remain and take up a little bit of CPU, CPU affinity of processes, as well as microarchitectural sources of variability such as caches, out of order execution and hyperthreading.

Current model and code

Currently I'm using the BEST model, implemented by the following code in Python using PyMC3 (heavily inspired by the linked document), in case it is of interest. The arguments are timings of the baseline version (baseline) and improved version (opt):

def statistical_analysis(baseline, opt):
    # Inspired by https://docs.pymc.io/notebooks/BEST.html
    y = pd.DataFrame(
            value=np.r_[baseline, opt],
            group=np.r_[['baseline']*len(baseline), ['opt']*len(opt)]

    μ_m = y.value.mean()
    μ_s = y.value.std()
    σ_low = µ_s/1000
    σ_high = µ_s*1000

    with pm.Model() as model:
        baseline_mean = pm.Normal('baseline_mean', mu=μ_m, sd=1000*μ_s)
        opt_mean = pm.Normal('opt_mean', mu=μ_m, sd=1000*μ_s)
        baseline_std = pm.Uniform('baseline_std', lower=µ_s/1000,
        opt_std = pm.Uniform('opt_std', lower=µ_s/1000, upper=1000*µ_s)
        ν = pm.Exponential('ν_minus_one', 1/29.) + 1
        λ_baseline = baseline_std**-2
        λ_opt = opt_std**-2

        dist_baseline = pm.StudentT('baseline', nu=ν, mu=baseline_mean,
                                    lam=λ_baseline, observed=baseline)
        dist_opt = pm.StudentT('opt', nu=ν, mu=opt_mean,
                               lam=λ_opt, observed=opt)

        diff_of_means = pm.Deterministic('difference of means',
                                         baseline_mean - opt_mean)
        ratio_of_means = pm.Deterministic('ratio of means',

        trace = pm.sample(draws=3000,tune=2000)

        baseline_hdi = az.hdi(trace['baseline_mean'])
        baseline_out = (baseline_hdi[0],

        opt_hdi = az.hdi(trace['opt_mean'])
        opt_out = (opt_hdi[0], trace['opt_mean'].mean(), opt_hdi[1])

        speedup_hdi = az.hdi(trace['ratio of means'])
        speedup = (speedup_hdi[0],
                   trace['ratio of means'].mean(),

        dif = trace['difference of means'] > 0
        prob = (dif > 0).sum()/len(dif)

    return (baseline_out, opt_out, speedup, prob)

The prob variable indicates how likely it is that a difference exists, and speedup includes the mean as well as 95% HDI for the ratio of execution time of the baseline version to the improved version. The remaining variables are the mean as well as 95% HDI of the execution time of baseline and improved versions.

Issues with the model

The BEST model assumes a Student t-distribution for the values of the execution time, but I have a hunch this is not an adequate modeling assumption.

Given a certain piece of code, one could in principle tally up every single instruction executed, and figure out exactly how fast an "undisturbed" CPU could run it, given the amount of execution resources like ALUs and load/store units, the latency of each instruction, etc. Therefore, there exists a minimum value, bounded by the CPU hardware capabilities, such that the code will never run faster than this. We cannot measure this minimum, though, because the measurements are contaminated by the sources of noise previously mentioned.

Thus, I'd like to think that my model should be the sum of a constant value (the minimum) and some distribution with positive values only, and probably a heavy tailed one, seeing as some outlier event may happen during the execution of the code (the system decides to update an app, or run a backup, or whatever).

Edit: some data

To give an idea of the kind of distribution that may be found in practice, I measured 5000 executions of the serial and a parallel version of the same code, for the same input data, and generated histograms for both, with 250 bins each. I'm not claiming this is necessarily representative, but it shows how inadequate the Student t-distribution is for this problem.

First, the serial version:

Histogram of serial version

And now for the parallel version:

enter image description here

The question

This leads me to the question:

What are some distributions that might be a good fit to this model?

  • 1
    $\begingroup$ These distributions superficially look gamma-distributed. Also, what if you log-transform the runtimes then try fitting a mixture of Gaussians? That might more closely match the sort of additive process that you were describing. $\endgroup$ Commented Oct 26, 2020 at 15:00
  • $\begingroup$ Since writing the question I've experimented with fitting gamma, lognormal and Weibull distributions. They're evidently a much better fit than a Student t-distribution, especially if I subtract the minimum of the data, but I can't achieve a fairly good fit -- and it seems related to the "hump" that is visible on the histograms, suggesting a possibility that the distribution is bimodal. Some new data that I collected shows this bimodal characteristic even more clearly. I may have to think a bit about this source of bimodality, as it's completely unexpected to me. $\endgroup$
    – swineone
    Commented Oct 26, 2020 at 15:10
  • $\begingroup$ By the way, something interesting that I also observed (coincidentally before your suggestion): the histogram of both the logarithm and the exponential of the data (as well as the original data itself) are very similar. This seems to suggest some sort of invariance to exponentials, similar to what we see in the calculation of derivatives. I have no idea if this means anything, but I thought I'd mention this to see if it helps someone figure something out. $\endgroup$
    – swineone
    Commented Oct 26, 2020 at 15:15
  • $\begingroup$ You may want to consider the Erlang distribution - see e.g. blog.newrelic.com/engineering/… $\endgroup$
    – MotiNK
    Commented Oct 26, 2020 at 16:40
  • $\begingroup$ I find this a bit difficult question for three reasons. 1 The two distributions are clearly different and it is a bit strange to be very accurate about the particular distribution. 2 When you do a paired comparison you may get even much more clear differences. Instead of comparing and analysing the two distributions, it is much easier to look at a single distribution for the difference. 3 I see several distributions being mentioned in the comments, but I'd prefer to analyze the software and see how it's performance is influenced, instead of guessing a distribution for a histogram. $\endgroup$ Commented Oct 27, 2020 at 16:05

1 Answer 1


Firstly, I don't think you really need a Bayesian approach here: you have lots of data and very little prior information. On the other hand, if you'd like to there's no harm, either, as long as your priors are sufficiently weakly informative (or informed by some reasonable prior information).

Secondly, the simplest thing to do is to log-transform the data. Clearly, run-time cannot be negative. Secondly, this might very well make a t-distribution be a pretty good approximation. As others have suggested a gamma distribution, Weibull distribution, exponential distribution or some other positive distribution may be sufficient.

Of course, you could also use non-parametric methods (that tends to be non-Bayesian), but since you don't seem to have any covariates etc., just a Wilcoxon test may be fine (and it comes with the Hodges-Lehmann estimate).

Finally, the results look so clear (if this is the real data), that you don't even need much statistics to tell one is better than the other. Or am I missing something?

  • $\begingroup$ While I haven't stopped to think about your answer with the focus it deserves, it appears to be a well thought off answer, and the only one, so I'll award the bounty. Thanks. $\endgroup$
    – swineone
    Commented Nov 2, 2020 at 22:48
  • 1
    $\begingroup$ As for the results looking so clear -- indeed in this case they do. Sometimes it's less clear cut. It's for these cases that I'm interested in the statistics. $\endgroup$
    – swineone
    Commented Nov 2, 2020 at 22:50
  • $\begingroup$ Did you look at the histogram of the log-transformed values? It's sort of hard to guess, but I'd guess they might quite possibly be close to normal (or t-distributed) and the lower limit (while it exists) may be for the purposes of an analysis be kind of irrelevant. $\endgroup$
    – Björn
    Commented Nov 3, 2020 at 0:05
  • 1
    $\begingroup$ I suppose, you could also, use Stan (or equivalently WinBugs, JAGS or your MCMC sampler of choice) and write out a full model that essentially says something like this in the model block: model{ log_min_time ~ student_t(4, 0,100); log_sd ~ normal(0,10); log_mean_extra_time ~ normal(0,10); log_time ~ student_t(4, log_sum_exp(log_min_time, log_mean_extra_time), log_sd);} (to be modified to allow different mean & SD by group; log_time is the logartihm of the observed times). $\endgroup$
    – Björn
    Commented Nov 3, 2020 at 0:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.