# Quantifying Performance Changes with Effect Size Confidence Intervals

TOMAS KALIBERA  RICHARD JONES  University of Kent

###### Abstract

Measuring performance and quantifying a performance change are core evaluation
techniques in programming language and systems research. Out of 122 recent
scientific papers published at PLDI, ASPLOS, ISMM, TOPLAS, and TACO, as many
as 65 included experimental evaluation that quantified a performance change
using a ratio of execution times. Unfortunately, few of these papers evaluated
their results with the level of rigour that has come to be expected in other
experimental sciences. The uncertainty of measured results was largely
ignored. Scarcely any of the papers mentioned uncertainty in the ratio of the
mean execution times, and most did not even mention uncertainty in the two
means themselves. Furthermore, most of the papers failed to address the non-
deterministic execution of computer programs (caused by factors such as memory
placement, for example), and none addressed non-deterministic compilation
(when a compiler creates different binaries from the same sources, which
differ in performance, for example again because of impact on memory
placement). It turns out that the statistical methods presented in the
computer systems performance evaluation literature for the design and summary
of experiments do not readily allow this either. This poses a hazard to the
repeatability, reproducibility and even validity of quantitative results.

Inspired by statistical methods used in other fields of science, and building
on results in statistics that did not make it to introductory textbooks, we
present a statistical model that allows us both to quantify uncertainty in the
ratio of (execution time) means and to design experiments with a rigorous
treatment of those multiple sources of non-determinism that might impact
measured performance. Better still, under our framework summaries can be as
simple as “system A is faster than system B by 5.5% $\pm$ 2.5%, with 95%
confidence”, a more natural statement than those derived from typical current
practice, which are often misinterpreted.

###### category:

D.2.8 Software Engineering Metrics

###### keywords:

Performance measures

###### category:

D.3.4 Programming Languages Processors

###### keywords:

Run-time environments

###### keywords:

statistical methods, random effects, effect size

††terms: Experimentation, Measurement, Performance

Author’s addresses: T. Kalibera (R.E. Jones), School of Computing, University
of Kent, Canterbury, CT2 7NF, UK.

_A preliminary version of a portion of this work was presented at the Third
European Performance Engineering Workshop._

##  1 Introduction

Quantification of performance change is a common task in experimental computer
science. Out of all 122 papers published in ASPLOS, ISMM, PLDI, TACO and
TOPLAS in 2011 up to 2nd August, 65 included empirical evaluation that
quantified performance change by giving the ratio of execution times (for
example, speed-ups of optimisations, or overheads of new techniques).
Quantification of performance change is also part of the software development
process (for example, detecting performance regressions), both open-source and
commercial.

Execution times for quantification are obtained by running benchmark
applications using the two systems to compare. Computer systems are becoming
more complex and increasingly resemble the direct products of nature observed
by physicists and natural scientists — an overwhelming number of factors
influence performance. Some are unknown, some are out of the experimenter’s
control, and some are non-deterministic. Repeated executions of the same
benchmark on the same system, even when the experimenter does everything
possible to enforce the same conditions, always report slightly different
execution times. This is even more pronounced if the same experiment is
repeated by an independent experimenter using different equipment. Thus, for
credible quantification of performance change, one needs to design experiments
and summarize their results in a repeatable and reproducible manner.
Experiment design and statistical inference are the fields of statistics that
address these issues, and are widely applied in physics, natural sciences and
social sciences.

Unfortunately, the practice in computer science lags behind. Of the 122 papers
at ASPLOS, ISMM, PLDI, TACO, and TOPLAS we looked at, 90 evaluated execution
time based on experiments. This includes ratio of execution times for two
systems and also the execution times for a single system. 71 of these 90
papers completely ignored the question of uncertainty in the measured times.
This lack of rigour makes repeatability difficult and undermines the validity
of the results. We hardly dare mention that simple rules for summarizing and
reporting uncertainties are taught to every student of physics and natural
sciences [Kirkup (1994)], are part of engineering practice [Taylor and Kuyatt
(1994)], and are based on statistics taught in introductory courses. Moreover,
advice on elementary experiment design and statistical inference has been
readily available to computer scientists in the form of textbooks aimed and
used for teaching [Jain (1991), Lilja (2000)]. A subset of the statistics of
the latter textbook has been advocated in the context of Java runtime
performance evaluation [Georges et al. (2007)] and compared with the practice
in scientific papers, with results not particularly flattering to our field.

Nevertheless, even though the best quantification method recommended so far in
our field ([Georges et al. (2007)] based largely on [Lilja (2000)]) is ahead
of current practice, it has a number of flaws. It exacts a high price for
statistical rigour in summarizing results: _we do not get a reliable estimate
of the metric that we are ultimately interested in, which is the ratio of
execution times_. Maybe this flaw on its own is a significant detractor from
wider adoption of the method in practice. Even worse, the rigour that the
method does provide is highly questionable, as the method is based on
statistical significance, a concept with a number of shortcomings, some of
which have been known for seven decades [Cohen (1994)]. While statistical
significance still dominates elementary statistical textbooks and is still
used by some researchers even in other fields, its limits have been well
described and arguments for its deprecation have been made in statistics
[Royall (1986)], psychology [Cohen (1994)], education [Coe (2002)], ecology,
medicine, bio-medicine, and biology [Nakagawa and Cuthill (2007)]. Critical
views of statistical significance have also been published in the context of
research in sociology, criminology, economics, marketing, chemistry, and
nursing (see fidlerps for references and details). Some journals explicitly
require alternative methods [Hill and Thompson (2005)]. And such methods are
available. As we demonstrate in this paper, the specialized statistical
literature [Fieller (1954), Davison and Hinkley (1997)] provides methods to
estimate the ratio of execution times we are interested in — inference of the
ratio of means is possible, and a so-called _effect size_ confidence interval
can be constructed for it.

Apart from the problems with the summary of results, georges07’s
quantification method, while ahead of the prevailing practice in the field,
lacks rigour in experimental design. It is based on a two-level 111Please note
that we use the term ‘level’ informally. In our text, it is _not_ what‘factor
level’ means in statistics. hierarchical experiment that repeats executions of
a benchmark, each consisting of repeated measurements. For example, to
evaluate an optimisation in a Java virtual machine, the method requires us to
run many invocations of the JVM and have each invocation run and measure many
iterations of a test application. Two levels are not enough. It is well known
that a large number of factors that influence performance are inherently non-
deterministic or need to be randomized by experimenter to avoid measurement
bias [Mytkowicz et al. (2009)]. The first group includes for example context
switches, hardware interrupts, memory placement due to virtual-to-physical
memory mapping [Kalibera et al. (2005)], randomized algorithms in compilation
[Kalibera and Tuma (2006)], or decisions of a just-in-time compiler on which
methods to compile. Non-deterministic compilation by itself means that
compilation needs to be repeated. The second group includes symbol names of
methods and variables in source code which impacts binary code layout [Gu et
al. (2004)], the size of the UNIX process environment or the linking order
[Mytkowicz et al. (2009)].

To efficiently run experiments with such as number of sources of non-
determinism, we need more than two levels. georges07’s quantification method,
however, does not allow this. It is, in fact, based on one-level only
summarization, which looks at the means from the lowest level (repeated
measurements) instead of the measurement themselves. There is a better way.
While not included in elementary statistical textbooks, there are statistical
methods which model so-called random effects [Charles E. McCulloch (2008)],
which can represent non-deterministic factors in our experiments. To our
knowledge, however, none of the existing models can be used directly. Thus, we
extend one such model to be more robust and to support an arbitrary number of
experiment levels. We also derive the optimum number of repetitions at each
level (with the best recommended method so far, there is no systematic
guidance to select the number of repetitions).

Our contributions in this work are:

  1. 1.

Statistical Inference for Quantification. A statistical model of a hierarchy
of random effects that contribute to fluctuations in performance. This allows
us to build a confidence interval for the mean performance metric (such as
execution time) within such a hierarchy. This is our first original
contribution to statistics.

  2. 2.

Experiment Planning for Quantification. After some initial measurements, the
model can also be used to guide the experiment. Repeating all experiments
blindly is often very expensive. However, the model will give the optimum
number of repetitions of experiments at different levels of the hierarchy that
are necessary to provide the narrowest confidence interval in available
experimentation time. This is our second original contribution to statistics.

  3. 3.

Asymptotic Parametric Quantification Method. A parametric approach to
constructing a confidence interval for the ratio of mean metric values, based
on Fieller’s theorem [Fieller (1954)], asymptotic normality of the mean, and
our statistical model. Hence, for example, the result of a study might be that
the change is 5.5% $\pm$ 2.5%, with 95% confidence. Such a statement provides
a clear but rigorous account of the magnitude of the performance change and
the uncertainty involved.

  4. 4.

Bootstrap Quantification Method. An alternative non-parametric approach to
constructing the same confidence interval using hierarchical random re-
sampling (a bootstrap method). This is an application of an existing
statistical method, used in our field much less than it deserves.

  5. 5.

Evaluation. Thorough evaluation of the two approaches on a set of benchmarks.
We estimate the true coverage of the intervals that can be obtained. We also
estimate a false alarm rate (the case where a change is detected although
there is none), given a threshold of how large changes need to be for us to
care. A particularly important benefit of our method is that it supports a
rigorous way to use such a threshold. This is our third original contribution.

Contributions 1, 2, 4, and partially 3 and 5 also apply to evaluations of a
single system. For contributions 1 and 2, we extend and correct a preliminary
version published in [Kalibera and Tuma (2006)], which in turn extends
[Kalibera et al. (2005)]. For item 5, we use our benchmark data from [Kalibera
and Tuma (2006)] as inputs for (new) statistical simulations.

##  2 BACKGROUND: Current Quantification Methods

In this section, we provide an overview of quantification methods currently
used, we describe their limitations, and we summarize the best-so-far
quantification methods recommended for overcoming these limitations.

###  2.1 Currently Used Methods

To map the current quantification methods used in computer science, we
analysed papers published at selected conferences (PLDI, ASPLOS, ISMM) and
journals (TOPLAS, TACO). We restricted our survey to papers published between
1st January and 2nd August, 2011. We thus cover all papers from the
conferences, issues 1–4 of TOPLAS, and issues 1–2 of TACO — 122 papers in
total.

As summarized in Table 2.1, nearly three quarters of the papers studied
included empirical performance evaluation that measured and reported execution
time. The majority of the execution time evaluations (over 70%) reported the
ratio of execution times as a measure of performance change (the remaining 30%
reported execution time for a single system). In total, over half of all
papers reported the ratio of execution times. We also noted papers where
execution time was important for the evaluation (counts in parentheses in the
table). More than half of the papers had empirical evaluations with execution
time as an important metric. Out of these papers where execution time was an
important metric, in nearly 90% it was the ratio of execution times that was
reported (only in 10% of cases it was the absolute execution time for a single
system). Nearly half of all papers had empirical evaluation with ratio of
execution times as an important metric. The ratio of execution times thus
seems to be a metric that people care about and report. Anecdotal evidence
confirms that this metric is also important to industrial developers.

In practice, where decisions have to be made based on the results of
performance quantification, there are certain minimum thresholds for relative
performance change above/below which it is considered of no practical
interest. These thresholds may depend on the evaluation context (a particular
system, regression test, optimisation, new feature). These thresholds are also
assumed and used in georges07. Of the papers in our survey, thresholds on
ratio of execution times are explicitly part of an auto-tuning algorithm by
mapreduce, where they have to be set (manually, based on experience) large
enough to cater for variabilities in the data. However, most papers that
quantify performance change report only the ratio of (mean) execution times,
and thus imply that this ratio alone, maybe along with an understanding of the
context of the study, is sufficient to assess importance of a performance
change. The use of ‘thresholds’ to assess importance is implied. We also have
confirmation from a senior member of a well-known industrial research
laboratory of the use of such thresholds in industry.

Current Practice of Performance Quantification  PLDI ASPLOS ISMM TOPLAS TACO
Total Number of Papers 55 32 13 13 9 122 Evaluated Execution Time 42 (28) 25
(20) 12 (11) 5 (2) 6 (6) 90 (67) Ignored Uncertainty In Measurement 39 (24) 18
(12) 5 (5) 4 (1) 5 (5) 71 (47) Evaluated Ratio of Execution Times 27 (25) 22
(19) 9 (8) 1 (1) 6 (6) 65 (59) Source:Scientific papers at selected venues
published in 2011, up to 2nd August. Counts are given in form ‘A (M)’, where A
is the number of all papers that reported the metric, and M the number of
papers where the metric was important or the main result in the evaluation.

Nearly 80% of the papers studied that evaluated execution time (both absolute
and relative) failed to mention anything about uncertainty in the figures they
reported. The numbers are not that much better if we only focus on papers
where execution time was an important metric: 70% of these papers failed to
mention uncertainty. Note that we were very generous in the classification —
if a paper reported that uncertainty was low, without any further details, we
gave it the benefit of the doubt and counted it as having taken uncertainty
into account. Ignoring uncertainty represents a severe lack of rigour in face
of the fact that execution time on today’s computer systems is always subject
to variation in performance. While majority of the papers studied ran
experiments on real hardware, it should be mentioned that some (particularly
from ASPLOS) were based on execution time measurements in a simulator. It is
certainly possible to have a simulator in which execution time is
deterministic, and we did not attempt to check the particular simulators used
in individual papers we analysed. However, if a simulator is deterministic, it
is not a realistic one [Alameldeen and Wood (2003)]. Thus, we have not
attempted to extract simulations from the summary — using a deterministic
simulator in our opinion falls into the category of ignoring uncertainty.

We used to hear arguments that, for very large speed-ups (say 2$\times$),
summarizing uncertainty is a waste of effort. This implies that experimenter
experience is that, in the given system, uncertainty would be much smaller.
However, in systems research much smaller performance changes are usually
reported. envsize found that out of 88 papers from ASPLOS 2008, PACT 2007,
PLDI 2007, CGO 2007 with experimental evaluation in a dedicated section, the
median of reported speed-ups was only 10%. In their work they show that the
measurement bias with SPEC CPU 2006 benchmarks can easily obfuscate this
speed-up when evaluating compiler optimisations. This means proper
randomization to avoid this bias would have lead to higher uncertainty than
10%, and hence proper handling of this uncertainty would be necessary (such a
randomization tool was later provided by stabilizer). A senior member of a
well-known industrial research lab told us that he would care only about
differences of 10% for new work, but look for differences of as little as 2-3%
or even 1% when looking for regressions (where he would look at uncertainty
too). Since such small performance changes are often of interest, ignoring
uncertainty in systems research represents a serious threat to validity. When
speed-ups are large compared to expected uncertainty, validity would not be
threatened, but ignoring it would still be a serious lack of rigour and in our
view not a property of good research.

Out of the 19 papers that touched on uncertainty in any way, most reported
uncertainty of single systems in the form of a standard deviation or of 95%
confidence intervals. It was exceptional for a paper to specify how the
confidence intervals were constructed; most were presumably based on the
normal distribution. While most papers reported the ratio of execution times,
they mostly showed uncertainty only for the means in isolation of each of the
two systems compared. Only three of the papers reported the uncertainty of the
ratios. However, it seems from the text of two of these that they did not take
uncertainty of both systems into account, but only normalized the uncertainty
of the new system against the mean of the old system. The third paper was
ours, and it took the uncertainty of both systems into account, but we have to
admit that our text did not say that. The remaining papers which address
uncertainty do so by saying that it was low (without giving numbers), and one
says that a proposed optimisation reduced it (without giving numbers).

The papers we studied gave extremely little information on the experimental
design used. From some papers, it is not even obvious if benchmarks were
executed more than once, and hence if the reported execution time is not just
a single measurement. Only very few papers mention that their experiments
repeat at two levels (i.e. execution of a benchmark and measurement within
that benchmark). However, all of those that do also seem to treat the data as
if it came from a one-level experiment. If executions were repeated, the
numbers of repetitions seemed to be arbitrary. Some experiments used an
adaptive number of iterations based on a heuristic — the benchmark code itself
calculates the standard deviation of the last few measurements taken and
terminates the experiment once the standard deviation is sufficiently small.

Some papers claim to use certain evaluation methodologies, such as those of
envsize, georges07 or georges08. However, these papers rarely made it clear
how the concrete experiments were performed. One paper claims to use the
method proposed by envsize to reduce variance, without any further details.
However, the cited paper argues for dealing with measurement bias through
randomization, which in fact increases the variance, but for very good
reasons. Some papers mention uncertainty due to the decisions of a just-in-
time compiler of which methods to compile. Confidence intervals or standard
deviation error bars are usually shown for a single system with one particular
configuration, even when multiple configurations were used. Uncertainty
quantified in this way is prone to bias.

##### Summary

The ratio of (mean) execution times is the major metric for quantification of
performance change. A threshold for the minimum/maximum ratio is used as a
measure of the practical importance of a difference. Uncertainty in the ratio
of execution times is almost never quantified, and when it is, only in an
unclear way. Sometimes separate uncertainties for the two systems (means) are
shown, but more commonly even these are ignored. Benchmarks are mostly
repeated only at one level, or exceptionally at two levels (executions and
measurements). Two-level experiments, if done at all, are evaluated as if they
were one-level only (we believe by merging all the measurements from the
different executions, but the texts are not explicit here).

###  2.2 What Is Wrong With Currently Used Methods

Any paper that reports execution time, yet ignores the uncertainty of the
final metric, be it absolute execution time of one system or the ratio of
execution times for two systems, is methodologically flawed. Moreover, failing
to address uncertainty makes the results unrepeatable. Repeating any
experiment will lead to slightly different results due to variability in
measurement. Without quantified uncertainty, we cannot know if we got a
‘slightly’ different result due to variability in measurement, or a very
different result due to an error. Note that a failure to address uncertainty
is not the only way to get unrepeatable results. A common problem is
measurement bias, when the experimenter fails to identify and then describe or
vary a factor that largely influences the results [Mytkowicz et al. (2009)].
In such a case, statistics even properly used for factors that were varied
cannot help and the summarized uncertainty would be misleadingnly small. The
responsibility for identifying these factors lies fully on the experimenter
and no statistics can help. But even when important factors are fixed and
described or varied properly, failing to address uncertainty is a potential
threat to validity when quantifying performance change, and this threat is
particularly high when the observed change is relatively small (e.g. 10%),
which it often is. When such quantification is used to show that an
optimisation provides a good speed-up, there is a risk that the ‘improvement’
lies within the range of uncertainty. A similar risk applies to some extent
when quantification is used to show that the overhead of a certain feature is
small — if the uncertainty of the new system is much larger than the old one,
the price may in practice be high even if the mean overhead is small.

Failure to report how many repetitions of a benchmark were executed makes
repeatability harder. Often, the reader is only left with suspicion that
insufficient repetitions were made and in particular that the researchers did
not repeat whole executions of benchmarks (that in turn repeat measurements
within a single execution). Sometimes, the reported number of measurements
used for averaging is suspiciously small (e.g. three or five). Failure to
repeat enough times, or indeed to repeat at all, increases the risk of
obtaining unrepeatable results and reaching incorrect conclusions.

In summary, our survey of published papers suggests that the current practice
of quantification of performance change in scientific papers is weak. We are
not the first to find this, but we were surprised at how prevalent this
practice is even in the most recent papers. The lack of reporting of
uncertainties while quantifying a performance change in Java runtimes research
was documented by georges07. envsize discovered that measurement bias was not
addressed in experimental computer systems performance evaluation; this is a
more general problem, but it also applies to the quantification of performance
change. These shortcomings have arisen despite well-known textbooks on
computer systems performance evaluation [Jain (1991), Lilja (2000)] having
described a number of errors common in performance evaluation.

###  2.3 Currently Proposed Methods

Experimental design and statistical inference are mature fields of
mathematical statistics with many results and described in a number of texts.
Basic recommendations on conducting experiments in physics and natural
sciences, measuring uncertainties, calculating propagation of uncertainties,
and reporting the results are available in textbooks such as [Kirkup (1994)].
More advanced rules of propagation of uncertainties can be found in a
recommendation by NIST [Taylor and Kuyatt (1994)]. High-level overview and
pointers to works on experiment design and the scientific method in general
can be found in [Wilson (1952)]. While physics does not deal with exactly the
same problems as experimental computer science, many ideas are general.
Moreover, basic methods of statistical inference and experimental design, in
the context of computer systems performance evaluation, have then been
presented in books [Jain (1991), Lilja (2000)]. These books have been and are
being used as textbooks for teaching performance evaluation courses. georges07
argued for wider adoption of a performance change quantification method based
on Lilja’s book, in the context of performance evaluation of Java runtimes.
They also added some practical recommendations for the number of repetitions
of benchmark measurements.

In this section, we summarize the current recommended best practice for
quantification of performance change, based on jain, lilja, and georges07.

####  2.3.1 Experiment Design

georges07 recommends repeating measurements within a single benchmark, as well
as repeating executions of the whole benchmark. The number of measurements
within a benchmark is controlled by a heuristic in the benchmark — as soon as
the standard deviation of the last few measurements is small enough, the
benchmark’s execution is finished. The number of benchmark executions is an
arbitrary choice of the experimenter, or, alternatively, additional executions
are adaptively added until confidence interval for the mean is sufficiently
narrow (its width is within 1% or 2% of the sample mean).

lilja and particularly jain also recommend experimental designs that cope with
various combinations of factors that can be fixed by the experimenter. These
designs apply to factors that can take only a finite, usually small, number of
values (i.e. a particular processor or benchmark). The effect of choosing a
particular value (say the processor) on performance is of interest to the
experimenter. Such effects are called _fixed effects_ in statistics (see
[Charles E. McCulloch (2008)] and [Searle et al. (1992)] for more detailed
discussion of fixed and random effects). The experimental designs then help in
planning which combinations of values of which factors to run to get
reasonable information given reasonable experimentation time. The methods
recommended by jain and lilja, however, do not address effects of factors that
can take a large number of values, usually at random and often not under our
control, for which the performance effect of a particular choice is not
interesting. Such a factor can be, say, a particular mapping of virtual to
physical pages chosen by the operating system. In statistics, these are called
_random effects_. In this work, we only address random effects. Our work could
be extended to include fixed effects, as there indeed are statistical models
that include both types of effects [Charles E. McCulloch (2008)], but we do
not attempt it here.

georges08 proposes a particular experiment design for performance evaluations
within managed environments with just-in-time compilation. The design assumes
generation of several compilation plans and then running both systems using
these plans. The motivation for this design, in contrast to just keeping the
compilation plans random (unpaired) for the two systems, is smaller
uncertainty within given experimentation time.

####  2.3.2 Summary of Experiment Results

A performance change is quantified using statistical significance [Georges et
al. (2007), Lilja (2000), Jain (1991)]. The output is a binary decision: “the
systems (may) have the same performance” or “(it is likely that) the systems
do not have the same performance”. This decision is based on the probability
that the actually observed difference or larger in (sample) means of the two
systems would occur if the (true) means were the same. This probability is
called the _p-value_ in statistical tests. The p-value is compared against a
pre-defined threshold called _significance level_ (i.e. 5%). If the p-value is
smaller, then the second decision, that the systems do not perform the same,
is chosen. Otherwise, the first decision, the default that they may have the
same performance, is chosen.

There are two recommended alternative methods that provide this kind of
quantification. The preferred one is a visual test using confidence intervals.
The other one is a statistical test. The visual test is as follows. We
construct two confidence intervals for the means of the two systems. We check
if they overlap. If they do not, the decision is that “it is likely that the
systems do not have the same performance”. If they do overlap, the method
recommended by lilja and georges07 finishes, concluding that “the systems may
have the same performance”. The method by jain adds another step — it falls
back to a statistical test if the intervals overlap only slightly, that is if
centre of neither interval lies within the other. The visual test is
preferred, because in addition to the binary decision, it is easy to give a
visual measure of how large the difference in the means actually is, compared
to their uncertainty. While this is a clear aid for an analyst, it lacks a
rigorous semantics (and in his dissertation, georgesphd recommends using the
statistical test instead).

The confidence interval used for the visual test has the same underlying
statistical backing as the statistical test. It is assumed that measurements
are independent identically distributed and follow the normal distribution
(the normality assumption is incorrectly omitted in many texts as we discuss
in greater detail in Section 3.2). jain and lilja suggest using Student’s
t-test for unequal variances (or the corresponding interval). georges07
additionally recommend (an interval based on) the z-test for sample sizes over
30. georges07 also suggest working with sample means of measurements from
individual executions, instead of with the measurements themselves. The
corresponding test/interval is hence only for these sample means. An
alternative solution that is sometimes used, but not part of their
recommendation, is to join all the measurements from the different executions
and treat them as coming from a one-level experiment.

lilja and particularly jain then describe analysis of variance (ANOVA) methods
to summarize experiments where multiple fixed factors are varied (i.e.
processor, memory size, etc.). georges07 summarizes also some of the ANOVA
methods from [Lilja (2000)], but does not recommend it for use in (Java)
performance evaluations as the outputs are too hard to understand and doing
all the measurements to provide the inputs is too time consuming. We do not
address these types of factors (‘fixed effects’) in this work, but our method
should be extensible to do so.

The summarisation method that georges08 recommends for the replay compilation
design leads to a test/interval based on paired t-test. The inputs are
differences of means obtained from the old and new systems, averaging over
multiple executions with a given compilation plan. We do not address paired
comparisons in this work, but most of our criticism of the recommended
quantification using the test/interval we mention later applies to it as well.

##  3 What Is Wrong With The Currently Proposed Methods

Even though the current recommended best practice for quantification of
performance change by [Jain (1991)], [Lilja (2000)], and [Georges et al.
(2007)] is ahead of the pitiful practice in the field, it has a number of
problems.

###  3.1 Experiment Design

Repetition of executions (and measurements) is not enough. It is necessary to
include all inherently random factors in the experiment, and hence to repeat
experiments at level(s) even higher than execution of the benchmark. For
example, randomized compilation necessitates both compilation and running the
experiments for different binaries. We have found this to be essential in our
earlier work with the GNU C++ Compiler when benchmarking CORBA middleware and
when benchmarking applications within the Mono platform [Kalibera and Tuma
(2006)]. The empirical evaluation in this paper confirms this need on the Mono
platform. Furthermore, we need to randomize and include additional (non-
random) factors that it does not make sense to fix. These factors can arise
from the operating system or the language environment. envsize provide
evidence that unexpected factors such as the size of UNIX process environment
or the linking order can impact performance. verbrugge observed a performance
impact of the names of identifiers in the source code. It is obvious that many
additional factors of this kind have not yet been discovered to have
significant performance impact, so they cannot be really fixed in experiments.
Sufficient randomization, called for by envsize, is hence necessary.
stabilizer provide a randomization tool applicable to applications studied in
[Mytkowicz et al. (2009)].

Note in particular that randomization of identifier names or linking order
leads in turn again to randomized compilation/building. Repeating compilation
or repeating at an even higher level can be very expensive, because large
systems take long to compile. Repeating too many times is a waste of
resources, but repeating too few times undermines repeatability and validity.
Hence, a good experiment design method should also allow the researcher to
derive an optimum number of repetitions at different levels, so as to make
efficient use of the time available for experimentation in order to reduce
uncertainty in the final result. The present recommended best practice does
not offer this.

###  3.2 Summarization of Experiment Results

The single most significant problem of the best recommended method is that _it
does not tell us what we want to know_. It does not give us the metric we are
ultimately interested in, a reliable estimate of the ratio of execution times.
Of course, we can use the recommended method and, in the case that the
(binary) conclusion is that the systems are likely to differ in performance,
we can in addition report the ratio of the sample means. But the recommended
method would not give us the uncertainty of this ratio, so we do not know how
much of the ratio is due to uncertainty. Since we do not report uncertainty of
the ratio, the results we provide are unrepeatable and may not be valid. This
problem also means that the recommended method does not allow us to compare
the ratio against a threshold (the minimal/maximal change that is important to
us).

The problem of not supporting comparison against a threshold is more
significant than it may seem at first — the recommended method looks for a
performance difference however small as long as it is unlikely to be by pure
chance. However, the larger the sample size is (the more measurements we
have), the more unlikely even a very small difference becomes. Hence, the
decision of the method is influenced by the number of measurements. In
practice this means that with a large sample size (and in our field it is easy
to generate very large samples), the decision will nearly always be “it is
likely that the systems do not have the same performance”, no matter how small
or large the difference actually is. The method then becomes of very little
use — it just adds an illusion of rigour to the results and, worse, only to
the results we are not interested in.

The fact that this does not tell us what we want, and that it often gives the
same answer only because the sample size is very large, is a fundamental
drawback to the statistical significance. This drawback has been known for the
last seven decades [Cohen (1994)]. It has been reiterated by researchers in
psychology [Cohen (1994)], education [Coe (2002)], and more recently in bio-
medicine and biology [Nakagawa and Cuthill (2007)]. This and other drawbacks
of the method have been brought up in criticism of statistical significance
also in sociology, criminology, economics, marketing, chemistry, and nursing
(see fidlerps for references and details). The drawbacks, however, seem to be
quite unknown to experimental computer scientists. Even worse, methods based
on statistical significance are notoriously hard to interpret. This may be
partially because they do not give us the answer to what we want to know, but
they also offer temptations, such as the belief that the p-value is actually
the probability that the systems have the same performance. The interpretation
of the results of statistical tests is so difficult that even statistical
textbooks sometimes get it wrong (examples are given by cohen). In a study by
oakes cited by cohen, 68 of 70 psychologists made an error in the
interpretation of a statistical test. Interpretation of the visual test using
confidence intervals seems simpler, but it is not much so. This is well
reflected by the wording used and stressed by lilja, and repeated by georges07
and georges08, to describe the positive outcome (the systems differ in
performance) of the visual test — “there is no evidence to suggest that there
is not a statistically significant difference”.

The wording is rather cryptic as it tries to be precise, which we do not find
it to be. In frequentist statistics which the texts rely on, the wording is
incorrect. The goal of the comparison is to learn something about the
difference in the _true_ but unknown means of two random variables. The true
means are either equal or different, and the statistical significance does
_not_ speak about this difference (it speaks about the difference of the
_sample_ means instead). Hence, the second part of the wording, “statistically
significant difference” is incorrect. The first part of the wording that
speaks about non-existence of evidence is incorrect as well. The test is about
verifying if the data we have measured are unlikely provided the true unknown
means are the same. Hence, we are checking if the data we have form an
evidence _against_ the assumed zero performance change. We can make no claims
about non-existence of evidence. And we are not looking for evidence _for_ the
true means being equal, but rather _against_ it.

lilja (also cited by georges07) explain the wording by stating that there is
always a certain probability, ‘$\alpha$’, that a large observed difference was
due to random fluctuations. This is correct, but does not justify the wording
in our opinion. Even more, the context in the texts suggests that in the
visual test of confidence interval overlap, this probability ($\alpha$) is the
significance of the confidence intervals used (i.e. $\alpha=0.05$ for 95%
confidence interval), which it is not. In fact, with the visual test this
probability is not known. If we use 95% confidence intervals, the probability
of such an error is not 5%, but under the normality assumption is below 1%
[Payton et al. (2003)]. This makes the visual test far more conservative than
it may seem, and hence its results are even harder to interpret, and more
likely to mislead. Also, one would need far more experimentation time to show
a performance change. Further information on the error of the visual test can
be found in overlapp,schenker.

With the (non-visual) statistical test, the theory tells us the probability of
erroneously concluding that compared systems do not have the same performance
when actually they do. The probability of this error is the threshold we are
comparing the p-value against (often 5%). For this reason, schenker (cited by
georgesphd) prefer statistical test over the visual test. Our understanding is
that the visual test is generally preferred, though, in our field, and that
both method share their key problems.

Yet another issue of visual and non-visual test is the use of parametric
methods on data that violate their assumptions. Computer performance
measurements cannot be assumed to be normally distributed. Often they are
multi-modal, with long-tails to the right. Deviations from normality may not
be fatal for the t-test/confidence interval though, and some practitioners in
other fields ignore them as well. The sample mean is asymptotically normal due
to the Central Limit Theorem. Many texts in our field omit the normality
assumption of the t-test or incorrectly state that the Central Limit Theorem
is enough to overcome it ([Lilja (2000), Jain (1991), Georges et al. (2007)]).
The Central Limit Theorem is _not_ enough, but some more involved studies show
that, under certain deviations from normality, the parametric methods work
well, even for reasonably small sample sizes (see for example cfirobust or
robustsim) given a reasonably high confidence (95% confidence intervals and
wider). Still, there is no general agreement that ignoring the normality
assumption is an acceptable practice, one should certainly make it clear if
such assumptions are ignored. Robust statistical methods do exist (see
robustas for summary and references), and there has been no study of how the
t-test/confidence interval is affected by violations from normality common in
computer performance data.

##  4 How to do things better

A good quantification method needs to provide an estimate of the ratio of
execution time means and an estimate of its uncertainty. These estimates
should be based on experiments that include all random factors to which the
measured system is subject, both factors naturally random and factors that we
randomize to avoid measurement bias. A statistical model of such an experiment
can then offer some way of planning the experiment: how many repetitions are
needed at each level.

###  4.1 Experiment Design with Random Effects

There are factors that influence performance, which are inherently random and
we cannot control them. To get valid results (avoid bias), we need to repeat
experiments at a high enough level to include all random factors. For example,
if our system is prone to randomized compilation which has an effect on
performance, we need to repeat compilations and measure multiple binaries.
With managed runtimes, the ‘binary’ may be the binary of the virtual machine
or the byte-code of the application. If we took an extreme position, it would
be all compiled code that influences performance observed by the benchmark. In
theory, we could start from the top level for each measurement, thus compiling
the binary, executing it only once and running only one iteration (say after
dropping the initial measurements that were prone to warm-up noise). This
approach would allow use of a one-level model. However, it would be an
extravagant waste of time if the variation in binaries had a far smaller
effect on performance than the variation in, say, executions or measurements.
In this case, intuitively it should be possible to do better by executing each
binary multiple times and reporting more than one measurement in an execution.

This situation can be modeled mathematically by a random effect models with
$n$-way classification. Such random effects model is based on a hierarchy of
$n$ random ways,222In statistical texts, one can say also ‘$n$ factors’, but
we refrain from doing so here to avoid confusion with the informal meaning of
a ‘factor’ used throughout the text. In our text, a ‘factor’ refers to a
technical/real cause that may impact performance and a ‘way’ refers to a
statistical model of (some of) such causes. which have effect on the random
distribution of the actual measurements. Hence, a model with $n$-way
classification corresponds to an $n+1$ level experiment. In particular, a
three-level experiment with repeated compilations, executions, and
measurements corresponds to a model in 2-way classification (measurements
influenced by execution and compilation). Random effects models are sometimes
also referred to as random effects ($n$-factor) ANOVA. References to
specialized literature on related models can be found in mcculloch, though we
have not found a model that would apply directly. Derivation of such a model
is one of the contributions of this work.

###  4.2 Quantification with Effect Size Confidence Interval

We propose to construct a confidence interval for the ratio of mean execution
times. This confidence interval will be a measure of uncertainty of the metric
that is of ultimate interest. Reported results can thus be for instance that
system A is 4%$\pm$1.5% faster than system B, with 95% confidence. Multiple
intervals would also lend themselves easily to graphical visualisation. The
method gives us a way to rigorously compare against a threshold. Say that we
only care about differences larger than 3%. Here, we would conclude that
performance of two systems is different if the upper bound of the confidence
interval is less than 0.97 or its lower bound greater than 1.03,

Note also that the method still allows a significance-based binary decision.
We conclude that the systems “may have the same performance” if the interval
includes 1. Otherwise, we would conclude that the systems are “likely not to
have the same performance”. The chance of erroneously concluding that the
systems differ in performance here is the confidence of the interval (5% for a
95% confidence interval). Hence, even for the significance based method only,
we have the advantage over the recommended practice that we know this error.
In summary, the confidence interval for the ratio of means subsumes the
currently best recommended practice.

Constructing a confidence interval for a metric that measures the difference
in two systems is a known concept in statistics. Such a metric is called the
_effect size_ and hence the _effect size confidence interval_. Effect size
confidence intervals have been proposed for quantification as a replacement
for significance methods in psychology [Cohen (1994)], education [Coe (2002)],
medicine, bio-medicine, and biology [Nakagawa and Cuthill (2007)]. The change
in statistical methods is not smooth, but evidence can be found that effect
size is already being used, at least to some extent, in medicine [Bland and
Altman (2000)] and psychology [Dunleavy et al. (2006)]. To support the change,
the use of effect size has sometimes been made an official requirement (eshill
list 23 journals mostly in psychology and education that require reporting
effect size, and the apa requires it in its Publication Manual).

There are different metrics used for measure the effect size in different
fields and one of them is the ratio of means. There are statistical methods
for interval construction, although they did not make it to introductory
textbooks (more information and references can be found in [Beyene and
Moineddin (2005), Schaarschmidt (2007), Dilba et al. (2007)]). From several
available methods, we have a closer look at two. One is based on statistical
simulation (bootstrap) [Davison and Hinkley (1997)]. The other is based on
Fieller’s theorem [Fieller (1954)], which gives a confidence interval for the
ratio of means of normally distributed variables. We show how both methods can
be applied in the experimental design we propose, and empirically evaluate
them.

###  4.3 Related Efforts

Some of the things we propose here have been proposed earlier in computer
science, but have not been widely adopted.

A bootstrap-based method for construction of confidence intervals is being
used in the Haskell community, supported by the Criterion benchmarking library
[O’Sullivan (2009)]. The library is for measuring performance of one system
only in single-level experiments. Apart from using robust methods it has also
the advantage of detecting outliers and auto-corellation of the data. The tool
is based on [Boyer (2008)].

Confidence interval for speed-up has been proposed by luo04 in the context of
processor simulation. The uncertainty does not come from random effects, but
rather from random sub-sampling – only randomly chosen execution intervals
from the whole application execution are simulated, the same intervals in both
of the systems to compare. The confidence interval proposed for this problem
comes from [Cochran (1977)], it is a parametric interval based on asymptotic
normality of not only the means, but also their ratio, which is only possible
for large samples. The method by fieller that we use in this work is more
general.

##  5 Outline of the New Method

The rest of our paper includes a detailed description of our quantification
method and its empirical evaluation based on real benchmark results. We offer
two alternative descriptions of our method. Section 6 provides a guide for
practitioners, while Section 7 is a variant for scientists with some
statistical background. Section 6 describes the method in practical terms and
with recommendations on how it can be used. Section 7 formulates the method in
statistical terms, gives the assumptions and necessary proofs, and discusses
further alternatives. We use the terminology common in the field.

In statistical terms, the core of the method is our statistical model of
random effects in $n$-way classification, which models performance in one
system. We apply it independently to both systems we aim to compare. For
clarity, the model is described for 2-way classification first (3 levels in
the experiment, Section 7.1.1), but later in the general form of $n$-way
classification (an arbitrary number of levels, Section 7.1.2). While we keep
referring to ‘execution time’ in the statistics to make presentation concrete,
the model is sufficiently general for any response variable, given the stated
assumptions. We then show how to construct the confidence interval for the
mean within this model, still for the mean of each system separately (Section
7.2). We give two alternative methods, a parametric method based on asymptotic
normality and a non-parametric method based on bootstrap. Later we show how to
construct the confidence interval for the ratio of the means of two systems,
again using a parametric and non-parametric method (Section 7.3). The non-
parametric bootstrap method is a natural extension of the method for a single
system. The parametric method uses Fieller’s theorem [Fieller (1954)].

##  6 For Practitioners: The New Method without Statistics

The technique we propose here comprises of the actual statistical method
(proven and precisely defined in the following section) and of practical
advice how to use it. We feel that it is important to stress that while the
statistical method is general, our advice may not be. The practicalities
depend on the context (the goal of the study, the criticality of the results
and the consequences of possible errors, the kind of computer system
evaluated, the runtime environment, infrastructure, and benchmarks). We
understand that having a single best technique that could be mechanically
followed in every study might be appealing, but we strongly believe that such
practice is an illusion, and we definitely do not provide one here.

###  6.1 Designing the Experiments

We describe our technique sequentially, but in practice it would be an
iterative process. At the end of this section, we provide a small worked
example.

##### Initial number of levels

When planning experiments, we first decide what the highest level of the
experiment should be. All factors influencing performance above that level
need to be controlled (and documented). The highest level might be compilation
(of the benchmark, runtime environment, or all). If compilation is found to be
fully deterministic, this may not be needed. Though note that even symbol
names have impact on performance [Gu et al. (2004)], so we may choose to
randomize these or otherwise randomize the building process to get rid of
other sources of measurement bias.

On the other hand, in practice, if repeating say compilation is very
expensive, and we know that it has relatively little impact on performance or
we find it out during the following process, we may choose not to repeat it.
The threshold we then use for comparisons, however, needs to be set
accordingly.

After deciding on the highest level, we continue by identifying lower levels.
We need to add an additional lower level to the experiment if the variability
(say in executions) is expected or known to be higher than in the level above
(say compilations), while the cost of repeating at the lower level is smaller
than at the higher level. If we create a level where it is in fact not needed,
we may detect this later in the following process. Also, we can measure the
costs of repetition in the following process as a side-effect (for example,
such cost is the time to build a new binary). In the formulae below, the
number of levels is $n+1$ (where $n$ is the number of ‘ways’ in the
statistical model, which is one less than the number of levels in the
experiment).

##### Numbers of repetitions

Having decided on initial levels of the experiment, we run the experiment a
few times to estimate variability at different levels. With our method, we
always use the same repetition counts at a particular level (thus all
executions have the same number of measurements, all binaries same number of
executions, etc.). We denote the repetition counts as $n_{i}$, where $1\leq
i\leq n+1$ is the level of the experiment ($n_{n+1}$ is the number at highest
level, say compilation).

For the initial experiment, we can set all these counts to the same value, say
30. Note that these counts do not include warming-up various components of the
experiment. Usually, we want to evaluate only measurements from a steady-state
of every execution. Say that we do so by discarding the first $c_{1}$
measurements of each execution before collecting the $n_{1}$ we need. $c_{1}$
must be found by the experimenter. Repetitions at higher levels, $c_{i}$,
$1<i\leq n$, also may have a non-zero cost. For instance, compilation has
notable cost as it takes time to compile. All costs are represented as counts
of measurements — the cost of compilation is the compilation time divided by
the (average) time for a measurement. Precision here is not crucial, as the
costs are only needed for optimisation.

To run the initial experiment, we need to have $c_{1}$ already and we use
$n_{i}=30$. During the experiment we collect measurements and also measure the
costs for repetition $c_{i}$ for $1<i\leq n$. Based on the costs and the
measurements, we decide on optimum counts of repetitions at all but the
highest level (the derivation is provided in Sections 7.4 and 7.5) :

$n_{1}=\left\lceil\,\sqrt{c_{1}\frac{T_{1}^{2}}{T_{2}^{2}}}\,\right\rceil,\qquad\forall
i,1<i\leq
n\,\,n_{i}=\left\lceil\,\sqrt{\frac{c_{i}}{c_{i-1}}\frac{T_{i}^{2}}{T_{i+1}^{2}}}\,\right\rceil.$

We obtain $T_{i}$, the unbiased estimator of the variance at level $i$,
through an iterative process. First, we calculate $S_{i}$, the biased
estimator of the variance at level $i$ (formulae below). Then, we start
calculating $T_{i}$ as (the derivation is provided in Section 7.6)

$\displaystyle T_{1}^{2}$ $\displaystyle=$ $\displaystyle S_{1}^{2},$
$\displaystyle\forall i,1<i\leq n+1\,\,T_{i}^{2}$ $\displaystyle=$
$\displaystyle S_{i}^{2}-\frac{S_{i-1}^{2}}{n_{i-1}}.$

If we should get $T_{i}^{2}\leq 0$ (or at least very small), then this level
of the experiment induces little variation so we can remove level $i$ from the
experiment. This is semantically equivalent of running the experiment all
again with fewer levels, but we obtain the same effect by dropping data from
the repetitions. We calculate $S_{i}$s as follows (the used notation is
detailed in Sections 7.1.2 and 7.1.3)

$S_{1}^{2}=\frac{1}{\prod_{k=2}^{n+1}n_{k}}\,\frac{1}{n_{1}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\dots\sum_{j_{1}=1}^{n_{1}}\left(Y_{j_{n+1}\dots
j_{1}}-\overline{Y}_{j_{n+1}\dots j_{2}\bullet}\right)^{2}$

for $i$, $2\leq i\leq n$

$S_{i}^{2}=\frac{1}{\prod_{k=i+1}^{n+1}n_{k}}\,\frac{1}{n_{i}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\dots\sum_{j_{i}=1}^{n_{i}}\left(\overline{Y}_{j_{n+1}\dots
j_{i}\underbrace{\bullet\dots\bullet}_{i-1}}-\overline{Y}_{j_{n+1}\dots
j_{i+1}\underbrace{\bullet\dots\bullet}_{i}}\right)^{2}$

and finally

$S_{n+1}^{2}=\frac{1}{n_{n+1}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\left(\overline{Y}_{j_{n+1}\underbrace{\bullet\dots\bullet}_{n}}-\overline{Y}_{\underbrace{\bullet\dots\bullet}_{n+1}}\right)^{2}$

$Y$ denotes the measurements, indexed by experiment levels (highest to
lowest). The bar over $Y$ denotes an arithmetic mean. The mean is always
calculated over all the indexes that are denoted by a bullet. While the
formulas for $S_{i}$s may seem involved in full generality, it is easy to see
how they are constructed. Suppose we have functions $M$ and $V$ that calculate
sample mean and variance of a vector,

$M(x_{1}\dots x_{k})=\frac{1}{k}\sum_{i=1}^{k}x_{i}\qquad V(x_{1}\dots
x_{k})=\frac{1}{k-1}\sum_{i=1}^{k}\left(x_{i}-M(x_{1}\dots x_{k})\right)^{2}$

Any statistical software would have these functions. For $n=2$, each of
$S^{2}_{1}$, $S^{2}_{2}$, and $S^{2}_{3}$ is created by three applications of
these functions to the data, where $V$ is applied once and $M$ two times, but
always in different order (see Table 6.1). Each application of $M$ or $V$
reduces the dimension of the data. For example, to calculate $S_{2}^{2}$ we
first apply $M$ on measurements from each execution, hence getting a two-
dimensional matrix (binaries $\times$ execution means) from a three-
dimensional matrix (binaries $\times$ executions $\times$ measurements).
Second, we apply $V$ on the execution means, getting a vector (one element per
binary) of variances of the execution means. Third, we apply $M$ on this
vector and get $S_{2}^{2}$.

Interpretation of $S_{i}^{2}$ in Three-level Experiment  Binaries Executions
Measurements Interpretation $S_{1}^{2}$ $M$ $M$ $\mathbf{V}$ Mean variance at
level 1 (due to measurements) $S_{2}^{2}$ $M$ $\mathbf{V}$ $M$ Mean variance
at level 2 (of execution means) $S_{3}^{2}$ $\mathbf{V}$ $M$ $M$ Variance at
level 3 (of binary means)

If we obtain a small number of repetitions at some level, we may decide to
remove that level to simplify the experiment. Indeed for practical purposes we
can also increase the number of repetitions at the lowest level (i.e.
measurements), if a measurement is very short. As there is no further overhead
for additional measurements at the lowest level (since we do not have to
repeat higher levels and nor perform further warm-up), it is cheap to get more
data useful for other things, such as spotting unusual behaviour. This method
would not work well for very small repetition numbers (say, below 5, but this
depends on experimenter’s judgment) — if at all possible, we should not use
such small counts. Note that in practice, to make the process simpler at the
cost of more expensive experiments in the end, we could assume $T_{i}=S_{i}$.
In that case we would not get problems with the negative numbers, but the
resulting numbers of repetitions at higher levels may be higher than needed.

Once we have decided on the final number of levels and repetition counts, we
re-run the experiments and continue as follows. In the following, $n_{i}$ and
$n$ refer to the new setting. In practice, the established setting would be
re-used for very similar experiments (same platform, benchmarks, etc).

##### Example

Let us assume that we are to dimension an experiment based on three-level data
($n=2$) shown in Table 6.1 (Raw Data). We have measurements from 3 binaries,
each executed 2 times, and reporting 2 measurements per each execution. In a
real experiment we will have more, but this is to demonstrate the method. To
calculate the $T^{2}_{i}$ we need the $S^{2}_{i}$, which we will calculate
iteratively as demonstrated in the table. First, we calculate the (matrices
of) execution means and variances ($\bullet\bullet M$ and $\bullet\bullet V$),
as recorded in the second and third table. Next, we calculate the column means
and variances of $\bullet\bullet M$, getting vectors $\bullet MM$ and $\bullet
VM$. Also we calculate $\bullet MV$. Finally, we get $S_{3}^{2}\doteq 3.6$
($VMM$), $S_{2}^{2}\doteq 2.6$ ($MVM$), and $S_{1}^{2}\doteq 16.5$ ($MMV$).
The grand mean, the mean of all measurements, is 6.5 ($MMM$). We now obtain
the $T_{i}$: $T_{1}^{2}=S_{1}^{2}\doteq 16.5$,
$T_{2}^{2}=S^{2}_{2}-S^{2}_{1}/2\doteq-5.7$, and $T_{3}^{2}\doteq 2.3$. As we
find that $T^{2}_{2}<0$, we can remove the second level of the experiment.
Hence, we only will run one execution per binary.

Example Data (Three-level)

| Raw Data
---
$\bullet\bullet\bullet$ | Binaries
Executions | 9 | 5 | 10 | 6 | 1 | 12
8 | 3 | 7 | 11 | 2 | 4
| Execution Means
---
$\bullet\bullet M$ | Binaries
Exec. | 7. | 0 | 8. | 0 | 6. | 5
5. | 5 | 9. | 0 | 3. | 0
| Execution Variances
---
$\bullet\bullet V$ | Binaries
Exec. | 8. | 0 | 8. | 0 | 60. | 5
12. | 5 | 8. | 0 | 2. | 0
| Binary Means
---
$\bullet MM$ | Binaries
| 6.3 | 8.5 | 4.8
| Variances of Exec. Means
---
$\bullet VM$ | Binaries
| 1.1 | 0.5 | 6.1
| Means of Exec. Variances
---
$\bullet MV$ | Binaries
| 10.3 | 8.0 | 31.3
| Variance of Binary Means
---
$VMM$ $S_{3}^{2}$
3.6
| Mean Variance of Exec. Means
---
$MVM$ $S_{2}^{2}$
2.6
| Mean of Measur. Variances
---
$MMV$ $S_{1}^{2}$
16.5
| Grand Mean
---
$MMM$
6.5

Example Data (Data from Table 6.1 Transformed to Two-level Data)

Raw Data
---
$\bullet\bullet$ | Binaries
Measur. | 9 | 10 | 1
8 | 7 | 2
| 5 | 6 | 12
| 3 | 11 | 4

Sample Means and Variances
---
$\bullet M$ | Binary Means | 6. | 3 | 8. | 5 | 4. | 8
$\bullet V$ | Binary Variances | 7. | 6 | 5. | 7 | 24. | 9
$VM$, $S^{2}_{2}$ | Variances of Binary Means | 3. | 6 |
$MV$, $S^{2}_{1}$ | Mean of Measurement Variances | 12. | 7 |
$MM$ | Grand Mean | 6. | 5 |

Let us also suppose that we observed that the system stabilises running a
single execution reliably after 10 measurements ($c_{1}=10$, suppose that
Table 6.1 already contains only the stable measurements). We need to find the
optimum repetition counts. As we decided to remove the second level of the
experiment, this optimum repetition count is just the number of measurements
per execution (the number of executions per binary will always be 1, and the
number of binaries we can increase any time to get more precise results). To
find the optimum we need to recalculate the $T^{2}_{i}$ and hence the
$S^{2}_{i}$ for two-level experiment ($n=1$). Applying the same method as
before, we get $S^{2}_{2}\doteq 3.6$ and $S^{2}_{1}\doteq 12.7$ (Table 6.1),
hence $T^{2}_{1}\doteq 12.7$ and $T^{2}_{2}\doteq 0.4$. And thus the optimum
number of measurements per binary, $n_{1}$, is

$n_{1}=\left\lceil\sqrt{c_{1}\frac{T^{2}_{1}}{T^{2}_{2}}}\right\rceil\doteq\left\lceil\sqrt{10\frac{12.7}{0.4}}\right\rceil=18$

We would thus run two-level experiment with $n_{1}=18$ measurements per
binary, executing each binary only once.

###  6.2 Summarizing the Results

##### Quantifying Performance of One System

We report performance of one system (in isolation) as the arithmetic mean of
all measurements $\overline{Y}$ (we omit the bullets in the notation for
simplicity). We estimate its uncertainty using a $(1-\alpha)$ confidence
interval (i.e. $\alpha=0.05$ gives a 95% confidence interval):

$\overline{Y}\pm
t_{1-\frac{\alpha}{2},\nu}\sqrt{\frac{S_{n+1}^{2}}{n_{n+1}}}=\overline{Y}\pm
t_{1-\frac{\alpha}{2},\nu}\sqrt{\frac{1}{n_{n+1}(n_{n+1}-1)}\sum_{j_{n+1}=1}^{n_{n+1}}\left(\overline{Y}_{j_{n+1}\underbrace{\bullet\dots\bullet}_{n}}-\overline{Y}_{\underbrace{\bullet\dots\bullet}_{n+1}}\right)^{2}}$

where $t_{1-\frac{\alpha}{2},\nu}$ is the $1-\frac{\alpha}{2}$-quantile of the
$t$-distribution with $\nu=n_{n+1}-1$ degrees of freedom. The mathematical
background is given in Section 7.2.2.

Alternatively, we can use statistical simulation (bootstrap) to calculate the
confidence interval. Say that we perform 1000 steps (or more if there is
time). Within each step, we use the real data to simulate a new experiment.
First, at the highest level, we randomly decide which iterations of the real
experiment to use (e.g. which binaries). We generate the same number of
iterations as the real experiment, but some of the real iterations can be used
multiple times while some not be used at all. We then apply this principle to
lower levels. In the end, we get the same number of measurements as in a real
experiment, with the same structure (repetition counts at each level and
number of levels), and we calculate a sample mean of all these measurements.
Hence, we get 1000 means, each originating from one step of the simulation.

We form a $(1-\alpha)$ confidence interval for the mean using $\alpha/2$ and
$1-\alpha/2$ sample quantiles of these means. If we have 1000 steps, we can do
this by ordering the simulated means and taking the 25th and 975th values.
Pseudo-code that illustrates this procedure is shown in Figure 1 on page 1 and
more details are given in Section 7.2.1.

##### Quantifying Performance Change

To quantify a performance change of a ‘new’ system over an ‘old’ system, we
use the same numbers of repetitions and levels for both. We report the ratio
of mean execution times of the two systems. As a measure of uncertainty, we
report a $(1-\alpha)$ confidence interval (again $\alpha=0.05$ gives a 95%
confidence interval) as follows. For detailed notation and derivation, see
Sections 7.3 and 7.3.2.

$\frac{\overline{{}^{O}Y}\cdot\overline{{}^{N}Y}\mp\sqrt{\left(\overline{{}^{O}Y}\cdot\overline{{}^{N}Y}\right)^{2}-\left(\left(\overline{{}^{O}Y}\right)^{2}-t_{\frac{\alpha}{2},\nu}^{2}\cdot\frac{{{}^{O}S}^{2}_{n+1}}{n_{n+1}}\right)\left(\left(\overline{{}^{N}Y}\right)^{2}-t_{\frac{\alpha}{2},\nu}^{2}\cdot\frac{{{}^{N}S}^{2}_{n+1}}{n_{n+1}}\right)}}{\overline{{}^{O}Y}^{2}-t_{\frac{\alpha}{2},\nu}^{2}\cdot{{}^{O}S}^{2}_{n+1}\cdot{n_{n+1}}^{-1}}$

The term $t_{\frac{\alpha}{2},\nu}$ denotes the $\frac{\alpha}{2}$-quantile of
the $t$-distribution with $\nu=n_{n+1}-1$ degrees of freedom. The left-super-
scripts ‘N’ and ‘O’ denote the new and the old system that we compare.

Alternatively, we can use statistical simulation (bootstrap) to calculate the
confidence interval. The algorithm is similar to the case of one system. We
again perform a number of steps (say 1000), each producing the metric of
interest from a simulated experiment, which now is the ratio of means. Hence,
in each step, we simulate measurements from both systems, calculate their
sample means, and then their ratio. When we have these 1000 ratios, we take
the sample quantiles for the confidence interval (i.e. 25th and 975th for a
95% interval). Pseudo-code that illustrates this procedure is shown in Figure
2 on page 2 (Section 7.3.1).

##### Example

Let us assume we have measurements from an old and a new system as shown in
Table 6.2. In practice, we would have more repetitions, but let us use a
simple example to demonstrate the method. Let us first show how to calculate
confidence interval for the mean of one system. We will demonstrate this on
the old system. We need to calculate the variance at highest level,
${{}^{O}S}^{2}_{3}$, but not the other variances, so there is less work to do
than when dimensioning the experiment. Variances at the highest level are
simply variances of means (of binaries, in our case). For the old system, the
mean for the first binary is 7.8, for the second 12.3 and for the third 11.5
(all numbers rounded to one decimal place). The grand mean is hence 10.5 and
variance of the means, ${{}^{O}S}^{2}_{3}$, is 5.8. The 95% confidence
interval for the mean of the old system is hence

$\overline{Y}\pm
t_{1-\frac{\alpha}{2},\nu}\sqrt{\frac{S_{n+1}^{2}}{n_{n+1}}}\doteq 10.5\pm
4.3\sqrt{\frac{5.8}{3}}\doteq 10.5\pm 6.0$

where 4.3 is a (rounded value of) the $0.975$ quantile of the $t$-distribution
with $2$ degrees of freedom. The confidence interval for the mean of the new
system would be obtained in the same way.

Example Raw Data (Three-level)

| Raw Data - Old System
---
$\bullet\bullet\bullet$ | Binaries
Executions | 9 | 11 | 16 | 13 | 15 | 7
5 | 6 | 12 | 8 | 10 | 14
| Raw Data - New System
---
$\bullet\bullet\bullet$ | Binaries
Executions | 10 | 12 | 9 | 1 | 8 | 5
6 | 7 | 11 | 4 | 3 | 2

Alternatively, we can use the bootstrap method to calculate the confidence
interval. Table 6.2 shows three bootstrap replicates of our measurements. In
practice, we would need much more. The first replicate uses binaries 1, 2, and
1 from the original old system (Table 6.2). For example, the third selected
binary, which corresponds to binary 1 of the old system of Table 6.2, takes
executions 1 and 1 from the original binary. For the first execution, it takes
measurements 1 and 1 from the original execution. For the second execution, it
takes measurements 2 and 1. The grand means of these three replicates are,
after rounding, $9.0$, $11.4$, and $10.2$. When sorted, we have $9.0$, $10.2$,
and $11.4$. Three replicates is certainly not enough, but if we had, say,
1000, we would take the 25th and 975th from such sorted sequence as the lower
and upper bound of a 95% confidence interval for the mean of the old system.

Bootstrap Replicates of Raw Data (Old System from Table 6.2)

| $\bullet\bullet\bullet$ | Binaries
---|---
Executions | 9 | 11 | 12 | 8 | 9 | 9
5 | 5 | 8 | 12 | 11 | 9
| Binaries
---
14 | 10 | 13 | 13 | 8 | 8
15 | 7 | 8 | 12 | 16 | 13
| Binaries
---
7 | 7 | 11 | 11 | 9 | 9
15 | 15 | 9 | 11 | 9 | 9

We are also interested in the 95% confidence interval for the ratio of the
mean execution times of the new and old system from Table 6.2, we proceed as
follows. We already know that ${{}^{O}S}^{2}_{3}\doteq 5.8$. Using the same
algorithm we calculate ${{}^{N}S}^{2}_{3}\doteq 4.6$. The grand means are
$\overline{{}^{O}Y}=10.5$ and $\overline{{}^{N}Y}=6.5$. From Student’s
$t$-distribution, $t^{2}_{\frac{\alpha}{2},\nu}=t^{2}_{\frac{0.05}{2},2}\doteq
18.5$. The confidence interval for the ratio of means is therefore

$\frac{10.5\cdot 6.5\mp\sqrt{(10.5\cdot
6.5)^{2}-\left(10.5^{2}-18.5\frac{5.8}{3}\right)\left(6.5^{2}-18.5\frac{4.6}{3}\right)}}{10.5^{2}-18.5\frac{5.8}{3}}\doteq\frac{68.3\mp
60.2}{74.5}.$

The confidence limits are thus $0.1$ and $1.7$ (90% performance improvement to
70% performance degradation). Such a wide interval would not be useful in
practice, but this example has used only very small repetition counts.

Alternatively, we could calculate the confidence interval for the ratio of
means of the two systems using the bootstrap method. The bootstrap replicates
would be created for both systems in the same way as we demonstrated it for
the old system. We would then calculate the ratio of means of these replicates
(the first of the new system over the first of the old system, the second of
the new over the second of the old, etc.). From these ratios, we would select
the respective quantiles.

##  7 For Scientists: The New Method and The Statistics Behind It

In this section, we formulate the method in statistical terms, give its
assumptions and provide proofs. We also discuss alternatives. In the
description of the statistics behind our method, we assume the execution time
of the operation of interest is a continuous random variable, $Y$. The range
of $Y$ is a subrange of real numbers:
$\exists\mathit{BCET},\mathit{WCET}\in\mathbb{R}:P(Y<\mathit{BCET})=P(Y>\mathit{WCET})=0$,
where $\mathit{BCET}$ is the best-case execution time and $\mathit{WCET}$ is
the worst-case execution time. We assume that the expectation of $Y$ exists
and that its variance is finite, $E(Y)=\mu$ and
$\mathop{var}\nolimits(Y)=\sigma^{2}$. The ‘operation of interest’ can be
anything that is the goal of the measurement, small or large, depending on the
benchmark we use. $Y$ models only the steady state duration of the operation,
and we only focus on steady state performance here.

###  7.1 Statistical Model for a Hierarchy of Random Effects

####  7.1.1 Two-way Classification

We first describe the model with three levels333In our text, the term ‘level’
always refers to levels in a benchmark experiment — it should not be confused
with its statistical meaning in texts that address fixed effects models. We
only have random effects. of hierarchy. For clarity of exposition, we shall
use our running example of random effects in compilation, execution, and
measurement. Later we show a general description for an arbitrary number of
levels. The intuition behind the model is simple: the times measured in a
single _execution_ are randomly distributed, with a mean that is also a random
variable (this notion is later formalised using conditional expectation). In
turn, the mean of these execution means in a _binary_ 444We use the term
‘binary’ to denote a single binary executable, that is, a product of
compilation. For the statistical model, it is just a factor. is a random
variable. And, finally, the mean of these binary means is an unknown constant,
the grand mean $\mu$ for the whole system we are interested in. We start with
a formalisation of the model and its properties.

We assume that measurements from a single execution are independent
identically distributed. Their mean differs for different executions, as the
measurements are influenced by random effects related to a particular
execution. This can be expressed using conditional expectation as
$E(Y|[\mu_{E}=m])=m$, where $\mu_{E}$ is a random variable which gives a mean
$m$ (a number), a mean of measurements of one particular execution. The
expression says that the random variable $E(Y|\mu_{E})$ has the value of $m$
when random variable $\mu_{E}$ has the value of $m$. We use a shorter notation
for this, $E(Y|\mu_{E})=\mu_{E}$. We assume that the variance of measurements
within execution is a constant
$\sigma_{E}^{2}=\mathop{var}\nolimits(Y|\mu_{E})$ (homoscedastic). This says
that the random variable $\mathop{var}\nolimits(Y|\mu_{E})$ has the value of
$\sigma_{E}^{2}$ no matter what is the value of $\mu_{E}$. We then assume that
the $\mu_{E}$ are independent identically distributed within a given binary,
but the mean of this distribution is again a random variable within a system:
$E(\mu_{E}|\mu_{B})=\mu_{B}$. We assume that the variances of execution means
within a binary are constant:
$\mathop{var}\nolimits(\mu_{E}|\mu_{B})=\sigma_{B}^{2}$. The binary mean
$\mu_{B}$ is a random variable for a given system, $E(\mu_{B})=\mu_{S}$, where
$\mu_{S}$ is a constant. We denote the variance of $\mu_{B}$ as
$\mathop{var}\nolimits(\mu_{B})=\sigma_{S}^{2}$, where $\sigma_{S}^{2}$ is a
constant.

It can be shown that $E(Y)=\mu_{S}\;(=\mu)$ and that
$\mathop{var}\nolimits(Y)=\sigma_{E}^{2}+\sigma_{B}^{2}+\sigma_{S}^{2}\;(=\sigma^{2})$.
We will do so later for the general case of $n$-way classification. In
summary, we have the following random effects model in two-way classification
(the distribution of the observed execution time is randomly influenced in two
ways, through execution and compilation):

$\begin{array}[]{rclrcl}E(Y)&=&\mu,&\mathop{var}\nolimits(Y)&=&\sigma^{2}=\sigma_{E}^{2}+\sigma_{B}^{2}+\sigma_{S}^{2}\\\\[5.69054pt]
E(Y|\mu_{E})&=&\mu_{E},&\mathop{var}\nolimits(Y|\mu_{E})&=&\sigma_{E}^{2}\\\\[2.84526pt]
E(\mu_{E}|\mu_{B})&=&\mu_{B},&\mathop{var}\nolimits(\mu_{E}|\mu_{B})&=&\sigma_{B}^{2}\\\\[2.84526pt]
E(\mu_{B})&=&\mu,&\mathop{var}\nolimits(\mu_{B})&=&\sigma_{S}^{2}\\\
\end{array}$

####  7.1.2 N-way Classification

In $n$-way classification, the measurements within an execution ($Y|\mu_{1}$)
are independent identically distributed with mean $\mu_{1}$ and variance
$\sigma_{1}^{2}$. $\mu_{1}$ is a random variable. $\sigma_{1}^{2}$ is a
constant:

$E\left(Y|\mu_{1}\right)=\mu_{1},\quad\mathop{var}\nolimits\left(Y|\mu_{1}\right)=\sigma^{2}_{1}.$
(1)

When $n$, the number of ways of the classification, is two or more, the mean
of $\mu_{1}|\mu_{2}$ is $\mu_{2}$, again a random variable. In general,

$\forall i,1\leq i\leq n-1,\quad
E\left(\mu_{i}|\mu_{i+1}\right)=\mu_{i+1},\quad\mathop{var}\nolimits\left(\mu_{i}|\mu_{i+1}\right)=\sigma^{2}_{i+1}.$
(2)

Finally, $\mu_{n}|\mu_{n+1}$ is a random variable with mean $\mu_{n+1}$, which
is a constant:

$E\left(\mu_{n}\right)=\mu_{n+1},\quad\mathop{var}\nolimits\left(\mu_{n}\right)=\sigma^{2}_{n+1}.$
(3)

######  Lemma 7.1 (Rule of Iterated Expectations)

If $X$ and $Y$ are random variables and the expectations exist,
$E\left[E(Y|X)\right]=E(Y)$. aos, Theorem 3.24, p. 55.

We will now show that $\mu_{n+1}=E(Y)$ (note that $E(Y)=\mu$ by definition):

$\displaystyle\mu=E(Y)$ $\displaystyle=^{(L\ref{lroi})}$ $\displaystyle
E\left[E\left(Y|\mu_{1}\right)\right]=^{(\ref{eM1})}E\left(\mu_{1}\right)$
$\displaystyle\forall i,1\leq i\leq n-1,\quad E\left(\mu_{i}\right)$
$\displaystyle=^{(L\ref{lroi})}$ $\displaystyle
E\left[E\left(\mu_{i}|\mu_{i+1}\right)\right]=^{(\ref{eMi})}E\left(\mu_{i+1}\right)$
$\displaystyle E\left(\mu_{n}\right)$ $\displaystyle=^{(\ref{eMn})}$
$\displaystyle\mu_{n+1}$

######  Lemma 7.2 (Property of conditional variance)

For random variables $X$ and $Y$,
$\mathop{var}\nolimits(Y)=E\left[\mathop{var}\nolimits(Y|X)\right]+\mathop{var}\nolimits\left[E(Y|X)\right]$.
aos, Theorem 3.27, p. 55.

We will now show that
$\sigma^{2}=\mathop{var}\nolimits(Y)=\sum_{i=1}^{n+1}\sigma_{i}^{2}$ (note
that $\mathop{var}\nolimits(Y)=\sigma^{2}$ by definition):

$\displaystyle\sigma^{2}=\mathop{var}\nolimits(Y)$
$\displaystyle=^{(\ref{lpcv})}$ $\displaystyle
E\left(\mathop{var}\nolimits\left(Y|\mu_{1}\right)\right)+\mathop{var}\nolimits\left(E\left(Y|\mu_{1}\right)\right)=^{(\ref{eM1})}E\left(\sigma_{1}^{2}\right)+\mathop{var}\nolimits\left(\mu_{1}\right)=$
$\displaystyle=$
$\displaystyle\sigma_{1}^{2}+\mathop{var}\nolimits\left(\mu_{1}\right)$
$\displaystyle\forall i,1\leq i\leq
n-1,\quad\mathop{var}\nolimits\left(\mu_{i}\right)$
$\displaystyle=^{(\ref{lpcv})}$ $\displaystyle
E\left(\mathop{var}\nolimits\left(\mu_{i}|\mu_{i+1}\right)\right)+\mathop{var}\nolimits\left(E\left(\mu_{i}|\mu_{i+1}\right)\right)=^{(\ref{eMi})}$
$\displaystyle=$ $\displaystyle
E\left(\sigma_{i+1}^{2}\right)+\mathop{var}\nolimits\left(\mu_{i+1}\right)=\sigma_{i+1}^{2}+\mathop{var}\nolimits\left(\mu_{i+1}\right)$
$\displaystyle\mathop{var}\nolimits\left(\mu_{n}\right)$
$\displaystyle=^{(\ref{eMn})}$ $\displaystyle\sigma^{2}_{n+1}.$

In summary, we have the following model in $n$-way classification:

$\begin{array}[]{rclrcl}E(Y)&=&\mu,&\mathop{var}\nolimits(Y)&=&\sigma^{2}=\sum_{i=1}^{n+1}\sigma_{i}^{2}\\\\[5.69054pt]
E(Y|\mu_{1})&=&\mu_{1},&\mathop{var}\nolimits(Y|\mu_{1})&=&\sigma_{1}^{2}\\\\[2.84526pt]
\forall i,1\leq i\leq n-1\quad
E(\mu_{i}|\mu_{i+1})&=&\mu_{i+1},&\mathop{var}\nolimits(\mu_{i}|\mu_{i+1})&=&\sigma_{i+1}^{2}\\\\[2.84526pt]
E\left(\mu_{n}\right)&=&\mu,&\mathop{var}\nolimits\left(\mu_{n}\right)&=&\sigma^{2}_{n+1}.\end{array}$

####  7.1.3 Properties of a Sample Mean with N-way Classification

We would like to estimate the unknown parameter of interest $\mu=E(Y)$ (for
example, mean execution time) based on (balanced555In ‘balanced’ experiments,
the numbers of repetitions at each level are constant, i.e. every binary is
executed the same number of times, every execution makes the same number of
measurements, etc.) measurements of $Y$. In this section, we will show that
the (sample) arithmetic mean is an unbiased estimator of $\mu$ and that it is
asymptotically normal. We will also derive the variance of this estimate, so
that we can later construct a confidence interval for $\mu$. Let $n_{i}$ be
the numbers of repetitions at each level of the experiment. With three levels,
we thus have a 2-way classification and $n_{3}$ is the number of binaries,
$n_{2}$ the number of executions per each binary, and $n_{1}$ number of steady
state measurements per each execution, and hence $n_{1}n_{2}n_{3}$ is the
total number of measurements. For $n$-way classification, we denote the sample
mean $\overline{Y}$ as

$\overline{Y}=\overline{Y}_{\underbrace{\bullet\dots\bullet}_{n+1}}=\frac{1}{\prod_{i=1}^{n+1}n_{i}}\left(\sum_{j_{n+1}=1}^{n_{n+1}}\sum_{j_{n}=1}^{n_{n}}\dots\sum_{j_{1}=1}^{n_{1}}Y_{j_{n+1}j_{n}\dots
j_{1}}\right)$ (4)

######  Lemma 7.3 (Lindeberg–Levy Central Limit Theorem)

Let $X_{1}\dots X_{n}$ be independent identically distributed with mean $\mu$
and finite positive variance $\sigma^{2}$. Then,
$\overline{X}_{\bullet}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ has an
_asymptotically_ normal distribution with mean $\mu$ and variance
$\sigma^{2}/n$, which we denote as $\overline{X}_{\bullet}\approx
N\left(\mu,\sigma^{2}/n\right).$

###### Lemma 7.4

Let $X_{1}\dots X_{n}$ be independent identically distributed normal
variables, $X_{i}\sim N\left(\mu_{i},\sigma^{2}_{i}\right)$. From the
properties of normal distribution [Wasserman (2004)], it follows that
$\overline{X}_{\bullet}$ has normal distribution with mean
$\overline{\mu}_{\bullet}$ and variance $\overline{\sigma^{2}}_{\bullet}/n$,
which we denote as $\overline{X}_{\bullet}\sim
N\left(\overline{\mu}_{\bullet},{\overline{\sigma^{2}}_{\bullet}}/n\right)$

###### Lemma 7.5

Let $f(t;\mu_{1},\sigma_{1})$,$f(t;\mu_{2},\sigma_{2})$ be density functions
of normal variables with means $\mu_{1}$, $\mu_{2}$ and variances
$\sigma_{1}^{2},\sigma_{2}^{2}$

$\int
f(\tau;\mu_{1},\sigma_{1})f(t-\tau;\mu_{2},\sigma_{2}))d\tau=f\left(t;\mu_{1}+\mu_{2},\sqrt{\sigma^{2}_{1}+\sigma^{2}_{2}}\right).$

In other words, a convolution of density functions of normal variables (the
left-hand side of the equation) is also a density function of a normal
variable. Moreover, the normal variable has mean $\mu_{1}+\mu_{2}$ and
variance $\sigma_{1}^{2}+\sigma_{2}^{2}$ (by a known property of the normal
distribution).

###### Lemma 7.6

Let $X$, $Y$ be random variables with expectations and finite variance, $X\sim
N\left(\mu_{X},\sigma^{2}_{X}\right)$ and $Y|\left[X=x\right]\sim
N\left(x,\sigma^{2}\right)$. Then, $Y\sim
N\left(\mu_{X},\sigma^{2}_{X}+\sigma^{2}\right)$.

_Proof._ Let $f$ be the probability density function of normal distribution
with mean $\mu$ and variance $\sigma^{2}$:

$f(x;\mu,\sigma)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right),\quad\mathrm{where}\;\exp(z)=e^{z}$

The density functions of $X$ and $Y|X$ from Lemma 7.6 are:

$f_{X}(x)=f(x;\mu_{X},\sigma_{X}),\quad f_{Y|X}(y|x)=f_{Y|x}(y)=f(y;x,\sigma)$

By the definition of conditional density:

$f_{Y,X}(y,x)=f_{Y|X}(y|x)\cdot f_{X}(x)$

It follows, that:

$\displaystyle f_{Y}(y)$ $\displaystyle=$ $\displaystyle\int
f_{Y,X}(y,x)\,dx=\int f_{Y|X}(y|x)f_{X}(x)\,dx$ $\displaystyle=$
$\displaystyle\int\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(y-x)^{2}}{2\sigma^{2}}\right)\cdot\frac{1}{\sigma_{X}\sqrt{2\pi}}\exp\left(-\frac{(x-\mu_{X})^{2}}{2\sigma_{X}^{2}}\right)dx$
$\displaystyle\left[\stackrel{{\scriptstyle\mathrm{substituting}}}{{u=x-\mu_{X}}}\right]$
$\displaystyle=$
$\displaystyle\int\frac{1}{\sigma\sqrt{2\pi}}\exp\left(-\frac{(y-\mu_{X}-u)^{2}}{2\sigma^{2}}\right)\cdot\frac{1}{\sigma_{X}\sqrt{2\pi}}\exp\left(-\frac{u^{2}}{2\sigma_{X}^{2}}\right)du$
$\displaystyle=$ $\displaystyle\int
f(y-u;\mu_{X},\sigma)f(u;0,\sigma_{X})\,du\quad\quad\quad$
$\displaystyle=^{(L\ref{lGauss})}$ $\displaystyle
f\left(y;\mu_{X},\sqrt{\sigma^{2}+\sigma^{2}_{X}}\right)$

Hence, $Y\sim N(\mu_{X},\sigma^{2}_{X}+\sigma^{2}).\quad\qed$

With these lemmas, we can infer asymptotic distributions of estimators in our
$n$-way model. Informally, the basic idea is simple. We get asymptotic normal
distributions for sample means of executions by the Central Limit Theorem
(Lemma 7.3), which propagates to the grand mean $\overline{Y}$ by the
properties of the normal distribution (Lemma 7.4). Similarly, the sample mean
of means of executions within a binary ($\mu_{1}|\mu_{2}$) is asymptotically
normal by the Central Limit Theorem, which propagates to the ‘mean of the
grand mean’ by properties of the normal distribution. The two are then joined
by convolution, with the help of Lemma 7.5.

Formally, the application of Lemma 7.3 on sample means for individual
executions is described as:

$\forall j_{n+1}\dots j_{2},\qquad\overline{Y}_{j_{n+1}\dots
j_{2}\bullet}\big{|}\left[{(\mu_{1})}_{j_{n+1}\dots j_{2}}=m\right]\approx
N\left(m,\frac{\sigma_{1}^{2}}{n_{1}}\right)$

The meaning is that for a particular representation $m$ of the respective
execution mean, we get a normally distributed sample mean with this
expectation $m$. The ranges for $j_{n+1}\dots j_{1}$ are as in (4). However,
for simplicity, we just will write

$\forall j_{n+1}\dots j_{2},\qquad\overline{Y}_{j_{n+1}\dots
j_{2}\bullet}\big{|}{(\mu_{1})}_{j_{n+1}\dots j_{2}}\approx
N\left({(\mu_{1})}_{j_{n+1}\dots j_{2}},\frac{\sigma_{1}^{2}}{n_{1}}\right)$

Further summarisation to a higher level (i.e. binaries) keeps the asymptotic
normal distribution by Lemma 7.4:

$\forall j_{n+1}\dots j_{3},\qquad\overline{Y}_{j_{n+1}\dots
j_{3}\bullet\bullet}\big{|}{\overline{(\mu_{1})}}_{j_{n+1}\dots
j_{3}\bullet}\approx N\left({\overline{(\mu_{1})}}_{j_{n+1}\dots
j_{3}\bullet},\frac{\sigma_{1}^{2}}{n_{1}n_{2}}\right)$

We now keep applying Lemma 7.4 for all $i$, $1<i\leq n$:

$\forall j_{n+1}\dots j_{i+1},\qquad\overline{Y}_{j_{n+1}\dots
j_{i+1}\underbrace{\bullet\dots\bullet}_{i}}\big{|}{\overline{(\mu_{1})}}_{j_{n+1}\dots
j_{i+1}\underbrace{\bullet\dots\bullet}_{i-1}}\approx
N\left({\overline{(\mu_{1})}}_{j_{n+1}\dots
j_{i+1}\underbrace{\bullet\dots\bullet}_{i-1}},\frac{\sigma_{1}^{2}}{\prod_{k=1}^{i}n_{i}}\right)$

finally getting for $i=n+1$ (again by Lemma 7.4)

$\overline{Y}_{\underbrace{\bullet\dots\bullet}_{n+1}}\big{|}{\overline{(\mu_{1})}}_{\underbrace{\bullet\dots\bullet}_{n}}\approx
N\left({\overline{(\mu_{1})}}_{\underbrace{\bullet\dots\bullet}_{n}},\frac{\sigma_{1}^{2}}{\prod_{k=1}^{n+1}n_{k}}\right)$

which could be written simply as

$\overline{Y}\big{|}\overline{\mu_{1}}\approx
N\left(\overline{\mu_{1}},\frac{\sigma_{1}^{2}}{\prod_{k=1}^{n+1}n_{k}}\right)$
(5)

Now we need to derive the distribution of $\overline{\mu_{1}}$. This can be
done quite similarly to the distribution of $\overline{Y}|\overline{\mu_{1}}$.
By applying the Central Limit Theorem (Lemma 7.3) on (unknown) means of
executions, we get an asymptotically normal distribution of the sample mean of
these means:

$\forall j_{n+1}\dots j_{3},\qquad{\overline{(\mu_{1})}}_{j_{n+1}\dots
j_{3}\bullet}\big{|}{(\mu_{2})}_{j_{n+1}\dots j_{3}}\approx
N\left({(\mu_{2})}_{j_{n+1}\dots j_{3}},\frac{\sigma_{2}^{2}}{n_{2}}\right)$

The ranges for $j_{n+1}\dots j_{1}$ are as in (4). Hence, in a 2-way
classification model, sample execution means $\overline{\mu_{1}}$ have two
sums (over binaries and over executions) and sample binary means
$\overline{\mu_{2}}$ have a single sum (over binaries). By applying Lemma 7.4,
we get

$\forall j_{n+1}\dots j_{4},\qquad{\overline{(\mu_{1})}}_{j_{n+1}\dots
j_{4}\bullet\bullet}\big{|}{\overline{(\mu_{2})}}_{j_{n+1}\dots
j_{4}\bullet}\approx N\left({\overline{(\mu_{2})}}_{j_{n+1}\dots
j_{4}\bullet},\frac{\sigma_{2}^{2}}{n_{2}n_{3}}\right)$

Now we keep applying Lemma 7.4 for all $i$, $1<i<n$:

$\forall j_{n+1}\dots j_{i+2},\,\,{\overline{(\mu_{1})}}_{j_{n+1}\dots
j_{i+2}\underbrace{\bullet\dots\bullet}_{i}}\big{|}{\overline{(\mu_{2})}}_{j_{n+1}\dots
j_{i+2}\underbrace{\bullet\dots\bullet}_{i-1}}\approx
N\left({\overline{(\mu_{2})}}_{j_{n+1}\dots
j_{i+2}\underbrace{\bullet\dots\bullet}_{i-1}},\frac{\sigma_{2}^{2}}{\prod_{k=2}^{i+1}n_{k}}\right)$

Finally, for $i=n$, by Lemma 7.4 we get

${\overline{(\mu_{1})}}_{\underbrace{\bullet\dots\bullet}_{n}}\big{|}{\overline{(\mu_{2})}}_{\underbrace{\bullet\dots\bullet}_{n-1}}\approx
N\left({\overline{(\mu_{2})}}_{\underbrace{\bullet\dots\bullet}_{n-1}},\frac{\sigma_{2}^{2}}{\prod_{k=2}^{n+1}n_{k}}\right)$

which can be written as

${\overline{\mu_{1}}}|{\overline{\mu_{2}}}\approx
N\left({\overline{\mu_{2}}},\frac{\sigma_{2}^{2}}{\prod_{k=2}^{n+1}n_{k}}\right)$
(6)

The same procedure for $\mu_{2}|\mu_{3}$ gives

${\overline{\mu_{2}}}|{\overline{\mu_{3}}}\approx
N\left({\overline{\mu_{3}}},\frac{\sigma_{3}^{2}}{\prod_{k=3}^{n+1}n_{k}}\right)$

For the general case of $\mu_{i}|\mu_{i+1}$, $1\leq i\leq n-1$, we then get

${\overline{\mu_{i}}}|{\overline{\mu_{i+1}}}\approx
N\left({\overline{\mu_{i+1}}},\frac{\sigma_{i+1}^{2}}{\prod_{k={i+1}}^{n+1}n_{k}}\right)$
(7)

The means at the highest level of non-determinism (the means of binaries in
the case of 2-way classification), the $\mu_{n}$, come from a single (non-
conditional) distribution with mean $\mu$ and variance $\sigma_{n+1}^{2}$. By
the Central Limit Theorem, Lemma 7.3, their sample mean is thus normally
distributed:

$\overline{\mu_{n}}\approx N\left(\mu,\frac{\sigma_{n+1}^{2}}{n_{n+1}}\right)$
(8)

By this we have the last missing bit for Lemma 7.6. Now, by using Lemma 7.6 on
(7) and (8), we get

$\overline{\mu_{n-1}}\approx
N\left(\mu,\frac{\sigma_{n+1}^{2}}{n_{n+1}}+\frac{\sigma_{n}^{2}}{n_{n}n_{n+1}}\right)$

and then by further applications of the lemma and (7) we get

$\overline{\mu_{1}}\approx
N\left(\mu,\sum_{i=2}^{n+1}\frac{\sigma_{i}^{2}}{\prod_{k=i}^{n+1}n_{k}}\right)$
(9)

Now we can apply Lemma 7.4 on (5) and (9), by which we get

$\overline{Y}\approx
N\left(\mu,\sum_{i=1}^{n+1}\frac{\sigma_{i}^{2}}{\prod_{k=i}^{n+1}n_{k}}\right)$
(10)

Hence, the sample arithmetic mean $Y$ is an unbiased estimator of $\mu$ with
$n$-way classification, is normally distributed, and we have an expression for
its variance.

###  7.2 Confidence Interval for the Mean of One System

Within the model described in the Section 7.1, we can construct a confidence
interval for the mean execution time ($\mu$) of a single system. In Section
7.3, we extend this to the confidence interval for the ratio of means of two
systems.

We show how to construct the interval using two alternative methods. The
bootstrap method is intuitively simple and works for additional metrics, such
as the median, as well as the mean. The parametric method based on asymptotic
normality of the model (shown in Section 7.1), only works for the mean, but
has the potential to provide narrower intervals for larger sample sizes. We
compare the two methods empirically in Section 8.

####  7.2.1 Bootstrap Confidence Interval for One System

Let us assume that we have $n+1$ levels of hierarchy in our experiment and
that we decided to repeat $n_{n+1},n_{n},\dots n_{1}$ times at each level.
When we have completed benchmarking, we would thus have
$n_{n+1}n_{n}\cdot\dots\cdot n_{1}$ measurements, which we denote as

$Y_{j_{n+1}j_{n}\dots j_{1}}\qquad\left(\forall i,1\leq i\leq n+1,\quad 1\leq
j_{i}\leq n_{i}\right)$

We assume that the measurements from a single execution, that is
$Y_{j_{n+1}\dots j_{2}\bullet}$ for any fixed $j_{n+1}\dots j_{2}$, are
independent identically distributed, and we assume independence of means at
higher levels as described in the previous section.

Very informally, the core idea behind bootstrap is to simulate (many)
experiments based on the real data and then calculate what we want on the
simulated data, in our case to calculate the confidence interval. So instead
of one realization of $\overline{Y}$ which comes from our real data, we can
have many realizations, generated by simulation, which is much faster than
real experimentation. From these simulated realizations, we construct the
confidence interval by selecting appropriate quantiles. Detailed discussion
and formulation of different bootstrap methods can be found in davison.

The simulation of the new realizations of $\overline{Y}$ in each iteration
randomly selects a subset of real data
($Y_{\underbrace{\bullet\dots\bullet}_{n+1}})$ with replacement. This means
that the size of the subset will be the same as of the original set
($n_{n+1}\cdot\dots\cdot n_{1}$), but some data points may not be present,
while others can be presented multiple times. In each iteration, the
simulation calculates a sample arithmetic mean of the subset. After all
iterations finish, we estimate the 0.025 and 0.975 quantiles for a 95%
confidence interval (if we have 1000 iterations, the quantiles can be
estimated as 25th and 975th ordered values, although different estimators for
quantiles exist). Pseudo-code for the described method is shown in Figure 1.
For a 95% interval, we select $\alpha=0.05$.

Out of the wide variety of bootstrap methods, the key decisions here were how
to do resampling, and how to construct the confidence interval given the
simulated means. For resampling, depending on the underlying distribution and
the statistic of interest, it is sometimes better to resample with replacement
only at higher levels of the hierarchy, but keep the lower levels intact [Ren
et al. (2010), Davison and Hinkley (1997)]. A naïve alternative is then also
to ignore the structure of the model and resample at random from all
measurements. We cover these alternatives in our evaluation later. The method
of resampling shown in Figure 1, replacement at all levels, seems to work best
(or at least not worse than others), and is a common default.

For the construction of the confidence interval given the simulated means, we
use the percentile method. This method is sensitive to non-symmetrical
distributions of the statistic, but it should not be much of a problem here as
we have shown earlier that $\overline{Y}$ is asymptotically normal. Still, one
can easily plug in alternative bootstrap methods for confidence interval
construction [Davison and Hinkley (1997)]. Statistical software packages, such
as R, implement plenty of different methods.

⬇

Input: $n$, $(n_{1},\dots,n_{n})$, $Y_{j_{n+1}j_{n}\dots j_{1}}$ where
$\forall i,\,1\leq j_{i}\leq n_{i}$, $nIterations=1000$, $\alpha=0.05$

Output: $lower$, $upper$

Uses: mean(x) ... arithmetic average

Uses: quantile(probability, x) ... select a sample quantile

Uses: resample (replacement, x) ... random resampling

$simulatedMeans$ = new vector[ $nIterations$ ]

foreach $iteration$ in 1..$nIterations$ {

$simulatedMeasurements$ = new vector[ $n_{n+1}\cdot\dots\cdot n_{1}$ ]

foreach $j_{n+1}$ in resample( 1..$n_{n+1}$, replacement = yes ) {

foreach $j_{n}$ in resample( 1..$n_{n}$, replacement = yes ) {

$\dots$

foreach $j_{1}$ in resample( 1..$n_{1}$, replacement = yes ) {

append $Y_{j_{n+1}j_{n}\dots j_{1}}$ to $simulatedMeasurements$

}

$\dots$

}

}

$simulatedMeans$[ $iteration$ ] = mean( $simulatedMeasurements$ )

}

$lower$ = quantile( probability = $\alpha/2$, $simulatedMeans$ )

$upper$ = quantile( probability = $1-\alpha/2$, $simulatedMeans$ )

Figure 1: Bootstrap Confidence Interval for One System.

####  7.2.2 Asymptotic Confidence Interval for One System

Alternatively to bootstrap, we can construct a confidence interval for the
mean with $n$-way classification using the asymptotic normality of
$\overline{Y}$ (equation 10 of Section 7.1.3). For this we need to estimate
the unknown variance of the sample mean. It is easier to do this directly than
by estimating the individual variances $\sigma_{i}^{2}$, $1\leq i\leq n+1$.

We use the following estimator:

$S_{n+1}^{2}=\frac{1}{n_{n+1}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\left(\overline{Y}_{j_{n+1}\underbrace{\bullet\dots\bullet}_{n}}-\overline{Y}_{\underbrace{\bullet\dots\bullet}_{n+1}}\right)^{2}$

We show in Section 7.6 that $S_{n+1}^{2}/n_{n+1}$ is an unbiased estimator of
the variance of the sample mean, that is

$E\left(\frac{S_{n+1}^{2}}{n_{n+1}}\right)=\sum_{i=1}^{n+1}\frac{\sigma_{i}^{2}}{\prod_{k=i}^{n+1}n_{k}}$

Note this also means that as long as we have a hierarchical experiment ($n\geq
1$), $S_{n+1}^{2}$ is _not_ an unbiased estimator of $\sigma_{n+1}^{2}$.
Relying on asymptotic normality even after the unknown variance of the sample
mean is replaced by its estimate, we get an asymptotic ($1-\alpha$) confidence
interval for $\mu$:

$\overline{Y}\pm
u_{1-\frac{\alpha}{2}}\sqrt{\frac{S_{n+1}^{2}}{n_{n+1}}}=\overline{Y}\pm
u_{1-\frac{\alpha}{2}}\sqrt{\frac{1}{n_{n+1}(n_{n+1}-1)}\sum_{j_{n+1}=1}^{n_{n+1}}\left(\overline{Y}_{j_{n+1}\underbrace{\bullet\dots\bullet}_{n}}-\overline{Y}_{\underbrace{\bullet\dots\bullet}_{n+1}}\right)^{2}}$
(11)

We do not claim that $\overline{Y}$ has a $t$-distribution, as we do not
assume normality of $Y|\mu_{1}$, $\mu_{1}|\mu_{2}$, …, $\mu_{n-1}|\mu_{n}$,
and $\mu_{n}$. A one-way model that makes such normality assumptions can be
found in mcculloch, including the respective confidence interval, which uses
the $t$-distribution with $n_{n+1}-1$ degrees of freedom. For large numbers of
degrees of freedom, the $t$ distribution converges to normal, so the choice is
not important. For small number of degrees of freedom, say smaller than 30 and
definitely smaller than 20, the confidence intervals become wider with the
$t$-distribution than with the Normal distribution. This means that under the
normality assumptions, one should definitely use the $t$ distribution,
otherwise the interval would be too narrow (its coverage will be smaller than
the projected $1-\alpha$). In the practice, a larger than projected coverage
is usually regarded as better than a smaller one, so it makes sense to use the
$t$ distribution anyway.

###  7.3 Confidence Interval for Ratio of Means

In this section, we show how the method for constructing the confidence
interval for the mean of one system (Section 7.2) can be extended to a
confidence interval of the ratio of means of two systems. We will refer to
these systems as ‘old’ and ‘new’, and make the same set of assumptions for the
two systems as we did for the single system so far (independence and identical
distributions at multiple levels). We denote the corresponding random
variables for execution time as ${}^{O}Y$ and ${}^{N}Y$, and the means as
${}^{O}\mu=E(^{O}Y)$ and ${}^{N}\mu=E(^{N}Y)$. Thus, we now have a bivariate
distribution ($F$) of ${}^{ON}Y=(^{O}Y,^{N}Y)\sim F(x_{O},x_{N})$. We assume
that ${}^{O}Y$ and ${}^{N}Y$ are independent. The variances
${}^{N}\sigma_{i}^{2}$, ${}^{O}\sigma_{i}^{2}$ and expectations in the two
systems can differ. Even the distributions may differ. The parameter of
interest is now

$\theta=t(F)=\frac{{}^{N}\mu}{{}^{O}\mu}=\frac{\int
x_{N}\,dF(x_{O},x_{N})}{\int x_{O}\,dF(x_{O},x_{N})}$

We estimate $\theta$ using measurements of the two systems. Again, we use
balanced measurements and the same repetition counts ($n_{1},\dots,n_{n+1}$)
for both systems. It can be shown that $T$,

$T=t(\widehat{F})=\frac{\overline{{}^{N}Y}}{\overline{{}^{O}Y}}$

is an unbiased (plugin) estimator for $\theta$ (davison, Example 2.2). Hence,
we can estimate $\theta$ using the ratio of arithmetic averages of the two
systems, and we will do so both in the bootstrap and the asymptotic parametric
method.

We estimate the unknown variances $\mathop{var}\nolimits(^{O}Y)$ and
$\mathop{var}\nolimits(^{N}Y)$ by ${{}^{O}S}^{2}_{n+1}$ and
${{}^{N}S}^{2}_{n+1}$ from Section 7.1.

####  7.3.1 Bootstrap Confidence Interval for Ratio of Means

To construct a bootstrap interval for the ratio of means $\theta$, we need to
simulate many realisations of its estimator $T$. $\alpha/2$\- and
$1-\alpha/2$\- quantiles of these realisations form a $(1-\alpha)$ percentile
bootstrap confidence interval for $\theta$. The pseudo-code is shown in Figure
2. Discussion of more elaborate bootstrap methods for construction of the
ratio of means can be found in luxfra.

⬇

Input: $n$, $(n_{1},\dots,n_{n})$, ${}^{O,N}Y_{j_{n+1}j_{n}\dots j_{1}}$ where
$\forall i\,1\leq j_{i}\leq n_{i}$, $nIterations=1000$, $\alpha=0.05$

Output: $lower$, $upper$

Uses: mean(x) ... arithmetic average

Uses: quantile(probability, x) ... select a sample quantile

Uses: resample (replacement, x) ... random resampling

function simulateMean( $oldnew$ ) {

$simulatedMeasurements$ = new vector[ $n_{n+1}\cdot\dots\cdot n_{1}$ ]

foreach $j_{n+1}$ in resample( 1..$n_{n+1}$, replacement = yes ) {

foreach $j_{n}$ in resample( 1..$n_{n}$, replacement = yes ) {

$\dots$

foreach $j_{1}$ in resample( 1..$n_{1}$, replacement = yes ) {

append ${}^{oldnew}Y_{j_{n+1}j_{n}\dots j_{1}}$ to $simulatedMeasurements$

}

$\dots$

}

}

return mean($simulatedMeasurements$)

}

$simulatedRatios$ = new vector[ $nIterations$ ]

foreach $iteration$ in 1..$nIterations$ {

$simulatedRatios$[ $iteration$ ] = $\frac{\texttt{simulateMean( oldnew = 'N'
)}}{\texttt{simulateMean( oldnew = 'O' )}}$

}

$lower$ = quantile( probability = $\alpha/2$, $simulatedRatios$ )

$upper$ = quantile( probability = $1-\alpha/2$, $simulatedRatios$ )

Figure 2: Bootstrap Confidence Interval for Ratio of Means of Two Systems.

####  7.3.2 Asymptotic Confidence Interval for Ratio of Means

We can construct an asymptotic interval for $\theta$ using a theorem by
fieller, which gives confidence limits for the ratio of means of two normally
distributed variables.

######  Lemma 7.7 (Fieller’s Theorem [Fieller (1954)])

Let $X,Y$ be normally distributed random variables, not necessarily
independent. Let $x,y$ be unbiased estimates of the means $E(X),E(Y)$. Let
$v_{xx},v_{yy}$ be variances of $x$ and $y$ (note, not of $X$ and $Y$). Let
$v_{xy}$ be sample covariance of $x$ and $y$. Then, the confidence limits for
$E(X)/E(Y)$ are:

$\alpha_{1},\alpha_{2}=\frac{\left(xy-t^{2}v_{xy}\right)\mp\sqrt{\left(xy-t^{2}v_{xy}\right)^{2}-\left(x^{2}-t^{2}v_{xx}\right)\left(y^{2}-t^{2}v_{yy}\right)}}{x^{2}-t^{2}v_{xx}}$

where $\alpha_{1}$ is the lower limit ($-$ sign) and $\alpha_{2}$ is the upper
limit ($+$ sign). $t$ is the critical value for the two-tail $t$-distribution,
that is, for a 95% interval it is the 0.025th quantile of the respective
$t$-distribution. The interval will be bounded and non-trivial, that is
$-\infty<\alpha_{1}<\alpha_{2}<\infty$ for $0<t^{2}<x^{2}/v_{xx}$. Details on
other cases can be found in fieller and luxfra.

Note that the requirement of $t^{2}<x^{2}/v_{xx}$ from the theorem is
intuitive. It states that the estimator of the denominator (the sample mean of
the old system) should be statistically significantly different from zero. In
our case with sample means of execution times, this condition is hardly ever
violated. Nevertheless, it can still happen due to the statistical nature of
the condition and hence the condition needs to be checked. As we detail later
in our evaluation, we have seen the condition violated for a sample size of
two in our experiments. The requirement that $t^{2}>0$ is also easily met in
practice, as we are not interested in intervals for confidence close to zero.
In our case, we also assume independence of the two variables, which implies
$v_{xy}=0$. The $(1-\alpha)$ confidence interval for $\theta$ is hence:

$\frac{\overline{{}^{O}Y}\cdot\overline{{}^{N}Y}\mp\sqrt{\left(\overline{{}^{O}Y}\cdot\overline{{}^{N}Y}\right)^{2}-\left(\left(\overline{{}^{O}Y}\right)^{2}-t_{\frac{\alpha}{2},\nu}^{2}\cdot\frac{{{}^{O}S}^{2}_{n+1}}{n_{n+1}}\right)\left(\left(\overline{{}^{N}Y}\right)^{2}-t_{\frac{\alpha}{2},\nu}^{2}\cdot\frac{{{}^{N}S}^{2}_{n+1}}{n_{n+1}}\right)}}{\overline{{}^{O}Y}^{2}-t_{\frac{\alpha}{2},\nu}^{2}\cdot{{}^{O}S}^{2}_{n+1}/{n_{n+1}}}$
(12)

where $t_{\frac{\alpha}{2},\nu}$ is the $\frac{\alpha}{2}$\- quantile of the
$t$-distribution with $\nu=n_{n+1}-1$ degrees of freedom. The Fieller’s result
can be extended for unbalanced experiments, where the calculation of the
numbers of degrees of freedom becomes more complicated. Details can be found
in mratios and reportfieller. In our case of a hierarchical experiment, this
number could be modified to match the number of degrees of freedom in our
variance estimates. Note that Fieller’s method assumes normality of the two
variables, that is of ${}^{O}Y$ and ${}^{N}Y$, which will be violated in
practice — execution time is not normally distributed. Relying only on the
asymptotic normality of the mean and looking only for an asymptotic interval,
we would use the quantiles of the standard Normal distribution instead of the
$t$-distribution in the interval. To be conservative, we would in practice use
the $t$-distribution anyway, as in the case of the interval for a single
system. The robustness of Fieller’s method to deviations from normality is
touched in [von Luxburg and Franz (2009)].

For alternative methods for construction of the confidence interval for the
ratio of means we refer the reader to ratiocficmp and luxfra. A common
alternative is based on the delta method.

###  7.4 Experiment Planning

To obtain an unbiased result and realistic error bars (confidence interval),
one has to make sure that all the randomness of the system is encapsulated by
the experiment. For example, if the system of interest has non-deterministic
compilation and different binaries also differ in performance, one has to run
multiple binaries within that experiment. The cost of building a new binary
can be high, e.g. for the Mono platform, which we use for validation of our
method, the compilation time was about 20 minutes. If the fluctuations in
performance due to non-determinism in compilation are small ($\sigma_{S}$ is
small) compared to, say, the one due to non-determinism in execution
($\sigma_{B}$), which is the usual case in our experiments, time allocated to
experimentation would be better spent by running the existing binaries
multiple times rather than building new binaries. We can formalise such trade-
offs and, based on initial experiments that repeat a few times at all levels
to estimate the variances, we can calculate the optimum number of repetitions
at each level that results in the most precise result for a given
experimentation time.

Let us assume the random effects model in $n$-way classification of Section
7.1. We have run initial experiments to obtain estimates of the variances
$\sigma^{2}_{1},\sigma^{2}_{2},\dots,\sigma^{2}_{n+1}$. We will discuss later
which estimators can be used. We also assume we know the cost of a new
repetition at the (non-top) levels, $c_{1},c_{2},\dots,c_{n}$, so that the
total cost of experimentation is

$c=(c_{n}+\dots(c_{3}+(c_{2}+(c_{1}+n_{1})n_{2})n_{3}\dots)n_{n})n_{n+1}$ (13)

The unit of the costs ($c$, $c_{i}$) is the number of (the lowest-level)
measurements that could be obtained in that time. The practical meaning of the
costs and their use have been described informally in Section 6.1. With three
levels only, the total cost would be

$c=(c_{2}+(c_{1}+n_{1})n_{2})n_{3}$

where $n_{1}$ is the number of steady state measurements, $c_{1}$ is the
number of warm-up measurements of an execution that are not included into
summarization (cost of a new execution), $n_{2}$ is the number of executions,
$c_{2}$ is the number of measurements that could be done in the time needed to
build one binary, and $n_{3}$ is the number of binaries.

The costs $c_{i}$ can be calculated during the initial experiments that have
to be executed to estimate the variances $\sigma^{2}_{i}$. For the purpose of
experiment planning, precise estimates of the costs are not needed — the
planning is a back-of-the-envelope calculation. For instance, we can use
averages from the initial experiments. The optimisation problem is to find the
$n_{1},\dots,n_{n}$ that minimize $f$,

$f(n_{1},n_{2},\dots,n_{n+1})=\sum_{i=1}^{n+1}\frac{\sigma_{i}^{2}}{\prod_{k=i}^{n+1}n_{k}}$

Function $f$ is a measure of the precision of the result, as it is a measure
of the width of the confidence interval for the mean — (equation 11 in Section
7.2). In this section, we will show that the optimal numbers of repetitions
$n_{1},\dots,n_{n}$ are

$n_{1}=\sqrt{c_{1}\frac{\sigma_{1}^{2}}{\sigma_{2}^{2}}},\qquad\forall
i,1<i\leq
n,\,\,n_{i}=\sqrt{\frac{c_{i}}{c_{i-1}}\frac{\sigma_{i}^{2}}{\sigma_{i+1}^{2}}}$
(14)

The optimal number of repetitions thus only depends on two adjacent levels of
the experiment, no matter what the total number of levels is. For example, in
a 2-way model, we have

$n_{1}=\sqrt{c_{1}\frac{\sigma^{2}_{1}}{\sigma^{2}_{2}}}\quad\textrm{and}\quad
n_{2}=\sqrt{\frac{c_{2}}{c_{1}}\frac{\sigma^{2}_{2}}{\sigma^{2}_{3}}}\\\ $

The number of samples $n_{1}$ should be large when the cost for warm-up
$c_{1}$ is high and/or when the variance in samples $\sigma^{2}_{1}$ is higher
than the variance in executions $\sigma^{2}_{2}$. The number of executions
$n_{2}$ should be large when the cost of a build $c_{2}$ is larger than the
cost of a run $c_{1}$ and/or when the variance in executions $\sigma^{2}_{2}$
is larger than the variance in binaries $\sigma^{2}_{3}$.

###  7.5 Derivation of the Optimalisation Formula

For the proof, we will use the following notation

$\displaystyle\forall i,1<i\leq n+1,\,\,s_{i}$ $\displaystyle=$
$\displaystyle\sigma_{i}^{2}$ $\displaystyle\forall i,1<i\leq n+1,\,\,p_{i}$
$\displaystyle=$ $\displaystyle\prod_{k=i}^{n+1}n_{k}$ $\displaystyle k_{1}$
$\displaystyle=$ $\displaystyle n_{1},\qquad\forall i,1<i\leq
n,\,\,k_{i}=(c_{i-1}+k_{i-1})n_{i}$

Note that $k_{i}$ is a recurrent definition of the cost of the experiment $c$.
In the new notation, we need to find $n_{1},\dots,n_{n+1}$ that minimize $f$
under the condition that $g=0$ for a given budget for the experiment $c$ (in
fact we only care about $n_{1},\dots,n_{n}$):

$\displaystyle f$ $\displaystyle=$ $\displaystyle
f(n_{1},\dots,n_{n+1})=\sum_{i=1}^{n+1}\frac{s_{i}}{p_{i}(n_{i},\dots,n_{n+1})}=\sum_{i=1}^{n+1}\frac{s_{i}}{p_{i}}$
$\displaystyle g$ $\displaystyle=$ $\displaystyle
g(n_{1},\dots,n_{n+1})=k_{n+1}(n_{1},\dots,n_{n+1})-c=k_{n+1}-c$ (15)

Note that we sometimes omit the formal arguments of functions for readability.
By the Lagrange Multiplier Theorem, optimum values can only be among solutions
of the system of equations:

$\displaystyle\forall i,1\leq i\leq n+1,\,\,\frac{\partial f}{\partial
n_{i}}(n_{1},\dots,n_{n+1})+\lambda\frac{\partial g}{\partial
n_{i}}(n_{1},\dots,n_{n+1})$ $\displaystyle=$ $\displaystyle 0$ $\displaystyle
g(n_{1},\dots,n_{n+1})$ $\displaystyle=$ $\displaystyle 0$ (16)

The partial derivatives are expressed in our notation as follows:

$\frac{\partial g}{\partial n_{1}}=n_{2}n_{3}\cdot\dots\cdot
n_{n+1}=p_{2}=\frac{k_{1}}{n_{1}}p_{2},\quad\frac{\partial g}{\partial
n_{i}}=^{(*)}\frac{k_{i}}{n_{i}}p_{i+1}$ $\frac{\partial f}{\partial
n_{1}}=-\frac{s_{1}}{n_{1}p_{1}},\qquad\frac{\partial f}{\partial
n_{2}}=-\frac{s_{1}}{n_{2}p_{1}}-\frac{s_{2}}{n_{2}p_{2}},\qquad\frac{\partial
f}{\partial n_{i}}=-\sum_{k=1}^{i}\frac{s_{k}}{n_{i}p_{k}}$

The derivation of $\frac{\partial f}{\partial n_{i}}$ is straightforward, but
the derivation of the marked $\frac{\partial g}{\partial n_{i}}$ deserves some
explanation. From the definitions of $k_{i}$ and $g$ we have that

$\displaystyle\frac{\partial g}{\partial n_{i}}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
n_{i}}(k_{n+1}-c)=\frac{\partial}{\partial
n_{i}}k_{n+1}=\frac{\partial}{\partial
n_{i}}(c_{n}+k_{n})n_{n+1}=\frac{\partial}{\partial
n_{i}}(c_{n}+(c_{n-1}+k_{n-1})n_{n}\,)n_{n+1}$ $\displaystyle=$
$\displaystyle\\!\\!\dots=\frac{\partial}{\partial
n_{i}}(c_{n}+(c_{n-1}+\dots(c_{i-1}+k_{i-1})n_{i}\dots\,)n_{n}\,)n_{n+1}$
$\displaystyle=$ $\displaystyle\left(\frac{\partial}{\partial
n_{i}}(c_{i-1}+k_{i-1})p_{i}\right)+0=(k_{i-1}+c_{i-1})p_{i+1}=p_{i+1}\frac{k_{i}}{n_{i}}.$

We will solve the system (7.5) by substitutions and induction, starting from
partial derivatives for $\frac{\partial}{\partial n_{1}}$. We first express
$n_{i}^{2}$ as follows:

$\displaystyle\frac{\partial f}{\partial n_{1}}+\lambda\frac{\partial
g}{\partial n_{1}}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\lambda\frac{k_{1}}{n_{1}}p_{2}$ $\displaystyle=$
$\displaystyle\frac{s_{1}}{n_{1}p_{1}}$ $\displaystyle\lambda k_{1}p_{1}p_{2}$
$\displaystyle=$ $\displaystyle s_{1}\qquad|\textrm{expand }k_{1},p_{1}$
$\displaystyle\lambda p^{2}_{2}n^{2}_{1}$ $\displaystyle=$ $\displaystyle
s_{1}$ $\displaystyle n^{2}_{1}$ $\displaystyle=$
$\displaystyle\frac{s_{1}}{\lambda p^{2}_{2}}=\frac{s_{1}}{\lambda
n^{2}_{2}p^{2}_{3}},\,\,n^{2}_{2}=\frac{s_{1}}{\lambda
n^{2}_{1}p^{2}_{3}},\,\,n^{2}_{3}=\frac{s_{1}}{\lambda
n^{2}_{1}n^{2}_{2}p^{2}_{4}},\ $ $\displaystyle n^{2}_{i}$ $\displaystyle=$
$\displaystyle\frac{s_{1}}{\lambda n^{2}_{1}n^{2}_{2}\dots
n^{2}_{i-1}p^{2}_{i+1}}.$ (17)

Now we can get a solution for $n_{1}$ from the equation for
$\frac{\partial}{\partial n_{2}}$

$\displaystyle\frac{\partial f}{\partial n_{2}}+\lambda\frac{\partial
g}{\partial n_{2}}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\lambda\frac{k_{2}}{n_{2}}p_{3}$ $\displaystyle=$
$\displaystyle\frac{s_{1}}{n_{2}p_{1}}+\frac{s_{2}}{n_{2}p_{2}}$ (18)
$\displaystyle\lambda k_{2}p_{1}p_{3}$ $\displaystyle=$ $\displaystyle
s_{1}+s_{2}n_{1}\qquad|\textrm{expand }k_{2},p_{1}$
$\displaystyle\lambda(c_{1}+n_{1})n_{1}n^{2}_{2}p^{2}_{3}$ $\displaystyle=$
$\displaystyle s_{1}+s_{2}n_{1}\qquad|\textrm{substitute }n^{2}_{2}\textrm{
using (\ref{eNI2})}\qquad$ (19) $\displaystyle(c_{1}+n_{1})s_{1}$
$\displaystyle=$ $\displaystyle s_{1}n_{1}+s_{2}n^{2}_{1}$ (20) $\displaystyle
n_{1}$ $\displaystyle=$ $\displaystyle\sqrt{\frac{s_{1}c_{1}}{s_{2}}}$ (21)

which is the initial case of (14) which we need to prove. We follow by
induction over $i$. Assuming (14) holds for $n_{1},\dots,n_{i}$, we prove it
for $n_{i+1}$:

$\displaystyle\frac{\partial f}{\partial n_{i+2}}+\lambda\frac{\partial
g}{\partial n_{i+2}}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\lambda\frac{k_{i+2}}{n_{i+2}}p_{i+3}$ $\displaystyle=$
$\displaystyle\frac{1}{n_{i+2}}\sum_{k=1}^{i+2}\frac{s_{k}}{p_{k}}=\frac{1}{n_{i+2}}\left(\sum_{k=1}^{i}\frac{s_{k}}{p_{k}}+\frac{s_{i+1}}{p_{i+1}}+\frac{s_{i+2}}{p_{i+2}}\right)$
$\displaystyle\lambda k_{i+2}\cdot p_{i+3}$ $\displaystyle=^{(*)}$
$\displaystyle\lambda
k_{i}p_{i+1}+\frac{s_{i+1}}{p_{i+1}}+\frac{s_{i+2}}{p_{i+2}}\qquad|\textrm{expand
}k_{i+2}$ $\displaystyle\lambda(c_{i+1}+(c_{i}+k_{i})n_{i+1})n_{i+2}$
$\displaystyle\cdot$ $\displaystyle p^{2}_{i+3}n_{i+2}n_{i+1}=\lambda
k_{i}p^{2}_{i+1}+s_{i+1}+s_{i+2}n_{i+1}$

The transformation marked by $(*)$ is a substitution for the same equation,
but for $\frac{\partial}{\partial n_{i}}$. By expanding the left hand side and
substituting $n^{2}_{i+2}$ using (7.5), we get

$\displaystyle c_{i+1}\frac{s_{1}}{n^{2}_{1}\cdot\dots\cdot
n^{2}_{i}n_{i+1}}+c_{i}\frac{s_{1}}{n^{2}_{1}\cdot\dots\cdot n^{2}_{i}}$
$\displaystyle=$ $\displaystyle s_{i+1}+s_{i+2}n_{i+1}$ $\displaystyle
c_{i+1}s_{1}+c_{i}s_{1}n_{i+1}$ $\displaystyle=$ $\displaystyle
s_{i+1}\,n^{2}_{1}\cdot\dots\cdot
n^{2}_{i}\,n_{i+1}+s_{i+2}\,n^{2}_{1}\cdot\dots\cdot n^{2}_{i}\,n^{2}_{i+1}$

It is easy to simplify $n^{2}_{1}\cdot\dots\cdot n^{2}_{i}$ by the induction
assumption for (14),

$n^{2}_{1}\cdot\dots\cdot
n^{2}_{i}=\frac{c_{1}}{1}\frac{s_{1}}{s_{2}}\cdot\frac{c_{2}}{c_{1}}\frac{s_{2}}{s_{3}}\cdot\frac{c_{3}}{c_{2}}\frac{s_{3}}{s_{4}}\cdot\dots\cdot\frac{c_{i}}{c_{i-1}}\frac{s_{i}}{s_{i+1}}=\frac{c_{i}s_{1}}{s_{i+1}}$

and thus

$\displaystyle c_{i+1}s_{1}$ $\displaystyle=$ $\displaystyle
s_{i+2}s_{1}\frac{c_{i}}{s_{i+1}}n^{2}_{i+1}$ $\displaystyle n_{i+1}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{c_{i+1}s_{1}s_{i+1}}{c_{i}s_{1}s_{i+2}}}=\sqrt{\frac{c_{i+1}}{c_{i}}\frac{s_{i+1}}{s_{i+2}}}$

which finishes the induction step. We have solutions of (7.5) for

$n_{1}=\sqrt{\frac{s_{1}c_{1}}{s_{2}}},\quad\forall i,2\leq i\leq
n,\,\,n_{i}=\sqrt{\frac{c_{i}s_{i}}{c_{i-1}s_{i+1}}}$

Technically, we could find a solution for $n_{n+1}$ based on the budget $c$
and solutions for $n_{1},\dots,n_{n}$, but we are not really interested in
$n_{n+1}$. Still, our solution of (7.5) need not be the minimum of (7.5). To
verify that it actually is, we proceed as follows. First, we reduce the
optimisation problem to a single function $h$, which we get by eliminating
$n_{n+1}$ from $f$ using $g$. From (7.5), since $g=0$, we get

$(c_{n}+k_{n})n_{n+1}-c=0\quad\text{hence}\quad n_{n+1}=\frac{c}{c_{n}+k_{n}}$

and

$\displaystyle f$ $\displaystyle=$ $\displaystyle
f(n_{1},\dots,n_{n+1})=\frac{1}{n_{n+1}}\left(s_{n+1}+\sum_{i=1}^{n}\frac{s_{i}}{n_{i}\cdot\dots\cdot
n_{n}}\right)$

For $h$ we have that

$\displaystyle h$ $\displaystyle=$ $\displaystyle
f(n_{1},\dots,n_{n})=\frac{c_{n}+k_{n}}{c}\left(s_{n+1}+\sum_{i=1}^{n}\frac{s_{i}}{n_{i}\cdot\dots\cdot
n_{n}}\right)=$ (22) $\displaystyle=$
$\displaystyle\frac{1}{c}\left(c_{n}+n_{1}\cdot\dots\cdot
n_{n}+\sum_{i=1}^{n-1}c_{i}\cdot n_{i+1}\cdot\dots\cdot
n_{n}\right)\cdot\left(s_{n+1}+\sum_{i=1}^{n}\frac{s_{i}}{n_{i}\cdot\dots\cdot
n_{n}}\right)$ (23)

Our optimisation problem of (7.5) is equivalent to minimizing $h$. From (22)
we know that $h$ has the form of

$h(n_{1},\dots,n_{n})=\alpha+\sum\beta n^{*}+\sum\gamma\frac{1}{n^{*}}$ (24)

where $\alpha,\beta,\gamma$ are non-negative real constants and $n^{*}$ are
products of combinations of variables from $n_{1},\dots,n_{n}$. For 1-way
classification, we have

$h(n_{1})=\alpha+\beta
n_{1}+\gamma\frac{1}{n_{1}}=\frac{(c_{1}+n_{1})(s_{1}+s_{2}n_{1})}{n_{1}c}$

for 2-way classification we have

$\displaystyle h(n_{1},n_{2})$ $\displaystyle=$
$\displaystyle\alpha+\beta_{1}n_{1}+\beta_{2}n_{2}+\beta_{3}n_{1}n_{2}+\gamma_{1}\frac{1}{n_{1}}+\gamma_{2}\frac{1}{n_{2}}+\gamma_{3}\frac{1}{n_{1}n_{2}}=$
$\displaystyle=$
$\displaystyle\frac{(c_{2}+n_{2}c_{1}+n_{1}n_{2})(s_{1}+s_{2}n_{1}+s_{3}n_{1}n_{2})}{n_{1}n_{2}c}$

From this general form it is easy to see that $h$ is continuous and positive
for $n_{i}>0$. Every element of the sum in (24) is always positive. $h$ does
not have a maximum, because it is unbounded for very small values of any and
all of its arguments, as well as for very large values. Each term of (24) is
convex, hence the sum of terms, $h$, is convex, so $h$ does not have an
inflection point, but has a global minimum.666We could have found the
solutions directly through partial differentiation of $h$. We have done this
for several values of $n$ using a symbolic algebra system, but we chose to use
the Lagrange multiplier method for the general case, as the computation is
simpler.

###  7.6 Estimating Unknown Variances

Estimating the unknown variances $\sigma_{1}^{2},\dots,\sigma_{n+1}^{2}$ is
harder than it may seem. Let us first define these ‘naïve’ estimators

$\displaystyle S_{1}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\prod_{k=2}^{n+1}n_{k}}\,\frac{1}{n_{1}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\dots\sum_{j_{1}=1}^{n_{1}}\left(Y_{j_{n+1}\dots
j_{1}}-\overline{Y}_{j_{n+1}\dots j_{2}\bullet}\right)^{2}$ $\displaystyle
S_{2}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\prod_{k=3}^{n+1}n_{k}}\,\frac{1}{n_{2}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\dots\sum_{j_{2}=1}^{n_{2}}\left(\overline{Y}_{j_{n+1}\dots
j_{2}\bullet}-\overline{Y}_{j_{n+1}\dots j_{3}\bullet\bullet}\right)^{2}$

for the general case of $i,2\leq i\leq n$

$S_{i}^{2}=\frac{1}{\prod_{k=i+1}^{n+1}n_{k}}\,\frac{1}{n_{i}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\dots\sum_{j_{i}=1}^{n_{i}}\left(\overline{Y}_{j_{n+1}\dots
j_{i}\underbrace{\bullet\dots\bullet}_{i-1}}-\overline{Y}_{j_{n+1}\dots
j_{i+1}\underbrace{\bullet\dots\bullet}_{i}}\right)^{2}$

and finally the already defined

$S_{n+1}^{2}=\frac{1}{n_{n+1}-1}\sum_{j_{n+1}=1}^{n_{n+1}}\left(\overline{Y}_{j_{n+1}\underbrace{\bullet\dots\bullet}_{n}}-\overline{Y}_{\underbrace{\bullet\dots\bullet}_{n+1}}\right)^{2}$

These estimators are, apart from $S_{1}^{2}$, not unbiased estimators for
$\sigma_{1}^{2},\dots,\sigma_{n+1}^{2}$. We will show later in this section
that

$\displaystyle E\left(S_{1}^{2}\right)$ $\displaystyle=$
$\displaystyle\sigma_{1}^{2}$ $\displaystyle E\left(S_{2}^{2}\right)$
$\displaystyle=$ $\displaystyle\sigma_{2}^{2}+\frac{\sigma_{1}^{2}}{n_{1}}$
$\displaystyle E\left(S_{3}^{2}\right)$ $\displaystyle=$
$\displaystyle\sigma_{3}^{2}+\frac{\sigma_{2}^{2}}{n_{2}}+\frac{\sigma_{1}^{2}}{n_{1}n_{2}}$

and for the general case of $i,2\leq i\leq n+1$

$E\left(S_{i}^{2}\right)=\sigma_{i}^{2}+\sum_{k=1}^{i-1}\frac{\sigma_{k}^{2}}{\prod_{l=k}^{i-1}n_{l}}=\sigma_{i}^{2}+E\left(\frac{S_{i-1}^{2}}{n_{i-1}}\right)$

Let us define estimators $T_{1},\dots,T_{n+1}$ as

$\displaystyle T_{1}^{2}$ $\displaystyle=$ $\displaystyle S_{1}^{2}$
$\displaystyle\forall i,1<i\leq n+1,\,\,T_{i}^{2}$ $\displaystyle=$
$\displaystyle S_{i}^{2}-\frac{S_{i-1}^{2}}{n_{i-1}}$

While $T_{i}^{2}$ are unbiased estimators of the unknown variances, they have
the issue that they may become negative for $i>1$. In our case we need the
estimators only to plan the experiment. Hence, we can simply iteratively
remove levels of the experiment where $T_{i}^{2}$ would be negative or zero.
Maximum likelihood (ML) and restricted maximum likelihood (REML) estimators
for one-way model with normality assumption can be found in mcculloch. Both
types have alternatives for $T_{i}^{2}>0$ and “$T_{i}^{2}<0$”, and both types
are biased. Estimators for a broad range of random and fixed effects models
can be found in [Searle et al. (1992)].

In the rest of this section, we will calculate the expectations of
$S^{2}_{i}$. We start with the lemmas we will use.

###### Lemma 7.8

Let $\overline{z}=\frac{1}{n}\sum_{i=1}^{n}z_{i}$. Then
$\sum_{i=1}^{n}\left(z_{i}-\overline{z}\right)^{2}=\left(\sum_{i=1}^{n}z^{2}_{i}\right)-n\left(\overline{z}\right)^{2}.$

###### Lemma 7.9

Let $z_{i}$, $i=1,..,n$ be random variables with constant expectation,
$Ez_{i}=Ez_{1}<\infty$. The following holds given the variances and
expectations exist:

$E\left(\sum_{i=1}^{n}\left(z_{i}-\overline{z}\right)^{2}\right)=\left(\sum_{i=1}^{n}\mathop{var}\nolimits(z_{i})\right)-n\mathop{var}\nolimits(\overline{z}).$

_Proof._ From Lemma 7.8 and the linearity of expectation we have that

$E\left(\sum_{i=1}^{n}\left(z_{i}-\overline{z}\right)^{2}\right)=^{(\ref{ll1})}E\left(\left(\sum_{i=1}^{n}z^{2}_{i}\right)-n\left(\overline{z}\right)^{2}\right)=\left(\sum_{i=1}^{n}E\left(z^{2}_{i}\right)\right)-nE\left(\left(\overline{z}\right)^{2}\right).$

From the well known property of variance,
$\mathop{var}\nolimits(X)=E(X^{2})-(EX)^{2}$, we have
$E(X^{2})=\mathop{var}\nolimits(X)+(EX)^{2}$, and hence

$\displaystyle\left(\sum_{i=1}^{n}E\left(z^{2}_{i}\right)\right)-nE\left(\left(\overline{z}\right)^{2}\right)$
$\displaystyle=$
$\displaystyle\left(\sum_{i=1}^{n}\mathop{var}\nolimits(z_{i})+(Ez_{i})^{2}\right)-n\left(\mathop{var}\nolimits\left(\overline{z}\right)+\left(E\overline{z}\right)^{2}\right)$
$\displaystyle=$
$\displaystyle\left(\sum_{i=1}^{n}\mathop{var}\nolimits(z_{i})\right)-n\mathop{var}\nolimits(\overline{z})\qed$

###### Lemma 7.10

Let $z_{i}$ , $i=1,..,n$ be random variables with constant variance,
$\mathop{var}\nolimits(z_{i})=\mathop{var}\nolimits(z_{1})<\infty$, and
constant covariance,
$\mathop{cov}\nolimits(z_{i},z_{j})=\mathop{cov}\nolimits(z_{1},z_{2})<\infty$
for $i\neq j$. Then

$\mathop{var}\nolimits\left(\sum_{i=1}^{n}z_{i}\right)=n\mathop{var}\nolimits(z_{1})+(n^{2}-n)\mathop{cov}\nolimits(z_{1},z_{2})$

and

$E\left(\frac{1}{n-1}\sum_{i=1}^{n}\left(z_{i}-\overline{z}\right)^{2}\right)=\mathop{var}\nolimits(z_{1})-\mathop{cov}\nolimits(z_{1},z_{2})$

Note that the first part of the lemma implies that

$\mathop{var}\nolimits(\overline{z})=\frac{\mathop{var}\nolimits(z_{1})}{n}+\frac{n-1}{n}\mathop{cov}\nolimits(z_{1},z_{2})$

_Proof._ The first part follows immediately from the property of the variance
of a sum of random variables:

$\mathop{var}\nolimits\left(\sum_{i=1}^{n}z_{i}\right)=\left(\sum_{i=1}^{n}\mathop{var}\nolimits(z_{i})\right)+2\sum_{i<j}^{n}\mathop{cov}\nolimits(z_{i},z_{j})=n\mathop{var}\nolimits(z_{1})+(n^{2}-n)\mathop{cov}\nolimits(z_{1},z_{2}).$

The second part follows from Lemma 7.9 and from the first part:

$\displaystyle
E\left(\frac{1}{n-1}\sum_{i=1}^{n}\left(z_{i}-\overline{z}\right)^{2}\right)$
$\displaystyle=$
${}^{(L\ref{ll2})}\frac{1}{n-1}\left(\sum_{i=1}^{n}\mathop{var}\nolimits(z_{i})\right)-\frac{n}{n-1}\mathop{var}\nolimits(\overline{z})$
$\displaystyle=$
$\displaystyle\frac{n}{n-1}\mathop{var}\nolimits(z_{1})-\frac{1}{n(n-1)}\mathop{var}\nolimits\left(\sum_{i=1}^{n}z_{i}\right)$
$\displaystyle=$
$\displaystyle\frac{n}{n-1}\mathop{var}\nolimits(z_{1})-\frac{1}{n-1}\mathop{var}\nolimits(z_{1})-\frac{n(n-1)}{n(n-1)}\mathop{cov}\nolimits(z_{1},z_{2})$
$\displaystyle=$
$\displaystyle\mathop{var}\nolimits(z_{1})-\mathop{cov}\nolimits(z_{1},z_{2})\qed$

The plan is to use Lemma 7.10 to calculate, in turn, the expectation of each
variance estimator $S^{2}_{i}$. We will use it for all $1\leq i\leq n+1$,
although $i=1$ and $i=n+1$ are special cases which could be solved in an
easier way. In each step, we will calculate as a side effect the variance of
$\sum_{i=1}^{n}z_{i}$, which will be handy for calculation of the term
“$\mathop{var}\nolimits(z_{1})$” in the next step. For this plan, we also need
the covariances of “$z_{1},z_{2}$”. We will first calculate these covariances
using the following three lemmas.

###### Lemma 7.11

For random variables $X_{i}$, $i=1,..,m$ and $Y_{j}$, $j=1.,,.n$ it follows
from a known property of covariance that

$\mathop{cov}\nolimits\left(\overline{X},\overline{Y}\right)=\frac{1}{mn}\sum_{i=1}^{m}\sum_{j=1}^{n}\mathop{cov}\nolimits(X_{i},Y_{j}).$

Note that if in addition $\mathop{cov}\nolimits(X_{i},Y_{j})$ is a constant,
then
$\mathop{cov}\nolimits(\overline{X},\overline{Y})=\mathop{cov}\nolimits(X_{1},Y_{1})$.

######  Lemma 7.12 (Law of Total Covariance)

For random variables $X$,$Y$,$Z$:
$\mathop{cov}\nolimits(X,Y)=E(\mathop{cov}\nolimits(X,Y|Z))+\mathop{cov}\nolimits(E(X|Z),E(Y|Z))$.

This is a known property which can be shown using the law of total expectation
and the definition of covariance.

######  Lemma 7.13 (Tower Property of Conditional Expectation)

For random variables $X$,$G$,$H$: $E(X|G)=E(E(X|G,H)|G).$

The proof is simple, but depends on “low-level” formal definition of
conditional expectation using $\sigma-$ algebras that we avoid here. Details
can be found in [Steele (2001)].

The covariance at the lowest level (measurements in one execution) can be
calculated as follows. For $k\neq l$

$\displaystyle\mathop{cov}\nolimits\left(Y_{j_{n+1}\dots
j_{2}k},Y_{j_{n+1}\dots j_{2}l}\right)$ $\displaystyle=$
${}^{(L\ref{lltc})}E\left(\mathop{cov}\nolimits\left(Y_{j_{n+1}\dots
j_{2}k},Y_{j_{n+1}\dots j_{2}l}|\left(\mu_{1}\right)_{j_{n+1}\dots
j_{2}}\right)\right)+$ $\displaystyle+$