##  0.1 RQ 1: Does influence testing behavior?

### Test Creation

intends to incentivize writing unit tests, and appears to achieve this: There
was only one participant in the group who did not write JUnit tests at all,
but six participants of the group which is significant ($p=0.039$) according
to the exact Fisher test [bower2003use]. In contrast, participants of the
group used the main method for manual testing more frequently: Only 2/17
participants in the group used the main method for testing at some point,
compared to 7/16 of the group. Notifications show some effect, but less than
the achievements: 2/8 participants wrote no unit tests, and 2/8 used main
testing which can be seen as an indication that the participants of the group
are slightly more motivated than the group, but less than the group. Focusing
on achievements seems to lead to the most motivation [8370133] to test: All
participants of the group wrote unit tests, and only one of them used main
testing, which suggests that achievements make developers test more than
needed to complete the study task.

[t]0.33  [t]0.33  [t]0.33  [t]0.33  [t]0.33  [t]0.33

Figure 1: Number of tests created over time Figure 2: Number of test
executions over time Figure 3: Number of test executions with coverage over
time Figure 4: Number of debug uses over time Figure 5: Number of passed tests
over time Figure 6: Number of levels over time Figure 7: Differences between ,
, and groups over time (see sec:procedure)

fig:teststime shows the JUnit tests written over time: Both the and groups
start early to write tests; the first participant of the group wrote a test
after a minute, and for the group after two minutes. From 22 minutes onwards
the number of tests in the group increases much faster than for the other
groups (except ), and from minute 30 on the improvement over the group is
significant (the confidence intervals do not overlap.) The group initially
behaved similarly to the group, but stagnated towards the end and stayed
between both the and groups, confirming that notifications contribute to, but
do not fully explain, the improvements. As expected, the group started sooner
and wrote significantly more tests than any other group between minutes 31 and
50, at which point the group caught up.

### Test Executions

[t]0.245  justification=centering [t]0.245  justification=centering [t]0.245
justification=centering [t]0.245  justification=centering [t]0.245
justification=centering [t]0.245  justification=centering [t]0.245
justification=centering [t]0.245  justification=centering

Figure 8: Number of test executions
during the experiment Figure 9: Number of test executions
with coverage in IntelliJ Figure 10: Number of times the IntelliJ
debug mode was used Figure 11: Number of tests written
by the participants Figure 12: Line coverage measured
by JaCoCo of the final test suites Figure 13: Branch coverage measured
by JaCoCo of the final test suites Figure 14: Mutation score measured
by PITest of the final test suites Figure 15: Number of failed tests accord-
ing to the golden test suite Figure 16: Differences between , , and groups

fig:testexecutions shows the number of test executions, with a mean of 22.8
for the group and 10.5 for the control group. The ($p<0.01$), ($p<0.01$) and
($p<0.01$) groups all executed significantly more tests than the group. In the
group, two values with more than 1.5 million and 250,000 test executions
(achieved using IntelliJ’s feature of repeating test executions a configurable
number of times) were omitted from the graph to make it more readable. These
two data points demonstrate the possible exploitation of when the participants
try to get as many achievements as possible, which may distract the developers
from their actual goal. Consequently, it seems important not to make the use
of mandatory for developers, but rather to let them decide on their own how to
use achievements as was done in the group. In addition, fig:testexecutionstime
illustrates the test executions over time together with the 84.6 % confidence
intervals. A significant difference between the and groups can be observed
starting from minute 27. The and groups behave similarly to each other, while
the group shows a significant difference to the group in the end. Overall,
participants with run their tests more often.

### Coverage Measurement

fig:testexecutionscoverage compares the number of times tests were executed
with coverage collection activated: Only one participant in the group used
this functionality of IntelliJ 10 times, resulting in an average number of
coverage executions of 0.6 overall participants. On the other hand, the group
used the coverage report with a mean of 4.7 executions per participant.
fig:testexecutionscoveragetime shows the coverage executions over time,
demonstrating that the first participant of the group used the coverage report
after 13 minutes, whereas in the group only after 35 minutes. This also means
that the difference is significant from the beginning on ($p<0.01$), when the
first participant executed their tests with coverage collection activated. The
group behaves similarly to the group, even though there is a dedicated
notification that suggests running tests with coverage. However, this is only
shown if tests are run without coverage. The group as expected measures
coverage the most during the last third except the very end. Overall, the
participants of the and groups seem to care more about coverage.

### Debug Executions

fig:debugmode shows how often tests were executed using the debug mode: On
average participants of the group used debug executions 2.6 times, and of the
group 1.8 times, although the difference is not significant ($p=0.82$). This
can also be seen in fig:debugtime, since the confidence intervals of all
groups except for the group overlap. We can see a significant difference
halfway through the experiment, but the overall difference is not significant
($p=0.074$). We assume that debugging will be more common during maintenance
and bug fixing, whereas our experiment task was to write new code.
Furthermore, since the class under test is quite simple, the participants may
not have had to use the debug mode to locate bugs found by failed tests.

RQ 1influences the testing behavior of the participants: They use JUnit tests
rather than main testing, create more tests, run their tests significantly
more often, and use coverage reports more frequently.

##  0.2 RQ 2: Does influence resulting test suites?

### Number of tests

fig:testnumber compares the resulting test suites in terms of the number of
JUnit tests they consist of. The participants of the group wrote a mean of
11.5 tests, while the group wrote 4.3 tests and the group 12.5. According to
the exact Wilcoxon-Mann-Whitney test, the ($p<0.01$) and ($p=0.036$) groups
wrote significantly more tests than the group. One can also see that the group
has more tests in their test suites (mean 7.38) than the group (mean 4.31),
which again confirms that notifications influence developers to write more
tests.

### Code coverage

fig:linecoverage compares the line coverage between the four groups, showing
that the ($p<0.01$) and ($p<0.01$) groups have significantly higher line
coverage than the group. On average, the and groups achieve 89 % and 48 % line
coverage, respectively. In the group, there is one outlier without any
coverage; this participant did not write any unit tests but only tested
manually. When asked about this peculiar behavior, the participant stated to
have rushed to complete and leave the study as quickly as possible. Results
are similar for branch coverage (fig:branchcoverage), which also shows a
significant difference between the control, ($p<0.01$) and ($p<0.01$) groups,
with a mean branch coverage of 85 % and 41 % for the and groups. The group has
no significantly different line ($p=0.44$) or branch coverage ($p=0.64$)
compared to the group, which is plausible since notifications only suggest
measuring coverage, but not maximizing it. Overall, the significantly higher
coverage for both the and groups confirms the influence of gamification on
test suite quality as measured with code coverage.

### Mutation score

The mean mutation score (fig:mutationscore) of the group is 84 %, while the
group achieves a mean score of 44 %, which is significantly lower ($p<0.01$).
The group shows no significant differences ($p=0.76$) compared to the group,
again confirming the influence of gamification; the group even achieves
slightly higher mutation scores ($p=0.01$) than the group.

RQ 2has a significant influence on the resulting test suites, increasing the
number of tests as well as the code coverage and mutation scores.

##  0.3 RQ 3: Do achievement levels reflect differences in test suites and
activities?

[b]0.48  justification=centering [b]0.48  lll Variable & $\bm{r^{2}}$ p-value
Line Coverage 0.21 ¡ 0.001
Branch Coverage 0.23 ¡ 0.001
Mutation Score 0.22 ¡ 0.001
Number of Tests 0.39 ¡ 0.001
justification=centering

Figure 17: The number of achievement levels reached Figure 18: The Pearson
correlations based on the number of levels Figure 19: The number of
achievement levels and their Pearson correlations with different test suite
metrics

RQ1 demonstrated that participants of the group engage most with testing
activities, followed by the group, then the group, and finally the group, who
motivated the least with testing. fig:levels compares how many achievement
levels the participants of the four groups reached, and confirms this ranking
of the groups with a mean of 2.6 for the group, 4 for the group, 6.1 for the
group, and 15.4 for the group; fig:levelstime shows that this holds throughout
the entire duration of the experiment. All groups reached significantly more
levels than the group ($p<0.01$, $p<0.01$, $p<0.01$). Thus, developers who
have been awarded achievements indeed are more committed to testing. To see
whether developers with more achievements and higher achievement levels also
produce overall better test suites, tab:corr shows the Pearson rank
correlations between test suite metrics (RQ2) and achievement levels. There is
a significant moderate correlation between the number of tests and
achievements and a significant weak positive correlation for all other
metrics. This is a strong affirmation of the gamification approach, as it
shows that acquiring achievements is indeed related to producing better test
suites.

RQ 3Higher achievement levels are indicators of better motivation with testing
and better resulting test suites.

##  0.4 RQ 4: Does influence the functionality of resulting code?

fig:passedteststime shows the number of passed tests of the golden test suite
during the experiment, where notably the and groups exhibit passing tests
earlier compared to the other groups. Particularly, the group reached a
plateau where most tests passed at minute 43, indicating that solving
achievements leads to faster identification and resolution of program faults
as well as faster implementation of new functionality.

At the end of the experiment, the average number of failing tests of the
golden test suite (fig:testsfailed) is 0.9 for the group, and 1.6 for the
group, with no statistically significant ($p=0.45$) difference. The and groups
also show similar results to the other groups, although participants in the
notifications group have the highest average number of failing tests (2.25).
The lack of significant differences can be attributed to several factors: (1)
the experiment duration was long enough for participants to complete the
implementation even with inferior testing behavior, (2) the golden test suite
consisted of only six test cases, and (3) the two methods to be implemented
allowed for a limited number of distinct possible bugs.

RQ 4The achievements of lead to earlier implementation of functionality,
although we observed no significant differences in functionality at the end of
the experiment.

##  0.5 RQ 5: Does influence the developer experience?

Figure 20: General survey responses Figure 21: Answers regarding the plugin

According to the exit survey (fig:exitall), the target class was well chosen.
All groups claim that they had enough time to implement and test the class and
that it was easy to understand and implement. Notable and significant
differences can be observed on whether participants believe to have tested
their code actively: All participants of the group claim to have actively
tested their code, while only 55 % of the group did this; this difference is
significant ($p=0.03$). Participants of the group are also more confident in
that they produced a good test suite than participants of the group, although
the difference is not significant ($p=0.16$). Participants of the group, who
invested more time into testing, believe even more often to have produced a
good test suite, significantly so when compared to the ($p=0.01$) and
($p=0.025$) groups, but not against the group ($p=0.14$). Throughout the
survey, we generally observe favorable responses by the group, although not
statistically significant for the other cases. Interestingly, members of the
group are least confident in that they produced a good implementation, which
indicates a disadvantage of their focus on the achievements rather than the
programming task.

Considering the experiences of the and groups, who both used achievements
(fig:exitplugin), more than half the participants perceived that improved
their programming behavior, and more than 60 % were driven by notifications,
including both encouraging and progress notifications. Most of the
participants understood the individual achievements and how to achieve them.
Nearly 60 % would like to use in their everyday work. More than 75 % of the
participants claim to have been motivated to improve testing, which indicates
that the plugin served its purpose.

RQ 5Users of perceived the task as easier, are more confident about their
results, and confirm that the achievements motivated them to apply more
testing.


# Sequential inductive prediction intervals

Benny Avelin  Combient Mix (Silo AI)
Mäster Samuelsgatan 56, 111 21 Stockholm, Sweden<EMAIL_ADDRESS>
(Date: August 31, 2024)

###### Abstract.

In this paper we explore the concept of sequential inductive prediction
intervals using theory from sequential testing. We furthermore introduce a
3-parameter PAC definition of prediction intervals that allows us via
simulation to achieve almost sharp bounds with high probability.

##  1\. Coverage and prediction intervals

If we let $X$ correspond to a prediction, then a $1-\alpha$ prediction
interval is in the simplest case nothing other than a set $C$ such that the
measure of that set is exactly $1-\alpha$, i.e.

$\displaystyle\mathbb{P}(X\in C)=1-\alpha.$

The set $C$ is unknown, and we would like to construct it based on
observations. Lets us consider a sequence of i.i.d. random variables
$X_{1},\ldots,X_{n}$ then say we have a way to determine
$C(X_{1},\ldots,X_{n-1})$ as a function of the first $n-1$ observations, then
the requirement of being a prediction interval is usually

$\displaystyle\mathbb{P}(X_{n}\in C(X_{1},\ldots,X_{n-1}))=1-\alpha.$ (1.1)

In this simple setting we could guarantee exactly the $1-\alpha$ level if we
take $C=[\widehat{q}_{\alpha/2},\widehat{q}_{1-\alpha/2}]$, where
$\widehat{q}_{\delta}$ is the empirical $\delta$ quantile, with this choice
and if $(n-1)\alpha/2\in\mathbb{N}$ then 1.1 holds. This is called a marginal
coverage guarantee.

In the context of regression the prediction interval statement becomes for a
sequence of i.i.d. observations $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$

$\displaystyle\mathbb{P}(Y_{n}\in
C((X_{1},Y_{1}),\ldots,(X_{n-1},Y_{n-1}),X_{n}))=1-\alpha.$

To achieve this, assume we have been given some data and say that we have
found a function $f(X)$ approximating the regression function
$\mathbb{E}[Y\mid X]$, then if we define the residual $Z_{i}=Y-f(X_{i})$ we
can form the prediction interval as $f(X_{n})+C$ where $C$ is the same
quantile based interval as in the previous example, but this time based on
$Z_{1},\ldots,Z_{n-1}$.

In the framework of conformal prediction [7, 4, 3], one would construct a
well-designed non-conformity score. The non-conformity score is in the formal
mathematical sense a function from the bag of previous observations and the
new observation to a real number. A simple example would be to take what
empirical quantile the new observation belongs to if we add it to the bag of
previous observations.

An extension of the idea of marginal guarantee is that of a sequential
guarantee. Consider again a sequence of random variables $Z_{1},\ldots,Z_{n}$
and let’s say we have a sequence of prediction intervals
$C_{i}=C(Z_{1},\ldots,Z_{i})$, the transductive prediction problem ([1]) is to
guarantee that

$\displaystyle\mathbb{P}(Z_{i}\in C_{i-1})=\alpha$

for all $i$ and that the sequence of errors

$\displaystyle E_{i}=\mathbb{I}_{Z_{i}\in C_{i-1}}$

is an i.i.d. sequence of random variables. With a simple randomization trick
(to avoid the discreteness of the problem) we can essentially use the
empirical quantiles here, see [6].

The interesting thing with respect to this paper, is that we can make this
guarantee, but only for the next sample. If we instead would like to guarantee
that

$\displaystyle\mathbb{P}(\mathbb{P}(Z_{i}\in C(Z_{1},\ldots,Z_{i-1})\mid
Z_{1},\ldots,Z_{i-1})\geq 1-\alpha)\geq 1-\delta$ (1.2)

this is denoted as the conditional coverage and is phrased as a 2-parameter
PAC (Probably Approximately Correct, [5]) problem, essentially introduced in
[6]. What is interesting here, is that if we fix the value of $i$ then this
denotes that with probability greater than $1-\delta$ we can, given the
calibration data $Z_{1},\ldots,Z_{i-1}$, construct an interval that has at
least the correct coverage probability. We interpret this to mean that the
prediction interval holds for an arbitrary number of future observations,
which is the so-called inductive setting (see [6]). That is, we want to
construct an interval which with high probability can be trusted.

Let us consider a scenario that highlights the differences between the
transductive and the inductive form of prediction intervals. In the first we
create a sequence of intervals for which the next sample will land with a
guaranteed probability, which is only true for one observation. This would be
useful in the online setting when you are trying to bet on the next upcoming
value in the context of say finance. This might seem subtle, but the fact that
the consecutive errors $E_{i}$ are i.i.d. does not guarantee that the
probability of the next sample being in the interval given the previous
observations is anything exact or above $1-\alpha$, if we have many samples it
will however be close with high probability.

In the second setting, i.e. when we talk about the inductive form, we are
trying to create an interval that we can trust for all eternity, or at least
for a fixed number of new samples which was considered in [6]. This is
different because we have instead phrased it as a decision problem. Namely, if
the procedure we design produces a prediction interval, then if we say that
this interval has a coverage of at least $1-\alpha$ we would be right
$1-\delta$ of the time. This is very similar to producing a confidence
interval, where the procedure of computing the interval and deciding that the
true parameter lies in the computed interval, we would be right $1-\delta$ of
the time.

In terms of selling a product where we either train the model once and sell it
or are intermittently updating the model with new data, we cannot control how
many times the customer will use the model for predictions. This means that
the transductive setting is not useful and this second style of producing
inductive intervals makes more sense.

The statement in 1.2 considers that we have chosen a sample size $i$ and
construct our prediction interval based on that sample size. This does not
leverage the fact that as you see more testing points you could tighten the
prediction interval to be closer to the true quantiles of the distribution.
This is intimately close to the idea of sequential testing, and by the duality
of tests and confidence intervals, the construction of confidence sequences
for parameters.

A confidence sequence is something that satisfies

$\displaystyle\mathbb{P}(\forall n:\text{parameter in set defined by
}(X_{1},\ldots,X_{n}))\geq 1-\delta.$

For confidence intervals this would be a sequence of confidence intervals
constructed such that the probability that our parameter is in the set is at
least $1-\delta$ for all $n$ uniformly. This means that we are allowed to peek
at our test at all times, the price we pay for this peeking is that the
intervals are necessarily bigger.

The same idea can be phrased for prediction, in this case we call $C_{n}$ the
prediction intervals and $X$ is then the true value

$\displaystyle\mathbb{P}\left(\forall n:\mathbb{P}(X\in C_{n})\geq
1-\alpha\right)\geq 1-\delta.$

That is, our sequence of prediction intervals will have coverage at least
$1-\alpha$ with probability greater than $1-\delta$. This statement would be
moot if we did not require $C_{n}$ to get tighter as $n\to\infty$,
specifically we would like $\mathbb{P}(X\in C_{n})\to 1-\alpha$ as
$n\to\infty$ with the best possible rate.

###  1.1. Confidence sequence for the quantile and prediction intervals

We will base our discussion on a recent result [2].

Let $X$ be a real valued random variable. Write $F(x):=\mathbb{P}(X\leq x)$ as
the distribution function, define the empirical distribution of an i.i.d.
sequence $X_{1},\ldots,X_{n}$ sampled from $F$ as
$\widehat{F}_{n}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}_{X_{i}\leq x}$. Define
the upper and lower quantile functions for $p\in[0,1]$,

$\displaystyle Q(p)$ $\displaystyle:=\sup\\{x:F(x)\leq p\\}$ $\displaystyle
Q^{-}(p)$ $\displaystyle:=\sup\\{x:F(x)<p\\},$

respectively. Denote $\widehat{Q}_{n}(p)$ and $\widehat{Q}_{n}^{-}(p)$ the
upper and lower quantile functions w.r.t. $\widehat{F}_{n}$.

Let $Q$ be the upper quantile function, let $\widehat{F}$ be the empirical
distribution function and $\widehat{Q}$ be the empirical quantile function.
Superscript $-$ denotes the corresponding lower quantile functions.

######  Theorem 1 ([2]).

Let $Z_{1},\ldots\overset{IID}{\sim}F$, and let $Q$ be the corresponding
quantile function. Then, for a fixed $\delta\in(0,1)$ the following holds

$\displaystyle\mathbb{P}(\forall i\in\mathbb{N},\forall
p\in(0,1):Q(p)\in[\widehat{Q}_{i}(p-g_{i}),\widehat{Q}_{i}^{-}(p+g_{i})])\geq
1-\delta,$

where

$\displaystyle
g_{i}=0.85\sqrt{\frac{1}{i}[\log(\log(ei))+0.8\log(1612/\delta)]}$

The good thing about the above is that the probability is over all $i$ and
over all $p$-quantiles simultaneously (this is because it is a consequence of
a DKW style inequality, i.e. a bound on the $L^{\infty}$ norm of
$F-\widehat{F}$). We can leverage this fact to get stronger control and obtain
also an upper bound on the coverage in 1.2. We will start with the lower bound
as it is easier. Choose a confidence level $\alpha$ and a failure rate
$\delta$ and note that Theorem 1 implies

$\displaystyle\mathbb{P}\left(\forall
i,\alpha\in(0,1):\widehat{Q}_{i}\left(\frac{\alpha}{2}-g_{i}\right)\leq
Q\left(\frac{\alpha}{2}\right);Q\left(1-\frac{\alpha}{2}\right)\leq\widehat{Q}_{i}^{-}\left(1-\frac{\alpha}{2}+g_{i}\right)\right)\geq
1-\delta,$ (1.3)

which implies that if we define the two intervals

$\displaystyle C^{l,\alpha}_{i}$
$\displaystyle=[\widehat{Q}_{i}(\alpha/2-g_{i}),\widehat{Q}_{i}^{-}(1-\alpha/2+g_{i})]$
$\displaystyle C^{u,\alpha}_{i}$
$\displaystyle=[\widehat{Q}_{i}^{-}(\alpha/2+g_{i}),\widehat{Q}_{i}(1-\alpha/2-g_{i})]$

then

$\displaystyle\mathbb{P}(\forall i,\alpha\in(0,1):\mathbb{P}(Z\in
C^{l,\alpha}_{i})\geq 1-\alpha)\geq 1-\delta,$ (1.4)

the upper bound can be obtained similarly, saying that

$\displaystyle\mathbb{P}(\forall i,\alpha\in(0,1):\mathbb{P}(Z\in
C^{u,\alpha}_{i})\leq 1-\alpha)\geq 1-\delta.$

Note that from 1.3 tells us that the upper and lower bound holds
simultaneously with probability greater than $1-\delta$, so we don’t have to
use the union bound to combine them. This provides good control over the
coverage probability, the issue is that we obtain two sets, one bigger with
the lower bound and one smaller with the upper bound. We can rephrase this for
a single interval, where we get an upper and lower bound, however we have to
choose if we want a deterministic lower bound and a random upper bound, or we
can choose a deterministic upper bound and a random lower bound.

Let us show how we do this, we begin by defining

$\displaystyle\alpha^{l}_{i}:=\sup\\{\widetilde{\alpha}\in[0,1]:\widehat{Q}^{-}_{i}(\widetilde{\alpha}/2+g_{i})<\widehat{Q}_{i}(\alpha/2-g_{i})\\}$

where the supremum is $0$ if the set is empty. We also define

$\displaystyle\alpha^{u}_{i}:=\sup\\{\widetilde{\alpha}\in[0,1]:\widehat{Q}_{i}(1-\widetilde{\alpha}/2-g_{i})>\widehat{Q}_{i}(1-\alpha/2+g_{i})\\}$

where the supremum is $0$ if the set is empty. Then define

$\displaystyle\alpha^{\ast}_{i}(\alpha)=\min(\alpha^{l}_{i}(\alpha),\alpha^{u}_{i}(\alpha)).$

From this construction we get that $C_{i}^{l,\alpha}\subset
C_{i}^{u,\alpha^{\ast}_{i}(\alpha)}$, which together with the fact that 1.3
holds for all $\alpha$, we get the following theorem

######  Theorem 2.

We have for the observable value $\alpha^{\ast}_{i}$ that

$\displaystyle\mathbb{P}(\forall
i,\alpha\in(0,1):1-\alpha^{\ast}_{i}(\alpha)\geq\mathbb{P}(Z\in
C^{l,\alpha}_{i})\geq 1-\alpha)\geq 1-\delta.$

This is a very useful estimate, as we can guarantee the lower bound with
probability $1-\delta$, the upper bound we can determine using the data.

It is also clear that studying the definition, we see that
$\alpha_{i}^{\ast}\to\alpha$ as $i\to\infty$ since the empirical quantiles
converge a.s. to the true quantile and $g_{i}\to 0$.

###  1.2. Conformal prediction and confidence sequences

Let us now show how to use confidence sequences in the context of an inductive
prediction setting using ideas from conformal prediction. Namely, assume that
we have a score function $S(x,y)$ which measures the error we make. Typical
choice is

$\displaystyle S(x,y)=\frac{|f(x)-y|}{\sigma(x)}$

where $f$ is the function making the prediction and $\sigma$ could for
instance be the square root of an estimator of the conditional variance. For a
specific example see the kernel variance estimator in Appendix B.

Now, for an i.i.d sequence of calibration points we will create the random
variable $Z_{i}=S(X_{i},Y_{i})$ and for this we can create the confidence
sequence

$\displaystyle\mathbb{P}(\forall i:\mathbb{P}(Z\in C_{i}(\alpha)\mid
Z_{1},\ldots,Z_{i})\geq 1-\alpha)\geq 1-\delta.$ (1.5)

This confidence sequence $C_{i}(\alpha)$ can be created using for instance
Theorem 2. Now note that $C_{i}(\alpha)=[L_{i}(\alpha),U_{i}(\alpha)]$ is an
interval, therefore if we define the score function to be the signed version
of $S$ above, i.e.

$\displaystyle S(x,y)=\frac{y-f(x)}{\sigma(x)}$

then we can define the scaled and shifted interval

$\displaystyle
D_{i}(\alpha)=[f(X)+\sigma(X)L_{i}(\alpha),f(X)+\sigma(X)U_{i}(\alpha)],$

and 1.5 becomes

$\displaystyle\mathbb{P}(\forall i:\mathbb{P}(Y\in
D_{i}(\alpha)\mid(X_{1},Y_{1}),\ldots,(X_{i},Y_{i}))\geq 1-\alpha)\geq
1-\delta$

The fact that our score function is signed allows us to treat problems in
where the conditional distribution of $Y\mid X$ is skew, see Fig. 6.

###  1.3. IID bounded residual

Let us see another application of these intervals in the context where we
cannot guarantee identical distribution of the residual. Assume that we have a
heteroscedastic problem, i.e.

$\displaystyle Y=f(X)+Z$

and assume that $Z$ is mean zero and independent of $X$. Then $f$ is the
regression function. Assume that we have used data to estimate $f$, call the
estimate $\widehat{f}$, and assume that we managed to achieve
$\|f-\widehat{f}\|_{L^{\infty}}<\varepsilon$, which can be achieved with high
probability if for instance $f$ is Lipschitz continuous. We are interested in
this situation to create prediction intervals around $Y$ using $\widehat{f}$,
but where we take deterministic $x$.

Consider now a sequence of given points $x_{1},\ldots,x_{n}$, then
$Y_{i}=f(x_{i})+Z_{i}$. The residual on the other hand is
$R_{i}=Y_{i}-\widehat{f}(x_{i})$, and since $\widehat{f}\neq f$ this is only a
sequence of independent random variables, it is not identically distributed.
We wish to create a prediction interval that contains $R_{i}$ with high
probability. Therefore, we introduce the following concept

######  Definition 1.1.

Let $X_{1},\ldots$ be an independent sequence of random variables and assume
that there exists constants $c_{1},\delta_{1},c_{2},\delta_{2}$ and an IID
sequence of random variables $Y_{1},\ldots$ (a coupling) such that

$\displaystyle c_{1}Y_{i}+\delta_{1}\leq X_{i}\leq c_{2}Y_{i}+\delta_{2}.$

Then we say that $X_{1},\ldots$ is an _IID bounded sequence_.

The above definition allows us to compare the distribution functions

######  Lemma 1.2.

Let $X_{1},\ldots$ be an IID bounded sequence with constants $c_{1}$,
$\delta_{1}$, $c_{2}$, $\delta_{2}$ and IID bounding sequence $Y_{1},\ldots$.
Then

$\displaystyle F_{X_{i}}(c_{1}x+\delta_{1})\leq F_{Y_{i}}(x)\leq
F_{X_{i}}(c_{2}x+\delta_{2})$

Define the empirical distribution functions as
$\widehat{F}_{Y_{n}}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}_{Y_{i}\leq x}$ and
$\widehat{F}_{X_{n}}(x)=\frac{1}{n}\sum_{i=1}^{n}\mathbb{I}_{X_{i}\leq x}$,
then

$\displaystyle\widehat{F}_{X_{n}}(c_{1}x+\delta_{1})\leq\widehat{F}_{Y_{n}}(x)\leq\widehat{F}_{X_{n}}(c_{2}x+\delta_{2})$

###### Proof.

It is immediate from the definition of IID bounded sequences, that

$\displaystyle F_{Y_{i}}(x)=\mathbb{P}(Y_{i}\leq
x)\leq\mathbb{P}((X_{i}-\delta_{2})/c_{2}\leq
x)=F_{X_{i}}(c_{2}x+\delta_{2}).$

The result for the empirical distribution function follows in the same way. ∎

######  Lemma 1.3.

Let the sequence $X_{i}$ be as in Lemma 1.2. Then define the quantile function
$Q_{X_{i}}(p)=\sup\\{x:F_{X_{i}}(x)\leq p\\}$ and the lower quantile function
$Q^{-}_{X_{i}}(p)=\sup\\{x:F_{X_{i}}(x)<p\\}$. We define $\widehat{Q}_{X_{n}}$
and $\widehat{Q}^{-}_{X_{n}}$ similarly. Then, for $\alpha\in(0,1)$ with
$n\alpha\in\mathbb{N}$, it holds

$\displaystyle
c_{1}\widehat{Q}_{Y_{n}}(\alpha)+\delta_{1}\leq\widehat{Q}_{X_{n}}(\alpha)\leq
c_{2}\widehat{Q}_{Y_{n}}(\alpha)+\delta_{2},$

and

$\displaystyle
c_{1}\widehat{Q}^{-}_{Y_{n}}(\alpha)+\delta_{1}\leq\widehat{Q}^{-}_{X_{n}}(\alpha)\leq
c_{2}\widehat{Q}^{-}_{Y_{n}}(\alpha)+\delta_{2}.$

###### Proof.

Let $n\alpha\in\mathbb{N}$ and $\alpha\in(0,1)$, then by Lemma 1.2

$\displaystyle\alpha+1/n=\widehat{F}_{Y_{n}}(\widehat{Q}_{Y_{n}}(\alpha))\leq\widehat{F}_{X_{n}}(c_{2}\widehat{Q}_{Y_{n}}(\alpha)+\delta_{2})$

this implies that
$c_{2}\widehat{Q}_{Y_{n}}(\alpha)+\delta_{2}\not\in\\{x:\widehat{F}_{X_{n}}(x)\leq\alpha\\}$,
and as such we obtain

$\displaystyle\widehat{Q}_{X_{n}}(\alpha)\leq
c_{2}\widehat{Q}_{Y_{n}}(\alpha)+\delta_{2}.$

On the other hand we have

$\displaystyle\alpha+1/n=\widehat{F}_{X_{n}}(\widehat{Q}_{X_{n}}(\alpha))\leq\widehat{F}_{Y_{n}}\left(\frac{\widehat{Q}_{X_{n}}(\alpha)-\delta_{1}}{c_{1}}\right)$

which in the same way implies that

$\displaystyle
c_{1}\widehat{Q}_{Y_{n}}(\alpha)+\delta_{1}\leq\widehat{Q}_{X_{n}}(\alpha).$

Let us now consider the lower quantile function, we first note that

$\displaystyle\alpha=\widehat{F}_{Y_{n}}(\widehat{Q}^{-}_{Y_{n}}(\alpha))\leq\widehat{F}_{X_{n}}(c_{2}\widehat{Q}^{-}_{Y_{n}}(\alpha)+\delta_{2})$

which again implies that
$c_{2}\widehat{Q}^{-}_{Y_{n}}(\alpha)+\delta_{2}\not\in\\{x:\widehat{F}_{X_{n}}(x)<\alpha\\}$
hence

$\displaystyle\widehat{Q}^{-}_{X_{n}}(\alpha)\leq
c_{2}\widehat{Q}^{-}_{Y_{n}}(\alpha)+\delta_{2},$

the lower bound now follows in the same way as above. ∎

Now, for IID bounded sequences we can construct prediction intervals similarly
to IID sequences but at a cost.

######  Theorem 3.

Define the interval

$\displaystyle
C_{n}=\left[\frac{\widehat{Q}_{X_{n}}(\frac{\alpha}{2}-g_{i})-\delta_{2}}{c_{2}},\frac{\widehat{Q}_{Y_{n}}^{-}(1-\frac{\alpha}{2}+g_{i})-\delta_{1}}{c_{1}}\right]$

then

$\displaystyle\mathbb{P}(\mathbb{P}(Y_{n+1}\in C_{n}\mid
X_{1},\ldots,X_{n})\geq 1-\alpha)\geq 1-\delta.$

Correspondingly, by enlarging $C_{n}=(l,r)$ as follows

$\displaystyle\widetilde{C}_{n}=\left(c_{1}l+\delta_{1},c_{2}r+\delta_{2}\right)$

we obtain a similar statement for $X_{n+1}$

$\displaystyle\mathbb{P}(\mathbb{P}(X_{n+1}\in\widetilde{C}_{n}\mid
X_{1},\ldots,X_{n})\geq 1-\alpha)\geq 1-\delta.$

###### Proof.

By Theorem 2 we have the prediction interval

$\displaystyle[\widehat{Q}_{Y_{n}}(\alpha/2-g_{i}),\widehat{Q}_{Y_{n}}^{-}(1-\alpha/2+g_{i})].$

This interval will with probability $\delta$ have coverage greater than
$1-\alpha$. Now from Lemma 1.3 we see that we can cover the prediction
interval above as follows

$\displaystyle\left[\widehat{Q}_{Y_{n}}(\frac{\alpha}{2}-g_{i}),\widehat{Q}_{Y_{n}}^{-}(1-\frac{\alpha}{2}+g_{i})\right]\subset\left[\frac{\widehat{Q}_{X_{n}}(\frac{\alpha}{2}-g_{i})-\delta_{2}}{c_{2}},\frac{\widehat{Q}_{Y_{n}}^{-}(1-\frac{\alpha}{2}+g_{i})-\delta_{1}}{c_{1}}\right].$

To get the last statement, we simply apply Lemma 1.3 again. ∎

######  Example 1.4.

For the day index $t$ we sell $X_{t}$ number of items, and
$\mathbb{E}[X_{t}]=\lambda(t)$. We assume that we have used previous years
data to estimate the rate function $\lambda(t)$ which we can with high
probability do with uniform error $\varepsilon$.

Our observations are thus

$\displaystyle X_{t}=\lambda(t)+Z_{t}$

where $Z_{t}$ are i.i.d.

Our estimated rate function is $\widehat{\lambda}(t)$ and we know that
$|\lambda(t)-\widehat{\lambda}(t)|\leq\varepsilon$. Let us construct the
sequence of errors defined as

$\displaystyle
E_{t}=X_{t}-\widehat{\lambda}(t)=Z_{t}+\lambda(t)-\widehat{\lambda}(t)=:Z_{t}+\varepsilon_{t},$

with $|\varepsilon_{t}|\leq\varepsilon$.

We are now exactly in the setting where we can apply Theorem 3.

##  2\. Sharper bounds

In this section we will use simulation in order to sharpen the bound in
Theorem 1. This is especially the case for smaller values of $t$, as can be
seen in Fig. 1, where we have the correct shape but slightly inflated values.
Instead of trying to come up with a smarter way to estimate the probabilities
involved we can simulate (at least for finite values).

Figure 1. Some trajectories of the quantile process and the time-uniform bound
provided in [2].

To begin, we recall the notation of [2] and define

$\displaystyle S_{t}=\sum_{i=1}^{t}(\mathbb{I}_{X_{i}\leq Q(p)}-p)$ (2.1)

where $X_{i}\sim F$ is an i.i.d. sequence and $Q(p)$ is the $p$ quantile
w.r.t. $F$. We note that by definition $\widehat{F}_{t}(Q(p))=S_{t}/t$, as
such, if we can find a function $-p<f(t,p,\delta)<0$ such that

$\displaystyle\mathbb{P}(\exists t:\,S_{t}/t<f(t,p,\delta))\leq\delta$ (2.2)

it would imply by the definition of the quantile function, that

$\displaystyle\mathbb{P}(\exists
t:\,Q(p)<\widehat{Q}_{t}(p+f(t,p,\delta)))\leq\delta.$ (2.3)

Thus, if we wish to have a lower bound on a quantile (think small quantile),
we can study the Bernoulli process $S_{t}$. In order to explain the simulation
setup we begin by requiring that we have an initial waiting time $m$ and a
cutoff-time $M$ (for the simulation), and consider

$\mathbb{P}(\exists t:\,S_{t}/t<f(t,p,\delta))\leq\mathbb{P}(\exists
t\in[m,M]:\,S_{t}/t<f(t,p,\delta))\\\ +\mathbb{P}(\exists
t>M:\,S_{t}/t<f(t,p,\delta))$

The first term we will handle via simulation and the second term will be
analytically bounded using the following quantile specific version of Theorem
1, which can be found in [2]:

######  Theorem 4.

For any $p,\delta\in(0,1)$, $s,\eta>1$, and $M\geq 1$, we have

$\displaystyle\mathbb{P}(\exists t\geq
M:Q(p)<\widehat{Q}_{t}(p-g(t,M,p,\delta)))\leq\delta$

with

$\displaystyle
g(t,M,p,\delta)=\sqrt{k_{1}^{2}p(1-p)tl(t)+k_{2}^{2}c_{p}^{2}l^{2}(t)}+c_{p}k_{2}l(t)$

where

$\displaystyle\begin{cases}l(t)=s\log(\log(\eta
t/M))+\log(2\zeta(s)/(\delta\log^{s}\eta))\\\
k_{1}=(\eta^{1/4}+\eta^{-1/4})/\sqrt{2}\\\ k_{2}=(\sqrt{\eta}+1)/2\\\
c_{p}=(1-2p)/3,\end{cases}$

where $\zeta$ is the Riemann zeta function.

Figure 2. Histogram over $c_{Z}$.

We now start by defining the random variable

$\displaystyle Z=\min_{t\in[m,M]}\frac{S_{t}}{tg(t,p,\delta/2)}$

where the function $g$ is given in Theorem 4 with some fixed parameters, we
have chosen the same as in [2], namely $\eta=2.04$ and $s=1.4$. The goal of
the simulation is to estimate the $\delta/2$ quantile of $Z$ (see Fig. 2). If
we consider an empirical quantile $\widehat{Q}_{t,Z}(\delta/2)$, formed using
$t$ i.i.d. copies of $Z$ which we denote $\mathbf{Z}_{t}$, then using Lemma
A.1 we get that for any $\widetilde{\delta}\in(0,1)$,

$\displaystyle\mathbb{P}(\mathbb{P}(\widetilde{Z}\leq\widehat{Q}_{t,\mathbf{Z}_{t}}(\delta/2-\epsilon(t,\delta/2,\widetilde{\delta}))\mid\mathbf{Z}_{t})\leq\delta/2)\geq
1-\widetilde{\delta}$

where $\widetilde{Z}$ is an independent copy of $Z$. Translating this back to
the quantile process we obtain

$\displaystyle\mathbb{P}\bigg{[}\mathbb{P}\big{(}\exists
t\in[m,M]:\,\widehat{F}_{t}(Q(p))<\widehat{Q}_{t,\mathbf{Z}_{t}}\bigg{(}\frac{\delta}{2}-\epsilon\big{(}t,\frac{\delta}{2},\widetilde{\delta}\big{)}\bigg{)}g\big{(}t,p,\frac{\delta}{2}\big{)}\mid\mathbf{Z}_{t}\big{)}\leq\frac{\delta}{2}\bigg{]}\geq
1-\widetilde{\delta}.$

Thus, the introduction of the simulation step makes the quantile confidence
statement of a $2$-parameter PAC type, as described for the prediction
interval. If we later employ this in the same way as we did before but for the
prediction interval we will have a $3$-parameter PAC statement. We can see how
this sharpens the bounds in Fig. 3.

Figure 3. Quantile process with a simulated quantile compared to the bound
obtained by Theorem 4.

If we denote
$c_{Z}=-\widehat{Q}_{t,\mathbf{Z}_{t}}(\delta/2-\epsilon(t,\delta/2,\widetilde{\delta}))$
(think of $c_{Z}>0$) and $g_{t}=g(t,p,\delta/2)$, then the above can be
written more compactly as

$\displaystyle\mathbb{P}(\mathbb{P}(\exists
t\in[m,M]:\,S_{t}/t<-c_{Z}g_{t}\mid\mathbf{Z}_{t})\leq\delta/2)\geq
1-\widetilde{\delta}.$

Define now the $\mathbf{Z}_{t}$ dependent event

$\displaystyle
A=\\{Q_{Z}(\delta/2)\geq\widehat{Q}_{t,\mathbf{Z}_{t}}(\delta/2-\epsilon(t,\delta/2,\widetilde{\delta}))\\},$

then, on $A$ we have

$\displaystyle\mathbb{P}(\exists
t\in[m,M]:\,S_{t}/t<-c_{Z}g_{t})\leq\delta/2.$

Using the union bound we can estimate on the event $A$

$\mathbb{P}(\exists t\geq m:S_{t}/t<-\mathbb{I}_{t\leq
M}c_{Z}g_{t}-\mathbb{I}_{t>M}g(t,M,p,\delta/2))\\\ \leq\mathbb{P}(\exists
t\in[m,M]:S_{t}/t<-c_{Z}g_{t})+\mathbb{P}(\exists
t>M:S_{t}/t<-g(t,M,p,\delta/2))\leq\delta$

######  Example 2.1.

For $p=0.05,\delta=0.05,m=10,\widetilde{\delta}=0.0000001$, and $M=10000$ we
obtain that $c_{Z}$ is roughly $0.6$ and that parameters $\eta,s$ can be
chosen such that $g(M,M,p,\delta/2)$ is essentially the same size as
$c_{Z}g_{M}$.

Assembling everything we finally get,

$\displaystyle\mathbb{P}(\mathbb{P}(\exists t\geq
m:\widehat{F}_{t}(Q(p))<p-\mathbb{I}_{t\leq
M}c_{Z}g_{t}-\mathbb{I}_{t>M}g(t,M,p,\delta/2)\mid\mathbf{Z}_{t})\leq\delta)\geq
1-\widetilde{\delta}$ (2.4)

We can compare the bounds obtained in Fig. 4 for our particular choices.

Figure 4. Comparison of the different bounds.

###  2.1. $3$-parameter PAC for prediction intervals

Combining the arguments leading to 2.4 and 1.4 we can redefine the event $A$
to be (where we have defined $Z$ w.r.t. $p=\alpha/2$)

$\displaystyle
A=\\{Q_{Z}(\delta/4)\geq\widehat{Q}_{t,\mathbf{Z}_{t}}(\delta/4-g(t,\delta/4,\widetilde{\delta}))\\},$

then the event $A$ happens with probability $1-\widetilde{\delta}$. On the
event $A$ we then obtain that our prediction interval can be taken to be

$\displaystyle C^{l,\alpha}_{i}$
$\displaystyle=[\widehat{Q}_{i}(\alpha/2-|c_{Z}|g_{i}),\widehat{Q}_{i}^{-}(1-\alpha/2+|c_{Z}|g_{i})]$
(2.5)

and on the event $A$ we have the following version of 1.4

$\displaystyle\mathbb{P}(\forall i:\mathbb{P}(X\in C^{l,\alpha}_{i})\geq
1-\alpha)\geq 1-\delta.$ (2.6)

This constitutes a $3$-parameter PAC definition with parameters
$(\widetilde{\delta},\delta,\alpha)$, where the first can be taken as small as
we wish at the cost of our simulation.

##  3\. Application examples

In this section we will apply our $3$-parameter PAC formulation for prediction
intervals in conjunction with what is developed in Section 1.2. The first
example is the following regression problem

$\displaystyle r(x):=\mathbb{E}[Y\mid X=x]=\sin(3x)$
$\displaystyle\sigma(x):=\sqrt{\mathbb{V}[Y\mid X=x]}=|\sin(10x)|.$

We assume that

$\displaystyle Y=r(X)+\sigma(X)\epsilon,\quad\epsilon\sim N(0,1),$

and that $X\sim\text{unif}([0,1])$. We estimate $r$ using the Nadaraya-Watson
estimator with Gaussian kernel of bandwidth $1/500$ and estimate $\sigma(x)$
using again a Nadaraya-Watson estimator as in Appendix B. We then construct
the score function

$\displaystyle S=\frac{Y-r(X)}{\sigma(X)}$

and construct a prediction interval for $S$, then this is used to construct a
prediction interval for $Y$. We can see the result in Fig. 5.

Figure 5. Simulation of the regression problem. The solid lines are the upper
and lower limit of the prediction interval using the $2$-PAC definition and
the dotted lines represent the intervals with the $3$-PAC definition for
increasing number of samples. The dashed lines represent the true regression
function as well as the true conditional quantiles.

Since our method of constructing prediction intervals is asymmetric we can
have a skew distribution for $\epsilon$, in Fig. 6 we see the corresponding
result when $\epsilon$ is assumed to be a standardized log-Normal.

We can see an example of the estimated coverage as we keep adding data in Fig.
7

Figure 6. Simulation of the skew regression problem. The solid lines are the
upper and lower limit of the prediction interval using the $2$-PAC definition
and the dotted lines represent the intervals with the $3$-PAC definition for
increasing number of samples. The dashed lines represent the true regression
function as well as the true conditional quantiles. Figure 7. Estimated
coverage of the log-Normal noise using the 3-parameter PAC in Section 2.1.

If we instead of relying on the 3-parameter PAC we use Theorem 2, we get also
the data-dependent upper bound on the true coverage. A resulting simulation
can be found in Fig. 8.

Figure 8. Estimated coverage using Theorem 2.

##  Appendix A Helper lemmas

######  Theorem 5 (Bennett’s inequality).

Let $X_{1},\ldots,X_{n}$ be i.i.d. random variables with finite variance such
that $\mathbb{P}(X_{i}\leq b)=1$ with mean zero. Let and
$\sigma^{2}=\mathbb{V}[X_{i}]$. Then for any $\epsilon>0$,

$\displaystyle\mathbb{P}(\overline{X}_{n}\geq\epsilon)\leq\exp\left(-\frac{n\sigma^{2}}{b^{2}}h\left(\frac{b\epsilon}{\sigma^{2}}\right)\right)$

where $h(u)=(1+u)\log(1+u)-u$ for $u>0$.

######  Lemma A.1.

Let $Z_{1},\ldots$ be an i.i.d. sequence of random variables, let $Q$ and $F$
be the quantile function and distribution function respectively, then

$\displaystyle\mathbb{P}(Q(p)\leq\widehat{Q}_{t}(p-\epsilon))\leq\exp\left(-tp(1-p)h\left(\frac{\epsilon}{p(1-p)}\right)\right)$

where $h(u)=(1+u)\log(1+u)-u$ for $u>0$. Thus for a given $\delta>0$ we can
find $\epsilon(t,p,\delta)$ such that

$\displaystyle\mathbb{P}(Q(p)\leq\widehat{Q}_{t}(p-\epsilon(t,p,\delta)))\leq\delta$

###### Proof.

First note that by the definition of the quantile function, we have

$\displaystyle\mathbb{P}(Q(p)\leq Q_{t}(p-\epsilon))=$
$\displaystyle\mathbb{P}(\widehat{F}_{t}(Q(p))\leq\widehat{F}_{t}(Q_{t}(p-\epsilon)))$
$\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\frac{1}{t}\sum_{i=1}^{t}\mathbb{I}_{\\{Z_{i}\leq
Q(p)\\}}\leq p-\epsilon\right)$

Now $X_{i}=\mathbb{I}_{\\{Z_{i}\leq Q(p)\\}}$ is an i.i.d. sequence of
Bernoulli random variables with probability $p$, thus an application of
Theorem 5 gives the result. ∎

##  Appendix B Nadaraya-Watson for variance

The basis for the classical Nadaraya-Watson estimator is to first estimate the
joint density of the pair $(X,Y)\sim F_{XY}$, where $X,Y\in\mathbb{R}$, using
a Kernel density estimator. The Nadaraya-Watson estimator is then nothing but
the regression function $\mathbb{E}[Y\mid X]$, where the expectation is taken
with respect to the estimated density. Similarly we can estimate the
conditional variance, namely we can estimate $\mathbb{V}[Y\mid X]$ by
computing it using the estimated density. Here we provide the calculation for
reference:

We begin by selecting a kernel $K_{h}:\mathbb{R}\to\mathbb{R}_{+}$ such that
$\int K_{h}(x)dx=1$, $\int xK_{h}(x)dx=0$ and $\int x^{2}K_{h}(x)dx=h^{2}$,
i.e. a non-negative function, a classical choice is the Gaussian. Having now
an i.i.d. sequence of observations $Z_{i}=(X_{i},Y_{i})$, for $i=1,\ldots,n$
we construct the kernel density estimator as

$\displaystyle\widehat{\rho}_{n}(z)=\frac{1}{n}\sum_{i=1}^{n}K_{h}(x-X_{i})K_{h}(y-Y_{i}).$

From this definition we compute

$\displaystyle\mathbb{V}_{\widehat{\rho}_{n}}[Y\mid X]=\int
y^{2}\frac{\widehat{\rho}_{n}(z)}{\int\widehat{\rho}(x,y)dy}dy-\mathbb{E}_{\widehat{\rho}_{n}}[Y|X]^{2}.$

From which we see that we first need to compute

$\displaystyle\int\widehat{\rho}(x,y)dy$
$\displaystyle=\int\frac{1}{n}\sum_{i=1}^{n}K_{h}(x-X_{i})K_{h}(y-Y_{i})dy$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}K_{h}(x-X_{i})$

Furthermore,

$\displaystyle\int
y^{2}\frac{\widehat{\rho}_{n}(z)}{\int\widehat{\rho}(x,y)dy}dy$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\int
y^{2}\frac{K_{h}(x-X_{i})K_{h}(y-Y_{i})}{\frac{1}{n}\sum_{j=1}^{n}K_{h}(x-X_{j})}dy$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{K_{h}(x-X_{i})}{\frac{1}{n}\sum_{j=1}^{n}K_{h}(x-X_{j})}\int
y^{2}K_{h}(y-Y_{i})dy$
$\displaystyle=\frac{\sum_{i=1}^{n}K_{h}(x-X_{i})(h^{2}+Y_{i}^{2})}{\sum_{j=1}^{n}K_{h}(x-X_{j})}$

thus we conclude that the Nadaraya-Watson kernel estimator of variance is

$\displaystyle\mathbb{V}_{\widehat{\rho}_{n}}[Y\mid
X]=\frac{\sum_{i=1}^{n}K_{h}(x-X_{i})(h^{2}+Y_{i}^{2})}{\sum_{j=1}^{n}K_{h}(x-X_{j})}-\left(\frac{\sum_{i=1}^{n}K_{h}(x-X_{i})Y_{i}}{\sum_{j=1}^{n}K_{h}(x-X_{j})}\right)^{2}.$

## References

  * [1] A. Gammerman, V. Vovk and V. Vapnik  “Learning by Transduction”  In _Proceedings of the Fourteenth Conference on Uncertainty in Artificial Intelligence_ , UAI’98  Madison, Wisconsin: Morgan Kaufmann Publishers Inc., 1998, pp. 148–155
  * [2] Steven R Howard and Aaditya Ramdas  “Sequential estimation of quantiles with applications to A/B testing and best-arm identification”  In _Bernoulli_ 28.3 Bernoulli Society for Mathematical StatisticsProbability, 2022, pp. 1704–1728
  * [3] Jing Lei and Larry Wasserman  “Distribution-free prediction bands for non-parametric regression”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ 76.1 Oxford University Press, 2014, pp. 71–96
  * [4] Harris Papadopoulos, Kostas Proedrou, Volodya Vovk and Alex Gammerman  “Inductive confidence machines for regression”  In _Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19–23, 2002 Proceedings 13_ , 2002, pp. 345–356  Springer
  * [5] Leslie G Valiant  “A theory of the learnable”  In _Communications of the ACM_ 27.11 ACM New York, NY, USA, 1984, pp. 1134–1142
  * [6] Vladimir Vovk  “Conditional validity of inductive conformal predictors”  In _Asian conference on machine learning_ , 2012, pp. 475–490  PMLR
  * [7] Vladimir Vovk, Alexander Gammerman and Glenn Shafer  “Algorithmic learning in a random world”  Springer, 2005


[c,d,e]Harvey B. Meyer

# Deep inelastic scattering off quark-gluon plasma
and its photon emissivity

Marco Cè  Tim Harris  Arianna Toniato  Csaba Török

###### Abstract

The photon emissivity of quark-gluon plasma probes the interactions in the
medium and differs qualitatively between a weakly coupled and a strongly
coupled plasma in the soft-photon regime. The photon emissivity is given by
the product of kinematic factors and a spectral function associated with the
two-point correlator of the electromagnetic current at lightlike kinematics. A
certain Euclidean correlator at imaginary spatial momentum can be calculated
in lattice QCD and is given by an integral over the relevant spectral function
at lightlike kinematics. I present a first exploratory lattice calculation of
this correlator. Secondly, I show how Euclidean correlators at imaginary
spatial momenta can also be used to probe the regime of deep inelastic
scattering off quark-gluon plasma, which reveals its parton distribution
function.

Preprint number: CERN-TH-2021-167

##  1 Introduction

In this presentation I will be concerned with the vector spectral functions of
strong-interaction matter at non-vanishing temperature, which can be defined
as (in diag$(+,-,-,-)$ metric)

$\rho^{\mu\nu}(q)=\int d^{4}x\,e^{iq\cdot
x}\,\frac{1}{Z}\sum_{n}e^{-E_{n}/T}\langle
n|[j^{\mu}(x),j^{\nu}(0)]|n\rangle,\qquad Z=\sum_{n}e^{-E_{n}/T}.$ (1)

Here $Z$ is the canonical partition function, and $j^{\mu}$ is a conserved
vector current, the case of the electromagnetic current
$j^{\mu}=\frac{2}{3}\bar{u}\gamma^{\mu}u-\frac{1}{3}\bar{d}\gamma^{\mu}d-\frac{1}{3}\bar{s}\gamma^{\mu}s+\dots$
being particularly important for the phenomenology of heavy-ion collisions [1,
2] and the physics of the early universe (see e.g. [3]), since it is via this
current that the medium interacts with photons. Therefore we will have this
case in mind in the following.

The components of $\rho^{\mu\nu}(q)$ can be parametrized by two independent
kinematic variables, which can be chosen to be the virtuality $q^{2}$ and the
photon energy $q^{0}$ in the rest frame of the thermal medium. At vanishing
virtuality $q^{2}=0$, the spectral functions describe the rate at which the
medium emits photons (see Eq. (10) below), and for positive virtuality, they
describe the production rate of lepton pairs via a timelike photon [4]. At
spacelike virtualities, the vector spectral functions measure the ability of
the medium to convert the energy stored in external electromagnetic fields
into heat [5]. A conceivable way to create electromagnetic fields with
frequencies on the order of a GeV is to scatter a lepton of energy $E$ on the
medium. The process is illustrated in Fig. 1. In the one-photon exchange
approximation, the cross-section for scattering off the medium of volume
$L^{3}$ reads

$\frac{d^{2}\sigma}{dE^{\prime}d\Omega}=\frac{e^{4}\,L^{3}}{8\pi^{2}Q^{4}}\left(\frac{E^{\prime}}{E}\right)\ell_{\mu\nu}W_{>}^{\mu\nu}(u,q),$
(2)

in the rest frame of the medium, with $E^{\prime}$ the final-state energy of
the lepton and $\Omega$ the solid angle of its outgoing momentum relative to
the incident momentum [6]. Here $\ell_{\mu\nu}$ is the (exactly known)
leptonic tensor and $W_{>}^{\mu\nu}$ the ‘hadronic’ tensor, names that we
borrow from the literature on deep-inelastic scattering on the nucleon. The
hadronic tensor is defined as

$W_{>}^{\mu\nu}(u,q)=\frac{1}{4\pi\,Z}\sum_{n}e^{-\beta E_{n}}\int
d^{4}x\;e^{iq\cdot x}\,\langle n|j^{\mu}(x)\,j^{\nu}(0)|n\rangle,$ (3)

with $u$ the four-velocity of the medium, and is related by the Kubo–Martin-
Schwinger relation to the spectral functions in the rest frame of the medium
[7],

$W_{>}^{\mu\nu}(u,q)=\frac{1}{4\pi(1-e^{-\beta
q^{0}})}\,\rho^{\mu\nu}(q^{0},\vec{q}),\qquad u=(1,\vec{0}).$ (4)

Particularly for studying the deep-inelastic scattering (DIS) limit, a
suitable tensor decomposition is

$\displaystyle W_{>}^{\mu\nu}(u,q)$ $\displaystyle=$ $\displaystyle
F_{1}(u\cdot q,Q^{2})\Big{(}-g^{\mu\nu}+\frac{q^{\mu}q^{\nu}}{q^{2}}\Big{)}$
(5) $\displaystyle+\frac{T}{u\cdot q}F_{2}(u\cdot
q,Q^{2})\Big{(}u^{\mu}-(u\cdot
q)\frac{q^{\mu}}{q^{2}}\Big{)}\Big{(}u^{\nu}-(u\cdot
q)\frac{q^{\nu}}{q^{2}}\Big{)},$

with $F_{1}$, $F_{2}$ being the structure functions of the medium. We
introduce the Bjorken variable

$x=\frac{Q^{2}}{2Tq^{0}},$ (6)

in which the temperature $T$ represents the energy scale of the medium used
here to make $x$ dimensionless and $Q^{2}=-q^{2}$. Note that $x$ is not
bounded from above by unity in the present context.

Figure 1: Scattering of a lepton on quark-gluon plasma at global thermal
equilibrium moving with a four-velocity $u$ in the one-photon exchange
approximation. The picture shows the interaction occurring with a cubic fluid
cell, producing an unobserved QCD final-state $X$.

##  2  Probing the structure functions with Euclidean correlation functions

Figure 2:  The $x_{3}$-integrand to obtain $H_{E}^{T}(\omega_{n},0)$ based on
Eq. (11) for the isovector current with quark charges ${\cal Q}_{u}=-{\cal
Q}_{d}=1/\sqrt{2}$, after integrating over $d{\vec{x}}_{\perp}\equiv
dx_{0}dx_{1}dx_{2}$. The area under both displayed correlators would be
$T^{2}/2$ for free massless quarks [9]. In general, the non-negativity of
$\sigma^{T}$ implies that the area can be no smaller for the $n=2$ than for
the $n=1$ correlator. The calculation is performed on a $24\times 96^{3}$
lattice QCD ensemble with two dynamical flavours of light quarks and $T\simeq
254\,$MeV, the X7 ensemble described in [5]. The displayed data is based on
1331 configurations, with 64 point sources used per configuration. At the
largest $x_{3}$, the effect of the periodicity of this variable can be seen in
the $n=1$ correlator.

We focus here on the spatially transverse channel represented by the structure
function $F_{1}$ and define the corresponding spectral function

$\sigma^{T}(q^{0},Q^{2})=8\pi(1-e^{-\beta q^{0}})F_{1}(q^{0},Q^{2}),$ (7)

viewed as a function of the photon energy and virtuality. Choose now a
Matsubara frequency $\omega_{n}=2\pi nT$ and a particular value of
$Q^{2}\leq\omega_{n}^{2}$. In [6], it is shown that the $q^{0}$-dependence of
this spectral function can be probed via the Euclidean correlator

$\displaystyle H_{E}^{T}(\omega_{n};Q^{2})$ $\displaystyle=$
$\displaystyle-2\int_{0}^{\beta}dx_{0}\int
d^{3}x\;e^{i\omega_{n}x_{0}+x_{3}\sqrt{\omega_{n}^{2}-Q^{2}}}\;\big{\langle}j_{1}(x)\,j_{1}(0)\big{\rangle}.$
(8)

This correlator contains a logarithmic ultraviolet divergence111An exception
to this statement is the case of $Q^{2}=0$, when the theory is regularized in
a Lorentz-invariant way., which cancels in the difference on the left-hand
side of Eq. (9). Indeed this difference would vanish in the vacuum, since in
the latter state $H_{E}^{T}$ only depends on $Q^{2}$. In order to reach these
conclusions, it is sufficient to Taylor-expand around $x=0$ the weight
function with which the position-space vector correlator is multiplied, taking
into account that this correlator is even under the coordinate transformation
$x_{0}\to-x_{0}$ as well as under $x_{3}\to-x_{3}$. For $\omega_{n}^{2}\geq
Q^{2}$, note the ‘imaginary spatial momentum’ involved in computing the
correlator $H_{E}^{T}$. This correlator admits the fixed-virtuality dispersive
representation

$H_{E}^{T}(\omega_{n};Q^{2})-H_{E}^{T}(\omega_{r};Q^{2})=\int_{0}^{\infty}\frac{dq^{0}}{\pi}\>q^{0}\>\sigma^{T}(q^{0},Q^{2})\Big{[}\frac{1}{(q^{0})^{2}+\omega_{n}^{2}}-\frac{1}{(q^{0})^{2}+\omega_{r}^{2}}\Big{]}\>.$
(9)

Since $\sigma^{T}(q^{0},Q^{2}=0)\sim(q^{0})^{1/2}$ at weak coupling [8], one
subtraction is expected to be sufficient to guarantee convergence of the
dispersion relation. In practice, one has the option to compute the left-hand
side of Eq. (9) using at long distances
$e^{x_{3}\sqrt{\hat{\omega}_{n}^{2}-Q^{2}}}$ with
$\hat{\omega}_{n}^{2}=\omega_{n}^{2}-{\rm O}(a^{2})$ slightly reduced in
magnitude, in order to avoid convergence issues in the infrared at non-zero
lattice spacing.

Of special interest is the case of zero virtuality, since the differential
photon emissivity of the medium is given by

$\frac{d\Gamma_{\gamma}}{dq^{0}}=\frac{\alpha}{\pi}\,\frac{q^{0}}{e^{\beta
q^{0}}-1}\;\sigma^{T}(q^{0},Q^{2}=0).$ (10)

At $Q^{2}=0$, $H_{E}^{T}$ is ultraviolet-finite in a regularisation respecting
Lorentz symmetry. On the lattice however, it is necessary to perform a
subtraction in order to obtain the correct continuum limit for the correlator
$H_{E}^{T}(\omega_{n})$ [9]. For that purpose, we introduce here the estimator

$H_{E}^{T}(\omega_{n},0)=-2\int_{0}^{\beta}dx_{0}\,\int
d^{3}x\,(e^{i\omega_{n}x_{0}}-e^{i\omega_{n}x_{2}})\,e^{\omega_{n}x_{3}}\;\big{\langle}j_{1}(x)j_{1}(0)\big{\rangle}.$
(11)

Note that the subtraction term ($e^{i\omega_{n}x_{2}}$) vanishes in the
continuum thermal theory, but on the lattice removes an ultraviolet divergence
associated with the lack of Lorentz symmetry at finite lattice spacing. This
is easiest to see from the fact that estimator (11) of $H_{E}^{T}(\omega_{n})$
automatically vanishes on the $(T=0)$ infinite lattice. With estimator (11),
no cumbersome subtraction of vacuum correlators is needed, which would require
costly simulations. The integrand of Eq. (11) in a lattice QCD calculation at
$T\approx 250\,$MeV is illustrated in Fig. 2. Similar to the comment below Eq.
(9), the exponential weight factor $e^{\omega_{n}x_{3}}$ can be replaced by a
factor with a slightly (O($a^{2}$)) reduced exponent, for instance by
substituting $\omega_{n}\to\frac{2}{a}\sin\frac{a\omega_{n}}{2}$. Given that
$H_{E}^{T}(0,0)=0$, we can write a simplied dispersion relation for the light-
like case,

$H_{E}^{T}(\omega_{n},0)=-\frac{\omega_{n}^{2}}{\pi}\int_{0}^{\infty}\frac{dq^{0}}{q^{0}}\;\frac{\sigma^{T}(q^{0})}{(q^{0})^{2}+\omega_{n}^{2}}\;.$
(12)

##  3 Thermal medium structure functions in the DIS regime

Figure 3:  The structure function $F_{1}(x,Q^{2})$ of the free plasma for
different values of $Q^{2}$, together with the Bjorken limit Eq. (14). Here we
have set $N_{c}=3$ and $\sum_{f}{\cal Q}_{f}^{2}=1$.

As in the case of the nucleon, one can show that the structure functions admit
a partonic interpretation in the DIS limit [6]. With $u$ being the four-
velocity of the fluid, let $f_{f}(\xi)\,d\xi$ represent the number of partons
of type $f$ in the fluid cell that carry momentum $\xi\>T\>u$. This number is
of order the cell volume, $L^{3}$. In the DIS limit, the structure function
reads

$F_{1}(u\cdot q,Q^{2})\stackrel{{\scriptstyle
Q^{2}\to\infty}}{{\longrightarrow}}\frac{1}{4L^{3}T}\sum_{f}{\cal
Q}_{f}^{2}\;f_{f}(x),$ (13)

with $x=Q^{2}/(2Tu\cdot q)$ kept fixed. In words, $4\,F_{1}\cdot dx$ is the
square-electric-charge weighted number of partons carrying a momentum $xT$
times the fluid four-velocity $u^{\mu}$ per unit transverse area in a slab of
fluid which in its rest frame has thickness $1/T$ in the longitudinal
direction.

In the theory of free quarks, in which the spectral functions can be computed
analytically [10], one finds in the DIS limit

$\lim_{Q^{2}\to\infty}F_{1}(x,Q^{2})=\frac{(\sum_{f}\mathcal{Q}_{f}^{2})N_{c}T^{2}}{4\pi^{2}}x\log(1+e^{-x/2}),$
(14)

thus showing that the parton distribution is proportional to
$x\log(1+e^{-x/2})$, and normalized such that the partons altogether carry the
entire momentum of the aforementioned fluid cell as given by its rest-frame
enthalpy, consistently with the expectation from ideal hydrodynamics [6]. The
approach to the DIS limit as a function of $x$ is illustrated in the free-
quark theory in figure 3. One can see that it is only in the limit
$Q^{2}\to\infty$ that the structure function acquires the normalizable form of
a probability distribution in the variable $x$. This represents a qualitative
difference with respect to ordinary DIS on the nucleon.

How can the DIS regime of the structure functions be addressed using lattice
QCD? Using the twist expansion of the product of two vector currents, one
finds the moment sum rules obeyed by the structure functions [6],

$\int_{0}^{\infty}dx\>x^{n-1}\>[F_{1}(x,Q^{2})]_{\textrm{leading-
twist}}=\frac{1}{2}\sum_{f,j}{\mathcal{Q}_{f}^{2}}M_{fj}(Q,\tilde{\mu})\langle
O_{nj}\rangle\>,\quad n=2,4,\dots\>,$ (15)

where the $M_{fj}$ are the coefficients parametrizing the operator mixing and
the $O_{nj}^{\mu_{1}\dots\mu_{n}}$ are the usual twist-two operators of mass
dimension $(n+2)$, whose thermal expectation values are parametrized as
follows,

$\langle O_{nj}^{\mu_{1}\dots\mu_{n}}\rangle=T^{n}[u^{\mu_{1}}\dots
u^{\mu_{n}}-\mathrm{traces}]\langle O_{nj}\rangle\>.$ (16)

The index $j$ runs over the quark flavours and the gluon operator. The right-
hand side of the moment sum rule is computable on the lattice for the lowest
few values of $n$. A difference with ordinary DIS is however that the leading-
twist component of $F_{1}$ must be taken before calculating the $x^{n-1}$
moment in Eq. (15), as the integral otherwise does not converge, as can be
seen from Fig. 3. Alternatively to computing the lowest $x$-moments of the
structure functions, the approach to the DIS limit can be studied by computing
the correlators $H_{E}^{T}(\omega_{n};Q^{2})$ for a sequence of Matsubara
frequencies $\omega_{n}$. Performing the change of variables
$x=Q^{2}/(2Tq^{0})$ in Eq. (9), one obtains

$H_{E}^{T}(\omega_{n};Q^{2})-H_{E}^{T}(\omega_{r};Q^{2})=\int_{0}^{\infty}\frac{dx}{\pi}\>x\>\hat{\sigma}^{T}(x,Q^{2})\frac{a_{r}^{2}-a_{n}^{2}}{(1+a_{n}^{2}x^{2})(1+a_{r}^{2}x^{2})}\>,$
(17)

with $a_{n}=2T\omega_{n}/Q^{2}$. Thus if $a_{n}$ and $a_{r}$ are kept fixed as
$Q^{2}$ is varied, the only $Q^{2}$ dependence on the left-hand side comes
from the spectral function
$\hat{\sigma}^{T}(x,Q^{2})\equiv\sigma^{T}(q^{0},Q^{2})$. The extent to which
it becomes independent of $Q^{2}$ in the Bjorken limit can thus be probed in
this way.

##  4 Conclusion

Dispersion relations at fixed spacelike virtuality, rather than at fixed
spatial momentum, open up new perspectives on the thermal spectral functions.
Here we have shown that certain moments of the spectrum of emitted photons can
be computed in lattice QCD without solving an inverse problem. The general
lesson is that when an analytic continuation in one variable of a Euclidean
correlation function depending on several kinematic variables is performed,
which numerically amounts to solving an inverse problem, it plays an important
role which of these kinematic variables are kept fixed.

We have also seen that moments of the in-medium structure functions can be
computed in lattice QCD at fixed spacelike virtuality. In this kinematic
regime, the structure functions describe the scattering of a lepton on the
medium, and, at sufficiently high virtuality in the Bjorken regime, should
reveal its quark-gluon constituents. Such a regime is however difficult to
reach with the current simulation techniques, due to the real exponential
$e^{x_{3}\sqrt{\omega_{n}^{2}-Q^{2}}}$ emphasizing the large distances in the
spatial directions. One must also remember that even in the theory of free
quarks, the Bjorken scaling is only reached after ‘higher-twist’
(O($1/Q^{2}$)) contributions are sufficiently suppressed.

## Acknowledgments

This work was supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program through Grant
Agreement No. 771971-SIMDAMA, as well as by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) through the Cluster
of Excellence “Precision Physics, Fundamental Interactions and Structure of
Matter” (PRISMA+ EXC 2118/1) funded by the DFG within the German Excellence
strategy (Project ID 39083149). The work of M.C. is supported by the European
Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie Grant Agreement No. 843134-multiQCD. T.H. is supported by UK
STFC CG ST/P000630/1. The generation of gauge configurations as well as the
computation of correlators was performed on the Clover and Himster2 platforms
at Helmholtz-Institut Mainz and on Mogon II at Johannes Gutenberg University
Mainz. The authors gratefully acknowledge the Gauss Centre for Supercomputing
e.V. (www.gauss-centre.eu) for funding project IMAMOM by providing computing
time through the John von Neumann Institute for Computing (NIC) on the GCS
Supercomputer JUWELS [11] at Jülich Supercomputing Centre (JSC).

## References

  * [1] P. Braun-Munzinger, V. Koch, T. Schaefer, and J. Stachel, Properties of hot and dense matter from relativistic heavy ion collisions, Phys. Rept. 621 (2016) 76–126, [arXiv:1510.0044].
  * [2] G. David, Direct real photons in relativistic heavy ion collisions, Rept. Prog. Phys. 83 (2020), no. 4 046301, [arXiv:1907.0889].
  * [3] T. Asaka, M. Laine, and M. Shaposhnikov, On the hadronic contribution to sterile neutrino production, JHEP 06 (2006) 053, [hep-ph/0605209].
  * [4] L. D. McLerran and T. Toimela, Photon and Dilepton Emission from the Quark - Gluon Plasma: Some General Considerations, Phys.Rev. D31 (1985) 545.
  * [5] M. Cè, T. Harris, H. B. Meyer, A. Steinberg, and A. Toniato, The rate of photon production in the quark-gluon plasma from lattice QCD, Phys. Rev. D 102 (2020), no. 9 091501, [arXiv:2001.0336].
  * [6] M. Cè, T. Harris, H. B. Meyer, and A. Toniato, Deep inelastic scattering on the quark-gluon plasma, JHEP 03 (2021) 035, [arXiv:2012.0752].
  * [7] H. B. Meyer, Transport Properties of the Quark-Gluon Plasma: A Lattice QCD Perspective, Eur.Phys.J. A47 (2011) 86, [arXiv:1104.3708].
  * [8] S. Caron-Huot, P. Kovtun, G. D. Moore, A. Starinets, and L. G. Yaffe, Photon and dilepton production in supersymmetric Yang-Mills plasma, JHEP 0612 (2006) 015, [hep-th/0607237].
  * [9] H. B. Meyer, Euclidean correlators at imaginary spatial momentum and their relation to the thermal photon emission rate, Eur. Phys. J. A54 (2018), no. 11 192, [arXiv:1807.0078].
  * [10] M. Laine, NLO thermal dilepton rate at non-zero momentum, JHEP 1311 (2013) 120, [arXiv:1310.0164].
  * [11] Jülich Supercomputing Centre, JUWELS: Modular Tier-0/1 Supercomputer at the Jülich Supercomputing Centre, Journal of large-scale research facilities 5 (2019), no. A135.