KwInputInput KwOutputOutput toto

# Adaptive Sampling for Estimating Distributions: A Bayesian Upper Confidence
Bound Approach

Dhruva Kartik<EMAIL_ADDRESS>
Neeraj Sood<EMAIL_ADDRESS>
Urbashi Mitra<EMAIL_ADDRESS>
University of Southern California  Los Angeles  Tara Javidi<EMAIL_ADDRESS>
University of California  San Diego

###### Abstract

The problem of adaptive sampling for estimating probability mass functions
(pmf) uniformly well is considered. Performance of the sampling strategy is
measured in terms of the worst-case mean squared error. A Bayesian variant of
the existing upper confidence bound (UCB) based approaches is proposed. It is
shown analytically that the performance of this Bayesian variant is no worse
than the existing approaches. The posterior distribution on the pmfs in the
Bayesian setting allows for a tighter computation of upper confidence bounds
which leads to significant performance gains in practice. Using this approach,
adaptive sampling protocols are proposed for estimating SARS-CoV-2
seroprevalence in various groups such as location and ethnicity. The
effectiveness of this strategy is discussed using data obtained from a
seroprevalence survey in Los Angeles county.

###### keywords:

Adaptive sampling, Distribution estimation, Bayesian methods, Upper Confidence
Bound

##  1 Introduction

We frequently encounter scenarios where we need to estimate a finite
collection of pmfs from their samples. We have a limited sample budget and we
would like to adaptively choose which pmf to sample from so that we have
uniformly good estimates of all the pmfs. The _goodness_ of our estimate can
be evaluated using various distribution distance metrics. In this paper, we
will focus on the mean squared error. A concrete application that fits this
setting is that of estimating SARS-CoV-2 incidence rates in various geographic
regions. Incidence rate, the fraction of _currently_ infected individuals, is
a crucial metric for assessing the safety of a given region. Incidence rate
cannot be estimated in an unbiased way simply based on the reported cases, it
needs to be estimated by conducting randomized tests in each region of
interest. Typically, such randomized surveys can be conducted by a central
authority which has a limited testing budget. A natural question that arises
in this scenario is that of allocating the available tests to various regions
in an efficient manner. We may also be interested in other parameters such as
SARS-CoV-2 seroprevalence which is an indicator of an individual’s potential
immunity against the virus, and other kinds of groups like ethnicity, age etc.
in addition to the geographic location. Other applications of this setting
include dynamics estimation in Markov Decision Processes (MDPs) and text
compression which are discussed in greater detail in [Shekhar et
al.(2020)Shekhar, Ghavamzadeh, and Javidi].

When we are interested in minimizing the mean squared error for every pmf, an
efficient sampling strategy would be to let the number of tests be
proportional to the _variance_ of the associated pmf. However, the key problem
is that we do not know the variance of the pmf. Therefore, we need to adapt
our allocation as we learn the distribution over time. An upper confidence
bound (UCB) based approached was proposed in [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi]. The main idea in this approach is to first form an
upper bound on the variance of each pmf. As a function of this upper bound and
the number of samples acquired so far for each pmf, the next sample is chosen
accordingly. It is shown that the number of samples obtained in this manner
would approximately be the same as that of an _oracle_ that allocates samples
using the knowledge of the variance of the underlying distributions. While
this strategy satisfies many such interesting theoretical properties, we
observed that it tends be more conservative than necessary in practice.

Our main goal in this paper is to devise a less conservative approach while
ensuring the same performance guarantees as in [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi]. Our contributions towards this goal are (i) we
assume a Dirichlet prior on the pmfs and compute aforementioned upper
confidence bounds using the posterior belief on the pmfs. This, combined with
a stronger set of inequalities gives us noticeably tighter bounds on the
variance of the pmfs; (ii) we analytically prove that this Bayesian approach
for computing confidence bounds is no worse than the one in [Shekhar et
al.(2020)Shekhar, Ghavamzadeh, and Javidi]; and (iii) we employ these methods
for estimating SARS-CoV-2 seroprevalence in various categories using data
obtained by a randomized survey [Sood et al.(2020)Sood, Simon, Ebner, Eichner,
Reynolds, Bendavid, and Bhattacharya] conducted in Los Angeles county. Our
analysis of the survey data illustrates the importance of adaptive sampling in
order to obtain better estimates. We would like to emphasize that the Bayesian
approach used herein to compute upper confidence bounds and the associated
proof methodologies can also be adapted to other UCB based approaches that are
widely used in online learning [Efroni et al.(2020)Efroni, Mannor, and
Pirotta, Rosenberg and Mansour(2019)] and dynamics estimation [Tarbouriech and
Lazaric(2019)].

#### Related Work

Our work is closely related to [Shekhar et al.(2020)Shekhar, Ghavamzadeh, and
Javidi]. We adopt their sampling strategy and many of the proof methodologies.
The key difference between our work and [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi] is our Bayesian approach for computing the upper
confidence bounds on the variances. Our approach for computing these bounds
has a significant impact on the performance in practice, especially in the
non-asymptotic regime. Multiple distance metrics have been considered in
[Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi], whereas we restrict
our attention to the mean squared error. It is however possible to extend our
methodology to the other distance metrics as well. When the pmfs of interest
have a binary support, which is the case in the problem of seroprevalence
estimation, then the problem of estimating distributions reduces to the
problem of estimating means of random variables. Substantial amount of work
has been done in the area of estimating the means of random variables
[Carpentier et al.(2011)Carpentier, Lazaric, Ghavamzadeh, Munos, and Auer]. In
[Carpentier et al.(2011)Carpentier, Lazaric, Ghavamzadeh, Munos, and Auer],
two approaches for computing upper confidence bounds were proposed, one using
the Chernoff-Hoeffding bound and the other using the empirical Bernstein
bound. The Chernoff-Hoeffding bound approach is very close to that of [Shekhar
et al.(2020)Shekhar, Ghavamzadeh, and Javidi]. We observe that the performance
of our Bayesian approach is better than both the approaches in [Carpentier et
al.(2011)Carpentier, Lazaric, Ghavamzadeh, Munos, and Auer]. However, the
regret bounds in [Carpentier et al.(2011)Carpentier, Lazaric, Ghavamzadeh,
Munos, and Auer] for the Bernstein approach are tighter than ours.

#### Notation

In general, subscripts denote time indices unless stated otherwise. For a
sampling strategy $g$ and a collection of pmfs $p$, we use
${\mathbb{P}}^{g}_{p}[\cdot]$ and ${\mathbb{E}}^{g}_{p}[\cdot]$ to indicate
respectively that the probability and expectation depend on the choice of $g$,
and that they are conditioned on the model $p$. We denote the indicator
function associated with an event $\mathcal{E}$ with
$\mathbbm{1}_{\mathcal{E}}$.

##  2 Problem Formulation

We are interested in estimating $K$ probability mass functions (pmfs) over a
finite space $\\{1,\dots,L\\}$ from their samples. Let the $k$-th pmf be
denoted by ${p}^{(k)}\doteq\left({p}^{(k,1)},\dots,{p}^{(k,L)}\right)$ and let
$p\doteq\left[{p}^{(1)},\dots,{p}^{(K)}\right]$. For convenience, we will
refer to each of these pmfs as an _arm_.

The total number of samples that can be obtained is $N$. At any given time, we
can adaptively choose an arm to obtain a sample. The $n$-th arm sampled is
denoted by $U_{n}$. Let the outcome of the $n$-th sample be denoted by
$Y_{n}$. We have $Y_{n}\sim{p}^{(U_{n})}.$ Let the history of sampled arms and
their corresponding outcomes be denoted by
$I_{n}\doteq\\{U_{1},Y_{1},\dots,U_{n-1},Y_{n-1}\\}$. Let the strategy used to
collect samples at time $n$ be $g_{n}$. In other words, $U_{n}=g_{n}(I_{n})$.
After collecting $n$ samples, we form an estimate ${{\hat{p}}}^{(k)}_{n}$ of
each pmf ${p}^{(k)}$, which in our case, is the empirical distribution. The
mean squared error between the true distribution ${p}^{(k)}$ and the estimate
is denoted by

$\displaystyle{\mathcal{L}}^{(k)}_{N}({p}^{(k)},g)={\mathbb{E}}_{p}^{g}\left[\mathrm{MSE}\left({p}^{(k)},{{\hat{p}}}^{(k)}_{N}\right)\right].$
(1)

######  Remark 2.1.

Note that we can scale the loss ${\mathcal{L}}^{(k)}$ by an appropriate weight
${w}^{(k)}$. Our sampling strategy and proofs can be easily adapted to this
weighted loss.

As discussed earlier, our sampling strategy will involve the computation of
upper confidence bounds on certain parameters. When Bayesian methods are used
to derive upper confidence bounds, it is generally not possible to provide
performance guarantees for every instance $p$ of the underlying pmfs. In such
cases, the loss $\mathcal{L}$ is _locally averaged_ over a small region
containing the true pmf $p$. We will therefore consider a loss that is locally
averaged with respect to $p$ using a distribution $\varrho$. We refer to
[Brown et al.(2001)Brown, Cai, and DasGupta, Bayarri and Berger(2004)] for a
detailed discussion on the interpretation of local averaging.

######  Remark 2.2.

Note that the distribution $\varrho$ used for local averaging is not the same
as the prior on $p$. We will use a Bayesian interpretation of the pmfs $p$
only to compute tighter upper confidence bounds. However, our formulation and
performance guarantees are frequentist in nature since the distribution
$\varrho$ is supported on an arbitrarily small region around $p$.

Let $\Delta$ denote the $L-1$ dimensional simplex and let
$\varrho={\varrho}^{(1)}\times\dots\times{\varrho}^{(K)}$ be a product
distribution on $\Delta^{K}$. Our goal is to minimize the following cost
associated with a sampling strategy $g$ which is given by

$\displaystyle
J_{N}(g,\varrho)\doteq\max_{k}\left\\{{\mathbb{E}}^{{\varrho}}\left[{\mathcal{L}}^{(k)}_{N}({\pi}^{(k)},g)\right]\right\\},$
(2)

where ${\pi}\sim{\varrho}$. The distribution ${\varrho}^{(k)}$ has the
following form. For any Borel measurable set $\mathcal{A}\subseteq\Delta$,

$\displaystyle{\varrho}^{(k)}(\mathcal{A})=\mu\left(\mathcal{A}\cap\eta({p}^{(k)})\right)/\eta,$
(3)

where $\mu$ is the uniform distribution over $\Delta$ and $\eta({p}^{(k)})$ is
an $\ell_{2}$-ball around ${p}^{(k)}$, such that
$\mu\left(\eta({p}^{(k)})\right)=\eta$. Clearly, the cost in (2) approaches
the cost in [Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi, Carpentier
et al.(2011)Carpentier, Lazaric, Ghavamzadeh, Munos, and Auer] when
$\eta\rightarrow 0$ and thus, it is desirable to select small $\eta$. However,
the performance guarantees we provide in this paper become weaker as $\eta$
becomes small. We will discuss the effects of the choice of $\eta$ in
subsequent sections.

#### The Oracle

Let there be an oracle that knows the probability distributions $p$. Using
this information, the oracle allocates a fixed number of samples for each arm.
If the number of samples allocated to arm $k$ is ${T}^{(k)}$, the mean-squared
error associated with arm $k$’s estimate is given by

$\displaystyle{\mathbb{E}}_{p}^{g}\left[\mathrm{MSE}\left({p}^{(k)},{{\hat{p}}}^{(k)}_{N}\right)\right]$
$\displaystyle=\frac{\sum_{l}{p}^{(k,l)}(1-{p}^{(k,l)})}{{T}^{(k)}}$ (4)
$\displaystyle\doteq\frac{{c}^{(k)}}{{T}^{(k)}}\doteq\varphi({c}^{(k)},{T}^{(k)}),$
(5)

where $\varphi$ is referred to as the _tracking function_ and ${c}^{(k)}$ is
referred to as the _tracking parameter_ [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi]. The oracle solves the following optimization problem
to obtain an allocation.

$\displaystyle\min_{{T}^{(1)},\dots,{T}^{(K)}}\;\max_{k}\varphi({c}^{(k)},{T}^{(k)})$
$\displaystyle\mathrm{s.t.}\;\sum_{k}{T}^{(k)}=N.$ (P1)

While Problem (P1) is a combinatorial optimization problem, an approximate
solution can easily be obtained by making a convex relation. The non-integer
solution and the optimal value for this problem can be expressed in closed
form as

$\displaystyle{T}^{(k)}_{*}=\frac{{c}^{(k)}N}{\sum_{i}{c}^{(i)}};$
$\displaystyle\varphi^{*}(p,N)=\frac{\sum{c}^{(k)}}{N}.$ (6)

The regret with respect to the oracle that knows the underlying distributions
is

$\displaystyle\mathcal{R}_{N}(g,\varrho)\doteq
J_{N}(g,\varrho)-{\mathbb{E}}^{{\varrho}}\left[\varphi^{*}(p,N)\right].$ (7)

##  3 Definitions and Framework

We will first define some important quantities and compute an upper confidence
bound on the tracking parameter ${c}^{(k)}$. This bound will be used in our
sampling strategy described in Section 4.

###  3.1 Prior Belief

Let us assume that each of the pmfs ${p}^{(k)}$ is independently drawn from a
Dirichlet distribution [Kotz et al.(2004)Kotz, Balakrishnan, and Johnson] with
parameters
${\alpha}^{(k)}_{1}\doteq\left({\alpha}^{(k,1)}_{1},\dots,{\alpha}^{(k,L)}_{1}\right)$.
In other words,

$\displaystyle({p}^{(1)},\dots,{p}^{(K)})\sim{\mathscr{D}}({\alpha}^{(1)}_{1})\times\dots\times{\mathscr{D}}({\alpha}^{(K)}_{1}).$
(8)

This prior distribution111Note that this is distinct from the distribution
$\varrho$ that was used for local averaging. will be denoted by $\rho_{1}$. We
will use the uniform prior which is equivalent to setting
${\alpha}^{(k,l)}_{1}=1$ for every $k,l$. Thus, $\rho_{1}$ is the same as the
distribution $\mu$ used in (41).

######  Remark 3.1.

Note that the prior $\rho_{1}$ is just something we assume. In a frequentist
setting with local averaging, we are not provided with any prior information
on the pmfs $p$. Since the distribution $\varrho$ is defined using the uniform
distribution $\mu$ in (41), the uniform prior $\rho_{1}=\mu$ happens to be a
convenient choice for computing our confidence bounds and subsequently
analyzing them. In some cases however, we do have some prior information on
$p$. A discussion on how to incorporate certain kinds of priors is in Appendix
D.

###  3.2 Posterior Update

When the prior on the pmfs is a factored Dirichlet distribution as in (8), the
posterior distribution on the pmfs after collecting $n$ samples is
$\prod_{k}{\mathscr{D}}({\alpha}^{(k)}_{n+1})$ [Kotz et al.(2004)Kotz,
Balakrishnan, and Johnson], where

$\displaystyle{\alpha}^{(k,l)}_{n+1}$
$\displaystyle={\alpha}^{(k,l)}_{n}+\mathbbm{1}_{(U_{n}=k)}\mathbbm{1}_{(Y_{n}=l)}.$
(9)

Further, since the pmf ${p}^{(k)}$ is Dirichlet distributed, the probability
${p}^{(k,l)}$ is distributed according to the Beta distribution. More
precisely, conditioned on information $I_{n}$, we have

$\displaystyle{p}^{(k,l)}\sim\mathrm{Beta}\left({\alpha}^{(k,l)}_{n},{\alpha}^{(k,0)}_{n}-{\alpha}^{(k,l)}_{n}\right),$
(10)

where ${\alpha}^{(k,0)}_{n}\doteq\sum_{l=1}^{L}{\alpha}^{(k,l)}_{n}$.

###  3.3 Upper Confidence Bounds

Let the number of times arm $k$ has been sampled until time $n$ be

$\displaystyle{T}^{(k)}_{n}=\sum_{\tau=1}^{n-1}{U}^{(k)}_{\tau}.$ (11)

Based on the posterior distribution of ${p}^{(k,l)}$, we can compute a lower
bound ${a}^{(k,l)}_{n}$ and an upper bound ${b}^{(k,l)}_{n}$ on ${p}^{(k,l)}$
such that

$\displaystyle{\mathbb{P}}^{\rho_{1}}[{p}^{(k,l)}>{b}^{(k,l)}_{n}\mid I_{n}]$
$\displaystyle\leq\delta_{n}/2$ (12)
$\displaystyle{\mathbb{P}}^{\rho_{1}}[{p}^{(k,l)}<{a}^{(k,l)}_{n}\mid I_{n}]$
$\displaystyle\leq\delta_{n}/2,$ (13)

where $\delta_{n}=\frac{\delta}{KLn(1+\log N)}$ and $\delta=\eta N^{-5/2}$.
Let us define an interval
${E}^{(k,l)}_{n}\doteq[{a}^{(k,l)}_{n},{b}^{(k,l)}_{n}]\cap{E}^{(k,l)}_{n-1}$.
This ensures that the interval ${E}^{(k,l)}_{n}$ is non-increasing in $n$.
Note that we can compute this interval based on the inverse cdf of the
posterior distribution, which is a Beta distribution.

For arm $k$, let ${q}^{(k)}_{n}$ and ${u}^{(k)}_{n}$ respectively be the
optimal solution and the optimum value of the following quadratic program

$\displaystyle\max_{q\in\Delta}\;\sum_{l}q_{l}(1-q_{l})$
$\displaystyle\mathrm{s.t.}\;q_{l}\in{E}^{(k,l)}_{n}\;\forall l.$ (14)

Notice that since ${E}^{(k,l)}_{n}$ is non-increasing in $n$, the upper bound
${u}^{(k)}_{n}$ is also non-increasing in $n$. Let us define the event
$\mathcal{E}$ as

$\displaystyle\mathcal{E}\doteq\left\\{{p}^{(k,l)}\in{E}^{(k,l)}_{n}\;\forall
k,l,n\right\\}.$ (15)

Note that because of the way the upper bound ${u}^{(k,l)}_{n}$ is defined in
(14), we have ${p}^{(k,l)}(1-{p}^{(k,l)})\leq{u}^{(k,l)}_{n}$ under the event
$\mathcal{E}$. Further, we have ${\mathbb{P}}^{\rho_{1}}[\mathcal{E}]\geq
1-\delta$ because of the union bound $\sum_{n}\delta_{n}\leq\delta$.

##  4 Strategies and Performance Bounds

The main idea is simply to sample the arm that has maximum
$\varphi({u}^{(k)}_{n},{T}^{(k)}_{n})$ at time $n$. The precise strategy is
stated as Algorithm 4.

Bayesian UCB n = 1 →N k = 1 →K Compute ${u}^{(k)}_{n}$ Assign
$U_{n}=\operatorname*{arg\,max}_{k}\varphi({u}^{(k)}_{n},{T}^{(k)}_{n})$
Sample arm $U_{n}$ to obtain $Y_{n}$ and update the posterior according to (9)
${{\hat{p}}}^{(k)}_{N}$

###  4.1 Regret Bounds

Using Algorithm 4, we can achieve nearly the same mean squared error as the
oracle would have achieved with the knowledge of underlying distributions. The
gap between the oracle and Algorithm 4 is characterized by the regret bound
below.

######  Theorem 4.1.

Using the sampling strategy in Algorithm 4, we have

$\displaystyle\mathcal{R}_{N}(g,\varrho)\leq\mathcal{O}\left(\frac{\ln(\eta)}{N^{\frac{3}{2}}}\right).$
(16)

Our methodology for deriving the regret bound in Theorem 4.1 is adapted from
that in [Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi]. There are two
key distinctions from the approach in [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi]: (i) the gap ${e}^{(k)}_{n}$ between the upper bound
${u}^{(k)}_{n}$ and the true tracking parameter ${c}^{(k)}$, and (ii) the use
of local averaging with respect to the distribution $\varrho$. These
distinctions arise due to the Bayesian nature of the confidence bounds in
(14). In [Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi], the gap
${e}^{(k)}_{n}$ is bounded by a simple closed form expression (see Lemma 3 in
[Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi]) which is then used to
derive the regret bound. Characterizing this gap in our Bayesian setting is
non-trivial due to the relatively complicated construction of the upper bound
in (14). In the following two lemmas, we first derive a bound on
${e}^{(k)}_{n}$.

######  Lemma 4.2.

The length of the interval ${E}^{(k,l)}_{n}$ satisfies

$\displaystyle\mathrm{len}\left({E}^{(k,l)}_{n}\right)\leq\sqrt{\frac{2\ln{\frac{2}{\delta_{n}}}}{({\alpha}^{(k,0)}_{n}+1)}}\leq\sqrt{\frac{2\ln{\frac{2}{\delta_{n}}}}{({T}^{(k)}_{n}+1)}}.$
(17)

######  Proof 4.3.

See Appendix A.

######  Lemma 4.4.

Under event $\mathcal{E}$, we have

$\displaystyle{e}^{(k)}_{n}\doteq{u}^{(k)}_{n}-{c}^{(k)}\leq\sqrt{\frac{8\ln{\frac{2}{\delta_{n}}}}{({T}^{(k)}_{n}+1)}}.$
(18)

######  Proof 4.5.

See Appendix B.

Note that the bound here is as tight as the one in Lemma 3 of [Shekhar et
al.(2020)Shekhar, Ghavamzadeh, and Javidi]. These lemmas address the first
aforementioned distinction from [Shekhar et al.(2020)Shekhar, Ghavamzadeh, and
Javidi]. Using Lemma 4.4, we can follow the steps of the proof of Theorem 1 in
[Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi] to complete the proof
of Theorem 4.1. We refer the reader to Appendix C for the full proof with
details related to local averaging.

##  5 Numerical Results

We consider an experimental setup with two arms ($K=2$) and a binary support
($L=2$). The pmfs associated with arms 1 and 2 are $[0.99,0.01]$ and
$[0.7,0.3]$. The total number of samples $N=2500$ and the value of $\eta=1/N$.
Figure 1 depicts the regret $\mathcal{R}_{N}$ for the sampling strategy in
[Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi] and our sampling
strategy in Algorithm 4. We will refer to the former strategy as UCB and the
latter as Bayesian UCB. There is a clear improvement in performance in terms
of regret. We also plot the average number of samples acquired by both
strategies in Figure 1. We observe that the number of samples acquired by our
Bayesian UCB approach is closer to the oracle allocation.

[][b] [][b] [][b]

Figure 1: The plot on the left depicts the regret $\mathcal{R}_{N}$ for the
strategy in Algorithm 4 vs. the strategy in [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi]. The plot in the middle depicts the upper confidence
bounds ${u}^{(k)}_{n}$. The plot on the right depicts the expected number of
samples obtained and includes the empirical Bernstein approach in [Carpentier
et al.(2011)Carpentier, Lazaric, Ghavamzadeh, Munos, and Auer] as well.

Recall that the main difference between the UCB and the Bayesian UCB
strategies is the way the upper bound ${u}^{(k,l)}_{n}$ is computed. In both
these approaches, we observe that for every component, the value of the
tracking function ${u}^{(k)}_{N+1}/{T}^{(k)}_{N+1}$ is nearly the same for all
arms $k$.

######  Proposition 5.1.

For every arm $k$, we have

$\displaystyle\frac{\sum_{k}{u}^{(k)}_{N+1}}{N}\times\frac{1}{f_{N}}\leq\frac{{u}^{(k)}_{N+1}}{{T}^{(k)}_{N+1}}$
$\displaystyle\leq\frac{\sum_{k}{u}^{(k)}_{N+1}}{N}\times f_{N},$ (19)

where $f_{N}=1+o(1)$. Therefore, ${T}^{(k)}_{N+1}\propto{u}^{(k)}_{N+1}$,
approximately.

In Figure 1, we plot a typical realization of the upper bounds ${u}^{(k)}_{n}$
as a function of time $n$. Notice that the bounds computed in our approach are
significantly tighter. This tightness can be attributed to two factors: (i)
confidence intervals on the probabilities ${p}^{(k,l)}$ computed using the
inverse cdf of the posterior distribution, and (ii) the quadratic program for
computing ${u}^{(k)}_{n}$ in (14). The bounds computed using the method in
[Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi] do capture the
_difference_ between the tracking parameters for arms 1 and 2. However,
according to Proposition 5.1, it is the _ratio_ between the upper bounds that
matters. The ratio between the upper bounds in [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi] is close to 1 even for $N=2500$. While this ratio
will approach ${c}^{(1)}/{c}^{(2)}$ eventually, it is quite slow in practice.
While our Bayesian approach provides tighter bounds, there is a computational
cost associated with it, especially because of the diminishing probability
$\delta_{n}$. For the binary support case, there are simpler closed form
approximations [Bayarri and Berger(2004), Brown et al.(2001)Brown, Cai, and
DasGupta] that can be used to compute ${a}^{(k,l)}_{n}$ and ${b}^{(k,l)}_{n}$.

We propose to use our sampling strategy for tracking SARS-CoV-2 seroprevalence
in populations classified on the basis of location, ethnicity, age etc. Each
category of interest can be modeled as an arm with a Bernoulli pmf, where the
probability ${p}^{(k,2)}$ represents the positivity associated with category
$k$. We use data from a randomized seroprevalence survey conducted in Los
Angeles county [Sood et al.(2020)Sood, Simon, Ebner, Eichner, Reynolds,
Bendavid, and Bhattacharya]. The plots in Figures 2 and 2 depict the gap
between optimal sample allocation and the random allocation strategy used by
the survey. This gap illustrates the importance of adaptive randomized
sampling for obtaining uniform estimates of seroprevalence. Without careful
sampling, we can have bad estimates for certain groups that may be
particularly vulnerable.

[][t] [][t]

Figure 2: Population classified on the basis of ethnicity and location,
respectively in plots 2 and 2. The plots depict the number of sampled obtained
in each category under each strategy. Here, _Actual_ refers to the number of
samples obtained by the survey [Sood et al.(2020)Sood, Simon, Ebner, Eichner,
Reynolds, Bendavid, and Bhattacharya]

######  Remark 5.2.

For estimating SARS-CoV-2 seroprevalence, there are some practical aspects
that require a slightly different formulation and sampling strategy. We
provide a heuristic approach for addressing these issues and the analysis of
this approach is a problem for future work.

##  6 Conclusions and Future Work

For the problem of uniformly estimating pmfs in a mean squared sense, a UCB
sampling approach is proposed. Unlike the state-of-the-art approach in
[Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi], a Dirichlet prior is
assumed on the pmfs and the upper confidence bounds are computed based on the
posterior. It is analytically shown that the performance of this Bayesian
approach can be now worse than the approach in [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi], and has noticeably better performance in practice.
We discuss the potential application of the sampling approach for SARS-CoV-2
seroprevalence and incidence rate estimation.

There are two potential directions for future work. One is to incorporate
temporal variation in the underlying pmfs. This is particularly important in
studying SARS-CoV-2 incidence rates. The other is to incorporate a prior on
the collection pmfs so that the correlation between them can be exploited to
improve our estimates and the sampling strategy. For both these directions, a
Bayesian model can be more helpful.

##  Appendix A Proof of Lemma 4.2

A Beta distribution with parameters $\alpha,\beta$ is sub-Gaussian with
parameter $\frac{1}{4(\alpha+\beta+1)}$ [Marchal et al.(2017)Marchal, Arbel,
et al.]. And because of (10) we can conclude using the Chernoff bound
[Ross(2014)] that

$\displaystyle{\mathbb{P}}^{\rho_{n}}\left[{p}^{(k,l)}>\frac{{\alpha}^{(k,l)}_{n}}{{\alpha}^{(k,0)}_{n}}+\sqrt{\frac{\ln\frac{2}{\delta_{n}}}{2({\alpha}^{(k,0)}_{n}+1)}}\right]$
$\displaystyle\leq\delta_{n}$ (20)
$\displaystyle{\mathbb{P}}^{\rho_{n}}\left[{p}^{(k,l)}<\frac{{\alpha}^{(k,l)}_{n}}{{\alpha}^{(k,0)}_{n}}-\sqrt{\frac{\ln\frac{2}{\delta_{n}}}{2({\alpha}^{(k,0)}_{n}+1)}}\right]$
$\displaystyle\leq\delta_{n}.$ (21)

Since ${a}^{(k,l)}_{n}$ and ${b}^{(k,l)}_{n}$ are selected using the inverse
cdf of the posterior distribution, we have

$\displaystyle{b}^{(k,l)}_{n}-{a}^{(k,l)}_{n}\leq
2\sqrt{\frac{\ln\frac{2}{\delta_{n}}}{2({\alpha}^{(k,0)}_{n}+1)}}.$ (22)

The second inequality in the lemma simply follows from the fact the
${\alpha}^{(k,0)}_{n}\geq{T}^{(k)}_{n}$.

##  Appendix B Proof of Lemma 4.4

Under event $\mathcal{E}$, we have

$\displaystyle|{u}^{(k)}_{n}-{c}^{(k)}|$ (23)
$\displaystyle=\left|\sum_{l}\left(({p}^{(k,l)})^{2}-({q}^{(k,l)}_{n})^{2}-{p}^{(k,l)}+{q}^{(k,l)}_{n}\right)\right|$
(24) $\displaystyle\stackrel{{\scriptstyle
a}}{{\leq}}\sum_{l}\left(\left|({p}^{(k,l)})^{2}-({q}^{(k,l)}_{n})^{2}\right|\right)$
(25)
$\displaystyle=\sum_{l}\left|{p}^{(k,l)}-{q}^{(k,l)}_{n}\right|\left({p}^{(k,l)}+{q}^{(k,l)}_{n}\right)$
(26)
$\displaystyle\leq\sum_{l}\mathrm{len}\left({E}^{(k,l)}_{n}\right)\left({p}^{(k,l)}+{q}^{(k,l)}_{n}\right)$
(27) $\displaystyle\stackrel{{\scriptstyle
b}}{{\leq}}2\sqrt{\frac{\ln{\frac{2}{\delta_{n}}}}{2({T}^{(k)}_{n}+1)}}\sum_{l}\left({p}^{(k,l)}+{q}^{(k,l)}_{n}\right)$
(28)
$\displaystyle=4\sqrt{\frac{\ln{\frac{2}{\delta_{n}}}}{2({T}^{(k)}_{n}+1)}}.$
(29)

Here, inequality in $(a)$ is because of the triangle inequality and the
inequality in $(b)$ is due to Lemma 4.2.

##  Appendix C Proof of Theorem 4.1

We have

$\displaystyle\delta\geq{\mathbb{P}}^{\rho_{1}}[\mathcal{E}^{c}]$
$\displaystyle={\mathbb{E}}^{\rho_{1}}\left[\delta_{p}\right]$ (30)
$\displaystyle\geq{\mathbb{E}}^{\rho_{1}}\left[\delta_{p}\mathbbm{1}_{\eta(p)}(p)\right]$
(31) $\displaystyle={\mathbb{E}}^{\varrho}\left[\delta_{p}\right]\eta,$ (32)

where $\delta_{p}\doteq{\mathbb{P}}^{g}[\mathcal{E}^{c}\mid p]$. Note that the
last equality in the display above follows from a simple change of measure
argument. Since $\delta_{p}$ is the probability of the unfavorable event
$\mathcal{E}^{c}$ conditioned on the model $p$, we can obtain a regret bound
in terms of $\delta_{p}$ without local averaging using the same arguments as
in the proofs of Lemma 1 and Theorem 1 in [Shekhar et al.(2020)Shekhar,
Ghavamzadeh, and Javidi, Appendix D.3]. Thus, we have

$\displaystyle{\mathcal{L}}^{(k)}_{N}($
$\displaystyle{p}^{(k)},g)-\varphi^{*}(p,N)$ (33)
$\displaystyle\leq\frac{(K+5)LM}{(\lambda_{\mathrm{min}})^{2}N^{3/2}}+\frac{\delta_{p}}{(\lambda_{\mathrm{min}})^{2}N}\left(1+6M\sqrt{N}\right)+\frac{6(K-1)M^{2}}{(\lambda_{\mathrm{min}})^{3}N^{2}}+L\delta_{p},$
(34)

where

$\displaystyle M$
$\displaystyle\doteq\frac{\lambda_{\mathrm{max}}\sqrt{8\log(2/\delta_{N})}}{\lambda_{\mathrm{min}}C}=\mathcal{O}\left(\sqrt{\log(\eta
N)}\right)$ (35) $\displaystyle\lambda_{\mathrm{max}}$
$\displaystyle\doteq\min_{k,\pi\in\eta(p)}\frac{{c}^{(k)}}{\sum_{i}{c}^{(i)}};\quad\lambda_{\mathrm{max}}\doteq\max_{k,\pi\in\eta(p)}\frac{{c}^{(k)}}{\sum_{i}{c}^{(i)}}$
(36) $\displaystyle C$ $\displaystyle\doteq\sum_{i}{c}^{(i)},$ (37)

and ${c}^{(i)}=\sum_{l}{\pi}^{(i,l)}(1-{\pi}^{(i,l)})$. Using (32), we can
then locally average the regret as

$\displaystyle{\mathbb{E}}^{{\varrho}}\left[{\mathcal{L}}^{(k)}_{N}({\pi}^{(k)},g)-\varphi^{*}(p,N)\right]$
(38)
$\displaystyle\leq\frac{(K+5)LM}{(\lambda_{\mathrm{min}})^{2}N^{3/2}}+\frac{\delta}{(\eta\lambda_{\mathrm{min}})^{2}N}\left(1+6M\sqrt{N}\right)+\frac{6(K-1)M^{2}}{(\lambda_{\mathrm{min}})^{3}N^{2}}+L\frac{\delta}{\eta}$
(39) $\displaystyle=\mathcal{O}(\log(\eta)N^{-3/2}).$ (40)

##  Appendix D Incorporating Priors

In this section, we will consider a scenario where we are provided with a
prior $\rho_{1}$ on $p$. We will first modify the distribution $\varrho$ used
for local averaging as

$\displaystyle{\varrho}(\mathcal{A})=\rho_{1}\left(\mathcal{A}\cap\eta({p})\right)/\eta,$
(41)

for any Borel measurable set $\mathcal{A}\subseteq\Delta^{K}$. Here, $\eta(p)$
is a ball around $p$ with $\rho_{1}(\eta(p))=\eta$.

#### Dirichlet priors

When the prior $\rho_{1}$ is any factored Dirichlet distribution (_i.e_. not
necessarily uniform) of the form in (8), one can easily extend all the results
with the modified $\varrho$.

#### Intervals for Bernoulli parameters

Another scenario of interest is when $L=2$ and we are given that the Bernoulli
parameter associated with ${p}^{(k)}$ lies in some interval
$[{\gamma}^{(k)}_{l},{\gamma}^{(k)}_{u}]$. We can then assume that the
Bernoulli parameter ${p}^{(k)}$ is distributed uniformly and independently
over the interval $[{\gamma}^{(k)}_{l},{\gamma}^{(k)}_{u}]$. Since this prior
does not correspond to a Beta distribution, the posterior belief may not be a
Beta distribution. However, we still can easily find the cdf of the posteriors
in this case because the posterior on the parameters is just going to be a
truncated Beta distribution. Let ${F}^{(k)}_{n}$ be the cdf of the posterior
at time $n$ with a uniform prior and let ${\tilde{F}}^{(k)}_{n}$ be the
posterior with the truncated prior. Then one can show that

$\displaystyle{\tilde{F}}^{(k)}_{n}(x)=\frac{{F}^{(k)}_{n}(x)-{F}^{(k)}_{n}({\gamma}^{(k)}_{l})}{{F}^{(k)}_{n}({\gamma}^{(k)}_{u})-{F}^{(k)}_{n}({\gamma}^{(k)}_{l})}.$
(42)

The confidence intervals in (12) and (13) can then be computed using the
inverse of ${\tilde{F}}^{(k)}_{n}$.

##  Appendix E Practical Considerations: Overall Estimate and Batch Sampling

In seroprevalence estimation, we generally allocate samples in batches of size
$B$. Also, we are generally interested in estimating the positivity of each
category _as well as_ the positivity in the overall population. Let the
fraction of individuals of category $k$ in the overall population be
${w}^{(k)}$. Then the overall positivity $r$ and its estimate $\hat{r}$ are
given by

$\displaystyle r$ $\displaystyle=\sum_{k}{w}^{(k)}{p}^{(k)}$ (43)
$\displaystyle\hat{r}_{N}$
$\displaystyle=\sum_{k}{w}^{(k)}{\hat{p}}^{(k)}_{N}.$ (44)

The mean squared error between $r$ and $\hat{r}_{N}$ is given by

$\displaystyle\mathrm{MSE}(r,\hat{r}_{N})=\sum_{k}\frac{({w}^{(k)})^{2}{c}^{(k)}}{{T}^{(k)}}.$
(45)

If the mean squared error associated with $r$ is not considered, then it may
so happen that a tiny group (small $w_{k}$) with high positivity will be
allocated too many samples. The contribution of this small group to the
overall estimate would be small and thus, allocating too many samples to it
could compromise the quality of the overall estimate $r$. Therefore, we need
to determine an allocation that accounts for the quality of the overall
estimate $r$ as well.

A suitable way to formalize this notion is to pose the following constraints
on the oracle allocation ${T}^{(k)}$

$\displaystyle\frac{{c}^{(k)}}{{T}^{(k)}}\leq{\theta}^{(k)},\;k=1,\dots,K$
(C1)
$\displaystyle\sum_{k}\frac{({w}^{(k)})^{2}{c}^{(k)}}{{T}^{(k)}}\leq{\theta}^{(0)}$
(46) $\displaystyle\sum_{k}{T}^{(k)}\leq N$ (47) $\displaystyle{T}^{(k)}\geq
0,\;k=1,\dots,K,$ (48)

where ${\theta}^{(k)}$, $k=0,\dots,K$ are predetermined constants. A solution
to the above set of constraints can be obtained by solving the following
optimization problem

$\displaystyle\min_{{T}^{(1)},\dots,{T}^{(K)}}\quad$
$\displaystyle\max\left\\{\max_{k}\left\\{\frac{{c}^{(k)}}{{\theta}^{(k)}{T}^{(k)}}\right\\},\sum_{k}\frac{({w}^{(k)})^{2}{c}^{(k)}}{{\theta}^{(0)}{T}^{(k)}}\right\\}$
(C2) $\displaystyle\mathrm{s.t.}\quad$ $\displaystyle\sum_{k}{T}^{(k)}\leq N$
(49) $\displaystyle{T}^{(k)}\geq 0,\;k=1,\dots,K.$ (50)

The Problem (C1) is feasible if and only if the optimum value of the
optimization problem above is less than or equal to 1. In that case, the
solution to Problem (C2) is a solution to Problem (C1). Notice that Problem
(C1) is very similar to Problem (P1) except for the additional mean squared
error term associated with the overall estimate $r$. Because of this
distinction, it is not clear whether one can view Problem (C2) as a particular
instance of Problem (P1) with appropriate modifications and simply apply the
adaptive sampling strategy in Section 4 to the modified problem. Nonetheless,
we provide a similar heuristic sampling approach that tracks the oracle quite
well (See Figure 3). At each time time $n$, the heuristic is to allocate
samples within each batch according to the solution of the following
optimization problem

$\displaystyle\min_{{\tau}^{(1)},\dots,{\tau}^{(K)},\lambda}\quad$
$\displaystyle\lambda$ (C3) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\frac{u_{n}^{(k)}}{T_{n}^{(k)}+{\tau}^{(k)}}\leq{\theta}^{(k)}\lambda$
(51)
$\displaystyle\sum_{k}\frac{({w}^{(k)})^{2}u_{n}^{(k)}}{T_{n}^{(k)}+{\tau}^{(k)}}\leq{\theta}^{(0)}\lambda$
(52) $\displaystyle\sum_{k}{\tau}^{(k)}=B$ (53) $\displaystyle{\tau}^{(k)}\geq
0\;\forall k.$ (54)

Figure 3: This plot represents the number of samples collected by the
seroprevalence survey in [Sood et al.(2020)Sood, Simon, Ebner, Eichner,
Reynolds, Bendavid, and Bhattacharya], the oracle allocation in Problem (C1)
with appropriate constants ${\theta}^{(k)}$, and the allocation by the
heuristic (denoted by Adaptive). Notice that the heuristic tracks the oracle
closely. This plot tells us that in order get a good overall estimate, we
should allocate fewer samples to an underrepresented group like the Pacific
Islanders (${w}^{(k)}\approx 0.003$) than the number suggested by the oracle
(P1) (See 2). However, it also tells us that we can allocate substantially
more samples than the number in the survey ($\approx{w}^{(k)}N$) to this group
for a better estimate of their positivity without comprising the quality of
the overall estimate.

This research was supported, in part, by National Science Foundation under
Grant NSF CCF-1817200, CCF-1718560, CPS-1446901, Grant ONR N00014-15-1-2550,
and Grant ARO W911NF1910269.

## References

  * [Bayarri and Berger(2004)] M Jésus Bayarri and James O Berger.  The interplay of bayesian and frequentist analysis.  _Statistical Science_ , pages 58–80, 2004.
  * [Brown et al.(2001)Brown, Cai, and DasGupta] Lawrence D Brown, T Tony Cai, and Anirban DasGupta.  Interval estimation for a binomial proportion.  _Statistical science_ , pages 101–117, 2001.
  * [Carpentier et al.(2011)Carpentier, Lazaric, Ghavamzadeh, Munos, and Auer] Alexandra Carpentier, Alessandro Lazaric, Mohammad Ghavamzadeh, Rémi Munos, and Peter Auer.  Upper-confidence-bound algorithms for active learning in multi-armed bandits.  In _International Conference on Algorithmic Learning Theory_ , pages 189–203. Springer, 2011.
  * [Efroni et al.(2020)Efroni, Mannor, and Pirotta] Yonathan Efroni, Shie Mannor, and Matteo Pirotta.  Exploration-exploitation in constrained mdps.  _arXiv preprint arXiv:2003.02189_ , 2020.
  * [Kotz et al.(2004)Kotz, Balakrishnan, and Johnson] Samuel Kotz, Narayanaswamy Balakrishnan, and Norman L Johnson.  _Continuous multivariate distributions, Volume 1: Models and applications_.  John Wiley & Sons, 2004.
  * [Marchal et al.(2017)Marchal, Arbel, et al.] Olivier Marchal, Julyan Arbel, et al.  On the sub-gaussianity of the beta and dirichlet distributions.  _Electronic Communications in Probability_ , 22, 2017.
  * [Rosenberg and Mansour(2019)] Aviv Rosenberg and Yishay Mansour.  Online convex optimization in adversarial markov decision processes.  In _International Conference on Machine Learning_ , pages 5478–5486, 2019.
  * [Ross(2014)] Sheldon M Ross.  _Introduction to probability models_.  Academic press, 2014.
  * [Shekhar et al.(2020)Shekhar, Ghavamzadeh, and Javidi] Shubhanshu Shekhar, Mohammad Ghavamzadeh, and Tara Javidi.  Adaptive sampling for estimating probability distributions.  In _International Conference on Machine Learning_ , 2020.
  * [Sood et al.(2020)Sood, Simon, Ebner, Eichner, Reynolds, Bendavid, and Bhattacharya] Neeraj Sood, Paul Simon, Peggy Ebner, Daniel Eichner, Jeffrey Reynolds, Eran Bendavid, and Jay Bhattacharya.  Seroprevalence of sars-cov-2–specific antibodies among adults in los angeles county, california, on april 10-11, 2020.  _Jama_ , 2020.
  * [Tarbouriech and Lazaric(2019)] Jean Tarbouriech and Alessandro Lazaric.  Active exploration in markov decision processes.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 974–982, 2019.


∎

corrauthor _correspondence_<EMAIL_ADDRESS>alsoMEPHI _also at:_ NRNU
MEPhI, Moscow, Russia alsoLev _also at:_ Moscow Inst. of Physics and
Technology, Russia nowIKZ _present address:_ Leibniz-Institut für
Kristallzüchtung, Berlin, Germany alsoDubna _also at:_ Dubna State University,
Dubna, Russia nowMacolino _present address:_ LAL, CNRS/IN2P3, Université
Paris-Saclay, Orsay, France 11institutetext: INFN Laboratori Nazionali del
Gran Sasso and Gran Sasso Science Institute, Assergi, Italy 22institutetext:
INFN Laboratori Nazionali del Gran Sasso and Università degli Studi
dell’Aquila, L’Aquila, Italy 33institutetext: INFN Laboratori Nazionali del
Sud, Catania, Italy 44institutetext: Institute of Physics, Jagiellonian
University, Cracow, Poland 55institutetext: Institut für Kern- und
Teilchenphysik, Technische Universität Dresden, Dresden, Germany
66institutetext: Joint Institute for Nuclear Research, Dubna, Russia
77institutetext: European Commission, JRC-Geel, Geel, Belgium 88institutetext:
Max-Planck-Institut für Kernphysik, Heidelberg, Germany 99institutetext:
Department of Physics and Astronomy, University College London, London, UK
1010institutetext: Dipartimento di Fisica, Università Milano Bicocca, Milan,
Italy 1111institutetext: INFN Milano Bicocca, Milan, Italy 1212institutetext:
Dipartimento di Fisica, Università degli Studi di Milano and INFN Milano,
Milan, Italy 1313institutetext: Institute for Nuclear Research of the Russian
Academy of Sciences, Moscow, Russia 1414institutetext: Institute for
Theoretical and Experimental Physics, NRC “Kurchatov Institute”, Moscow,
Russia 1515institutetext: National Research Centre “Kurchatov Institute”,
Moscow, Russia 1616institutetext: Max-Planck-Institut für Physik, Munich,
Germany 1717institutetext: Physik Department, Technische Universität München,
Germany 1818institutetext: Dipartimento di Fisica e Astronomia, Università
degli Studi di Padova, Padua, Italy 1919institutetext: INFN Padova, Padua,
Italy 2020institutetext: Physikalisches Institut, Eberhard Karls Universität
Tübingen, Tübingen, Germany 2121institutetext: Physik-Institut, Universität
Zürich, Zurich, Switzerland

# Calibration of the Gerda experiment

The Gerda collaborationcorrauthor
M. AgostiniUCL,TUM  G. AraujoUZH  A.M. BakalyarovKU  M. BalataALNGS  I.
BarabanovINRM  L. BaudisUZH  C. BauerHD  E. BellottiMIBF,MIBINFN  S.
BelogurovITEP,INRM,alsoMEPHI  A. BettiniPDUNI,PDINFN  L. BezrukovINRM  V.
BiancacciPDUNI,PDINFN  E. BossioTUM  V. BotheHD  V. BrudaninJINR  R.
BrugneraPDUNI,PDINFN  A. CaldwellMPIP  C. CattadoriMIBINFN  A.
ChernogorovITEP,KU  T. ComellatoTUM  V. D’AndreaAQU  E.V. DemidovaITEP  N. Di
MarcoALNGS  E. DoroshkevichINRM  F. FischerMPIP  M. FominaJINR  A.
GangapshevINRM,HD  A. GarfagniniPDUNI,PDINFN  C. GoochMPIP  P. GrabmayrTUE  V.
GurentsovINRM  K. GusevJINR,KU,TUM  J. HakenmüllerHD  S. HemmerPDINFN  R.
HillerUZH  W. HofmannHD  J. HuangUZH  M. HultGEEL  L.V. InzhechikINRM,alsoLev
J. Janicskó CsáthyTUM,nowIKZ  J. JochumTUE  M. JunkerALNGS  V. KazalovINRM  Y.
KermaïdicHD  H. KhushbakhtTUE  T. KihmHD  I.V. KirpichnikovITEP  A.
KlimenkoHD,JINR,alsoDubna  R. KneißlMPIP  K.T. KnöpfleHD  O. KochetovJINR
V.N. KornoukhovITEP,INRM  P. KrauseTUM  V.V. KuzminovINRM  M. LaubensteinALNGS
M. LindnerHD  I. LippiPDINFN  A. LubashevskiyJINR  B. LubsandorzhievINRM  G.
LutterGEEL  C. MacolinoALNGS,nowMacolino  B. MajorovitsMPIP  W. ManeschgHD  L.
ManzanillasMPIP  M. MiloradovicUZH  R. MingazhevaUZH  M. MisiaszekCR  P.
MoseevINRM  Y. MüllerUZH  I. NemchenokJINR,alsoDubna  L. PandolaCAT  K.
PelczarGEEL,CR  L. PertoldiPDUNI,PDINFN  P. PiseriMILUINFN  A. PulliaMILUINFN
C. RansomUZH  L. RauscherTUE  S. RiboldiMILUINFN  N. RumyantsevaKU,JINR  C.
SadaPDUNI,PDINFN  F. SalamidaAQU  S. SchönertTUM  J. SchreinerHD  M. SchüttHD
A-K. SchützTUE  O. SchulzMPIP  M. SchwarzTUM  B. SchwingenheuerHD  O.
SelivanenkoINRM  E. ShevchikJINR  M. ShirchenkoJINR  L. ShtembariMPIP  H.
SimgenHD  A. SmolnikovHD,JINR  D. StukovKU  A.A. VasenkoITEP  A.
VeresnikovaINRM  C. VignoliALNGS  K. von SturmPDUNI,PDINFN  T. WesterDD  C.
WiesingerTUM  M. WojcikCR  E. YanovichINRM  B. ZatschlerDD  I. ZhitnikovJINR
S.V. ZhukovKU  D. ZinatulinaJINR  A. ZschockeTUE  A.J. ZsigmondMPIP  K.
ZuberDD  and G. ZuzelCR

###### Abstract

The GERmanium Detector Array (Gerda) collaboration searched for neutrinoless
double-$\beta$ decay in 76Ge with an array of about 40 high-purity
isotopically-enriched germanium detectors. The experimental signature of the
decay is a monoenergetic signal at $Q_{\beta\beta}$$=2039.061(7)$ keV in the
measured summed energy spectrum of the two emitted electrons. Both the energy
reconstruction and resolution of the germanium detectors are crucial to
separate a potential signal from various backgrounds, such as neutrino-
accompanied double-$\beta$ decays allowed by the Standard Model. The energy
resolution and stability were determined and monitored as a function of time
using data from regular 228Th calibrations. In this work, we describe the
calibration process and associated data analysis of the full Gerda dataset,
tailored to preserve the excellent resolution of the individual germanium
detectors when combining data over several years.

††journal: Eur. Phys. J. C

##  1 Introduction

Neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay is a hypothetical,
second-order weak interaction process in which a nucleus changes its charge
number by two units with the emission of two electrons but without
accompanying anti-neutrinos. This lepton-number violating process is only
permitted if neutrinos are massive Majorana fermions, i.e. if there is a
Majorana mass term in the Lagrangian of the underlying theory. Such a term
appears in many extensions of the Standard Model of particle physics and could
explain why neutrino masses are much smaller than those of all other fermions
Mohapatra:2006gs . The GERmanium Detector Array (Gerda) collaboration searched
for the $0\nu\beta\beta$ decay of the isotope 76Ge with a $Q$-value of
${Q_{\beta\beta}=2039.061(7)}$ keV Mount:2010zz  by operating high-purity
germanium (HPGe) detectors isotopically enriched to $>$86% in 76Ge, making
them also the potential source of $0\nu\beta\beta$ decay.

We used three types of enriched germanium detectors: 30 broad energy germanium
(BEGe) detectors, 7 coaxial detectors, and 5 newer inverted coaxial (IC)
detectors. The BEGe detectors are smaller (average 0.7 kg) but offer superior
energy resolution and pulse shape discrimination (PSD) properties compared to
the coaxial detectors Agostini:2019mwn , while the IC detectors provide energy
resolution and PSD properties similar to the BEGe detectors Domula:2017mei
but with a larger mass (average 1.9 kg) comparable to that of the coaxial
detectors (average 2.3 kg), allowing for the easier design of larger germanium
arrays.

The array of germanium detectors was immersed in a cryostat filled with 64 m3
of liquid argon (LAr). The top of the cryostat and the surrounding water tank
houses a clean room containing a glove box and lock system for deploying the
HPGe detectors and calibration sources. The entire setup was located
underground at the Laboratori Nazionali del Gran Sasso (LNGS) of INFN, Italy,
and is described in detail in Agostini:2017hit .

The first phase of the experiment was operated with 18 kg of coaxial detectors
(inherited from the Heidelberg-Moscow Gunther:1997ai  and Igex Morales:1998hu
collaborations) between November 2011 and September 2013. Phase II started in
December 2015, after 20 kg of BEGe detectors produced for the Gerda experiment
were added and the liquid argon volume around the detector array was
instrumented with photosensors as a veto against radioactivity
Agostini:2017hit . During an upgrade in mid-2018, referred to as the Phase II
upgrade, IC detectors with a total mass of 9.6 kg were added, and the LAr
instrumentation was upgraded. Phase II ended in November 2019. While the
calibration procedure of Phase I data has been discussed in Ackermann:2012xja
; Agostini:2013mzu , the focus of this paper is the calibration of the Phase
II data.

In all recent $0\nu\beta\beta$ decay experiments, the signature of the rare
nuclear transition is a monoenergetic peak in the measured energy spectrum of
the two electrons at $Q_{\beta\beta}$. Consequently, a crucial parameter to
distinguish a signal from the background is the energy estimator. The better
the energy resolution of the detectors, the narrower the signal energy region
effectively becomes, and an excess over the continuous background can be more
clearly identified. One strength of HPGe detectors is their unparalleled
energy resolution (typically $\sigma$/E$\sim$0.1% at $Q_{\beta\beta}$). It
permits the almost complete rejection of background events from regular two-
neutrino-accompanied double-$\beta$ decays Abramov:2019hhh , an otherwise
irreducible background in $0\nu\beta\beta$ decay searches Elliott:2002xe ;
Maneschg:2017mzu .

Given the central role of the energy observable, adequate measures must be
taken to accurately determine the energy scale and resolution, monitor their
stability over the full data acquisition period, and determine the relevant
uncertainties entering the statistical analysis for the $0\nu\beta\beta$ decay
search. In Sect. 2 we detail the calibration procedure, while in Sect. 3 we
discuss the analysis of the calibration data and the energy scale
determination, including the procedures employed to monitor and maintain the
stability of the HPGe detectors over time. In Sect. 4 we describe the
determination of the energy resolution for the $0\nu\beta\beta$ decay
analysis, and in Sect. 5 we provide an evaluation of the associated
uncertainties. In Sect. 6, we discuss the determination of the residual energy
bias and its uncertainty. In Sect. 7, we compare the results from calibration
data with those in the physics data (data used for the $0\nu\beta\beta$ decay
search) for the resolutions of the lines from decays of 40K and 42K. We close
in Sect. 8 with a summary and a discussion of our main results.

##  2 Energy calibration process

To perform the calibrations, we regularly exposed the HPGe detectors to three
custom-made low-neutron emission 228Th calibration sources Baudis:2015sba ,
each with an activity of about 10 kBq. These sources were stored within
shielding above the lock system, at a vertical distance of at least 8 m to the
HPGe detector array, during physics data acquisition. Since 228Th has a half-
life of 1.9 yr, the sources were replaced during Phase II to ensure a
sufficient level of activity.

During calibration runs of the HPGe detectors, the 228Th sources were lowered
into the LAr cryostat to reach the level of the detector array by three source
insertion systems tarkathesis ; baudis2013monte . Each of these deploy a
single source, placed on tantalum absorbers (h = 60 mm, $\diameter$ = 32 mm)
During calibration, each source was placed at three different heights to
expose the detector array more homogeneously, and data were acquired at each
location for up to 30 min. With this exposure, typically around $(1-3)\times
10^{3}$ events are observed in the prominent 208Tl $\gamma$ line at 2614.5 keV
in a BEGe detector, and $(0.6-1)\times 10^{4}$ events in a coaxial or IC
detector. Calibration data were acquired every 7-10 days with a total of 142
calibration runs used for the analysis of the Phase II data.

The triggering energy threshold for this acquisition during calibration
corresponds to $\sim$400 keV. This threshold was set to include the strong
$\gamma$ line of the 228Th spectrum at 583.2 keV while keeping the event rate
at a manageable level for the data acquisition system. The detector signals
were read out with charge sensitive amplifiers, and digitized by a 100 MHz
14-bit flash analog-to-digital converter (FADC). As for physics data
acquisition, for each trigger a 160 $\upmu$s long waveform is recorded at a
sampling rate of 25 MHz, centered around the trigger time and covering an
energy range up to $\sim$6 MeV. During calibration, every 2 s a test pulse was
injected into the amplifier electronics of all germanium detectors to monitor
the stability of their gains. Between successive calibration runs, i.e. during
physics data acquisition, test pulsers were injected every 20 s for the same
purpose.

Data from the FADCs are transformed into the MGDO (ROOT-based) format
Brun:1997pa ; Agostini:2011nf  and processed to analyze properties of the
recorded waveforms using the Gelatio software Agostini:2011xe , as is the case
for physics data. The energy is estimated from the amplitude of the waveform
after applying a digital filter which reduces the impact of noise and thus
improves the resolution. As a fast first estimate for monitoring and cross-
check purposes, a pseudo-Gaussian filter is applied to obtain an energy
estimator as part of the online analysis gatti1986processing . An improved
energy resolution for the final data analysis is achieved with a Zero Area
Cusp (ZAC) filter Agostini:2015pta , which removes the effect of low-frequency
noise. This filter is optimized offline for each calibration run and HPGe
detector to minimise the resolution of the highest energy $\gamma$ line in the
228Th spectrum valerio:thesis .

A set of heuristic event selection criteria is applied to ensure that events
recorded during calibration are of a physical origin, and to reduce pile-up
events. The underestimation of energy for these events cause low-energy tails
in the spectra of $\gamma$ lines and can bias the estimated energy resolution.
These selection criteria are based on the properties of the waveform, such as
the baseline stability and slope, trigger time, number of triggered events,
and rise time of the pulse. The probability of rejecting physical
interactions, estimated with events from the regularly injected test pulses,
is below 0.1% lazzarothesis .

##  3 Analysis of energy spectra

Nuclei of the 228Th isotope decay in a chain via $\alpha$ and $\beta$ decays
to the stable 208Pb with the emission of multiple monoenergetic $\gamma$ rays.
These result in sharp peaks in the recorded energy spectra, as shown in Fig. 1
in the combined spectra of each detector type. The pattern of observed peaks
is used to identify the $\gamma$ lines and thereby determine their energy and
resolution.

Figure 1: Combined energy spectrum for 228Th calibration data for all enriched
detectors of BEGe, coaxial, and IC type during Phase II after rebinning to 3
keV. The inset shows the fit to the 2.6 MeV line in the spectrum of the
detector GD91A before the Phase II upgrade with 0.3 keV binning, with the
components of the fit drawn separately (linear and step backgrounds are
combined). The energies of the nine peaks that typically contribute to the
formation of calibration curves are labelled.

The TSpectrum class of ROOT is used to find peak positions in the uncalibrated
spectrum, such that all peaks with amplitudes exceeding 1/20 of the amplitude
of the most prominent peak are found. This threshold was chosen to avoid the
detection of spurious peaks. The peak with the highest energy is identified as
the full energy peak (FEP) of the $\gamma$ ray from the decay of 208Tl, a
daughter of 228Th, at $E_{\mathrm{FEP}}=2614.5$ keV. A preliminary calibration
for the energy estimator $T$ is applied assuming a linear energy scale without
offset:

$\displaystyle E_{0}(T)=\frac{E_{\mathrm{FEP}}}{T_{\mathrm{FEP}}}\cdot T\,.$
(1)

A candidate peak is confirmed if its preliminary estimated energy is
consistent within 6 keV with the energy of a known line in the 228Th spectrum.
The 6 keV value permits the accurate identification of peaks while allowing
for some non-linearity of the energy scale. The known peaks correspond to
$\gamma$ rays from isotopes in the 228Th decay chain with energies above 500
keV and branching ratios above 0.3%, including the detector specific single
escape peak (SEP) at 2103.5 keV and double escape peak (DEP) at 1592.5 keV
resulting from the 2.6 MeV $\gamma$ ray of 208Tldecays. In the context of this
paper, without ambiguity, FEP, SEP, and DEP always refer to those of 208Tl.
The double peak due to the 511.0 keV annihilation line and 510.7 keV $\gamma$
line from 208Tl is excluded from the analysis, in particular since the
resolution of the annihilation peak is broadened due to the Doppler effect
coldwell2016experimental .

###  3.1 Peak fitting and calibration curves

To determine the position $\mu$ and energy resolution in terms of $\textrm{the
full width at half maximum (FWHM)}=2.35\cdot\sigma$ of the identified peaks,
fits are performed locally in an energy window of $10-20$ keV around the peak
position obtained from the preliminary calibration. These are configured
manually and separately for each peak to avoid interference from neighbouring
peaks.

Minimally, a Gaussian $g(E)$ is used to model the peak, and a linear function
$f_{\mathrm{lin}}(E)$ is used to model the background:

$\displaystyle g(E)$
$\displaystyle=\frac{n}{\sqrt{2\pi}\sigma}\exp\quantity[-\frac{(E-\mu)^{2}}{2\sigma^{2}}],$
(2) $\displaystyle f_{\mathrm{lin}}(E)$ $\displaystyle=a+b\cdot E,$ (3)

where $n$, $\mu$ and $\sigma$ are the intensity, position, and width of the
peak, and $a$ and $b$ give the intercept and slope of the linear function,
respectively.

For high statistics peaks (583.2 keV, 727.3 keV, 763.5 keV, 860.5 keV, and
2614.5 keV), the SEP, and the DEP, a step function is used to model the flat
backgrounds occurring only above or below the peak from multiple Compton
scatters:

$\displaystyle
f_{\mathrm{step}}(E)=\frac{d}{2}\mathrm{erfc}\left(\frac{E-\mu}{\sqrt{2}\sigma}\right),$
(4)

where $d$ is the height of the step function, and erfc denotes the
complementary error function.

For the high statistics peaks as defined above, a low-energy tail is
additionally used to model the effects of incomplete charge collection and the
residual presence of pile-up events:

$\displaystyle
h(E)=\frac{c}{2\beta}\exp\left(\frac{E-\mu}{\beta}+\frac{\sigma^{2}}{2\beta^{2}}\right)\mathrm{erfc}\left(\frac{E-\mu}{\sqrt{2}\sigma}+\frac{\sigma}{\sqrt{2}\beta}\right),$
(5)

where $\beta$ and $c$ are the height and slope of the tail, respectively. An
example of the FEP peak fit is shown in the inset of Fig. 1.

Peaks are excluded after the fit if any of the following heuristic rules are
fulfilled: (i) the estimated FWHM is above 11 keV or below 1.5 keV; (ii) the
peak maximum is lower than 2.5 times the linear component of the background or
lower than 10 counts; (iii) the fitting error on the FWHM is larger than the
FWHM itself. These rules are purely heuristic and designed to remove peaks
that cannot be fitted well, mainly due to low statistics.

Typically around 5-8 peaks per detector survive all selection criteria. The
FEP is always identified, since the peak identification algorithm requires it.
In >80% of cases the lines at 583.2 keV, 860.5 keV, 1592.5 keV (DEP) and
2103.5 keV (SEP), and in (15-60)% the lines at 727.3 keV, 785.4 keV, 1078.6
keV and 1620.7 keV are found. All other $\gamma$ lines are seen in $<$3% of
the spectra from individual detectors of a single calibration run, due to
insufficient statistics.

From the obtained peak positions, we determine the calibration curve which is
a function to convert the uncalibrated energy estimator $T$ into a physical
energy in keV. We plot the peak positions in terms of the uncalibrated energy
estimator $T$ of identified peaks against their physical energies $E$
according to literature values litvals , and then fit with a linear function

$\displaystyle E(T)=p_{0}+p_{1}\cdot T.$ (6)

In the rare case where only the FEP is found and successfully fitted, the
resulting calibration curve has an intercept of zero. Such cases only occured
during periods of instability for a detector which were excluded from data
analysis (see Sect. 3.3). For most detectors, the linear calibration curve
describes the peak positions within a few tenths of a keV, as shown in Sect.
6. Discussion of a quadratic correction to the energy scale introduced for
five detectors after the Phase II upgrade can be found in Sect. 3.2.

Typically, a calibration curve is used to calibrate the physics data following
a calibration run. However, if a calibration is taken after changes in the
experimental setup, or instabilities in the detector array, the calibration
curves may be applied retrospectively to additionally calibrate the period
between the changes or instabilities and the calibration run. The unstable
period itself would not be used for physics analysis.

###  3.2 Quadratic correction

After the Phase II upgrade, several detectors (the new IC detectors and one
coaxial detector) exhibited larger residuals in their calibration curves
compared to other detectors, up to 2.5 keV at 1.5 MeV. Whether these effects
are directly related to change in cable routing during the Phase II upgrade,
is unknown. These effects could be largely accounted for by the incorporation
of a quadratic correction to the calibration curves,

$\displaystyle E(T)=m_{0}+m_{1}\cdot E_{0}(T)+m_{2}\cdot E_{0}^{2}(T),$ (7)

where $E_{0}$ is the energy estimator after the application of the linear
calibration curve as described in Sect. 3.1.

The parameters $m_{0}$, $m_{1}$ and $m_{2}$ were determined by fitting the
residuals of each detector’s calibration curves per calibration run. One such
example is shown in Fig. 2. These parameters were observed to be stable over
time, with the exception of a single jump for three IC detectors following
hardware activities in the clean room. Therefore a single correction for each
detector’s stable period was applied, using the average parameters of all
calibration runs in that period. After the quadratic correction, the remaining
residuals were within a few tenths of a keV.

Figure 2: Fitting the residuals of the calibration curve with a quadratic
function, as shown for detector ANG2 for the calibration on 15th October 2018.

###  3.3 Detector performance and stability

To consistently combine data over an extended period of time while preserving
the excellent energy resolution of the HPGe detectors, it is vital to monitor
the stability of the energy scale between calibrations and exclude periods
with significant shifts and fluctuations which would contribute to the width
of the peaks. As previously mentioned, test pulses are regularly injected into
the readout electronics to monitor the stability of the data acquisition
system. Their signal magnitude corresponds to an energy of about 3 MeV.
Periods with significant jumps or drifts (>1 keV) in the amplitude of the test
pulses are excluded from data analysis and a calibration is performed once the
detector stabilizes. The corresponding loss of exposure is at the few-percent
level. The origin of these shifts is largely unknown, but may be caused by
temperature changes in the electronics.

Additionally, we monitor the stability of the FEP position in the calibration
spectrum. If the position of this line changes by more than 1 keV between
successive calibrations without an identifiable reason (maintenance, longer
breaks, specific incident), the data of the respective detector are discarded
from the analysis for that period of time. The corresponding exposure loss is
at the few-percent level. Smaller or temporary drifts may still affect the
obtainable energy resolution and are discussed as a systematic uncertainty in
Sect. 5.2.

Due to hardware changes, the detectors may experience changes in their energy
resolution and energy scale over longer periods of time. To more accurately
reflect the properties of a detector at a certain time, for the final Gerda
analysis Agostini:2020xta  we divide the full data acquisition period for each
detector into stable sub-periods called partitions. The stability is judged
based on two parameters: the FWHM at the FEP and the residual at SEP. The
former reflects the changes in the detector resolution, while the latter
catches the changes in the energy bias at the energy peak closest to
$Q_{\beta\beta}$ (see Sect. 6 for more discussions on the bias). After the
Phase II upgrade and cable rerouting, the resolutions improved for most of the
detectors. Therefore, for simplicity, we start a new partition for all
detectors after the upgrade. There are one to four partitions for each
detector. The majority of the detectors have only two partitions, split at the
time of the Phase II upgrade. An example of the partitions is shown in Fig. 3.

Figure 3: FWHM of the FEP as a function of time for detector GD76B, one of the
BEGe detectors. Each data point comes from one calibration run. The full data
acquisition period is divided into three partitions, shown in solid circle
(blue), triangle (green), and diamond (red), respectively. The partition shown
in triangles is due to the Phase II upgrade and coincident improvement in
resolution. The partition shown in diamonds is due to the jump in resolution
in January 2019 when a hardware change took place.

##  4 Energy resolutions from the combined calibration spectra

Depending on the specific physics analysis, we calculated the energy
resolution either by partition, described in the previous section, or by
detector type. For the $0\nu\beta\beta$ decay search reported in
Agostini:2019hzm , where the data before the Phase II upgrade were analyzed,
an effective resolution for each detector type was employed. For the final
$0\nu\beta\beta$ decay search of Gerda reported in Agostini:2020xta , where
all Gerda data were analyzed, we calculated a resolution for each partition, a
much more fine-grained approach. At the expense of increased complexity, the
partition approach improves the physics result by capturing the variations
among the detectors as well as the variation over time.

Since both methods are applicable for Gerda and any other experiment with a
modular detector setup, here, we discuss both approaches. While the detector
type approach was replaced in favor of partitioning for the final
$0\nu\beta\beta$ decay search, the former gives an overview of the overall
detector performance. For this reason all the illustrations and calibration
parameters are provided by detector type. For simplicity, we refer to a
collection of detectors of the same type as a dataset.

###  4.1 By partition

To obtain the $\gamma$ line resolutions for each detector partition, we first
produce combined calibration spectra. The energy spectra obtained from each
calibration run within one partition are first normalised to account for
differing statistics, and then weighted according to the time span for which
the corresponding calibration curves were used to calibrate physics data. The
resulting $\gamma$ peaks in a combined spectrum will be representative of the
average performance of that detector in that partition.

The peak identification and fit procedure described in Sect. 3 is then applied
to each combined calibration spectrum.

The SEP is broadened due to the known Doppler effect and is thus excluded
coldwell2016experimental . We also observe broadening in the DEP. This is
hypothesised to originate due to events occurring more frequently in the outer
regions of the detectors and thus being more susceptible to incomplete charge
collection hull1996charge . This line is therefore excluded as well.

The dependence of the $\gamma$ line resolutions on the calibrated energy $E$
is then fitted with the function Agostini:2015pta

$\displaystyle\sigma(E)=\sqrt{a+bE},$ (8)

where $a$ and $b$ are fit parameters. The former accounts for the
contributions from electronic noise, while the latter accounts for statistical
fluctuations in the number of charge carriers. The resolution at
$\text{Q}_{\beta\beta}$ is then given by using $E=Q_{\beta\beta}$ in Eq. 8.

The resultant FWHM resolutions at $\text{Q}_{\beta\beta}$ of the partitions
vary between 2.3 keV and 8.8 keV, as shown in Fig. 4. Values for each
partition can be found in chloethesis . Systematic errors are calculated via a
dedicated study as explained in Sect. 5.

Figure 4: Distribution of FWHM resolution at $\text{Q}_{\beta\beta}$ per
detector partition. The detector partitions with resolutions > 6 keV are due
to two coaxial detectors whose resolutions degraded after the Phase II
upgrade.

###  4.2 By detector type

The appropriate method for calculating effective resolutions by detector type
depends on the specific application.

####  4.2.1 Background modeling

For background modeling, energy dependent resolutions are required, i.e.
resolution curves. To calculate these for datasets, the procedure is similar
to that for the partitions, though weighting is now required to combine the
resolutions from different detector partitions. When data from multiple
partitions are combined by adding their energy spectra, Gaussian peaks in the
individual spectra combine to become a Gaussian mixture, namely the sum of
multiple Gaussian distributions with different centroids and resolutions. The
resolution of individual partitions in a dataset is stable within a factor of
1.7 for BEGe and IC detectors. For coaxial detectors there is a slightly
higher fluctuation, but still within a factor of three.The variation in
position of the centroid is much smaller than the energy resolution, typically
around 0.2 keV. Therefore the shape of a peak in the combined energy spectra
remains approximately Gaussian and can be characterized by an effective
resolution, computed from the resolution of individual partitions.

The variance of a Gaussian mixture is given by:

$\displaystyle\mathrm{\sigma}^{2}=\sum_{i}w_{i}\left(\sigma_{i}^{2}+\mu_{i}^{2}\right)-\left(\sum_{i}w_{j}\cdot\mu_{i}\right)^{2},$
(9)

where the sum goes over Gaussians with means and standard deviations $\mu_{i}$
and $\sigma_{i}$, with weights $w_{i}$, representing the relative contribution
to expected peak counts of individual Gaussians gaussianmixturethesis .

For a dataset comprised of individual partitions, these parameters stand for
the individual partitions’ resolution $\sigma_{i}$, and peak position
$\mu_{i}$, which can be different due to independent systematic effects on the
energy scale. The weights are the expected relative event count contribution
of individual partitions. Since peak counts are proportional to exposure
$\mathcal{E}_{i}=m_{i}\cdot t_{i}$, with individual detector’s active mass
$m_{i}$ and live time $t_{i}$, the relative exposure contribution is:

$\displaystyle w_{i}=\frac{\mathcal{E}_{i}}{\mathcal{E}},$ (10)

where $\mathcal{E}=\sum_{j}\mathcal{E}_{j}$ is the total exposure of the
dataset.

Since the biases in the energy scale are small, we can neglect the differences
in the peak positions. Eq. 9 therefore simplifies to:

$\displaystyle\sigma=\sqrt{\frac{1}{\mathcal{E}}\sum_{i}\mathcal{E}_{i}\,\sigma_{i}^{2}},$
(11)

with total error $\delta_{\sigma}$ from the statistical fitting errors of
individual partition resolutions $\delta_{\sigma_{i}}$:

$\displaystyle\delta_{\sigma}=\sqrt{\frac{1}{\mathcal{E}^{2}\sigma^{2}}\sum_{i}(\mathcal{E}_{i}\sigma_{i}\delta_{\sigma_{i}})^{2}},$
(12)

with negligible uncertainty in the weights.

For instance, the simplified model of the FEP is a Gaussian with a mean of
2614.5 keV and a width fixed to the effective resolution (see Eq. 11) of the
dataset. On the other hand, a Gaussian mixture model would consist of the sum
of a Gaussian for each partition, each with its own resolution and centroid.
Fig. 5 shows the Gaussian mixture and simplified signal models for the IC and
coaxial datasets. For the IC and BEGe datasets, the Gaussian mixture model is
very close to a Gaussian shape, as the centroid differences are small and the
partitions in each dataset have similar resolutions. The resolutions among the
coaxial detectors are more varied and thus using a Gaussian signal model may
be less appropriate.

Figure 5: Comparison of simplified Gaussian signal model (dashed blue) to the
more detailed Gaussian mixture signal model (solid black) of the FEP, for
datasets formed of the partitions of BEGe (left), coaxial (middle) and IC
(right) detectors. Red lines show Gaussian shaped peaks for individual
partitions, which have been rescaled by a factor of 20/5/1 for the
BEGe/Coax/IC detectors for visibility.

To calculate the effective resolution curves for each dataset, first the
$\gamma$ line resolutions are obtained for each of the partitions as in Sect.
4.1. For all $\gamma$ lines whose resolution was reliably determined for all
partitions in that dataset, an effective resolution of the dataset at that
energy is calculated using Eq. 11. All other lines which were missing in at
least one detector partition are excluded.

Once the effective resolutions for each energy and dataset have been
determined by weighting partition resolutions with Eq. 11, their energy
dependence is fitted with Eq. 8.

The obtained effective resolutions and functions of the three detector types
are shown in Fig. 6 and Tab. 1. The statistical errors are obtained from the
fit.

Detector type | $a$ | $b$
---|---|---
[keV2] | [$10^{-4}$ keV]
BEGe | 0.551(1) | 4.294(9)
coaxial | 0.985(2) | 10.73(2)
IC | 0.280(2) | 5.83(2)
Table 1: Parameters of resolution curves (Eq. 8) obtained for datasets of each
detector type. Figure 6: Effective resolution curves for BEGe (purple),
coaxial (blue) and IC (green) datasets. Open points indicate broadened lines
not used to form the resolution curves, namely the double- and single-escape
peaks of the 2.6 MeV line due to 208Tl decay. Square markers indicate the
exposure-weighted resolutions of the lines in physics data due to 40K (1460.8
keV) and 42K (1524.7 keV) decays.

####  4.2.2 $0\nu\beta\beta$ decay search

As mentioned before, in earlier Gerda $0\nu\beta\beta$ decay analyses such as
Agostini:2019hzm , partitioning was not performed, and data from multiple
detectors were combined to form a dataset for each detector type. In the case
of the Gerda Phase II data, very few events (in fact, only one) are observed
close to $\text{Q}_{\beta\beta}$, so using a signal model of a Gaussian with
an effective resolution as in Sect. 4.2.1 is not appropriate. Instead, a
simple weighted average of the partition resolutions at
$\text{Q}_{\beta\beta}$ gives the resolution expectation value to be
associated with events in the region of interest, i.e.:

$\displaystyle\sigma=\sum_{i}w_{i}\sigma_{i},$ (13)

where the sum goes over the partitions with resolutions $\sigma_{i}$ and
weights $w_{i}$. For the three detector types we obtain the resolutions at
$\text{Q}_{\beta\beta}$ as given in Tab. 2.

Detector type | Resolution at $\text{Q}_{\beta\beta}$ [keV]
---|---
BEGe | 2.8 $\pm$ 0.3
coaxial | 4.0 $\pm$ 1.3
IC | 2.9 $\pm$ 0.1
Table 2: FWHM resolutions at $\text{Q}_{\beta\beta}$ for datasets of each
detector type, reported as exposure-weighted averages. The uncertainty is
given by the standard deviation among the detector partitions.

##  5 Energy resolution uncertainty at $Q_{\beta\beta}$

The statistical uncertainty on the energy resolution decreases with rising
statistics over time, and is on the order of only a few eV. As such, the
uncertainty on the energy resolution is dominated by systematic effects. We
consider various sources of systematic uncertainty, given here in decreasing
order of their contribution: (i) resolution shifts over time; (ii) energy
scale shifts over time; (iii) choice of the resolution fitting function. Due
to the nature of these uncertainties, their magnitude will not decrease over
time, but could change if the detector setup or analysis methods change.

In the following sections, we explain how individual contributions to the
systematic uncertainty were determined (Sect. 5.1 to Sect. 5.3), and how they
are combined together to give a total uncertainty per partition (Sect. 5.4).

###  5.1 Resolution stability

We consider a systematic uncertainty estimated from the fluctuations in the
resolution obtained for each calibration over time. For each partition, we
calculate the standard deviation of the resolution at FEP,
$\sigma_{\mathrm{FEP}}$, among individual calibration runs in that partition.
Assuming that in Eq. 8, any systematic fluctuation of the energy resolution is
caused by the two correlated parameters changing proportionally, the energy
resolution uncertainty $\delta$ divided by the energy resolution $\sigma$ is
independent of energy. This is supported by the high degree of correlation
between the fit parameters $a$ and $b$ of Eq. 8 of $-0.81$ for the fitted
partition resolution curves. With this specific model, we can translate the
uncertainty at the FEP energy to $Q_{\beta\beta}$:

$\displaystyle\delta_{Q_{\beta\beta}}=\frac{\sigma_{Q_{\beta\beta}}}{\sigma_{\mathrm{FEP}}}\delta_{\mathrm{FEP}}.$
(14)

The mean value for this component across all partitions is 0.11 keV, with a
standard deviation of 0.06 keV.

###  5.2 Pulser stability

Once the energy scale has been determined via a calibration as described in
Sect. 3.1, the calibration curves are used until the next calibration. While
several parameters are monitored to ensure detector stability, fluctuations of
the energy scale can still deteriorate the resolution for physics data
compared to calibration data. Fluctuations on time scales smaller than the
typical calibration duration (1.5 h) are also present in the calibration data.
The effect from these short-term fluctuations will thus be included in the
calculated effective resolution. Fluctuations on larger time scales, up to
around one week, can, within the restraints of our data quality requirements,
contribute additionally to the resolution in physics data compared to the
resolution obtained from calibration data.

This additional contribution was estimated using the position of test pulser
events (see Sect. 3.3). Shifts in the test pulser positions averaged over 1.5
h, normalised by their statistical uncertainty, were analyzed. Were the
variation in energies due only to statistical fluctuations, these normalized
residuals would be distributed normally with a mean of 0 and a standard
deviation of 1. The observed deviation from this standard normal distribution
can be quantified as an additional contribution to the resolution, which is
typically on the order of 0.2 keV (1$\sigma$) or 0.6 keV (FWHM). As an
example, for a detector partition with a resolution of $\mathrm{FWHM}=3.0$
keV, the additional systematic uncertainty is given by:

$\displaystyle\delta_{\mathrm{sys.}}$
$\displaystyle=\sqrt{\mathrm{FWHM}^{2}+\mathrm{(0.6\,keV})^{2}}-\mathrm{FWHM}$
$\displaystyle=0.06\,\mathrm{keV}.$ (15)

The mean value for this component is 0.08 keV, with a standard deviation of
0.07 keV among partitions.

###  5.3 Choice of the resolution function

We used the square root of a linear function to model the resolution as a
function of energy (Eq. 8). While this choice is physically well-motivated,
including both statistical variations in the number of charge carriers, and
effects due to the electronics, there are some common alternatives. For
example, one could add a quadratic term under the square root to model the
effects of incomplete charge collection or integration,

$\displaystyle\sigma=\sqrt{a+bE+cE^{2}}.$ (16)

To estimate the variation of the resolution at $Q_{\beta\beta}$ for the
different choices of functions, the values obtained for the two discussed
choices are compared. Using the square root of linear (Eq. 8) and quadratic
(Eq. 16) functions, an average difference of 0.05 keV is obtained, with a
standard deviation of 0.05 keV among partitions.

###  5.4 Total resolution uncertainty by partition

The total resolution uncertainty is obtained by summing individual
contributions in quadrature, thereby assuming no correlations. The resultant
FWHM resolution uncertainties at $\text{Q}_{\beta\beta}$ of the partitions
vary between 0.04 keV and 0.37 keV, with a mean (standard deviation) of 0.13
(0.07) keV.

##  6 Energy bias and uncertainty

Due to the different assumptions and approximations in the calibration
procedure, slight biases in the energy scale may remain. Such biases may, for
example, be caused by the integral non-linearity of the FADCs Abgrall:2020jto
. Small non-linearities in the energy scale are for example neglected due to
the use of a linear calibration function. Therefore a peak from a $\gamma$ ray
with well defined energy might be displaced towards higher or lower energies.
Correspondingly, for each individual event, while its reconstructed energy
will fluctuate according to the resolution, it may also be systematically
displaced.

To evaluate the energy bias per partition near $Q_{\beta\beta}$, we look at
the residual at the SEP defined in Sect. 3.3 in the combined calibration
spectrum, since the SEP is very close to $Q_{\beta\beta}$ with a difference of
64.5 keV. The statistics is sufficient to reach a precision of
$\mathcal{O}$(0.01 keV) for the SEP position. The average bias is found to be
$-0.07$ keV, with a standard deviation of 0.29 keV among the partitions. Since
the $0\nu\beta\beta$ decay search is extremely sensitive to the energy of the
events close to $Q_{\beta\beta}$, in the final Gerda analysis Agostini:2020xta
, we correct for the energy bias of the events that fall into the energy range
considered for the $0\nu\beta\beta$ decay search (1930 keV to 2190 keV), by
adding the amount of bias to the calibrated event energy. This approach is
justified by studying the residuals at the 42K peak (1525 keV) and the DEP
(1592.5 keV), which are two closely located peaks with the former appearing in
the physics data Abramov:2019hhh  and the latter in the calibration data. The
relation between them is consistent with that between $Q_{\beta\beta}$ and the
SEP.

For the uncertainty of the bias, we use the residual fluctuations of the SEP
over time near $Q_{\beta\beta}$. We additionally include a systematic
uncertainty of 0.02 keV accounting for the potential difference between the
bias at SEP and that at $Q_{\beta\beta}$. It was estimated by performing a
linear interpolation between the residuals at the DEP and the SEP which are on
the two sides of $Q_{\beta\beta}$. In total, the average bias uncertainty is
0.17 keV.

##  7 Comparison to physics data

The two strongest $\gamma$ lines in our physics data spectrum are those due to
40K (1460.8 keV) and 42K (1524.7 keV) decays Agostini:2020xta . The measured
resolution of these peaks allows for a cross-check to the conclusions drawn
solely from calibration data. For every partition, the background energy
spectrum around each of these lines is fitted using a Gaussian for the signal
and a linear function for the background. The background rate was constrained
to be non-negative across the fitting window. Partitions with potassium peaks
with low counting statistics, i.e. those whose best-fit is compatible with
zero counts, are excluded from further analysis.

Given their proximity in energy, the extracted resolution for each of the two
lines is expected to coincide within 0.05 keV. Indeed, no significant
difference between the resolutions of the two peaks was found.

We compared the resolution obtained in the potassium lines with the one
predicted from the resolution curves extracted from the combined energy
spectra (see Sect. 4), as shown in Fig. 7. The systematic uncertainty for the
calibration resolution is calculated in the same way as described in Sect. 5.
The measured resolutions and predicted values from calibration data show a
high degree of correlation, with a Pearson correlation coefficient of 0.92,
and with 66% compatibility within one $\sigma$. Similar results are obtained
for the 40K line.

Figure 7: Resolution of the 1524.7 keV 42K line as measured from physics data
and extracted from calibration data, for each detector partition. The red line
shows the case of perfect agreement.

##  8 Conclusions

A reliable and stable energy scale is crucial to the search for
$0\nu\beta\beta$ decay of 76Ge performed with the Gerda experiment. The event
energies are reconstructed using the ZAC filter to minimize the effects of
low-frequency noise. To preserve the excellent energy resolution of the
germanium detectors when combining data over a long period of time, they are
calibrated weekly using 228Th sources. By identifying $\gamma$ peaks in the
recorded spectrum the energy scale and energy resolution can be determined.

For each calibration, the stability of the energy scale and resolution is
monitored via the 2.6 MeV FEP from 208Tl decays. Between successive
calibration runs the energy scale is monitored via test pulser events injected
into the readout electronics of the HPGe detectors. Data with short-term
instabilities are discarded from further analysis.

To more accurately reflect the properties of a detector at a certain time, we
have introduced the partitioning of the detectors’ data into stable sub-
periods. The stability is based on the long-term changes of the energy
resolution at the FEP and the residual at the SEP.

For each partition, a combined calibration analysis is performed to calculate
the energy resolution used for the $0\nu\beta\beta$ decay analysis. For this
purpose, calibration data in a partition are combined into a single spectrum.
The resolution curve is obtained by fitting a resolution model function to the
obtained resolutions of individual peaks in the combined spectrum. Among the
partitions, the calculated resolutions at $\text{Q}_{\beta\beta}$ range from
2.3 keV to 8.8 keV, with an exposure-weighted mean (standard deviation) of
3.0(0.8) keV.

Alternatively, effective resolution curves per detector type are calculated by
modeling the signal by a single Gaussian with a width according to the
standard deviation of a Gaussian mixture of the individual detector partition
contributions. Over Phase II we obtained exposure-weighted average resolutions
at $\text{Q}_{\beta\beta}$ for the BEGe/coaxial/IC detectors of ($2.8\pm 0.3$)
keV, ($4.0\pm 1.3$) keV, and ($2.9\pm 0.1$) keV respectively.

Dedicated studies were performed to determine the resolution systematic
uncertainties for the $0\nu\beta\beta$ decay analysis. Various sources of
systematic uncertainty on the resolution were considered: the fluctuations of
the resolution and energy scale over time, and the choice of resolution
function. The average total systematic uncertainty across all partitions is
0.13 keV.

The energy bias for the events near $Q_{\beta\beta}$ is estimated and
corrected based on the residual of the SEP. Among the partitions, the average
bias is -0.07 keV with a standard deviation of 0.29 keV. The average
uncertainty of these biases is 0.17 keV.

The energy scale, partitioning, resolutions, and energy biases discussed in
this paper are essential to the final search for $0\nu\beta\beta$ decay with
Gerda described in Agostini:2020xta . The success of the Gerda program in
reaching the world’s most stringent $0\nu\beta\beta$ decay half-life
constraint given by $T_{1/2}^{0\nu}>1.8\cdot 10^{26}$ yr at 90% C.L, was
achieved in part due to the excellent energy resolution offered by germanium
detectors and the analysis described in this work. This is an important step
towards Legend in developing the next generation of $0\nu\beta\beta$ decay
76Ge experiments Abgrall:2017syy .

##  9 Acknowledgements

The Gerda experiment is supported financially by the German Federal Ministry
for Education and Research (BMBF), the German Research Foundation (DFG), the
Italian Istituto Nazionale di Fisica Nucleare (INFN), the Max Planck Society
(MPG), the Polish National Science Centre (NCN), the Foundation for Polish
Science (TEAM/2016-2/17), the Russian Foundation for Basic Research, and the
Swiss National Science Foundation (SNF). This project has received
funding/support from the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreements No 690575 and No
674896. This work was supported by the Science and Technology Facilities
Council (ST/T004169/1). J. Huang and C. Ransom thank the UZH for the Postdoc
and Candoc Forschungskredit fellowships respectively. The institutions
acknowledge also internal financial support.

The Gerda collaboration thanks the directors and the staff of the LNGS for
their continuous strong support of the Gerda experiment.

## References

  * (1) R.N. Mohapatra, A.Y. Smirnov, Ann. Rev. Nucl. Part. Sci. 56, 569 (2006).  DOI 10.1146/annurev.nucl.56.080805.140534
  * (2) B.J. Mount, M. Redshaw, E.G. Myers, Phys. Rev. C 81, 032501 (2010).  DOI 10.1103/PhysRevC.81.032501
  * (3) M. Agostini et al. (GERDA Collaboration), Eur. Phys. J. C 79, 978 (2019).  DOI 10.1140/epjc/s10052-019-7353-8
  * (4) A. Domula, M. Hult, Y. Kermaïdic, G. Marissens, B. Schwingenheuer, T. Wester, K. Zuber, Nucl. Instrum. Meth. A 891, 106 (2018).  DOI 10.1016/j.nima.2018.02.056
  * (5) M. Agostini et al. (GERDA Collaboration), Eur. Phys. J. C78, 388 (2018).  DOI 10.1140/epjc/s10052-018-5812-2
  * (6) M. Gunther, et al. (Heidelberg-Moscow Collaboration), Phys. Rev. D 55, 54 (1997).  DOI 10.1103/PhysRevD.55.54
  * (7) A. Morales, Nucl. Phys. B Proc. Suppl. 77, 335 (1999).  DOI 10.1016/S0920-5632(99)00440-5
  * (8) K.H. Ackermann et al. (GERDA Collaboration), Eur. Phys. J. C 73, 2330 (2013).  DOI 10.1140/epjc/s10052-013-2330-0
  * (9) M. Agostini et al. (GERDA Collaboration), Phys. Rev. Lett. 111, 122503 (2013).  DOI 10.1103/PhysRevLett.111.122503
  * (10) M. Agostini et al. (GERDA Collaboration), JHEP 03, 139 (2020).  DOI 10.1007/JHEP03(2020)139
  * (11) S.R. Elliott, P. Vogel, Ann. Rev. Nucl. Part. Sci. 52, 115 (2002).  DOI 10.1146/annurev.nucl.52.050102.090641
  * (12) W. Maneschg, in _Prospects in Neutrino Physics_ (2017)
  * (13) L. Baudis, G. Benato, P. Carconi, C.M. Cattadori, P. De Felice, K. Eberhardt, R. Eichler, A. Petrucci, M. Tarka, M. Walter, JINST 10, P12005 (2015).  DOI 10.1088/1748-0221/10/12/P12005
  * (14) M. Tarka, Studies of Neutron Flux Suppression from a $\gamma$-ray Source and The GERDA Calibration System.  Ph.D. thesis, Universität Zürich (2012).  URL https://doi.org/10.5167/uzh-74790
  * (15) L. Baudis, A.D. Ferella, F. Froborg, M. Tarka, Nucl. Instrum. Meth. A 729, 557 (2013)
  * (16) R. Brun, F. Rademakers, Nucl. Instrum. Meth. A 389, 81 (1997).  DOI 10.1016/S0168-9002(97)00048-X
  * (17) M. Agostini, et al., J. Phys. Conf. Ser. 375, 042027 (2012).  DOI 10.1088/1742-6596/375/1/042027
  * (18) M. Agostini, L. Pandola, P. Zavarise, O. Volynets, JINST 6, P08013 (2011).  DOI 10.1088/1748-0221/6/08/P08013
  * (19) E. Gatti, P.F. Manfredi, La Rivista del Nuovo Cimento (1978-1999) 9, 1 (1986)
  * (20) M. Agostini et al. (GERDA Collaboration), Eur. Phys. J. C75, 255 (2015).  DOI 10.1140/epjc/s10052-015-3409-6
  * (21) V. D’Andrea, Improvement of Performance and Background Studies in GERDA Phase II.  Ph.D. thesis, Gran Sasso Science Institute (GSSI) (2017).  URL http://hdl.handle.net/20.500.12571/9641
  * (22) A. Lazzaro, Signal processing and event classification for a background free neutrinoless double beta decay search with the GERDA experiment.  Ph.D. thesis, Technische Universität München (2019).  URL mediatum.ub.tum.de/node?id=1507626
  * (23) R.L. Coldwell, G.P. Lasche, Journal of Radioanalytical and Nuclear Chemistry 307, 2509 (2016)
  * (24) I.N.E..E. Laboratory.  Gamma-ray spectrum catalog of isotopes.  URL {http://www.radiochemistry.org/periodictable/gamma_spectra/index.html}
  * (25) M. Agostini et al. (GERDA Collaboration), Phys. Rev. Lett. 125, 252502 (2020).  DOI 10.1103/PhysRevLett.125.252502
  * (26) M. Agostini et al. (GERDA Collaboration), Science 365, 1445 (2019).  DOI 10.1126/science.aav8613
  * (27) E.L. Hull, J. Xing, D.L. Friesel, R.H. Pehl, N.W. Madden, T.W. Raudorf, L.S. Varnell, _Charge Collection Physics in Semiconductor Detectors_ (1996)
  * (28) C. Ransom, Energy calibration for the GERDA and LEGEND-200 experiments.  Ph.D. thesis, Universität Zürich (2021)
  * (29) C. Améndola, Algebraic Statistics of Gaussian Mixtures.  Ph.D. thesis, Technische Universität Berlin (2017).  URL https://depositonce.tu-berlin.de/handle/11303/7284
  * (30) N. Abgrall et al. (MAJORANA Collaboration), IEEE Trans. Nucl. Sci. 68, 359 (2021).  DOI 10.1109/TNS.2020.3043671
  * (31) N. Abgrall et al. (LEGEND Collaboration), AIP Conf. Proc. 1894, 020027 (2017).  DOI 10.1063/1.5007652