to the presence of _noise_ in the photon stream. From Bose-Einstein statistics
the noise in the photon stream of a blackbody source, described by the
variance in the number of photons emitted per unit time $\sigma_{\gamma}^{2}$
can be related to the mean photon production per unit time $\mu_{\gamma}$,
photon energy $E_{\gamma}=h\nu$ and thermal energy of the source $E_{T}=kT$ by
[18]

$\sigma_{\gamma}^{2}=\mu_{\gamma}\left[\frac{e^{h\nu/kT}}{e^{h\nu/kT}-1}\right],$

where the quantity in square brackets is the so-called _boson factor_. In the
regime where the photon energy is significantly larger than the thermal energy
$h\nu\gg kT$, which corresponds to photon wavelengths of
$0.3-30\,\mu\mathrm{m}$ and temperatures $T<500\mathrm{K}$, the boson factor
is near unity such that $\sigma_{\gamma}^{2}\approx\mu_{\gamma}$ and the
number of photons $k$ emitted from the source per unit time is accurately
modeled by the Poisson mass function

$P(k)=\frac{\mu_{\gamma}^{k}e^{-\mu_{\gamma}}}{k!},\quad k\in\mathbb{N}_{0}.$

Light sources with these characteristics are referred to as _Poissonian_ and
are assumed in the photon transfer method.

####  3.1.2 Photon transfer $\gamma\to\mathrm{DN}$ and conversion gain

Detection of photons is done via electro-optical image sensors where the
pixels comprising the sensor convert, that is, _transfer_ photons $(\gamma)$
to electrons $(e\mathchar 45\relax)$. After exposing the sensor to a light
source for some _integration time_ , the packets of electrons collected by
each pixel are passed through the sensor’s circuitry whereby each charge
packet is first converted to a voltage $(V)$ and then digitized to produce a
_digital number_ $(\mathrm{DN})$ representing the intensity of the source at
each pixel. As one can see, several conversions, e.g. $\gamma\to e\mathchar
45\relax$, $e\mathchar 45\relax\to V$, and $V\to\mathrm{DN}$, take place in
the process of image formation. Treating each conversion process as a mapping
from one unit to another gives rise to the concept of the system _transfer
function_ $\mathcal{T}:\gamma\to\mathrm{DN}$ which represents the aggregate
mapping of photons to digital numbers. It is important to note that in general
each pixel will have its own unique transfer function and must be
characterized individually. For sensors where each pixel’s transfer function
can be assumed identical we say the sensor is _uniform_. Unless otherwise
specified, whenever we discuss the transfer function and measurement thereof,
we will assume it is for an individual pixel.

In most cases it is desirable to characterize a pixel in terms of electrons;
thus removing the dependence of photon wavelength. Luckily, for Silicon based
electro-optical sensors and photons in the visible band
$0.4-0.7\,\mu\mathrm{m}$ the conversion of photons to electrons, described by
the quantum yield $\eta$ is unity; therefore rendering the photon-to-digital
number transfer function $\mathcal{T}:\gamma\to\mathrm{DN}$ equivalent to the
electron-to-digital number transfer function $\mathcal{T}:e\mathchar
45\relax\to\mathrm{DN}$. The rather lofty goal of the photon transfer method
is to measure the transfer function solely by observing the pixel’s digital
output and then use the estimated transfer function to convert the digital
numbers back to a physical quantity of electrons produced in the pixel over
the integration time. The ability to reverse engineer the transfer function in
this way ultimately allows one to characterize the pixel in terms of key
imaging performance metrics like read noise and dynamic range all while
treating the sensor system as a black box.

Perhaps the simplest transfer function one could have is that represented by
the linear equation $\mathcal{T}(e\mathchar 45\relax)=e\mathchar 45\relax/g$
where $g$ is a constant with units $(e\mathchar 45\relax/\mathrm{DN})$. Pixels
that admit such a transfer function are naturally called _linear_ and the
constant $g$ is referred to as the _conversion gain_. In the special case of
linearity, an expression for $g$ is simple to derive. Assuming Poisson photon
statistics and unity quantum yield $\eta=1\,(\gamma/e\mathchar 45\relax)$, the
number of electrons collected by a pixel over the integration time is modeled
by $N\sim\operatorname{Poisson}(\mu_{e\mathchar 45\relax})$, where
$\mu_{e\mathchar 45\relax}=\mu_{\gamma}\operatorname{QE}_{\lambda}$ is the
mean number of electrons collected per integration time, $\mu_{\gamma}$ is the
corresponding mean number of incident photons per integration time, i.e. the
illumination level, and $\operatorname{QE}_{\lambda}$ is the wavelength
dependent quantum efficiency describing the probability that an incident
photon with wavelength $\lambda$ produces a free electron. In this context,
the transfer function can be viewed as a transformation of the random variable
$N$ so that if we let $P=\mathcal{T}(N)$ denote the random variable
representing the _photon induced_ digital output then
$\mu_{\mathrm{p}}\,(\mathrm{DN})=\mathsf{E}P=\mu_{e\mathchar 45\relax}/g$ and
$\sigma_{\mathrm{p}}^{2}\,(\mathrm{DN}^{2})=\mathsf{Var}P=\mu_{e\mathchar
45\relax}/g^{2}$; giving rise to the fundamental photon transfer relationship

$g=\frac{\mu_{\mathrm{p}}}{\sigma_{\mathrm{p}}^{2}}.$ (3.1)

The remarkable aspect of the model that led up to the relationship (3.1) is
that it provides a natural and simple way to estimate $g$ based solely on
observing the random output of a pixel exposed to incident illumination. If we
denote $\mathbf{P}=\\{P_{i}\\}_{i=1}^{n}$ as a sample of $n$ digital
observations generated by a pixel under illumination and
$T(\mathbf{P})=(\bar{P},\hat{P})$, where
$\bar{P}=\frac{1}{n}\sum_{i=1}^{n}P_{i}$ is the sample mean and
$\hat{P}=\frac{1}{n-1}\sum_{i=1}^{n}(P_{i}-\bar{P})^{2}$ as the sample
variance, then $T$ is an unbiased estimator of the parameter
$\theta=(\mu_{\mathrm{p}},\sigma_{\mathrm{p}}^{2})$ so that $g$ can be
estimated by

$G=\frac{\bar{P}}{\hat{P}}.$ (3.2)

While simple, the estimator (3.2) has an inherent weakness due to the fact
that the relationship $(\ref{eq:fundamental_g_relationship})$ does not account
for the background signal and noise produced by real pixels in the absence of
photon interaction. In some cases, given a sensor of sufficiently high
quality, it can be possible to expose the pixels to high levels of
illumination so that this background noise is dominated by photon noise and
$g\approx\mu_{\mathrm{p}}/\sigma_{\mathrm{p}}^{2}$ approximately holds. Such
sensors are said to achieve a _shot noise limited_ response and the
measurement of $g$ is said to be performed in the shot noise limited region of
the pixel’s dynamic range.

It is of no surprise that many sensors cannot achieve a shot noise-limited
response such that the background signal and noise produced by the pixels
cannot be ignored. To extend our model to include such cases we first let
$D\sim F_{D}$ represent the digital output of the pixels in the absence of
illumination (dark), which is distributed according to some distribution
$F_{D}$ with mean $\mu_{\mathrm{d}}\,(\mathrm{DN})=\mathsf{E}D$ and variance
$\sigma_{\mathrm{d}}^{2}\,(\mathrm{DN}^{2})=\mathsf{Var}D$. Assuming $D$ is
independent of $P$, the population mean and variance of the pixels digital
output under illumination becomes

$\displaystyle\mathsf{E}(P+D)\coloneqq\mu_{\mathrm{p}+\mathrm{d}}$
$\displaystyle=\mu_{\mathrm{p}}+\mu_{\mathrm{d}}=\mu_{e\mathchar
45\relax}/g+\mu_{\mathrm{d}}$
$\displaystyle\mathsf{Var}(P+D)\coloneqq\sigma_{\mathrm{p}+\mathrm{d}}^{2}$
$\displaystyle=\sigma_{\mathrm{p}}^{2}+\sigma_{\mathrm{d}}^{2}=\mu_{e\mathchar
45\relax}/g^{2}+\sigma_{\mathrm{d}}^{2},$

respectively, which leads us to the modify the gain relationship
$(\ref{eq:fundamental_g_relationship})$ as

$g=\frac{\mu_{\mathrm{p}+\mathrm{d}}-\mu_{\mathrm{d}}}{\sigma_{\mathrm{p}+\mathrm{d}}^{2}-\sigma_{\mathrm{d}}^{2}}.$
(3.3)

Unlike the one-sample, shot noise-limited estimator for $g$ given by (3.2),
one can see from (3.3) that the measurement of $g$ in a sub shot noise-limited
regime will require two separate samples captured under dark and illuminated
conditions. Denoting the sample under illumination by
$\mathbf{X}=\\{X_{i}\\}_{i=1}^{n_{1}}$ and the sample in the dark by
$\mathbf{Y}=\\{Y_{i}\\}_{i=1}^{n_{2}}$, the unknown parameter vector
$\theta=(\mu_{\mathrm{p}+\mathrm{d}},\sigma_{\mathrm{p}+\mathrm{d}}^{2},\mu_{\mathrm{d}},\sigma_{\mathrm{d}}^{2})$
can be estimated with the unbiased two-sample statistic
$T(\mathbf{X},\mathbf{Y})=(\bar{X},\hat{X},\bar{Y},\hat{Y})$ with $\bar{X}$
and $\hat{X}$ representing the sample mean and sample variance of $\mathbf{X}$
and likewise for the sample $\mathbf{Y}$. Substituting the components of $T$
directly into $(\ref{eq:general_g_relationship})$ gives an estimator for $g$
of the form

$G=\frac{\bar{X}-\bar{Y}}{\hat{X}-\hat{Y}}=\frac{\bar{P}}{\hat{P}},$ (3.4)

where $\bar{P}$ and $\hat{P}$ are the new estimators of the photon induced
mean $\mu_{\mathrm{p}}$ and variance $\sigma_{\mathrm{p}}^{2}$, respectively.

By accounting for the effects of dark signal and noise with
$(\bar{Y},\hat{Y})$, the estimator (3.4) does relieve the need for a pixel to
achieve a shot noise limited response. However, as was discussed in the
introduction of Chapter 1, the introduction of these dark corrections, in
particular that of $\hat{Y}$, lead to (3.4) exhibiting ill-behaved
characteristics in low illumination conditions. To combat this one may measure
$g$ with $(\ref{eq:sub_shot_noise_g_estimator})$ under high illumination
conditions where the estimator will be most well-behaved. That said, such a
procedure becomes invalid if the pixel admits a nonlinear transfer function.

Characterizing pixels with nonlinear transfer functions using only sample
statistics of the output signal is a much more challenging task. The inherent
complexity of characterizing nonlinear pixels comes from the fact that we can
no longer impose the very restrictive assumption that $g$ is a constant and
instead must treat it as some unknown function $g(\cdot)$ that varies with
illumination level. In years past, two methods for nonlinear characterization
have been proposed, namely, the nonlinear compensation (nlc) and nonlinear
estimation (nle) techniques [16, 21, 4]. In the nlc technique, the approach
taken is to not directly measure $g$ but rather measure analogous gains for
the first two moments of the photon induced signal $P$. What makes the nlc
method work is the observation that pixels with nonlinear transfer functions
typically exhibit a linear response at low illumination [15, 16]. As such, if
one assumes Poisson photon statistics and exposes the pixel to a sufficiently
low level of illumination $\mu_{\gamma}^{\ast}$, the mean number of electrons
collected in the pixel can be described by

$\mu_{e\mathchar
45\relax}^{\ast}=(\mu_{\mathrm{p}+\mathrm{d}}^{\ast}-\mu_{\mathrm{d}})\times
g^{\ast},$

where $\mu_{\mathrm{p}+\mathrm{d}}^{\ast}$ and $g^{\ast}$ represent the the
corresponding quantities at the illumination level $\mu_{\gamma}^{\ast}$. If
we then define the relative illumination level
$r=\mu_{\gamma}/\mu_{\gamma}^{\ast}$ one can subsequently define the
illuminated population parameters and mean electron signal as functions of $r$
whereby $\mu_{e\mathchar 45\relax}^{\ast}=\mu_{e\mathchar 45\relax}(1)$,
$\mu_{\mathrm{p}+\mathrm{d}}^{\ast}=\mu_{\mathrm{p}+\mathrm{d}}(1)$,
$\sigma_{\mathrm{p}+\mathrm{d}}^{2\ast}=\sigma_{\mathrm{p}+\mathrm{d}}^{2}(1)$.
Using this notation the expected mean number of electrons at relative
illumination level $r$ becomes

$\mu_{e\mathchar 45\relax}(r)=\mu_{e\mathchar 45\relax}^{\ast}\times r.$

Having an explicit expression for this quantity then allows us to define a
_signal gain_ and _noise gain_ in terms of relative illumination level as

$\displaystyle s(r)$ $\displaystyle=\frac{\mu_{e\mathchar
45\relax}(r)}{\mu_{\mathrm{p}+\mathrm{d}}(r)-\mu_{\mathrm{d}}}$ (3.5)
$\displaystyle n(r)$ $\displaystyle=\sqrt{\frac{\mu_{e\mathchar
45\relax}(r)}{\sigma_{\mathrm{p}+\mathrm{d}}^{2}(r)-\sigma_{\mathrm{d}}^{2}}},$

where both functions have units of $(e\mathchar 45\relax/\mathrm{DN})$.

In practice, nlc photon transfer characterization is performed by first
identifying the linear region of the pixel and measuring $g^{\ast}$. Then,
$s(r)$ and $n(r)$ are measured at several illumination levels. At each
illumination level, $r$ is recorded and the transfer functions are estimated
by according to the formulas in $(\ref{eq:signal_noise_transfer_functions})$
above. The estimated transfer functions are then fit by some curve to produce
signal and noise gains, which allow one to convert sample means and standard
deviations in units of digital numbers back into meaningful quantities of
electrons.

###  3.2 Noise model

Dark noise $D\sim F_{D}$ represents the sum of many different noise sources
present in the pixel and downstream circuitry and therefore is justifiably
modeled as normal. We further assume that individual observations $D_{i}$ of
dark noise are mutually independent such that the sample
$\mathbf{D}=(D_{1},\dots,D_{n})^{\mathsf{T}}$ is modeled by
$\mathbf{D}\sim\mathcal{N}(\mu_{\mathrm{d}}\mathbf{1}_{n},\sigma_{\mathrm{d}}^{2}\mathbf{I}_{n})$
where $\mathbf{1}_{n}\in\mathbb{R}^{n\times 1}$ denotes a column vector of
ones and $\mathbf{I}_{n}\in\mathbb{R}^{n\times n}$ is the $n\times n$ identity
matrix. Furthermore, photon induced noise $P$, although Poissonian by
assumption, quickly approaches a normal approximation for even small values of
$\mu_{\gamma}$111A typical rule of thumb is the normal approximation is useful
for $\mu_{\gamma}>30$.. Since the arrival of photons is independent of the
dark noise $D$ and individual observations $P_{i}$ are mutually independent,
the observed vector
$\mathbf{P}+\mathbf{D}=((P+D)_{1},\dots,(P+D)_{n})^{\mathsf{T}}$ is modeled as
$(\mathbf{P}+\mathbf{D})\sim\mathcal{N}(\mu_{\mathrm{p}+\mathrm{d}}\mathbf{1}_{n},\sigma_{\mathrm{p}+\mathrm{d}}^{2}\mathbf{I}_{n})$.
With these underlying assumptions we will let
$\mathbf{X}=(X_{1},\dots,X_{n_{1}})^{\mathsf{T}}$ be a sequences of $n_{1}$
digital observations of a pixel under some level of illumination and
$\mathbf{Y}=(Y_{1},\dots,Y_{n_{2}})^{\mathsf{T}}$ be a separate sample of
$n_{2}$ observations of the same pixel in the dark so that the joint vector is
modeled by

$\left(\mathbf{X}\atop\mathbf{Y}\right)\sim\mathcal{N}\left(\begin{pmatrix}\mu_{\mathrm{p}+\mathrm{d}}\mathbf{1}_{n_{1}}\\\
\mu_{\mathrm{d}}\mathbf{1}_{n_{2}}\end{pmatrix},\begin{pmatrix}\sigma_{\mathrm{p}+\mathrm{d}}^{2}\mathbf{I}_{n_{1}}&\mathbf{0}\\\
\mathbf{0}&\sigma_{\mathrm{d}}^{2}\mathbf{I}_{n_{2}}\end{pmatrix}\right).$
(3.6)

Under the proposed model, the two-sample statistic
$T(\mathbf{X},\mathbf{Y})=(\bar{X},\hat{X},\bar{Y},\hat{Y})^{\mathsf{T}}$,
where

$\displaystyle\bar{X}$ $\displaystyle=\frac{1}{n_{1}}\sum_{k=1}^{n_{1}}X_{k},$
$\displaystyle\hat{X}$
$\displaystyle=\frac{1}{n_{1}-1}\sum_{k=1}^{n_{1}}(X_{k}-\bar{X})^{2},$
$\displaystyle\bar{Y}$ $\displaystyle=\frac{1}{n_{2}}\sum_{k=1}^{n_{2}}Y_{k},$
$\displaystyle\hat{Y}$
$\displaystyle=\frac{1}{n_{2}-1}\sum_{k=1}^{n_{2}}(Y_{k}-\bar{Y})^{2},$

forms a vector of mutually independent components and constitutes a complete-
sufficient statistic for the unknown parameter vector
$\theta=(\mu_{\mathrm{p}+\mathrm{d}},\sigma_{\mathrm{p}+\mathrm{d}}^{2},\mu_{\mathrm{d}},\sigma_{\mathrm{d}}^{2})$.
Normal sampling theory gives the distributional results

$\displaystyle\bar{X}$
$\displaystyle\sim\mathcal{N}(\mu_{\mathrm{p}+\mathrm{d}},\sigma_{\mathrm{p}+\mathrm{d}}^{2}/n_{1}),$
$\displaystyle\hat{X}$
$\displaystyle\sim\mathcal{G}(\alpha_{1},\alpha_{1}/\sigma_{\mathrm{p}+\mathrm{d}}^{2}),$
$\displaystyle\bar{Y}$
$\displaystyle\sim\mathcal{N}(\mu_{\mathrm{d}},\sigma_{\mathrm{d}}^{2}/n_{2}),$
$\displaystyle\hat{Y}$
$\displaystyle\sim\mathcal{G}(\alpha_{2},\alpha_{2}/\sigma_{\mathrm{d}}^{2}),$

where $\alpha_{i}=(n_{i}-1)/2$. Since the shapes parameters $\alpha_{i}$ are
directly related to the sample sizes $n_{i}$ in this manner, we will use the
term _sample size_ to refer to both quantities when it is expedient to do so.
Further denoting the estimator for the photon induced mean as
$\bar{P}=\bar{X}-\bar{Y}$ also gives the sometimes useful result

$\bar{P}\sim\mathcal{N}(\mu_{\mathrm{p}},\sigma_{\mathrm{p}+\mathrm{d}}^{2}/n_{1}+\sigma_{\mathrm{d}}^{2}/n_{2}).$

We will further assume that the observed data $(\mathbf{X}^{\mathsf{T}}\
\mathbf{Y}^{\mathsf{T}})^{\mathsf{T}}$ is produced by a pixel with a linear
transfer function or equivalently is produced in the linear region of a pixel
with nonlinear transfer function. This assumption implies
$\mu_{\mathrm{p}+\mathrm{d}}=\mu_{\mathrm{d}}+\mu_{e\mathchar 45\relax}/g$ and
$\sigma_{\mathrm{p}+\mathrm{d}}^{2}=\sigma_{\mathrm{d}}^{2}+\mu_{e\mathchar
45\relax}/g^{2}$ where $g$ is the conversion gain given by
$(\ref{eq:general_g_relationship})$. Since $\mu_{e\mathchar 45\relax}\geq 0$
it follows that $\mu_{\mathrm{p}+\mathrm{d}}\geq\mu_{\mathrm{d}}>0$ and
$\sigma_{\mathrm{p}+\mathrm{d}}^{2}\geq\sigma_{\mathrm{d}}^{2}>0$ with
equality reached at zero illumination. Throughout the remaining sections, we
will also use the notation of Chapter 2 to denote the ratio of dark and
illuminated variances by
$\zeta=\sigma_{\mathrm{d}}^{2}/\sigma_{\mathrm{p}+\mathrm{d}}^{2}$ where
$\zeta\in[0,1]$, $\zeta=1$ is achieved at zero illumination, and $\zeta=0$ is
achieved at infinite illumination which is the shot noise limit222$\zeta=0$
can be treated as the mathematical definition of the shot noise limit. In
practice, the term _shot noise limited_ is a subjective term meant to describe
a pixel that can achieve a $\zeta$-value near zero before saturating..

###  3.3 The estimator $\mathscr{G}_{\nu}$

We are finally ready to present our new estimator for the photon transfer
conversion gain. Our first result shows that under the proposed noise model,
if an unbiased estimator for $g$ exists, it must have infinite variance on at
least a portion of the $\zeta$-domain.

######  Theorem 14 (Corollary of Theorem 1).

Under the normal model of pixel noise in $(\ref{eq:pixel_noise_model})$, if an
unbiased estimator of $g$ exists it must have infinite variance for at least
$\zeta\in(1/2,1]$.

###### Proof.

Assuming the model $(\ref{eq:pixel_noise_model})$,
$T(\mathbf{X},\mathbf{Y})=(\bar{P},\hat{X},\hat{Y})^{\mathsf{T}}$ is a
complete-sufficient statistic for the parameter
$\theta=(\mu_{\mathrm{p}},\sigma_{\mathrm{p}+\mathrm{d}}^{2},\sigma_{\mathrm{d}}^{2})^{\mathsf{T}}$.
Since the components of $T$ are mutually independent and
$\mathsf{E}\bar{P}=\mu_{\mathrm{p}}$, it follows that if unbiased estimator of
$g$ exists it will be of the form
$\mathscr{G}(T)=\bar{P}\times\mathscr{T}(\hat{X},\hat{Y})$ with
$\mathscr{T}(\hat{X},\hat{Y})$ denoting an unbiased estimator of
$(\sigma_{\mathrm{p}+\mathrm{d}}^{2}-\sigma_{\mathrm{d}}^{2})^{-1}$. Again
making use of the independence of $\bar{P}$ and $\mathscr{T}$ we then write

$\mathsf{Var}\mathscr{G}=(\mathsf{E}\bar{P}^{2})\mathsf{Var}\mathscr{T}+(\mathsf{E}\mathscr{T})^{2}\mathsf{Var}\bar{P}.$

But according to Theorem 1, if the estimator $\mathscr{T}$ exists it must have
infinite variance on at least $\zeta\in(1/2,1]$. Since $\mathscr{G}$ is an
umvue for $g$ and $\mathsf{Var}\mathscr{T}<\mathsf{Var}\mathscr{G}$ the
desired result immediately follows. ∎

Given that there is no unbiased estimator for $g$ that can achieve finite
variance on the entire $\zeta$-domain we again turn to biased estimation. The
following theorem presents a biased estimator based on the estimator
$\mathscr{T}_{\nu}$ derived in Chapter 2.

######  Proposition 5.

Under the assumed model we have
$\sigma_{\mathrm{p}+\mathrm{d}}^{2}=\sigma_{\mathrm{d}}^{2}/\zeta$,
$\sigma_{\mathrm{p}}^{2}=\sigma_{\mathrm{d}}^{2}\zeta^{-1}(1-\zeta)$, and
$\mu_{\mathrm{p}}=\sigma_{\mathrm{d}}^{2}\zeta^{-1}(1-\zeta)g$.

###### Proof.

The first two relationships are derived from
$\sigma_{\mathrm{p}+\mathrm{d}}^{2}=\sigma_{\mathrm{p}}^{2}+\sigma_{\mathrm{d}}^{2}$
and $\zeta=\sigma_{\mathrm{d}}^{2}/\sigma_{\mathrm{p}+\mathrm{d}}^{2}$. To
derive the relationship for $\mu_{\mathrm{p}}$ we use
$\mu_{\mathrm{p}}=\mu_{e\mathchar 45\relax}/g$ and
$\sigma_{\mathrm{p}}^{2}=\mu_{e\mathchar 45\relax}/g^{2}$ to write
$\mu_{\mathrm{p}}=\sigma_{\mathrm{p}}^{2}g$. ∎

######  Theorem 15.

Let $T=(\bar{P},\hat{X},\hat{Y})^{\mathsf{T}}$ and

$\mathscr{G}_{\nu}(T)=\bar{P}\times\mathscr{T}_{\nu}(\hat{X},\hat{Y}),\quad\nu>0.$

Then, $\mathscr{G}_{\nu}$ is the unique uniformly minimum variance unbiased
estimator of

$\mathsf{E}\mathscr{G}_{\nu}=(1-\zeta^{\nu})\,g,$

with

$\mathsf{Var}\mathscr{G}_{\nu}=\frac{\sigma_{\mathrm{d}}^{2}}{\zeta}\left(\frac{1}{n_{1}}+\frac{\zeta}{n_{2}}+\frac{(\sigma_{\mathrm{d}}g)^{2}}{\zeta}(1-\zeta)^{2}\right)\mathsf{E}\mathscr{T}_{\nu}^{2}-(1-\zeta^{\nu})^{2}g^{2}$

and $\mathsf{E}\mathscr{T}_{\nu}^{2}$ is given by Theorem 9 for
$\kappa_{1}=\sigma_{\mathrm{d}}^{2}/\zeta$ and $\alpha_{i}=(n_{i}-1)/2$.

###### Proof.

First, note that $\sigma_{\mathrm{p}+\mathrm{d}}^{2}>\sigma_{\mathrm{d}}^{2}$
for any nonzero illumination level; thus, by Remark 2 we must require $\nu>0$
to guarantee $\mathsf{ARB}\mathscr{T}_{\nu}<1$. Next, for the proposed normal
model $T(\mathbf{X},\mathbf{Y})=(\bar{X},\bar{Y},\hat{X},\hat{Y})$ is a
complete-sufficient statistic for the parameter
$\theta=(\mu_{\mathrm{p}+\mathrm{d}},\mu_{\mathrm{d}},\sigma_{\mathrm{p}+\mathrm{d}}^{2},\sigma_{\mathrm{d}}^{2})$;
thus by the Lehmann-Scheffé theorem $\mathscr{G}_{\nu}$ is the unique umvue of
its expected value. To evaluate the expected value, we use the independence of
$\bar{P}$ and $\mathscr{T}_{\nu}(\hat{X},\hat{Y})$ to write
$\mathsf{E}\mathscr{G}_{\nu}=(\mathsf{E}\bar{P})(\mathsf{E}\mathscr{T}_{\nu})$.
We know $\mathsf{E}\bar{P}=\mu_{\mathrm{p}}$ and $\mathsf{E}\mathscr{T}_{\nu}$
is given in Theorem 4. Using the formula for $g$ in
$(\ref{eq:general_g_relationship})$ then gives the desired result for
$\mathsf{E}\mathscr{G}_{\nu}$. To obtain the variance we again use
independence of $\bar{P}$ and $\mathscr{T}_{\nu}$ to write

$\mathsf{Var}\mathscr{G}_{\nu}=(\mathsf{Var}\bar{P}+(\mathsf{E}\bar{P})^{2})(\mathsf{E}\mathscr{T}_{\nu}^{2})-(1-\zeta^{\nu})^{2}g^{2}.$

Recalling
$\bar{P}\sim\mathcal{N}(\mu_{\mathrm{p}},\sigma_{\mathrm{p}+\mathrm{d}}^{2}/n_{1}+\sigma_{\mathrm{d}}^{2}/n_{2})$
and the relations of Proposition 5 then leads to the desired result. ∎

In the following analysis we will derive many results involving the
coefficient of variation of $\mathscr{G}_{\nu}$. Given the restrictions
$\nu>0$ and $\zeta\in(0,1)$ we know that $\mathsf{E}\mathscr{G}_{\nu}>0$ so
that $\mathsf{CV}\mathscr{G}_{\nu}=\mathsf{ACV}\mathscr{G}_{\nu}$. As such, to
maintain some level of congruence between the following results and those of
Chapter 2 we will use these quantities interchangabley when its convenient to
do so. We now show that $\mathsf{ARB}\mathscr{G}_{\nu}$ and
$\mathsf{ACV}\mathscr{G}_{\nu}$ are strictly increasing functions of $\zeta$.

######  Proposition 6.

$\mathsf{ARB}\mathscr{G}_{\nu}$ is a strictly increasing function of $\zeta$.
Furthermore, if $\nu>1$ then $\mathsf{ACV}\mathscr{G}_{\nu}$ is also a
strictly increasing function of $\zeta$.

###### Proof.

Given the definition of the absolute relative bias we have

$\mathsf{ARB}\mathscr{G}_{\nu}=\left\lvert\frac{(1-\zeta^{\nu})g-g}{g}\right\rvert=\zeta^{\nu}.$

Since $\nu>0$ it immediately follows that
$\partial_{\zeta}\mathsf{ARB}\mathscr{G}_{\nu}>0$ on $\zeta\in(0,1)$. Next, we
use the independence of $\bar{P}$ and $\mathscr{T}_{\nu}$ to write

$\mathsf{ACV}^{2}\mathscr{G}_{\nu}=\mathsf{CV}^{2}\mathscr{T}_{\nu}+(\mathsf{CV}^{2}\mathscr{T}_{\nu})(\mathsf{CV}^{2}\bar{P})+\mathsf{CV}^{2}\bar{P}$
(3.7)

so that

$\partial_{\zeta}\mathsf{ACV}^{2}\mathscr{G}_{\nu}=(\mathsf{CV}^{2}\bar{P}+1)(\partial_{\zeta}\mathsf{CV}^{2}\mathscr{T}_{\nu})+(\mathsf{CV}^{2}\mathscr{T}_{\nu}+1)(\partial_{\zeta}\mathsf{CV}^{2}\bar{P}).$
(3.8)

Using the relations in Proposition 5 we have

$\mathsf{CV}^{2}\bar{P}=\frac{1}{(\sigma_{\mathrm{d}}g)^{2}}\frac{\zeta}{(1-\zeta)^{2}}\left(\frac{1}{n_{1}}+\frac{\zeta}{n_{2}}\right),$
(3.9)

which is easily shown to be decreasing in $\zeta$ on $\zeta\in(0,1)$333In
$(\ref{eq:CV2P_zeta_form})$ we see the quantity $\sigma_{\mathrm{d}}g$ which
is the sensor dark noise in units of electrons. This quantity represents the
sum of read noise and dark current noise present in the pixel and is a
function of integration time.. Likewise, for $\nu>1$ we know from Lemma 27
that $\partial_{\zeta}\mathsf{ACV}^{2}\mathscr{T}_{\nu}>0$ which implies all
terms in $(\ref{eq:CV2Gv_zeta_derivative_expanded})$ are positive. Noting that
$\mathsf{ACV}\mathscr{G}_{\nu}$ is positive completes the proof. ∎

In the following Lemma we establish a property of $\mathscr{G}_{\nu}$ that
agrees with a long standing observation about estimators of $g$ in the
literature, that is, the dispersion of $\mathscr{G}_{\nu}$ is dominated by the
dispersion of the estimator for
$(\sigma_{\mathrm{p}+\mathrm{d}}^{2}-\sigma_{\mathrm{d}}^{2})^{-1}$ as the
illumination increases. To accomplish this we will consider the quantity

$\mathscr{E}=\frac{\mathsf{CV}^{2}\mathscr{T}_{\nu}}{\mathsf{CV}^{2}\mathscr{G}_{\nu}}.$

This ratio is useful for studying the dominance of
$\mathsf{Var}\mathscr{G}_{\nu}$ by the dispersion of $\mathscr{T}_{\nu}$ for a
couple reasons. First, using the definition of the coefficient of variation
and independence of $\bar{P}$ and $\mathscr{T}_{\nu}(\hat{X},\hat{Y})$ we are
able to write

$\mathsf{CV}^{2}\mathscr{G}_{\nu}=\mathsf{CV}^{2}\mathscr{T}_{\nu}\left(1+(1+\mathsf{CV}^{-2}\mathscr{T}_{\nu})\mathsf{CV}^{2}\bar{P}\right),$

and thus

$\mathscr{E}=\left[1+(1+\mathsf{CV}^{-2}\mathscr{T}_{\nu})\mathsf{CV}^{2}\bar{P}\right]^{-1}.$
(3.10)

Since the quantity inside the brackets is bounded below by one it follows that
$0<\mathscr{E}<1$ and so $\mathscr{E}$ gives an intuitive measure of
dominance. Second, if we define the random variable
$\mathscr{G}_{\nu}^{\ast}=\mu_{\mathrm{p}}\mathscr{T}_{\nu}(\hat{X},\hat{Y})$
which estimates $g$ when $\mu_{\mathrm{p}}$ is known, i.e. has zero variance,
then $\mathscr{E}$ is equivalent to

$\mathscr{E}=\frac{\mathsf{Var}\mathscr{G}_{\nu}^{\ast}}{\mathsf{Var}\mathscr{G}_{\nu}}.$

So in this context we see that $\mathscr{E}\approx 1$ corresponds to $\bar{P}$
behaving like a fixed constant in comparison to $\mathscr{T}_{\nu}$; hence,
$\mathsf{Var}\mathscr{G}_{\nu}$ is completely dominated by the dispersion of
$\mathscr{T}_{\nu}$.

######  Lemma 28.

As the illumination level increases $\mathscr{E}\searrow 1$.

###### Proof.

We have already shown that $\mathscr{E}<1$; thus, if it approaches one then it
must do so from below. All that is left is to show that $\mathscr{E}=1$ in the
limit of infinite illumination. Our starting point is the expression for
$\mathscr{E}$ in $(\ref{eq:E_ratio_expression})$, namely,

$\mathscr{E}=\left(1+(1+\mathsf{CV}^{-2}\mathscr{T}_{\nu})\mathsf{CV}^{2}\bar{P}\right)^{-1}.$

As the illumination level increases without bound, $\zeta\to 0$ and from Lemma
25

$\mathsf{CV}^{2}\mathscr{T}_{\nu}\to\begin{cases}\frac{2}{n_{1}-5},&n_{1}>5\\\
\infty,&n_{1}\leq 5.\end{cases}$

Hence, $1+\mathsf{CV}^{-2}\mathscr{T}_{\nu}$ approaches a positive and finite
constant in the limit. Next we turn our attention to the
$\mathsf{CV}^{2}\bar{P}$ term. Using $(\ref{eq:CV2P_zeta_form})$ it is
immediately obvious that $\mathsf{CV}^{2}\bar{P}\to 0$ as $\zeta\to 0$. Thus,
as the illumination increases without bound
$(1+\mathsf{CV}^{-2}\mathscr{T}_{\nu})\mathsf{CV}^{2}\bar{P}\to 0$ which
completes the proof. ∎

######  Theorem 16.

For the estimator $\mathscr{G}_{\nu}$ in Theorem 15

$\mathsf{ARB}\mathscr{G}_{\nu}=\mathsf{ARB}\mathscr{T}_{\nu}$

and as illumination increases

$\mathsf{ACV}\mathscr{G}_{\nu}\sim\mathsf{ACV}\mathscr{T}_{\nu}.$

In particular, if $\alpha_{1}+\alpha_{2}>2$, $\alpha_{1}>2\nu$, and $\nu>1$
then as illumination increases, $\zeta\nearrow 0$ and

$\mathsf{ACV}\mathscr{G}_{\nu}\sim\mathsf{ACV}\mathscr{T}_{\nu}\left(1+\frac{n_{1}-3}{4n_{1}}\frac{\zeta}{(\sigma_{\mathrm{d}}g)^{2}}+\mathcal{O}(\zeta^{2})\right).$

###### Proof.

The first two claims follow directly from Proposition 6 and Lemma 28,
respectively. For the last claim we assume $\nu>1$ and use
$(\ref{eq:gTilde_asym_form})$ to deduce

$\tilde{g}^{2}_{k,\ell}(\zeta,\nu)\sim\begin{cases}(k!)^{2}+\mathcal{O}(\zeta),&k=\ell\\\
\mathcal{O}(\zeta^{2}),&k>\ell.\end{cases}$

If in addition $\alpha_{1}+\alpha_{2}>2$ and $\alpha_{1}>2\nu$ then
$\mathsf{CV}\mathscr{T}_{\nu}$ is finite and it follows that as $\zeta\nearrow
0$

$\mathsf{CV}^{2}\mathscr{T}_{\nu}\sim\sum_{k=1}^{\infty}\frac{k!}{(\alpha_{1})_{k}}+\mathcal{O}(\zeta)=\frac{1}{\alpha_{1}-2}+\mathcal{O}(\zeta)$

and

$1+\mathsf{CV}^{-2}\mathscr{T}_{\nu}\sim\alpha_{1}-1+\mathcal{O}(\zeta)=\frac{n_{1}-3}{2}+\mathcal{O}(\zeta).$

Multiplying this result by

$\mathsf{CV}^{2}\bar{P}\sim\frac{1}{(\sigma_{\mathrm{d}}g)^{2}}\frac{\zeta}{n_{1}}+\mathcal{O}(\zeta^{2})$

and using $\sqrt{1+f(x)}\sim 1+\frac{1}{2}f(x)+\mathcal{O}(f^{2}(x))$ as
$f(x)\to 0$ then gives

$\mathscr{E}^{-1/2}=\frac{\mathsf{ACV}\mathscr{G}_{\nu}}{\mathsf{ACV}\mathscr{T}_{\nu}}\sim
1+\frac{n_{1}-3}{4n_{1}}\frac{\zeta}{(\sigma_{\mathrm{d}}g)^{2}}+\mathcal{O}(\zeta^{2}),$

which is the desired result. ∎

The main conclusion of Theorem 16 is that confidence intervals for
$\mathsf{ARB}\mathscr{T}_{\nu}$ are equal to those of
$\mathsf{ARB}\mathscr{G}_{\nu}$ and that for sufficiently high illumination
confidence intervals for $\mathsf{ACV}\mathscr{T}_{\nu}$ may be used to
approximate those of $\mathsf{ACV}\mathscr{G}_{\nu}$. That said, these results
do not give any indication as to how high the illumination level must be to
achieve a good approximation or how $\mathsf{ACV}\mathscr{G}_{\nu}$ compares
to $\mathsf{ACV}\mathscr{T}_{\nu}$ under low illumination when taking into
account other parameters like the dark noise $\sigma_{\mathrm{d}}$, conversion
gain $g$, and sample sizes. We will return to this problem later and for now
will be content with knowing that given a sufficiently high illumination level
$\mathsf{ACV}\mathscr{T}_{\nu}\approx\mathsf{ACV}\mathscr{G}_{\nu}$.

###  3.4 A demonstration of $g$-estimation with confidence intervals

With Theorem 16 at hand, we now bring together several of the results derived
thus far and demonstrate the process of estimating $g$ in the context of the
photon transfer method. Using the parameter values in Table 3.1, $N=10^{6}$
pseudo-random observations of $T=(\bar{P},\hat{X},\hat{Y})$ were generated.
The parameter values chosen represent what one might typically see in the
photon transfer method. Since the sample sizes used are large we use
$(\ref{eq:Tv_asym_hyper_form})$ to derive the _first order_ approximation

$\mathscr{G}_{\nu,1}(T)=\bar{P}\times\mathscr{T}_{\nu,1}(\hat{X},\hat{Y}),$

where

$\mathscr{T}_{\nu,1}(\hat{X},\hat{Y})=\frac{1}{\hat{Y}}\Biggl{(}\left(1-\frac{1}{\alpha_{1}}\right)\nu{F}\left({1,1+\nu\atop
2};1-\frac{\hat{X}}{\hat{Y}}\right)\\\
-\frac{1}{\alpha_{1}}(\nu)^{(2)}{F}\left({2,1+\nu\atop
3};1-\frac{\hat{X}}{\hat{Y}}\right)-\frac{\alpha_{1}+\alpha_{2}}{6\alpha_{1}\alpha_{2}}(\nu)^{(3)}{F}\left({3,1+\nu\atop
4};1-\frac{\hat{X}}{\hat{Y}}\right)\Biggr{)}$

and $\alpha_{i}=(n_{i}-1)/2$. Note that the presence of integers in the top
and bottom parameters of the hypergeometric terms mean they reduce to
elementary functions of the argument $1-\hat{X}/\hat{Y}$.

Computation of the confidence intervals was done by first computing
$Z_{\alpha}=\hat{Y}/\hat{X}\,F_{\alpha,2\alpha_{1},2\alpha_{2}}$ for each pair
of $(\hat{X},\hat{Y})$ and confidence level $\alpha=0.05$. The $95\%$
confidence intervals for $\mathsf{ARB}\mathscr{G}_{\nu}$ were then computed by
substituting the values for $Z_{\alpha}$ into Theorem 12. As for the
confidence intervals of $\mathsf{ACV}\mathscr{G}_{\nu}$, a quick computation
of $\mathscr{E}$ using the parameters in Table 3.1 shows that despite the
relatively low level of illumination used in the simulation
$\mathscr{E}=0.998685\dots$ so that we may approximate these intervals with
those for $\mathscr{T}_{\nu}$ in Theorem 13. The relatively large values used
for $\alpha_{1}$ and $\alpha_{2}$ means computation of these interval
estimates will have to be done with the series expansion for
$\mathsf{ACV}\mathscr{T}_{\nu}$ which introduces the practical problem of
truncation error. If we define

$\mathsf{ACV}_{n}\mathscr{T}_{\nu}=\left(\sum_{k=1}^{n}\sum_{\ell=0}^{k}\frac{\tilde{g}_{k,\ell}^{2}(z,\nu)}{(\alpha_{1})_{\ell}(\alpha_{2})_{k-\ell}\,\ell!\,(k-\ell)!}\right)^{1/2},$

then the absolute relative error incurred in using this truncated expansion is
given by

$R_{n,1}=\left\lvert\frac{\mathsf{ACV}_{n}\mathscr{T}_{\nu}}{\mathsf{ACV}\mathscr{T}_{\nu}}-1\right\rvert.$

Corollary 17 (see Appendix D) provides an upper bound, $R_{n,m,1}^{\star}$,
for $R_{n,1}$ so that one may compute the number of terms $n=K^{\star}$ that
guarantees $R_{n}<\epsilon$ according to

$K^{\star}=\min\\{K|R_{K,m,1}^{\star}\leq\epsilon\land K\in\mathbb{N}\\}.$

Note that the quantity $R_{n,m,1}^{\star}$ in Corollary 17 provides a single
upper bound for the entire set of $Z_{\alpha}$ values so that $K^{\star}$ is
the appropriate number of terms needed in computing all $10^{6}$ interval
estimates within the specified error tolerance. The approximate interval
estimates for $\mathsf{ACV}\mathscr{G}_{\nu}$ may then be computed by
substituting $\mathsf{ACV}_{K^{\star}}\mathscr{T}_{\nu}$ and $Z_{\alpha}$ into
Theorem 13.

Table 3.2 presents the results of the simulation. Using the simulated values
for $Z_{\alpha}$ and choosing $m=0$ and an absolute relative truncation error
tolerance $\epsilon=5\times 10^{-4}$ led to only needing the first
$K^{\star}=3$ terms in the series expansion for
$\mathsf{ACV}\mathscr{T}_{\nu}$. Substituting $K^{\star}$ into
$R_{n,0,1}^{\star}$ subsequently showed that using this number of terms
guaranteed the relative truncation error was less than $2.7\times 10^{-4}$ for
all $10^{6}$ values of $Z_{\alpha}$. Additionally, note the relative
difference between the exact and estimated values for
$\mathsf{E}\mathscr{G}_{\nu}$ and $\mathsf{Var}\mathscr{G}_{\nu}$ are

$R_{\mathsf{E}\mathscr{G}_{\nu}}=\left\lvert\frac{\widehat{\mathsf{E}}\mathscr{G}_{\nu,1}}{\mathsf{E}\mathscr{G}_{\nu}}-1\right\rvert\times
100\%=0.001\%$

and

$R_{\mathsf{Var}\mathscr{G}_{\nu}}=\left\lvert\frac{\widehat{\mathsf{Var}}\mathscr{G}_{\nu,1}}{\mathsf{Var}\mathscr{G}_{\nu}}-1\right\rvert\times
100\%=0.24\%,$

which indicates no significant error was incurred from using the asymptotic
approximation $\mathscr{G}_{\nu,1}$ in place of $\mathscr{G}_{\nu}$. Lastly,
we see that the estimated coverage probability for the $95\%$ confidence
intervals for $\mathsf{ARB}\mathscr{G}_{\nu}$ and
$\mathsf{ACV}\mathscr{G}_{\nu}$ agree with the target values, the latter of
which is due to the fact that $\mathscr{E}\approx 1$ for the chosen
parameters.

This exercise raises additional questions as no instruction is given on the
appropriate values for $\nu$, $n_{1}$ and $n_{2}$ to achieve desired values of
$\mathsf{ARB}\mathscr{G}_{\nu}$ and $\mathsf{ACV}\mathscr{G}_{\nu}$. In
particular, for real experiments, gathering large numbers of observations
takes time and so knowledge of _optimal_ sample sizes would be useful in
reducing the time required to gather data. We will take a closer look at this
problem in the following section.

Parameter | Equation | Value | Unit
---|---|---|---
$\alpha$ | $-$ | $0.05$ | $-$
$n_{1}$ | $-$ | $3001$ | $-$
$n_{2}$ | $-$ | $1501$ | $-$
$g$ | $-$ | $5$ | $e\mathchar 45\relax/\mathrm{DN}$
$\nu$ | $-$ | $\frac{1}{2}e^{\pi}$ | $-$
$\mu_{e\mathchar 45\relax}$ | $-$ | $150$ | $e\mathchar 45\relax$
$\mu_{\mathrm{d}}$ | $-$ | $10$ | $\mathrm{DN}$
$\sigma_{\mathrm{d}}^{2}$ | $-$ | $16$ | $\mathrm{DN}^{2}$
$\mu_{\mathrm{p}+\mathrm{d}}$ | $\mu_{\mathrm{d}}+\mu_{e\mathchar 45\relax}/g$ | $40$ | $\mathrm{DN}$
$\sigma_{\mathrm{p}+\mathrm{d}}^{2}$ | $\sigma_{\mathrm{d}}^{2}+\mu_{e\mathchar 45\relax}/g^{2}$ | $22$ | $\mathrm{DN}^{2}$
$\zeta$ | $\sigma_{\mathrm{d}}^{2}/\sigma_{\mathrm{p}+\mathrm{d}}^{2}$ | $\frac{8}{11}$ | $-$
$\mu_{\bar{P}}$ | $\mu_{\mathrm{p}+\mathrm{d}}-\mu_{\mathrm{d}}$ | $30$ | $\mathrm{DN}$
$\sigma_{\bar{P}}^{2}$ | $\sigma_{\mathrm{p}+\mathrm{d}}^{2}/n_{1}+\sigma_{\mathrm{d}}^{2}/n_{2}$ | $\frac{81038}{4504501}$ | $\mathrm{DN}^{2}$
$\alpha_{1}$ | $(n_{1}-1)/2$ | $1500$ | $-$
$\alpha_{2}$ | $(n_{2}-1)/2$ | $750$ | $-$
$\beta_{1}$ | $\alpha_{1}/\sigma_{\mathrm{p}+\mathrm{d}}^{2}$ | $\frac{750}{11}$ | $\mathrm{DN}^{-2}$
$\beta_{2}$ | $\alpha_{2}/\sigma_{\mathrm{d}}^{2}$ | $\frac{375}{8}$ | $\mathrm{DN}^{-2}$
Table 3.1: Simulation parameters. Quantity | Exact Value | Estimated Value | Unit
---|---|---|---
$K^{\star}$ | $3$ | $-$ | $-$
$R_{K^{\star},0,1}^{\star}$ | $0.02671$ | $-$ | $\%$
$\mathsf{E}\mathscr{G}_{\nu}$ | $4.87447$ | $4.87440$ | $e\mathchar 45\relax/\mathrm{DN}$
$\mathsf{Var}\mathscr{G}_{\nu}$ | $0.36671$ | $0.36759$ | $(e\mathchar 45\relax/\mathrm{DN})^{2}$
$\mathsf{P}(\operatorname{CI}_{\alpha}(\mathsf{ARB}\mathscr{T}_{\nu})\ni\mathsf{ARB}\mathscr{G}_{\nu})$ | $0.95$ | $0.95008$ | $-$
$\mathsf{P}(\operatorname{CI}_{\alpha}(\mathsf{ACV}\mathscr{T}_{\nu})\ni\mathsf{ACV}\mathscr{G}_{\nu})$ | $\sim 0.95$ | $0.94945$ | $-$
Table 3.2: Simulation results based on $N=10^{6}$ pseudo-random observations.

###  3.5 Optimal sample sizes

When discussing the notion of optimal sample sizes we must first begin with a
definition of what is optimal. Typically, experimenters wish to achieve the
smallest possible uncertainty with the fewest number of total observations.
This goal is manifested in the following definition.

######  Definition 11 (Optimal sample sizes).

Let $\mathbf{X}_{j}=(X_{i,j},\dots,X_{n_{j},j})$ with
$X_{i,j}\overset{\mathrm{iid}}{\sim}F_{j}(\theta_{j})$ denote a random sample
of size $n_{j}$ drawn from the distribution $F_{j}$ with parameters
$\theta_{j}$. Furthermore, let $T(\mathbf{X})$ denote an estimator based on
$M$ random samples $\mathbf{X}=(\mathbf{X}_{1},\dots,\mathbf{X}_{M})$ and
$\mathsf{ACV}T(n|\theta_{T})$ denote its absolute coefficient of variation as
a function of the sample sizes $n=(n_{1},\dots,n_{M})$ with parameters
$\theta_{T}=\cup_{j}\theta_{j}$. Then the optimal sample sizes are defined as
the vector-valued function
$n^{\text{opt}}(\mathrm{N},\theta_{T})=(n_{1}^{\text{opt}},\dots,n_{M}^{\text{opt}})(\mathrm{N},\theta_{T})$
that minimizes $\mathsf{ACV}T(n|\theta_{T})$ subject to
$\mathrm{N}=n_{1}+\cdots+n_{M}$ for some fixed $\mathrm{N}>0$.

In the case where there are two samples, so that $n=(n_{1},n_{2})$, this
definition of optimality requires solving the single-variable optimization
problem

$\displaystyle n_{2}^{\text{opt}}(\mathrm{N},\theta_{T})$
$\displaystyle=\operatorname*{arg\,inf}_{n_{2}\in[0,\mathrm{N}]}\mathsf{ACV}^{2}T(\mathrm{N}-n_{2},n_{2}|\theta_{T})$
(3.11a) $\displaystyle n_{1}^{\text{opt}}(\mathrm{N},\theta_{T})$
$\displaystyle=\mathrm{N}-n_{2}^{\text{opt}}(\mathrm{N},\theta_{T}),$ (3.11b)

for some total number of samples $\mathrm{N}$ and parameters $\theta_{T}$.
Note that while this formulation does give the optimal sample sizes according
to definition 11, it does not allow one to specify a predetermined value
(desired outcome) for $\mathsf{ACV}T$. Since it may be desirable to specify
the outcome of $\mathsf{ACV}T$ instead of $\mathrm{N}$ we can introduce the
constraint
$\mathsf{ACV}T(n_{1}^{\text{opt}},n_{2}^{\text{opt}}|\theta_{T})=\mathsf{acv}_{0}$
and reformulate (3.11) into the equivalent two-variable problem

$\displaystyle\inf_{n_{2}\in[0,\mathrm{N}]}$
$\displaystyle\mathsf{ACV}^{2}T(\mathrm{N}-n_{2},n_{2}|\theta_{T})\big{|}_{(\mathrm{N},n_{2})=(n_{1}^{\text{opt}}+n_{2}^{\text{opt}},n_{2}^{\text{opt}})}$
(3.12a)
$\displaystyle\mathsf{ACV}^{2}T(n_{1}^{\text{opt}},n_{2}^{\text{opt}}|\theta_{T})-\mathsf{acv}_{0}^{2}=0,$
(3.12b)

which can then be solved for $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ in
terms of $\mathsf{acv}_{0}$ and $\theta_{T}$.

Before proceeding, let’s consider what this system of equations says and how
we solve it. Three cases must be considered which we shall denote as:
_nondegenerate_ ($n_{1}^{\text{opt}},n_{2}^{\text{opt}}\neq 0$), _weakly
degenerate_ ($n_{1}^{\text{opt}}=0$ or $n_{2}^{\text{opt}}=0$ but not both),
and _degenerate_ ($n_{1}^{\text{opt}},n_{2}^{\text{opt}}=0$). For the
nondegenerate case we begin with equation (3.12a) by fixing the total number
of samples to $\mathrm{N}$ and then minimizing $\mathsf{ACV}^{2}T$ in $n_{2}$
subject to this constraint. This process gives the optimal sample size
$n_{2}^{\text{opt}}$ in terms of $\mathrm{N}$ and $\theta_{T}$444In fact, if
we were to stop here equation (3.12a) would be identical to (3.11a), which
upon substituting $\mathrm{N}=n_{1}^{\text{opt}}+n_{2}^{\text{opt}}$ yields an
equation representing the optimal relationship, a.k.a. the _optimality
relation_ , between the two sample sizes. Solving the optimality relation for
either $n_{1}^{\text{opt}}$ or $n_{2}^{\text{opt}}$ and then substituting into
(3.12b) we then obtain one of the optimal sample sizes parameterized in terms
of $\mathsf{acv}_{0}$ and $\theta_{T}$. This solution may then be used in
conjunction with the optimality relation to obtain the solution for the
remaining optimal sample size. In the case of weak degeneracy, the infimum of
$\mathsf{ACV}^{2}T(\mathrm{N}-n_{2},n_{2}|\theta_{T})$ occurs at the boundary
$n_{2}=0,\mathrm{N}$ and so our two-sample problem reduces to that of finding
a single optimal sample size via (3.12b). Lastly, the case of degeneracy
results when $\mathsf{ACV}T\to 0$ and so zero samples are needed to achieve
the desired outcome $\mathsf{ACV}T=\mathsf{acv}_{0}$. As will be seen in the
latter analysis, this situation can occur at boundary points of the parameter
space and thus be treated as a limit of the nondegenerate case.

####  3.5.1 Asymptotic properties at low illumination

One of the unanswered curiosities of Section 3.4 was how the quantity
$\mathscr{E}$ behaves at low illumination w.r.t. the parameters
$(\sigma_{\mathrm{d}},g,n_{1},n_{2})$. If we fix these parameters and consider
what only happens as the illumination decreases, i.e. as $\zeta\searrow 1$, we
see that $\mathsf{CV}\bar{P}\to\infty$ while $\mathsf{CV}\mathscr{T}_{\nu}$
approaches a constant and thus $\mathscr{E}\to 0$. From this observation one
might conclude that the only way to achieve $\mathscr{E}$-values near unity is
to perform measurements at high illumination; thus, extinguishing the
possibility of approximating confidence intervals and optimal sample sizes for
$\mathscr{G}_{\nu}$ with those of $\mathscr{T}_{\nu}$ under low illumination
conditions. However, if we instead consider $\mathscr{E}$ as a function of the
optimal samples sizes for $\mathscr{T}_{\nu}$ or $\mathscr{G}_{\nu}$, which
vary with the illumination level, the limiting behavior at low illumination
changes drastically. The following analysis takes a look at the behavior of
the optimal sample sizes for $\bar{P}$, $\mathscr{T}_{\nu}$, and
$\mathscr{G}_{\nu}$ under low illumination conditions. The main results are
found in Lemma 29, Theorem 17, and Corollary 14, which show for each estimator
that the optimal samples sizes are asymptotically equal and proportional to
$(1-\zeta)^{-2}$ as $\zeta\searrow 1$. We will then conclude this section by
using these asymptotic results in Corollary 15 to study how $\mathscr{E}$
behaves at low illumination when subject to optimal sample sizes.

######  Lemma 29.

Let $(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}})$ denote the optimal sample
sizes for $\bar{P}$ that also satisfy
$\mathsf{ACV}\bar{P}(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}}|\zeta,\sigma_{\mathrm{d}},g)=\mathsf{acv}_{0}$
for some fixed $\mathsf{acv}_{0}\in\mathbb{R}^{+}$. Then as illumination
decreases, $\zeta\searrow 1$, $n_{2}^{{\text{opt}}}/n_{1}^{{\text{opt}}}\to
1$, and $n_{i}^{{\text{opt}}}\sim C_{\bar{P}}(1-\zeta)^{-2}$ where
$C_{\bar{P}}=2((\sigma_{\mathrm{d}}g)^{2}\mathsf{acv}_{0}^{2})^{-1}$.

###### Proof.

For the sake of brevity, we begin by letting
$a=\frac{1}{(\sigma_{\mathrm{d}}g)^{2}}\frac{\zeta}{(1-\zeta)^{2}}$ and
writing

$\mathsf{ACV}^{2}\bar{P}(n_{1},n_{2}|\zeta,\sigma_{\mathrm{d}},g)=a\left(\frac{1}{n_{1}}+\frac{\zeta}{n_{2}}\right).$

If $\zeta\in(0,1)$, then $a\neq 0$ and
$\mathsf{ACV}^{2}\bar{P}(\mathrm{N}-n_{2},n_{2}|\theta_{\bar{P}})$ is smooth
and strictly convex in $n_{2}$ on $(0,\mathrm{N})$. Additionally, the limiting
value of $\mathsf{ACV}^{2}\bar{P}(\mathrm{N}-n_{2},n_{2}|\theta_{\bar{P}})$ at
the endpoints $n_{2}=0,\mathrm{N}$ is infinite and so we know our optimization
problem is nondegenerate. These observations lead us to conclude that
$(\ref{eq:optimal_SOEs})$ can be uniquely solved by

$\displaystyle\partial_{n_{2}}\mathsf{ACV}^{2}\bar{P}(\mathrm{N}-n_{2},n_{2}|\theta_{\bar{P}})\Big{|}_{(\mathrm{N},n_{2})=(n_{1}^{\text{opt}}+n_{2}^{\text{opt}},n_{2}^{\text{opt}})}$
$\displaystyle=0$
$\displaystyle\mathsf{ACV}^{2}\bar{P}(n_{1}^{\text{opt}},n_{2}^{\text{opt}}|\theta_{\bar{P}})-\mathsf{acv}_{0}^{2}$
$\displaystyle=0,$

which upon substituting appropriate values produces the system of equations

$\displaystyle(1)$ $\displaystyle
a\left(\frac{1}{(n_{1}^{\text{opt}})^{2}}-\frac{\zeta}{(n_{2}^{\text{opt}})^{2}}\right)$
$\displaystyle=0$ $\displaystyle(2)$ $\displaystyle
a\left(\frac{1}{n_{1}^{\text{opt}}}+\frac{\zeta}{n_{2}^{\text{opt}}}\right)-\mathsf{acv}_{0}^{2}$
$\displaystyle=0.$

Upon inspection, $(1)$ gives us the optimality relation
$n_{2}^{\text{opt}}/n_{1}^{\text{opt}}=\sqrt{\zeta}$ and so as $\zeta\searrow
1$ we have $n_{2}^{\text{opt}}/n_{1}^{\text{opt}}\to 1$ which is the first
claim. Substituting the optimality relation into $(2)$ and solving for
$n_{2}^{\text{opt}}$ further yields
$n_{2}^{\text{opt}}=a/\mathsf{acv}_{0}^{2}(\sqrt{\zeta}+\zeta)$. Making use of
this result and again calling on the optimality relation we have after
reintroducing $a$:

$\displaystyle n_{1}^{\text{opt}}$
$\displaystyle=\frac{1}{(\sigma_{\mathrm{d}}g)^{2}\mathsf{acv}_{0}^{2}}\frac{\zeta(1+\sqrt{\zeta})}{(1-\zeta)^{2}}$
$\displaystyle n_{2}^{\text{opt}}$
$\displaystyle=\frac{1}{(\sigma_{\mathrm{d}}g)^{2}\mathsf{acv}_{0}^{2}}\frac{\zeta(\sqrt{\zeta}+\zeta)}{(1-\zeta)^{2}}.$

Now as $\zeta\searrow 1$,
$\zeta(1+\sqrt{\zeta})\sim\zeta(\sqrt{\zeta}+\zeta)\sim
2+\mathcal{O}(1-\zeta)$; hence,

$n_{i}^{\text{opt}}\sim\frac{2}{(\sigma_{\mathrm{d}}g)^{2}\mathsf{acv}_{0}^{2}}(1-\zeta)^{-2}+\mathcal{O}((1-\zeta)^{-1}),\quad
i=1,2$

which proves the second claim. Recall that we assumed $\zeta\in(0,1)$ so as to
render our optimization problem nondegenerate. However, one can obtain the
optimal sample sample sizes at the endpoints $\zeta=1$ and $\zeta=0$ by
considering the appropriate limits, the latter of which gives us the
degenerate case $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})=(0,0)$. The proof is
now complete. ∎

In Lemma 29 we studied the optimal sample sizes for $\bar{P}$ that satisfied
Definition 11 and the constraint $\mathsf{ACV}\bar{P}=\mathsf{acv}_{0}$.
Unlike $\bar{P}$, the estimator $\mathscr{T}_{\nu}$ is biased so we will
investigate its optimal sample sizes w.r.t. a desired outcome for both
$\mathsf{ACV}\mathscr{T}_{\nu}$ and $\mathsf{ARB}\mathscr{T}_{\nu}$. We could
also choose $\mathsf{ARB}\mathscr{T}_{\nu}=\mathsf{arb}_{0}$ to be constant
but this may not always be desirable as will be seen in later sections. In
particular, notice that for any $\nu>0$,
$\mathsf{ARB}\mathscr{T}_{\nu}=\zeta^{\nu}\to 0$ in the shot noise limit
$\zeta\nearrow 0$. As such, forcing a bias on $\mathscr{T}_{\nu}$ in the shot
noise limit by requiring $\mathsf{ARB}\mathscr{T}_{\nu}=\mathsf{arb}_{0}$
could be viewed as counterproductive in certain contexts. These observations
steer us in the direction of specifying a constraint on
$\mathsf{ARB}\mathscr{T}_{\nu}$ that varies with illumination level and
vanishes in the shot noise limit. To help aid in our discussion we introduce
the notion of the _bias profile_.

######  Definition 12 (Bias profile).

Let $T\sim F_{T}(\theta_{T})$ be an estimator of $g(\theta_{T})\neq 0$ for
which $\mathsf{ARB}T=\rho(\theta_{T})$. Then we say $\rho(\theta_{T})$ is the
bias profile for $T$.

Given the simplicity of the expression for $\mathsf{ARB}\mathscr{T}_{\nu}$ we
can choose just about any bias profile we desire by setting
$\nu=\log_{\zeta}\rho$. With so many available options we will aim for
simplicity and proceed with the following.

######  Definition 13.

For $\zeta\in[0,1]$, $\mathsf{arb}_{0}\in(0,1)$, and $b\in\mathbb{R}_{0}^{+}$,
${\nu^{\smash{\dagger}}}\coloneqq\log_{\zeta}\mathsf{arb}_{0}+b$.

######  Proposition 7.

Let ${\nu^{\smash{\dagger}}}$ be as defined in Definition 13. Then, the bias
profile of $\mathscr{T}_{\nu^{\smash{\dagger}}}$ is
$\rho=\mathsf{arb}_{0}\zeta^{b}$.

This choice of bias profile is useful since it is simple and guarantees
$\mathsf{ARB}\mathscr{T}_{\nu^{\smash{\dagger}}}\leq\mathsf{arb}_{0}$
everywhere on $\zeta\in[0,1]$. Furthermore, the parameter $b$ provides the
flexibility to choose a constant bias profile $(b=0)$ or a profile that
vanishes to zero in the shot noise limit $(b>0)$. With this matter out of the
way we present Lemmas 30-32, which provide us with the necessary results for
studying the behavior of the optimal sample sizes of
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ under low-illumination conditions in
Theorem 17.

######  Definition 14 (Lower-incomplete gamma function).

For $\Re s>0$

$\gamma(s,z)\coloneqq\int_{0}^{z}t^{s-1}e^{-t}\,\mathrm{d}t.$

######  Lemma 30.

Let ${\nu^{\smash{\dagger}}}$ be as define in Definition 13. Then as
$\zeta\searrow 1$

$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})\sim\frac{\gamma(n+1,-\log\mathsf{arb}_{0})}{1-\mathsf{arb}_{0}}(1-\zeta)^{-n}+\mathcal{O}((1-\zeta)^{1-n}),$

where $\gamma(s,z)$ is the lower-incomplete gamma function of Definition 14.

###### Proof.

From Corollary 10 form $(\mathrm{iii})$ we have

$(1-\zeta)^{n}\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})=n!\frac{\zeta^{n-\omega}-\mathsf{arb}_{0}\zeta^{b}\sum_{k=0}^{n}\frac{(\log_{\zeta}\mathsf{arb}_{0}+b-n+\omega)_{k}}{k!}(1-\zeta)^{k}}{1-\mathsf{arb}_{0}\zeta^{b}}.$

As $\zeta\searrow 1$, one may use the generalized binomial theorem to deduce
$\zeta^{x}\sim 1+\mathcal{O}(1-\zeta)$. Furthermore, we also determine as
$\zeta\searrow 1$: $\log_{\zeta}\mathsf{arb}_{0}\to\infty$ so that
$(\log_{\zeta}\mathsf{arb}_{0}+b-n+\omega)_{k}\sim\log_{\zeta}^{k}\mathsf{arb}_{0}$
and
$[(1-\zeta)\log_{\zeta}\mathsf{arb}_{0}]^{k}\sim(-\log\mathsf{arb}_{0})^{k}(1+\mathcal{O}(1-\zeta))$.
After making the appropriate substitutions and collecting terms we find

$(1-\zeta)^{n}\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})\sim\Gamma(n+1)\frac{1-\mathsf{arb}_{0}\sum_{k=0}^{n}\frac{(-\log\mathsf{arb}_{0})^{k}}{k!}}{1-\mathsf{arb}_{0}}+\mathcal{O}(1-\zeta),$

which according to [27, Eq. $06.06.03.0009.01$] can be written in terms of the
lower incomplete gamma function. Solving for $\tilde{g}_{n,\omega}$ completes
the proof. ∎

######  Definition 15 (Polygamma function).

For $n\in\mathbb{N}_{0}$,

$\psi^{(n)}(z)\coloneqq\partial_{z}^{n+1}\log\Gamma(z).$

######  Lemma 31.

Let $(\alpha_{1},\alpha_{2})\in\mathbb{R}^{+}\times\mathbb{R}^{+}$ and

$f_{k}(\alpha_{1},\alpha_{2})=\sum_{\ell=0}^{k}\frac{1}{(\alpha_{1})_{\ell}(\alpha_{2})_{k-\ell}\ell!(k-\ell)!},\quad
k\in\mathbb{N}.$

Then for all $k$ and some constant $\mathrm{A}>0$, the optimal pair
$(\alpha_{1}^{\text{opt}},\alpha_{2}^{\text{opt}})$ that minimize $f_{k}$
subject to $\mathrm{A}=\alpha_{1}+\alpha_{2}$ is

$(\alpha_{1}^{\text{opt}},\alpha_{2}^{\text{opt}})=(\tfrac{1}{2}\mathrm{A},\tfrac{1}{2}\mathrm{A}).$

###### Proof.

Working with the properties of the Pochhammer symbol and Relation 12 we may
write $f_{k}$ in the equivalent form

$f_{k}(\alpha_{1},\alpha_{2})=\frac{1}{k!}\frac{(\alpha_{1}+\alpha_{2}-1+k)_{k}}{(\alpha_{1})_{k}(\alpha_{2})_{k}}.$

Substituting $\mathrm{A}=\alpha_{1}+\alpha_{2}$ and differentiating w.r.t.
$\alpha_{2}$ yields

$\partial_{\alpha_{2}}f_{k}(\mathrm{A}-\alpha_{2},\alpha_{2})=f_{k}(\mathrm{A}-\alpha_{2},\alpha_{2})\Psi_{k}^{(0)}(\mathrm{A}-\alpha_{2},\alpha_{2}),$

where

$\Psi_{k}^{(0)}(\alpha_{1},\alpha_{2})=\psi^{(0)}(\alpha_{1}+k)-\psi^{(0)}(\alpha_{1})-\psi^{(0)}(\alpha_{2}+k)+\psi^{(0)}(\alpha_{2})$

and $\psi^{(n)}(z)$ is the polygamma function of Definition 15. By repeated
application of the recurrence relation
$\psi^{(n)}(z+1)=\psi^{(n)}(z)+(-1)^{n}n!z^{-n-1}$ we deduce

$\Psi_{k}^{(0)}(\alpha_{2})=\sum_{m=0}^{k-1}\left(\frac{1}{\mathrm{A}-\alpha_{2}+m}-\frac{1}{\alpha_{2}+m}\right)\begin{cases}<0,&0<\alpha_{2}<\frac{1}{2}\mathrm{A}\\\
=0,&\alpha_{2}=\frac{1}{2}\mathrm{A}\\\
>0,&\frac{1}{2}\mathrm{A}<\alpha_{2}<\mathrm{A}.\end{cases}$

Further noting that $f_{k}(\mathrm{A}-\alpha_{2},\alpha_{2})$ is positive
proves it has a unique global minimum at $\alpha_{2}=\frac{1}{2}\mathrm{A}$,
which completes the proof. ∎

######  Definition 16 (Exponential integral function).

$\operatorname{Ei}(z)\coloneqq\int_{0}^{z}\frac{e^{t}-1}{t}\,\mathrm{d}t+\frac{1}{2}\left(\log
z-\log\left(\tfrac{1}{z}\right)\right)+\gamma,$

where $\gamma=0.577216\dots$ is the Euler-Mascheroni constant.

######  Lemma 32.

For $c,x\in\mathbb{R}$

$\int_{[0,x]^{2}}e^{-u+cuv-v}\,\mathrm{d}u\mathrm{d}v=\tfrac{1}{c}e^{-1/c}\left(\operatorname{Ei}\left(\tfrac{1}{c}\right)-2\operatorname{Ei}\left(\tfrac{1}{c}(1-cx)\right)+\operatorname{Ei}\left(\tfrac{1}{c}(1-cx)^{2}\right)\right),$

where $\operatorname{Ei}(z)$ is the exponential integral function of
Definition 16.

###### Proof.

Let $I$ denote the integral in question. Integrating w.r.t. $u$ and then
substituting $t=(cv-1)x$ yields

$I=\tfrac{1}{c}e^{-1/c}\int_{-x}^{x(cx-1)}e^{-t/(cx)}\frac{e^{t}-1}{t}\,\mathrm{d}t.$

With a bit of algebra we may further write

$I=\tfrac{1}{c}e^{-1/c}\left(\int_{-x}^{x(cx-1)}\frac{e^{(1-1/(cx))t}-1}{t}\,\mathrm{d}t-\int_{-x}^{x(cx-1)}\frac{e^{-1/(cx)t}-1}{t}\,\mathrm{d}t\right).$

Now substituting $s=(1-1/(cx))t$ and $s=-1/(cx)t$ into the first and second
integral, respectively, we have after some simplification

$I=\tfrac{1}{c}e^{-1/c}\left(\int_{0}^{1/c}-2\int_{0}^{(1-cx)/c}+\int_{0}^{(1-cx)^{2}/c}\right)\frac{e^{s}-1}{s}\,\mathrm{d}s.$

From Definition 16, if $z\in\mathbb{R}$ then
$\int_{0}^{z}\frac{1}{t}(e^{t}-1)\,\mathrm{d}t=\operatorname{Ei}(z)-\log|z|-\gamma$.
Substituting this result into $I$, all of the logarithmic and constant terms
cancel leaving us with the desired result. ∎

######  Theorem 17.

Let $(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}})$ denote the optimal sample
sizes for $\mathscr{T}_{\nu^{\smash{\dagger}}}$ that also satisfy
$\mathsf{ACV}\mathscr{T}_{\nu^{\smash{\dagger}}}(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}}|\zeta,{\nu^{\smash{\dagger}}})=\mathsf{acv}_{0}$
and
$\mathsf{ARB}\mathscr{T}_{\nu^{\smash{\dagger}}}=\mathsf{arb}_{0}\zeta^{b}$
for $\mathsf{acv}_{0}\in\mathbb{R}^{+}$, $\mathsf{arb}_{0}\in(0,1)$, and
$b\in\mathbb{R}_{0}^{+}$ fixed. Then as illumination decreases, $\zeta\searrow
1$, $n_{2}^{{\text{opt}}}/n_{1}^{{\text{opt}}}\to 1$, and
$n_{i}^{{\text{opt}}}\sim
C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}(1-\zeta)^{-2}$ where
$C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ is the solution to

$\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}(C_{\mathscr{T}_{\nu^{\smash{\dagger}}}},\mathsf{arb}_{0})-\mathsf{acv}_{0}^{2}=0$

with

$\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}=\frac{(C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}/4)e^{-C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}/4}}{(1-\mathsf{arb}_{0})^{2}}\biggl{(}\operatorname{Ei}\left(\tfrac{C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}{4}\right)-2\operatorname{Ei}\left(\log\mathsf{arb}_{0}+\tfrac{C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}{4}\right)\\\
+\operatorname{Ei}\left(\tfrac{(C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}/4+\log\mathsf{arb}_{0})^{2}}{C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}/4}\right)\biggr{)}-1$

and $\operatorname{Ei}(z)$ is the exponential integral function of Definition
16.

###### Proof.

First note that $\mathscr{T}_{\nu^{\smash{\dagger}}}$ satisfies the constraint
$\mathsf{ARB}\mathscr{T}_{\nu^{\smash{\dagger}}}=\mathsf{arb}_{0}\zeta^{b}$
independent of $n_{1}$ and $n_{2}$. Using the result of Lemma 30 and the
relationship $\alpha_{i}=(n_{i}-1)/2$ we deduce as $\zeta\searrow 1$

$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(n_{1},n_{2}|\zeta,\mathsf{arb}_{0},b)\sim\sum_{k=1}^{\infty}\frac{f_{k}\left(\tfrac{n_{1}-1}{2},\tfrac{n_{2}-1}{2}\right)}{(1-\zeta)^{2k}}(a_{k}^{2}(\mathsf{arb}_{0})+\mathcal{O}(1-\zeta)),$

where $f_{k}(\alpha_{1},\alpha_{2})$ is as defined in Lemma 31 and
$a_{k}(\mathsf{arb}_{0})=\gamma(k+1,-\log\mathsf{arb}_{0})/(1-\mathsf{arb}_{0})$.
By Lemma 31 we know that the optimal sample sizes
$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ that minimizes $f_{k}$ with
$\mathrm{N}=n_{1}+n_{2}$ fixed occurs when $n_{1}=n_{2}$; thus,
$n_{2}^{{\text{opt}}}/n_{1}^{{\text{opt}}}\to 1$ as $\zeta\searrow 1$ which is
our first claim. Letting
$n_{1}^{{\text{opt}}}=n_{2}^{{\text{opt}}}=n^{{\text{opt}}}$ we have

$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(n_{1}^{\text{opt}},n_{2}^{\text{opt}}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})\sim\sum_{k=1}^{\infty}\frac{f_{k}\left(\tfrac{n^{\text{opt}}-1}{2},\tfrac{n^{\text{opt}}-1}{2}\right)}{(1-\zeta)^{2k}}(a_{k}^{2}(\mathsf{arb}_{0})+\mathcal{O}(1-\zeta)),$

Since $(1-\zeta)^{-2k}$ becomes arbitrarily large in the neighborhood of
$\zeta=1$ the function $f_{k}$ must simultaneously become arbitrarily small as
to maintain a finite value of
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ in the limit. But
$f_{k}((n^{\text{opt}}-1)/2,(n^{\text{opt}}-1)/2)$ is a decreasing function of
$n^{\text{opt}}$ so it must be that $n^{\text{opt}}\to\infty$ as
$\zeta\searrow 1$. Making use of the asymptotic relation $(s)_{k}\sim
s^{k}(1+\mathcal{O}(s^{-1}))$ as $s\to\infty$ and binomial theorem we deduce
for $n^{\text{opt}}\to\infty$:

$f_{k}\left(\tfrac{n^{\text{opt}}-1}{2},\tfrac{n^{\text{opt}}-1}{2}\right)\sim\frac{1}{k!}\left(\frac{4}{n^{\text{opt}}}\right)^{k}+\mathcal{O}((n^{\text{opt}})^{-k-1}),$

which upon substitution into our expression for
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ gives

$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(n_{1}^{\text{opt}},n_{2}^{\text{opt}}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})\sim\sum_{k=1}^{\infty}\left(\frac{1}{k!}\left(\frac{4}{(1-\zeta)^{2}n^{\text{opt}}}\right)^{k}+\frac{\mathcal{O}((n^{\text{opt}})^{-k-1})}{(1-\zeta)^{2k}}\right)\\\
\cdots\times(a_{k}^{2}(\mathsf{arb}_{0})+\mathcal{O}(1-\zeta)).$

Given that we require
$\mathsf{ACV}\mathscr{T}_{\nu^{\smash{\dagger}}}=\mathsf{acv}_{0}<\infty$ in
the limit it follows that we must have $n^{\text{opt}}\sim
C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}(1-\zeta)^{-2}+\mathcal{O}((1-\zeta)^{-1})$
for some positive constant $C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ which
proves our second claim. Substituting this asymptotic form for
$n^{\text{opt}}$ into our expression for
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ we see that all error
terms are $\mathcal{O}(1-\zeta)$ so that upon passing to the limit we obtain

$\lim_{\begin{subarray}{c}\zeta\searrow 1\\\
\mathsf{ACV}\mathscr{T}_{\nu^{\smash{\dagger}}}=\mathsf{acv}_{0}\\\
\mathsf{ARB}\mathscr{T}_{\nu^{\smash{\dagger}}}=\mathsf{arb}_{0}\zeta^{b}\end{subarray}}\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(n_{1}^{\text{opt}},n_{2}^{\text{opt}}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=\sum_{k=1}^{\infty}\frac{\gamma^{2}(k+1,-\log\mathsf{arb}_{0})}{(1-\mathsf{arb}_{0})^{2}}\frac{\left(4/C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}\right)^{k}}{k!},$

where $C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ is the positive constant that
solves

$\sum_{k=1}^{\infty}\frac{\gamma^{2}(k+1,-\log\mathsf{arb}_{0})}{(1-\mathsf{arb}_{0})^{2}}\frac{\left(4/C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}\right)^{k}}{k!}=\mathsf{acv}_{0}^{2}.$
(3.13)

Letting $\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}$
denote the series on the l.h.s. of (3.13) we use the integral representation
of the lower incomplete gamma function in Definition 14 to write

$\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}=\lim_{n\to\infty}\int_{[0,-\log\mathsf{arb}_{0}]^{2}}\frac{e^{-(u+v)}}{(1-\mathsf{arb}_{0})^{2}}\sum_{k=1}^{n}\frac{\left(4uv/C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}\right)^{k}}{k!}\,\mathrm{d}u\mathrm{d}v.$

The integrand is the sum of nonnegative terms and is bounded above by its
limiting form as $n\to\infty$ which in turn is integrable on the domain of
integration; thus, by argument of dominated convergence we have

$\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}=\int_{[0,-\log\mathsf{arb}_{0}]^{2}}\frac{e^{-u+4uv/C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}-v}}{(1-\mathsf{arb}_{0})^{2}}\,\mathrm{d}u\mathrm{d}v-1,$

which is evaluated in terms of the exponential integral function with Lemma
32. The proof is now complete. ∎

Figure 3.1 plots $C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ on the unit square.
These values for $C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ were computed in
mathematica using a Newton-Raphson iteration with starting point

$C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}^{\ast}=\frac{2\gamma^{2}(2,-\log\mathsf{arb}_{0})}{\mathsf{acv}_{0}^{2}(1-\mathsf{arb}_{0})^{2}}\left(1+\left(1+2\mathsf{acv}_{0}^{2}(1-\mathsf{arb}_{0})^{2}\frac{\gamma^{2}(3,-\log\mathsf{arb}_{0})}{\gamma^{4}(2,-\log\mathsf{arb}_{0})}\right)^{1/2}\right),$

which is the solution to $(\ref{eq:CT_solution_series_form})$ using only the
$k=1,2$ terms.

Figure 3.1: Plot of
$C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}(\mathsf{arb}_{0},\mathsf{acv}_{0})$
on the unit square.

Since $\mathsf{ARB}\mathscr{G}_{\nu}=\mathsf{ARB}\mathscr{T}_{\nu}$ we may
combine the results of Lemma 29 and Theorem 17 to make the following claims
about the optimal sample sizes for $\mathscr{G}_{\nu^{\smash{\dagger}}}$ at
low illumination.

######  Corollary 14.

Let $(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}})$ denote the optimal sample
sizes for $\mathscr{G}_{\nu^{\smash{\dagger}}}$ that also satisfy
$\mathsf{ACV}\mathscr{G}_{\nu^{\smash{\dagger}}}(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}}|\zeta,{\nu^{\smash{\dagger}}},\sigma_{\mathrm{d}},g)=\mathsf{acv}_{0}$
and
$\mathsf{ARB}\mathscr{G}_{\nu^{\smash{\dagger}}}=\mathsf{arb}_{0}\zeta^{b}$
for $\mathsf{acv}_{0}\in\mathbb{R}^{+}$, $\mathsf{arb}_{0}\in(0,1)$, and
$b\in\mathbb{R}_{0}^{+}$ fixed. Then as illumination decreases, $\zeta\searrow
1$, $n_{2}^{{\text{opt}}}/n_{1}^{{\text{opt}}}\to 1$, and
$n_{i}^{{\text{opt}}}\sim
C_{\mathscr{G}_{\nu^{\smash{\dagger}}}}(1-\zeta)^{-2}$ where
$C_{\mathscr{G}_{\nu^{\smash{\dagger}}}}$ is the solution to

$\overline{\mathsf{ACV}^{2}\mathscr{G}_{\nu^{\smash{\dagger}}}}(C_{\mathscr{G}_{\nu^{\smash{\dagger}}}},\mathsf{arb}_{0},\sigma_{\mathrm{d}},g)-\mathsf{acv}_{0}^{2}=0$

with

$\overline{\mathsf{ACV}^{2}\mathscr{G}_{\nu^{\smash{\dagger}}}}=\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}(C_{\mathscr{G}_{\nu^{\smash{\dagger}}}},\mathsf{arb}_{0})+2\frac{\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}(C_{\mathscr{G}_{\nu^{\smash{\dagger}}}},\mathsf{arb}_{0})}{(\sigma_{\mathrm{d}}g)^{2}C_{\mathscr{G}_{\nu^{\smash{\dagger}}}}}+\frac{2}{(\sigma_{\mathrm{d}}g)^{2}C_{\mathscr{G}_{\nu^{\smash{\dagger}}}}}$

and $\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}$ given in
Theorem 17.

###### Proof.

Begin by writing

$\mathsf{ACV}^{2}\mathscr{G}_{\nu^{\smash{\dagger}}}=\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}+(\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}})(\mathsf{ACV}^{2}\bar{P})+\mathsf{ACV}^{2}\bar{P}.$

As was the case in Theorem 17, $\mathscr{G}_{\nu^{\smash{\dagger}}}$ satisfies
the constraint on $\mathsf{ARB}\mathscr{G}_{\nu}$ independent of $n_{1}$ and
$n_{2}$. Furthermore, from Lemma 29 and Theorem 17 we know in the limit
$\zeta\searrow 1$ that the optimal sample sizes for $\bar{P}$ and
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ are asymptotically equal and
proportional to $(1-\zeta)^{-2}$. Substituting
$n_{1}^{\text{opt}}=n_{2}^{\text{opt}}\sim
C_{\mathscr{G}_{\nu^{\smash{\dagger}}}}(1-\zeta)^{-2}$ and passing to the
limit gives the desired expression for
$\overline{\mathsf{ACV}^{2}\mathscr{G}_{\nu^{\smash{\dagger}}}}$. ∎

We are finally able to return to the questions surrounding the behavior of
$\mathscr{E}$ as a function of the optimal sample sizes at low illumination.

######  Corollary 15.

Let $T$ stand in place for $\mathscr{T}_{\nu^{\smash{\dagger}}}$ or
$\mathscr{G}_{\nu^{\smash{\dagger}}}$ and
$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ denote the optimal sample sizes for
$T$ such that $\mathsf{ACV}T=\mathsf{acv}_{0}$ and
$\mathsf{ARB}T=\mathsf{arb}_{0}\zeta^{b}$ with
$\mathsf{acv}_{0}\in\mathbb{R}^{+}$, $\mathsf{arb}_{0}\in(0,1)$, and
$b\in\mathbb{R}_{0}^{+}$ fixed. Furthermore, define $\mathscr{E}_{T}$ to be
the quantity $\mathscr{E}$ as a function of
$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ and
$\overline{\mathscr{E}_{T}}=\lim_{\zeta\searrow 1}\mathscr{E}_{T}$. Then,

$\overline{\mathscr{E}_{T}}=\left(1+(1+\overline{\mathsf{ACV}^{-2}\mathscr{T}_{\nu^{\smash{\dagger}}}}(C_{T},\mathsf{arb}_{0}))\frac{2}{(\sigma_{\mathrm{d}}g)^{2}C_{T}}\right)^{-1},$

where $\overline{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}$ is
given in Theorem 17. In particular, if
$T=\mathscr{T}_{\nu^{\smash{\dagger}}}$:

$\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}=\left(1+\frac{1+\mathsf{acv}_{0}^{-2}}{(\sigma_{\mathrm{d}}g)^{2}}\frac{\zeta}{(1-\zeta)^{2}}\left(\frac{1}{n_{1}^{\text{opt}}}+\frac{\zeta}{n_{2}^{\text{opt}}}\right)\right)^{-1}$

and

$\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}=\left(1+\frac{2(1+\mathsf{acv}_{0}^{-2})}{(\sigma_{\mathrm{d}}g)^{2}C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}\right)^{-1}.$

Corollary 15 shows us that when one substitutes the optimal sample sizes for
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ or $\mathscr{G}_{\nu^{\smash{\dagger}}}$
into $\mathscr{E}$, the limit $\lim_{\zeta\searrow
1}\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ is nonzero. Additionally,
for the special case where the optimal sample sizes for
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ are chosen, both
$\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ and
$\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}$ reduce to very
simple and easy to compute expressions. To get a sense for how
$\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}$ behaves w.r.t.
its parameters, Table 3.3 presents numerical values of this quantity for
select values of $\sigma_{\mathrm{d}}g$, $\mathsf{arb}_{0}$, and
$\mathsf{acv}_{0}$. One can see for the chosen parameter values that the dark
noise, denoted $\sigma_{\mathrm{d}}g$, plays a significant role in the
proximity of $\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}$ to
one. Indeed, for $\sigma_{\mathrm{d}}g$ greater than about $5\,e\mathchar
45\relax$ we observe
$\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}\approx 1$, which
implies that the optimal sample sizes for
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ and
$\mathscr{G}_{\nu^{\smash{\dagger}}}$ are nearly identical even at near zero
illumination.

Using the optimal sample sizes of $\mathscr{T}_{\nu^{\smash{\dagger}}}$ in
place of those for $\mathscr{G}_{\nu^{\smash{\dagger}}}$ is not just
convenient for producing a compact expression for $\mathscr{E}$. When
measuring $g$ in an actual experiment, one cannot know the optimal sample
sizes with certainty because they are based on unknown parameters. For
$\mathscr{T}_{\nu^{\smash{\dagger}}}$, the optimal sample sizes are a function
of $\zeta=\sigma_{\mathrm{d}}^{2}/\sigma_{\mathrm{p}+\mathrm{d}}^{2}$, which
is itself a function of two unknown parameters that must be estimated. This
problem is compounded when considering optimal sample sizes for
$\mathscr{G}_{\nu^{\smash{\dagger}}}$ which are a function of
$(\sigma_{\mathrm{d}}^{2},g,\zeta)$ and thus require estimating four
parameters, i.e.
$(\mu_{\mathrm{d}},\mu_{\mathrm{p}+\mathrm{d}},\sigma_{\mathrm{d}},\sigma_{\mathrm{p}+\mathrm{d}})$.
The two additional degrees of freedom render estimating the optimal sample
sizes for $\mathscr{G}_{\nu^{\smash{\dagger}}}$ in a real experiment
impractical due to the amount of uncertainty introduced in the process of
estimation. Instead, it is far better to measure $g$ in a region where
$\mathscr{E}\approx 1$ so that the optimal sample sizes for
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ can be substituted. As has been
demonstrated here, this situation is possible even when the illumination level
is near zero small provided a sufficiently large dark noise.

| $\sigma_{\mathrm{d}}g=1\,e\mathchar 45\relax$ | $\sigma_{\mathrm{d}}g=2\,e\mathchar 45\relax$
---|---|---
$\mathsf{arb}_{0}/\mathsf{acv}_{0}$ | $0.01$ | $0.02$ | $0.05$ | $0.01$ | $0.02$ | $0.05$
$0.01$ | $0.645$ | $0.645$ | $0.646$ | $0.879$ | $0.879$ | $0.879$
$0.02$ | $0.629$ | $0.629$ | $0.629$ | $0.871$ | $0.871$ | $0.872$
$0.05$ | $0.587$ | $0.587$ | $0.587$ | $0.850$ | $0.850$ | $0.850$
| $\sigma_{\mathrm{d}}g=5\,e\mathchar 45\relax$ | $\sigma_{\mathrm{d}}g=10\,e\mathchar 45\relax$
$\mathsf{arb}_{0}/\mathsf{acv}_{0}$ | $0.01$ | $0.02$ | $0.05$ | $0.01$ | $0.02$ | $0.05$
$0.01$ | $0.978$ | $0.978$ | $0.979$ | $0.995$ | $0.995$ | $0.995$
$0.02$ | $0.977$ | $0.977$ | $0.977$ | $0.994$ | $0.994$ | $0.994$
$0.05$ | $0.973$ | $0.973$ | $0.973$ | $0.993$ | $0.993$ | $0.993$
Table 3.3: $\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}$ for
select values of $\sigma_{\mathrm{d}}g$, $\mathsf{arb}_{0}$, and
$\mathsf{acv}_{0}$.

####  3.5.2 Computation of optimal sample sizes for
$\mathscr{T}_{\nu^{\smash{\dagger}}}$

Given that the optimal samples sizes for $\mathscr{T}_{\nu^{\smash{\dagger}}}$
are most desirable for experimentation, this section is dedicated to further
studying their characteristics as well as computing them. Our approach will be
organized as follows. We will begin by using the analysis in Section 2.4.1 to
derive a few results pertaining to zeros and special values of
$\tilde{g}_{n,\omega}(z,\nu)$ evaluated at
$(z,\nu)=(\zeta,{\nu^{\smash{\dagger}}})$. These results will then allow us to
write the system of equations needed to solve for the optimal sample sizes
$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ as well as provide exact solutions
at two points of interest. From here we will derive explicit approximations
for the optimal sample sizes and use our previous findings to determine the
usefulness of these approximations. Finally, the derived approximations will
be used as the starting point in numerical methods to compute
$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$, which will be demonstrated for some
sample values.

The following proposition defines some notation needed to discuss the zeros of
$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})$ in Lemma 33.

######  Proposition 8.

Let ${\nu^{\smash{\dagger}}}(\zeta)=\log_{\zeta}\mathsf{arb}_{0}+b$ be a
function of $\zeta$ as defined in Definition 13,
$\mathbb{N}_{b}\coloneqq\\{n\in\mathbb{N}_{0}:n\geq b\\}$,
$\zeta_{b,n}\coloneqq\sqrt[n-b]{\mathsf{arb}_{0}}\mathds{1}_{b<n}$, and
$Z_{b}\coloneqq\\{\zeta_{b,n}:n\in\mathbb{N}_{b}\\}$. Then,

  * $(1)$

${\nu^{\smash{\dagger}}}(\zeta)$ is strictly increasing from $[0,1]$ onto
$[b,\infty)$,

  * $(2)$

${\nu^{\smash{\dagger}}}[Z_{b}]=\mathbb{N}_{b}$.

######  Lemma 33 (Zeros of
$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})$).

$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})=0\iff(n,n-\omega)\in\mathbb{N}^{2}\land(\zeta=0\lor(\zeta\in
Z_{b}\setminus\\{0\\}\land n-\omega\geq{\nu^{\smash{\dagger}}})).$

###### Proof.

The proof follows from Lemma 17 which states

$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})=0\iff A\lor B,$

with $A=(n\in\mathbb{N}\land
n>\omega\land{\nu^{\smash{\dagger}}}\in\mathbb{R}_{0}^{+}\land z=0)$ and
$B=(n-\omega-{\nu^{\smash{\dagger}}}\in\mathbb{N}_{0}\land{\nu^{\smash{\dagger}}}\neq
0)$. Starting with $A$ notice that
${\nu^{\smash{\dagger}}}\in\mathbb{R}_{0}^{+}$ always holds so that we may
equivalently write $A=((n,n-\omega)\in\mathbb{N}^{2}\land\zeta=0)$. Next,
consider $B$ and observe that a necessary condition for
${\nu^{\smash{\dagger}}}=0$ is $\zeta=0$. Since the condition $\zeta=0$ is
already covered by $A$ we may write
$B=(n-\omega-{\nu^{\smash{\dagger}}}\in\mathbb{N}_{0}\land\zeta\neq 0)$.
Furthermore, $n-\omega-{\nu^{\smash{\dagger}}}\in\mathbb{N}_{0}$ requires
${\nu^{\smash{\dagger}}}$ to be an integer and
$n-\omega\geq{\nu^{\smash{\dagger}}}$. But ${\nu^{\smash{\dagger}}}$ can only
be an integer if $\zeta\in Z_{b}$ so that
$B=(n-\omega\geq{\nu^{\smash{\dagger}}}\land\zeta\in Z_{b}\setminus\\{0\\})$.
But the conditions specified by $B$ implicitly require $n-\omega\geq
1\implies(n,n-\omega)\in\mathbb{N}^{2}$ so that this requirement can be
factored out of $A$ and $B$ leaving us with the desired result. ∎

######  Corollary 16.

For all $(\zeta,{\nu^{\smash{\dagger}}})\in[0,1]\times[b,\infty)$

$n=\omega\implies\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})\neq 0.$

Furthermore,

$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})=0,\
\forall(n,n-\omega)\in\mathbb{N}^{2}\iff\zeta\in\\{0,\zeta_{b,1}\\}.$

###### Proof.

The first claim immediately follows from Lemma 33 since $n=\omega\implies
n-\omega=0\notin\mathbb{N}$. Furthermore, the conditions of Lemma 33 are
always satisfied if $(n,n-\omega)\in\mathbb{N}^{2}$ and
$\zeta=0\lor{\nu^{\smash{\dagger}}}=1$, where the latter holds if and only if
$\zeta\in\\{0,\zeta_{b,1}\\}$. ∎

Corollary 16 will turn out to be very important in determining the behavior of
optimal sample sizes of $\mathscr{T}_{\nu^{\smash{\dagger}}}$. In particular,
the results of this corollary indicate that $\zeta=0$ and $\zeta=\zeta_{b,1}$
are the only values of $\zeta$ for which
$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})=0$ for all $n>\omega$. We
shall see that these two points of interest corresponds to the only
$\zeta$-values that render the optimization problem for the optimal sample
sizes of $\mathscr{T}_{\nu^{\smash{\dagger}}}$ weakly degenerate. Before
proceeding, we need one more result, which gives us the expression for
$\tilde{g}_{n,n}(\zeta,{\nu^{\smash{\dagger}}})$ at $\zeta=0,\zeta_{b,1}$.

######  Lemma 34.

Let $L\in\\{0,\zeta_{b,1}\\}$. Then $\lim_{\zeta\nearrow
L}\tilde{g}_{n,n}(\zeta,{\nu^{\smash{\dagger}}})=n!$ for all
$n\in\mathbb{N}_{0}$.

###### Proof.

From Corollary 10 form $(\mathrm{iii})$ we have for all $n\in\mathbb{N}_{0}$:

$\tilde{g}_{n,n}(\zeta,{\nu^{\smash{\dagger}}})=n!\frac{1-\mathsf{arb}_{0}\zeta^{b}\sum_{k=0}^{n}\frac{({\nu^{\smash{\dagger}}})_{k}}{k!}(1-\zeta)^{k}}{(1-\zeta)^{n}(1-\mathsf{arb}_{0}\zeta^{b})}.$

Each case is now examined separately.

  1. $(1)$

$\zeta\nearrow 0$.

###### Proof.

If $b>0$ then $\zeta\nearrow 0\implies\zeta^{b}\nearrow 0$ and so
$\tilde{g}_{n,n}(\zeta,{\nu^{\smash{\dagger}}}_{b>0})\to n!$. Likewise, if
$b=0$ then $\zeta^{b}=1$ while
$({\nu^{\smash{\dagger}}})_{k}\to\mathds{1}_{k=0}$ and so again we find
$\tilde{g}_{n,n}(\zeta,{\nu^{\smash{\dagger}}}_{b=0})\to n!$ as expected. ∎

  2. $(2)$

$\zeta\nearrow\zeta_{b,1}$.

###### Proof.

If $b\geq 1$ then $\zeta_{b,1}=0$, which was already covered in part $(1)$.
Now suppose $0\leq b\leq 1$. As $\zeta\nearrow\zeta_{b,1}$:
${\nu^{\smash{\dagger}}}\to 1$ and $\mathsf{arb}_{0}\zeta^{b}\to\zeta_{b,1}$;
hence

$\lim_{\zeta\nearrow\zeta_{b,1}}\tilde{g}_{n,n}(\zeta,{\nu^{\smash{\dagger}}}_{b\in[0,1]})=n!\frac{1-\zeta_{b,1}\sum_{k=0}^{n}(1-\zeta_{b,1})^{k}}{(1-\zeta_{b,1})^{n+1}}=n!,$

which completes the proof. ∎

∎

Our next task is to derive the system of equations needed for computing the
optimal sample sizes. Since $\mathsf{ACV}\mathscr{T}_{\nu^{\smash{\dagger}}}$
is naturally parameterized in terms of the shape parameters $\alpha_{1}$ and
$\alpha_{2}$ and not the sample sizes $n_{1}$ and $n_{2}$ our discussion will
focus on finding the optimal shape parameters according to Definition 11,
knowing that we may convert these optimal shape parameters to their
corresponding sample sizes through the relation
$n_{i}^{\text{opt}}=2\alpha_{i}^{\text{opt}}+1$. As was done in Lemma 29 we
will accomplish this task by letting $\mathrm{A}=\alpha_{1}+\alpha_{2}$ be
constant and then showing that
$\mathsf{ACV}^{2}\mathscr{T}_{\nu}(\mathrm{A}-\alpha_{2},\alpha_{2})$ is:
$(1)$ twice continuously differentiable and strictly convex on
$\alpha_{2}\in(0,\mathrm{A})$ and $(2)$ is infinite at the endpoints
$\alpha_{2}=0,\mathrm{A}$ with the exception of two points of interest. These
findings will then be used to show that the minimization in
$(\ref{eq:optimal_SOEs})$ can be uniquely solved by equating a derivative with
zero. The following result presents and intermediate step in demonstrating
convexity of
$\mathsf{ACV}^{2}\mathscr{T}_{\nu}(\mathrm{A}-\alpha_{2},\alpha_{2})$.

######  Lemma 35.

Let $k\in\mathbb{N}$, $\mathrm{A}>0$, $\zeta\in[0,1]$, and
${\nu^{\smash{\dagger}}}$ be as defined in Definition 13. Then,

$f_{k}(\alpha_{2}|\mathrm{A},\zeta,{\nu^{\smash{\dagger}}}):=\sum_{\ell=0}^{k}\frac{\tilde{g}_{k,\ell}^{2}(\zeta,{\nu^{\smash{\dagger}}})}{(\mathrm{A}-\alpha_{2})_{\ell}(\alpha_{2})_{k-\ell}\ell!(k-\ell)!}$

is positive and strictly convex on $\alpha_{2}\in(0,\mathrm{A})$. Furthermore,

$\lim_{\alpha_{2}\nearrow
0}f_{k}(\alpha_{2}|\cdot)=\begin{cases}k!/(\mathrm{A})_{k},&\zeta\in\\{0,\zeta_{b,1}\\}\\\
\infty,&\text{otherwise}\end{cases}$

and

$\lim_{\alpha_{2}\searrow\mathrm{A}}f_{k}(\alpha_{2}|\cdot)=\infty.$

###### Proof.

The claim of positivity comes from the fact that the summand of $f_{k}$ is
nonnegative and from Corollary 16 which implies
$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})\neq 0$ for at least one
$\ell\in\\{0,\dots,k\\}$. For the claim of convexity we first write

$\partial_{\alpha_{2}}^{2}\frac{1}{(\mathrm{A}-\alpha_{2})_{\ell}(\alpha_{2})_{k-\ell}}=\frac{[\Psi^{(0)}_{k,\ell}(\mathrm{A}-\alpha_{2},\alpha_{2})]^{2}+\Psi^{(1)}_{k,\ell}(\mathrm{A}-\alpha_{2},\alpha_{2})}{(\mathrm{A}-\alpha_{2})_{\ell}(\alpha_{2})_{k-\ell}},$

where

$\displaystyle\Psi^{(0)}_{k,\ell}(\alpha_{1},\alpha_{2})$
$\displaystyle=\psi^{(0)}(\alpha_{1})-\psi^{(0)}(\alpha_{1}+\ell)-\psi^{(0)}(\alpha_{2})+\psi^{(0)}(\alpha_{2}+k-\ell)$
$\displaystyle\Psi^{(1)}_{k,\ell}(\alpha_{1},\alpha_{2})$
$\displaystyle=\psi^{(1)}(\alpha_{1})-\psi^{(1)}(\alpha_{1}+\ell)+\psi^{(1)}(\alpha_{2})-\psi^{(1)}(\alpha_{2}+k-\ell).$

Certainly, $[\Psi^{(0)}_{k,\ell}(\mathrm{A}-\alpha_{2},\alpha_{2})]^{2}\geq 0$
for all $k$, $\ell$, and $\alpha_{2}\in(0,\mathrm{A})$. Furthermore,

$\Psi^{(1)}_{k,\ell}(\mathrm{A}-\alpha_{2},\alpha_{2})=\sum_{j=0}^{\ell-1}\frac{1}{(\mathrm{A}-\alpha_{2}+j)^{2}}+\sum_{j=0}^{k-\ell-1}\frac{1}{(\alpha_{2}+j)^{2}}.$

Since $(k,\ell)\in(\mathbb{N},\mathbb{N}_{0})$ with $\ell\leq k$, it follows
that at least one of these sums will be nonzero such that
$\Psi^{(1)}_{k,\ell}(\mathrm{A}-\alpha_{2},\alpha_{2})>0$ and
$[(\mathrm{A}-\alpha_{2})_{\ell}(\alpha_{2})_{k-\ell}]^{-1}$ is a strictly
convex function of $\alpha_{2}$ on $(0,\mathrm{A})$. Combining this
observation with the fact that there is always at least one
$\ell\in\\{0,\dots,k\\}$ for which
$\tilde{g}_{k,\ell}(\zeta,{\nu^{\smash{\dagger}}})\neq 0$ further implies that
$f_{k}$ is a finite sum of at least one positive strictly convex function with
the rest of the terms being either being zero or also positive and strictly
convex; hence, $f_{k}$ is positive and strictly convex on
$\alpha_{2}\in(0,\mathrm{A})$.

Now turning to the limits notice that if $\zeta\in\\{0,\zeta_{b,1}\\}$ then by
Corollary 16 and Lemma 34 we have
$f_{k}(\alpha_{2}|\cdot)=\tilde{g}_{k,k}^{2}(\zeta,{\nu^{\smash{\dagger}}})/((\mathrm{A}-\alpha_{2})_{k}k!)\to
k!/(\mathrm{A})_{k}$ as $\alpha_{2}\nearrow 0$. For the complementary case
$\zeta\in[0,1]\setminus\\{0,\zeta_{b,1}\\}$ we may use the the asymptotic
expansion
$\frac{1}{(x)_{m}}\sim\frac{1}{\Gamma(m)}(\frac{1}{x}-H_{m-1}+\mathcal{O}(x))$
for $x\to 0$ and Lemma 34 to further deduce

$f_{k}(\alpha_{2}|\cdot)\sim\frac{k!}{(\mathrm{A})_{k}}+\mathcal{O}(\alpha_{2}^{-1})\to\infty$

as $\alpha_{2}\nearrow 0$. For the other limit $\alpha_{2}\searrow\mathrm{A}$,
Corollary 16 tells us that the $\ell=0$ term in $f_{k}$ is always positive;
hence, as $\alpha_{2}\searrow\mathrm{A}$

$f_{k}(\alpha_{2}|\cdot)\sim\mathcal{O}((\mathrm{A}-\alpha_{2})^{-1})\to\infty,$

which completes the proof. ∎

######  Theorem 18 (Optimal sample sizes of
$\mathscr{T}_{\nu^{\smash{\dagger}}}$).

Let
$(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}})(\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})$
with
$\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}}=(\zeta,\mathsf{acv}_{0},\mathsf{arb}_{0},b)$
denote the optimal sample sizes for $\mathscr{T}_{\nu^{\smash{\dagger}}}$ as a
function of $\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ that also satisfy
$\mathsf{ACV}\mathscr{T}_{\nu^{\smash{\dagger}}}(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=\mathsf{acv}_{0}$
and
$\mathsf{ARB}\mathscr{T}_{\nu^{\smash{\dagger}}}=\mathsf{arb}_{0}\zeta^{b}$
for $\mathsf{acv}_{0}\in\mathbb{R}^{+}$, $\mathsf{arb}_{0}\in(0,1)$, and
$b\in\mathbb{R}_{0}^{+}$ fixed. Then, for $\zeta\in\\{0,\zeta_{b,1}\\}$:

$(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}})(\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=(2\mathsf{acv}_{0}^{-2}+5,1).$

Likewise, for all remaining $\zeta\in[0,1]\setminus\\{0,\zeta_{b,1}\\}$:

$(n_{1}^{{\text{opt}}},n_{2}^{{\text{opt}}})(\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=(2\alpha_{1}^{{\text{opt}}}+1,2\alpha_{2}^{{\text{opt}}}+1)(\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}}),$

where $\alpha_{1}^{{\text{opt}}}$ and $\alpha_{2}^{{\text{opt}}}$ satisfy

$\displaystyle\partial_{\alpha_{2}}\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})\Big{|}_{(\mathrm{A},\alpha_{2})=(\alpha_{1}^{\text{opt}}+\alpha_{2}^{\text{opt}},\alpha_{2}^{\text{opt}})}$
$\displaystyle=0$ (3.14)
$\displaystyle\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\alpha_{1}^{\text{opt}},\alpha_{2}^{\text{opt}}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})-\mathsf{acv}_{0}^{2}$
$\displaystyle=0.$

###### Proof.

We begin with the special case $\zeta\in\\{0,\zeta_{b,1}\\}$ and call upon
Corollary 16 and Lemma 34 to find

$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\alpha_{1},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=\sum_{k=1}^{\infty}\frac{(1)_{k}(1)_{k}}{(\alpha_{1})_{k}k!}=\frac{1}{\alpha_{1}-2}.$

Since $\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ is independent of
$\alpha_{2}$ it follows that the optimization problem (3.14) is weakly
degenerate so that $\alpha_{2}^{\text{opt}}=0$ and $\alpha_{1}^{\text{opt}}$
satisfies

$\frac{1}{\alpha_{1}^{\text{opt}}-2}=\mathsf{acv}_{0}^{2}.$

Solving for $\alpha_{1}^{\text{opt}}$ and then using the relation
$n_{i}^{\text{opt}}=2n_{i}^{\text{opt}}+1$ then gives the desired result for
the special case $\zeta\in\\{0,\zeta_{b,1}\\}$.

Now consider the complementary case
$\zeta\in[0,1]\setminus\\{0,\zeta_{b,1}\\}$. Upon inspection of the double
integral representation of $\mathsf{E}\mathscr{T}_{\nu}^{2}$ in Theorem 9 and
assuming appropriate values of $(\zeta,{\nu^{\smash{\dagger}}})$ to guarantee
convergence, we see that
$\mathsf{E}\mathscr{T}_{\nu^{\smash{\dagger}}}^{2}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})$
is a smooth function of $\alpha_{2}$ and hence so is
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})$.
Additionally, upon writing

$\mathsf{ACV}^{2}\mathscr{T}_{\nu}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=\sum_{k=1}^{\infty}f_{k}(\alpha_{2}|\mathrm{A},\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}}),$

we conclude from Lemma 35 that
$\mathsf{ACV}^{2}\mathscr{T}_{\nu}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})$
is: $1)$ positive and strictly convex on $\alpha_{2}\in(0,\mathrm{A})$ and
$2)$ infinite at the endpoints $\alpha_{2}=0,\mathrm{A}$. As such, it follows
that for every set of parameters
$\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ with
$\zeta\notin\\{0,\zeta_{b,1}\\}$ that
$\alpha_{2}^{\text{opt}}(\mathrm{A},\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=\operatorname*{arg\,inf}_{\alpha_{2}\in(0,\mathrm{A})}\mathsf{ACV}^{2}\mathscr{T}_{\nu}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})$
is unique and corresponds to

$\partial_{\alpha_{2}}\mathsf{ACV}^{2}\mathscr{T}_{\nu}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})_{\alpha_{2}=\alpha_{2}^{\text{opt}}}=0.$

This is turn implies the optimization problem for finding
$(\alpha_{1}^{\text{opt}},\alpha_{2}^{\text{opt}})$ is nondegenerate and is
given by (3.14). The proof is now complete. ∎

Taking a small detour, we now relate the optimal sample sizes for
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ in Theorem 18 back to a classical result
given by James Janesick in his book _Photon Transfer_.

######  Remark 3 (Optimal sample sizes in the shot noise limit).

Using the expression for $\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$
in Corollary 15 as well as the results from Theorem 18 it is straightforward
to deduce

$\lim_{\zeta\nearrow 0}\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}=1,$

which shows that the optimal sample sizes for
$\mathscr{T}_{\nu^{\smash{\dagger}}}$ equal those for
$\mathscr{G}_{\nu^{\smash{\dagger}}}$ in the shot noise limit. If one wants to
obtain a $1\%$ relative uncertainty for the measurement of $g$ in the shot
noise limit, i.e. $\lim_{\zeta\nearrow
0}\mathsf{ACV}\mathscr{G}_{\nu^{\smash{\dagger}}}=0.01$, then according to
Theorem 18 the optimal total number of samples needed is
$\mathrm{N}^{\text{opt}}=20,006$. This exact result agrees with Janesick’s
estimate $\mathrm{N}^{\text{opt}}\approx 20,000$ of [16, pgs. 79–81].

Up to this point we have made a significant amount of progress in studying the
optimal sample sizes of $\mathscr{T}_{\nu^{\smash{\dagger}}}$. In particular,
we have established uniqueness of these optimal sample sizes, derived their
exact expressions at the two points of interest $\zeta\in\\{0,\zeta_{b,1}\\}$,
determined a system of equations for evaluating them at all remaining
$\zeta\in[0,1]\setminus\\{0,\zeta_{b,1}\\}$, studied their asymptotic behavior
in the low illumination limit $\zeta\searrow 1$, and derived the expression
$\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$, which allow us to
determine when they give a good approximation to the optimal sample sizes of
$\mathscr{G}_{\nu^{\smash{\dagger}}}$. With so much theoretical understanding
of these optimal sample sizes we’re in a good place to begin discussing how we
actually compute them.

It is of no surprise that the system of equations for evaluating the optimal
sample sizes given in (3.14) lacks a closed-from solution and so we must turn
to numerical methods to solve it. This of course introduces the practical
problem of implementing root finding algorithms. We already know the bounds of
the parameters $\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ and so the only
thing left to implement (3.14) is a good approximation for
$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ to use as a starting point. Because
$\mathsf{ACV}\mathscr{T}_{\nu^{\smash{\dagger}}}$ is naturally parameterized
in terms of $(\alpha_{1},\alpha_{2})$ we will seek an approximate solution for
the optimal pair $(\alpha_{1}^{\text{opt}},\alpha_{2}^{\text{opt}})$ and then
transform it into $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ via the relation
$n_{i}^{\text{opt}}=2\alpha_{i}^{\text{opt}}+1$.

To begin deriving our approximation we call on the series expansion given in
Corollary 12 to write
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\alpha_{1},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})\approx\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}^{\ast}(\alpha_{1},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})$,
where

$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}^{\ast}(\alpha_{1},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=\frac{1}{\alpha_{1}}\tilde{g}_{1,1}^{2}(\zeta,{\nu^{\smash{\dagger}}})+\frac{1}{\alpha_{1}\alpha_{2}}\tilde{g}_{2,1}^{2}(\zeta,{\nu^{\smash{\dagger}}})+\frac{1}{\alpha_{2}}\tilde{g}_{1,0}^{2}(\zeta,{\nu^{\smash{\dagger}}})$
(3.15)

and $\tilde{g}_{k,\ell}(z,\nu)$ given by Corollary 10 $(\mathrm{i})$. Upon
substituting $\alpha_{1}=\mathrm{A}-\alpha_{2}$ into this approximation we see
that it possesses the same useful properties of the exact form of
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})$,
namely, positivity, smoothness, strictly convexity on
$\alpha_{2}\in(0,\mathrm{A})$, as well as the same limiting properties at the
endpoints $\alpha_{2}=0,\mathrm{A}$. For the sake of brevity we let
$a=\tilde{g}_{1,1}(\zeta,{\nu^{\smash{\dagger}}})$,
$b=\tilde{g}_{2,1}(\zeta,{\nu^{\smash{\dagger}}})$,
$c=\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}})$, and $d=\mathsf{acv}_{0}$
and then substitute
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}^{\ast}$ into (3.14)
yielding the system of equations

$\displaystyle(1)$
$\displaystyle\frac{a^{2}}{\alpha_{1}^{2}}+\frac{b^{2}}{\alpha_{1}^{2}\alpha_{2}}-\frac{b^{2}}{\alpha_{1}\alpha_{2}^{2}}-\frac{c^{2}}{\alpha_{2}^{2}}$
$\displaystyle=0$ $\displaystyle(2)$
$\displaystyle\frac{a^{2}}{\alpha_{1}}+\frac{b^{2}}{\alpha_{1}\alpha_{2}}+\frac{c^{2}}{\alpha_{2}}-d^{2}$
$\displaystyle=0.$

To put this system into a more useful form we perform the transformations
$(1)^{\prime}=(2)/\alpha_{1}-(1)$ and $(2)^{\prime}=(2)/\alpha_{2}+(1)$ to get
the equivalent set of equations

$\displaystyle(1)^{\prime}$ $\displaystyle\alpha_{1}$
$\displaystyle=\frac{1}{c^{2}}\left(d^{2}\alpha_{2}^{2}-c^{2}\alpha_{2}-b^{2}\right)$
$\displaystyle(2)^{\prime}$ $\displaystyle\alpha_{2}$
$\displaystyle=\frac{1}{a^{2}}\left(d^{2}\alpha_{1}^{2}-a^{2}\alpha_{1}-b^{2}\right).$

Now substituting the r.h.s. of $(1)^{\prime}$ into $(2)^{\prime}$ and the
r.h.s. of $(2)^{\prime}$ into $(1)^{\prime}$ separates the variables
$\alpha_{1}$ and $\alpha_{2}$ into two quartic equations, each of which yields
four roots for a total of sixteen possible solutions to our problem.
Substituting all sixteen potential solutions back into our original system we
find only two work, namely,

$(\alpha_{1},\alpha_{2})=\left(\frac{a^{2}\pm\sqrt{a^{2}c^{2}+b^{2}d^{2}}}{d^{2}},\frac{c^{2}\pm\sqrt{a^{2}c^{2}+b^{2}d^{2}}}{d^{2}}\right).$
(3.16)

To obtain the final solution we must determine which sign in front of the root
to take. Taking the positive sign clearly yields a nonnegative solution for
$\alpha_{1}$ and $\alpha_{2}$ so to verify this is the correct choice we must
show that choosing the negative sign yields a nonsensical solution. The
following two lemmas show that if we choose the negative sign then the
solution for $\alpha_{2}$ is nonpositive and thus is not the desired choice.
Since the choice of sign for both $\alpha_{1}$ and $\alpha_{2}$ must be the
same it follows that the positive sign is indeed the correct choice.

######  Lemma 36.

Let ${\nu^{\smash{\dagger}}}$ be as defined in Definition 13. Then for all
$\zeta\in(0,1)$:
$\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}})\geq-\frac{1}{2}$.

###### Proof.

We begin with the explicit form

$\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}})=\frac{\zeta}{1-\zeta}-\frac{(\log_{\zeta}\mathsf{arb}_{0}+b)\mathsf{arb}_{0}\zeta^{b}}{1-\mathsf{arb}_{0}\zeta^{b}}.$

Our approach will be to put a lower bound on
$\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}})$ by successively minimizing it
w.r.t. each of its variables. Differentiating w.r.t. $b$ we find

$\partial_{b}\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}})=\frac{x(x-\log
x-1)}{(1-x)^{2}}\Big{|}_{x=\mathsf{arb}_{0}\zeta^{b}}.$

For the specified parameter restriction we have
$0<\mathsf{arb}_{0}\zeta^{b}<\mathsf{arb}_{0}<1\implies x\in(0,1)$; thus, it
is straightforward to show that the above expression in $x$ is positive. It
follows that

$\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}})\geq\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}}_{b=0})=\frac{\zeta}{1-\zeta}-\frac{(\log_{\zeta}\mathsf{arb}_{0})\mathsf{arb}_{0}}{1-\mathsf{arb}_{0}}.$

Now differentiating w.r.t. $\mathsf{arb}_{0}$ we have

$\partial_{\mathsf{arb}_{0}}\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}}_{b=0})=\frac{1}{\log\zeta}\frac{\mathsf{arb}_{0}-\log\mathsf{arb}_{0}-1}{(1-\mathsf{arb}_{0})^{2}}.$

The quantity involving $\mathsf{arb}_{0}$ is positive on
$\mathsf{arb}_{0}\in(0,1)$ while $\log\zeta<0$ on $\zeta\in(0,1)$. Therefore,
$\partial_{\mathsf{arb}_{0}}\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}}_{b=0})<0$
and

$\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}}_{b=0})\geq\lim_{\mathsf{arb}_{0}\searrow
1}\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}}_{b=0})=\frac{\zeta}{1-\zeta}+\frac{1}{\log\zeta}.$

The resulting function of $\zeta$ is strictly decreasing in $\zeta$. Taking
the limit $\zeta\searrow 1$ subsequently gives

$\tilde{g}_{1,0}(\zeta,{\nu^{\smash{\dagger}}}_{b=0})\geq\lim_{\zeta\searrow
1}\frac{\zeta}{1-\zeta}+\frac{1}{\log\zeta}=-\frac{1}{2}.$

The proof is now complete. ∎

######  Lemma 37.

Choosing the negative sign in $(\ref{eq:aprx_system_potential_solutions})$
leads to $\alpha_{2}\leq 0$.

###### Proof.

For brevity we will use $\tilde{g}_{n,\omega}$ to denote
$\tilde{g}_{n,\omega}(\zeta,{\nu^{\smash{\dagger}}})$. We need to show

$\tilde{g}_{1,0}^{2}-\sqrt{\tilde{g}_{1,0}^{2}\tilde{g}_{1,1}^{2}+\tilde{g}_{2,0}^{2}\mathsf{acv}_{0}^{2}}\leq
0.$

Since all quantities under the radical are positive it follows that

$\tilde{g}_{1,0}^{2}-\sqrt{\tilde{g}_{1,0}^{2}\tilde{g}_{1,1}^{2}+\tilde{g}_{2,0}^{2}\mathsf{acv}_{0}^{2}}\leq\tilde{g}_{1,0}^{2}-\sqrt{\tilde{g}_{1,0}^{2}\tilde{g}_{1,1}^{2}}.$

From Lemma 13 we have $\tilde{g}_{1,1}=\tilde{g}_{1,0}+1$; hence,

$\tilde{g}_{1,0}^{2}-\sqrt{\tilde{g}_{1,0}^{2}\tilde{g}_{1,1}^{2}+\tilde{g}_{2,0}^{2}\mathsf{acv}_{0}^{2}}\leq\tilde{g}_{1,0}^{2}-\sqrt{\tilde{g}_{1,0}^{4}+\tilde{g}_{1,0}^{2}(2\tilde{g}_{1,0}+1)}.$

But now recall from Lemma 36 that
$\tilde{g}_{1,0}\geq-\frac{1}{2}\implies\tilde{g}_{1,0}^{2}(2\tilde{g}_{1,0}+1)\geq
0$. So it follows that

$\tilde{g}_{1,0}^{2}-\sqrt{\tilde{g}_{1,0}^{2}\tilde{g}_{1,1}^{2}+\tilde{g}_{2,0}^{2}\mathsf{acv}_{0}^{2}}\leq\tilde{g}_{1,0}^{2}-\sqrt{\tilde{g}_{1,0}^{4}}=0,$

which completes the proof. ∎

Now knowing that we must choose the positive sign in
$(\ref{eq:aprx_system_potential_solutions})$ we have after reintroducing $a$,
$b$, $c$, $d$, and $n_{i}=2\alpha_{i}+1$:

$\displaystyle n_{1}$
$\displaystyle=\frac{2}{\mathsf{acv}_{0}^{2}}\left(\tilde{g}_{1,1}^{2}+\sqrt{\tilde{g}_{1,1}^{2}\tilde{g}_{1,0}^{2}+\tilde{g}_{2,1}^{2}\mathsf{acv}_{0}^{2}}\right)+1$
$\displaystyle n_{2}$
$\displaystyle=\frac{2}{\mathsf{acv}_{0}^{2}}\left(\tilde{g}_{1,0}^{2}+\sqrt{\tilde{g}_{1,1}^{2}\tilde{g}_{1,0}^{2}+\tilde{g}_{2,1}^{2}\mathsf{acv}_{0}^{2}}\right)+1,$

which is the desired solution for the approximate optimal sample sizes of
$\mathscr{T}_{\nu^{\smash{\dagger}}}$. So how good are these approximations?
First note that with the results of Corollary 16 and Lemma 34 we have

$\lim_{\zeta\nearrow
L}(n_{1},n_{2})(\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})=(2\mathsf{acv}_{0}^{-2}+1,1),\quad
L\in\\{0,\zeta_{b,1}\\}.$

Comparing this with Theorem 18 we see that our approximate solutions can be
made exact at $\zeta=\\{0,\zeta_{b,1}\\}$ by instead using the slightly
modified solution

$\displaystyle n_{1}^{\ast}$
$\displaystyle=\frac{2}{\mathsf{acv}_{0}^{2}}\left(\tilde{g}_{1,1}^{2}+\sqrt{\tilde{g}_{1,1}^{2}\tilde{g}_{1,0}^{2}+\tilde{g}_{2,1}^{2}\mathsf{acv}_{0}^{2}}\right)+5$
(3.17) $\displaystyle n_{2}^{\ast}$
$\displaystyle=\frac{2}{\mathsf{acv}_{0}^{2}}\left(\tilde{g}_{1,0}^{2}+\sqrt{\tilde{g}_{1,1}^{2}\tilde{g}_{1,0}^{2}+\tilde{g}_{2,1}^{2}\mathsf{acv}_{0}^{2}}\right)+1.$

At the opposite end of the $\zeta$-domain, we may use Lemma 30 to deduce as
$\zeta\searrow 1$: $n_{2}^{\ast}/n_{1}^{\ast}\to 1$ and $n_{i}^{\ast}\sim
C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}^{\ast}(1-\zeta)^{-2}$, where

$C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}^{\ast}=\frac{2\gamma^{2}(2,-\log\mathsf{arb}_{0})}{\mathsf{acv}_{0}^{2}(1-\mathsf{arb}_{0})^{2}}\left(1+\left(1+\mathsf{acv}_{0}^{2}(1-\mathsf{arb}_{0})^{2}\frac{\gamma^{2}(3,-\log\mathsf{arb}_{0})}{\gamma^{4}(2,-\log\mathsf{arb}_{0})}\right)^{1/2}\right).$

Table 3.4 tabulates the ratio
$C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}^{\ast}/C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$
for a variety of values for $\mathsf{arb}_{0}$ and $\mathsf{acv}_{0}$. One can
see that the ratio is near unity for these values which indicates our
approximate solutions for $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ given by
(3.17) are not only exact at $\zeta=\\{0,\zeta_{b,1}\\}$ but are also
exceptionally accurate for small $\mathsf{arb}_{0}$ and $\mathsf{acv}_{0}$ in
the neighborhood of $\zeta=1$.

$\mathsf{arb}_{0}/\mathsf{acv}_{0}$ | $0.01$ | $0.02$ | $0.05$ | $0.10$ | $0.20$
---|---|---|---|---|---
$0.01$ | $0.9999$ | $0.9999$ | $0.9999$ | $0.9999$ | $0.9999$
$0.02$ | $0.9997$ | $0.9997$ | $0.9997$ | $0.9997$ | $0.9998$
$0.05$ | $0.9978$ | $0.9980$ | $0.9982$ | $0.9984$ | $0.9985$
$0.10$ | $0.9913$ | $0.9919$ | $0.9927$ | $0.9934$ | $0.9942$
$0.20$ | $0.9655$ | $0.9679$ | $0.9714$ | $0.9743$ | $0.9771$
Table 3.4:
$C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}^{\ast}/C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$
for various values of $\mathsf{arb}_{0}$ and $\mathsf{acv}_{0}$.

With a suitable approximation to the optimal sample sizes we may finally turn
to implementing a numerical routine for computing them. Numerically solving
the system (3.14) ultimately requires a method for computing
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ and its partial
derivative w.r.t. $\alpha_{2}$. The fact that the approximation for
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ in (3.15) lead to
excellent approximations for the optimal sample sizes across the entire
$\zeta$-domain tells us that we may compute the optimal sample sizes by
substituting

$\mathsf{ACV}_{n}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}=\sum_{k=1}^{n}\sum_{\ell=0}^{k}\frac{\tilde{g}_{k,\ell}^{2}(\zeta,{\nu^{\smash{\dagger}}})}{(\alpha_{1})_{\ell}(\alpha_{2})_{k-\ell}\ell!(k-\ell)!}.$
(3.18)

into (3.14) for some appropriate value of $n$. To determine an appropriate
value for $n$ we consider the absolute relative error

$R_{n,2}=\left\lvert\frac{\mathsf{ACV}_{n}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}-1\right\rvert,$

which according to Corollary 17 is bounded above by $R_{n,2}\leq
R_{n,m,2}^{\ast}$ with

$R_{n,m,2}^{\ast}=\left\lvert\left(1+\frac{E_{n,m}^{\ast}}{\mathsf{ACV}_{n}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}\right)^{-1}-1\right\rvert$

and $E_{n,m}^{\ast}$ given in Theorem 20. Before proceeding, the following
result gives us an exact value for $R_{n,2}$ at $\zeta\in\\{0,\zeta_{b,1}\\}$.

######  Lemma 38.

For $\zeta\in\\{0,\zeta_{b,1}\\}$

$R_{n,2}|_{(n_{1},n_{2})=(n_{1}^{\text{opt}},n_{2}^{\text{opt}})}=\frac{(n+1)!}{(\mathsf{acv}_{0}^{-2}+2)_{n}}.$

###### Proof.

The proof follows from noting that at $\zeta=0,\zeta_{b,1}$:

$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}=\frac{1}{\alpha_{1}-2},$
$\mathsf{ACV}_{n}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}=\frac{1}{\alpha_{1}-2}\left(1-\frac{(n+1)!}{(\alpha_{1})_{n}}\right),$

and $\alpha_{1}^{\text{opt}}=\mathsf{acv}_{0}^{-2}+2$. ∎

The sharpness of the upper bound $R_{n,m,2}^{\ast}$ improves with increasing
$m$ but the rate at which this happens significantly decreases for $m>n+1$;
thus, it seems reasonable to use $R_{n,n+1}^{\ast}$ as the upper bound in our
numerical routine. To get a sense of how $R_{n,n+1,2}^{\ast}$ behaves, we
evaluated it at $(n_{1},n_{2})=(n_{1}^{\ast},n_{2}^{\ast})$ for
$\mathsf{arb}_{0}=0.5$, $\mathsf{acv}_{0}=0.05$, and various values of $b$ and
$n$ as plotted in Figure 3.2. We observe that as $n$ increases, the upper
bound decreases rapidly except for values of $\zeta$ in the neighborhood of
$\zeta=0,\zeta_{b,1},1$. However, notice that for the exact value of $R_{n,2}$
given by Lemma 38: $R_{1,2}=4.98\times 10^{-3}$, $R_{3,2}=3.67\times 10^{-7}$,
$R_{5,2}=6.69\times 10^{-11}$, $R_{7,2}=2.26\times 10^{-14}$. Since the
approximations $n_{1}^{\ast}$ and $n_{2}^{\ast}$ are exact as $\zeta\to L$,
$L\in\\{0,\zeta_{b,1}\\}$, we can conclude that the high upper bound for
$R_{n,2}$ near these values of $\zeta$ is the result of a deficiency of
$R_{n,n+1,2}^{\ast}$ and not that the relative error isn’t decrease rapidly
with $n$. As such, estimates for the appropriate value of $n$ in (3.18)
determined using $R_{n,n+1,2}^{\ast}$ will be, in some cases, significantly
overestimated near these points.

Figure 3.2: Absolute relative error bound $R_{n,n+1,2}^{\ast}$ versus $\zeta$
for $\mathsf{arb}_{0}=0.5$, $\mathsf{acv}_{0}=0.05$, $b=0,1$, and $n=1,3,5,7$.

Given these observations, Algorithm 1 presents a simple procedure for
computing the optimal sample sizes for $\mathscr{T}_{\nu^{\smash{\dagger}}}$.

Algorithm 1 Optimal sample size pseudo-code.

1:procedure OptimalSamples($\zeta$, $\mathsf{arb}_{0}$, $b$,
$\mathsf{acv}_{0}$, $\epsilon$, $n_{\max}$)

2: if $\zeta\in\\{0,\zeta_{b,1}\\}$ then

3: $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})=(2\mathsf{acv}_{0}^{-2}+5,1)$;

4: return $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$;

5: else

6: $n=1$;

7: while $R_{n,n+1,2}^{\ast}(\zeta,n_{1}^{\ast},n_{2}^{\ast})>\epsilon$ and
$n<n_{\max}$ do

8: $n=n+1$;

9: Substitute $\mathsf{ACV}_{n}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ in
(3.14);

10: Solve for $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ using starting point
$(n_{1}^{\ast},n_{2}^{\ast})$;

11: if
$R_{n,n+1,2}^{\ast}(\zeta,n_{1}^{\text{opt}},n_{2}^{\text{opt}})\leq\epsilon$
or $n=n_{\max}$ then

12: return $(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$;

13: else

14: $n=n+1$;

15: goto line 9;

To demonstrate its effectiveness, Algorithm 1 was implemented in mathematica
using the FindRoot[] function for Line 10 to numerically solve the system
(3.14). The procedure was executed for $\zeta\in[0,1)$,
$\mathsf{arb}_{0}=0.5$, $\mathsf{acv}_{0}=0.05$, and four different values of
$b$ ranging from zero to one. Figure 3.3 plots the computed sample sizes along
with $\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ evaluated at a dark
noise of $\sigma_{\mathrm{d}}g=1e\mathchar 45\relax$. The values chosen for
$\mathsf{arb}_{0}$ and $\sigma_{\mathrm{d}}g$, while not typical, were chosen
to exaggerate the features of the optimal sample size curves discussed in the
preceding analysis.

Focusing first on the curve for $n_{2}^{\text{opt}}$ we observe for values of
$b<1$ that the curve has a cusp at $\zeta=0$ and $\zeta=\zeta_{b,1}$ where it
takes on a value of one. As, $\zeta\searrow 1$ this curve approaches the
$n_{1}^{\text{opt}}$ curve; both of which approach
$C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}(1-\zeta)^{-2}$ in the limit. As for
the $\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ curve, we see that it
equals one at $\zeta=0$ and
$\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}=0.15847\dots$ at
$\zeta=1$. Furthermore, the cusp present in the $n_{2}^{\text{opt}}$ curve is
also present in the curve for
$\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$, which highlights the
reason why choosing values of $b$ less than one may yield undesirable results,
in particular, unacceptably low values of
$\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$.

As for the computational details of the routine, it was decided to use
$\epsilon=5\times 10^{-7}$ so that
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ was approximated to at
least seven significant digits and $n_{\max}=20$ to avoid issues of numerical
stability. Figure 3.4 plots the number of terms, $n$, and upper bound on the
absolute relative error, $R_{n,n+1,2}^{\ast}$, as a function of $\zeta$
obtained from the routine for the case $b=0$. From the figure we observe that
$R_{n,n+1,2}^{\ast}$ exhibits a sharp increase near the points
$\zeta\in\\{0,\zeta_{b,1}\\}$; resulting in larger values for $n$. Only a few
$\zeta$-values near $\zeta_{b,1}$ hit the bound $n=20$, which is the same
behavior observed for all other values of $b$ in Figure 3.3. To check the
quality of the solution at these problematic points we note that in the
neighborhood of $\zeta=\zeta_{b,1}$ that ${\nu^{\smash{\dagger}}}$ and
$n_{2}^{\text{opt}}$ are quite small while $n_{1}^{\text{opt}}$ is large.
Recalling the double integral representation of
$\mathsf{E}\mathscr{T}_{\nu^{\smash{\dagger}}}^{2}$ we may then approximate
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ with

$\mathsf{ACV}_{n}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}=\int_{0}^{1}\int_{0}^{1}\sum_{k=0}^{n}\frac{(1+{\nu^{\smash{\dagger}}})_{k}^{2}}{(\alpha_{1})_{k}k!}\frac{\left((1-x)(1-y)\right)^{k}}{\left((1-(1-\zeta)x)(1-(1-\zeta)y)\right)^{k}}\\\
\times{F}\left({1-{\nu^{\smash{\dagger}}},1-{\nu^{\smash{\dagger}}}\atop\alpha_{2}};(1-x)(1-y)\right)f_{XY}(x,y)\,\mathrm{d}x\mathrm{d}y-1,$
(3.19)

where

$f_{XY}(x,y)=\frac{\left((1-(1-\zeta)x)(1-(1-\zeta)y)\right)^{-{\nu^{\smash{\dagger}}}-1}}{F(1,1+{\nu^{\smash{\dagger}}};2;1-\zeta)^{2}}$

and the error $E_{n}$ such that
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}=\mathsf{ACV}_{n}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}+E_{n}$
is positive and bounded above by

$E_{n}^{\ast}={F}\left({1+{\nu^{\smash{\dagger}}},1+{\nu^{\smash{\dagger}}}\atop\alpha_{1}};1\right)_{\\!\\!n+1}^{\\!c}\int_{0}^{1}\int_{0}^{1}\frac{\left((1-x)(1-y)\right)^{n+1}}{\left((1-(1-\zeta)x)(1-(1-\zeta)y)\right)^{n+1}}\\\
\times{F}\left({1-{\nu^{\smash{\dagger}}},1-{\nu^{\smash{\dagger}}}\atop\alpha_{2}};(1-x)(1-y)\right)f_{XY}(x,y)\,\mathrm{d}x\mathrm{d}y$

with $F(\alpha,\beta;\gamma;z)_{n}^{c}$ given in Proposition 9. Looking back
as the numerical data we observe that one of the problematic points in Figure
3.4 occurs at $\zeta=0.499443$. Using this value for $\zeta$ and the optimal
sample sizes computed from our routine we find

$\left\lvert\frac{\mathsf{ACV}_{20}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}{\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}}-1\right\rvert\leq
6.21\times 10^{-15},$

which shows that the approximation
$\mathsf{ACV}_{20}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ given by (3.19)
approximates the exact value
$\mathsf{ACV}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ to approximately
fourteen significant digits at this point of interest. Substituting
$\mathsf{ACV}_{20}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}$ into (3.14) we then
find

$\displaystyle\partial_{\alpha_{2}}\mathsf{ACV}_{20}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\mathrm{A}-\alpha_{2},\alpha_{2}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})\Big{|}_{(\mathrm{A},\alpha_{2})=(\alpha_{1}^{\text{opt}}+\alpha_{2}^{\text{opt}},\alpha_{2}^{\text{opt}})}$
$\displaystyle=-6.84\times 10^{-9}$
$\displaystyle\mathsf{ACV}_{20}^{2}\mathscr{T}_{\nu^{\smash{\dagger}}}(\alpha_{1}^{\text{opt}},\alpha_{2}^{\text{opt}}|\theta_{\mathscr{T}_{\nu^{\smash{\dagger}}}})-\mathsf{acv}_{0}^{2}$
$\displaystyle=9.48\times 10^{-10},$

which shows no significant deviation from the desired solution $(0,0)$. This
confirms that the solution to the optimal sample sizes obtained in our routine
near $\zeta=\zeta_{b,1}$ are not in significant error as the upper bound
$R_{n,n+1,2}^{\ast}$ might suggest.

Figure 3.3: $\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ and
$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})$ versus $\zeta$ for
$\mathsf{arb}_{0}=0.5$, $\mathsf{acv}_{0}=0.05$, $b=\\{0,0.369,0.613,1\\}$,
and $\sigma_{\mathrm{d}}g=1e\mathchar 45\relax$. Values of $\zeta_{b,1}$ and
$\overline{\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}}$ are also
indicated on each plot. Figure 3.4: $R_{n^{\ast},n^{\ast}+1,2}^{\ast}$ (top)
and $n^{\ast}$ (bottom) versus $\zeta$ for $\mathsf{arb}_{0}=0.5$, $b=0$, and
$\mathsf{acv}_{0}=0.05$.

###  3.6 Pixel-level conversion gain estimation: An introduction

Estimating conversion gain with the pt method is typically performed under the
assumption that each pixel in the sensor array exhibits identical
characteristics. When this assumption holds we say the sensor is _uniform_ and
treat each pixel as a source of i.i.d. random noise. Such an assumption is
very convenient for performing pt characterization as we can leverage the
large number of pixels in a sensor array to obtain low-variance estimates of
critical population values. Take for example the estimate of population dark
variance $\sigma_{\mathrm{d}}^{2}$. When the assumption of uniformity holds,
one can estimate this quantity from just two images captured under dark
conditions. Denoting $(D_{1},D_{2})$ as these two dark images each with $R$
rows of pixels and $C$ columns of pixels we first compute the pixelwise
difference $\Delta D=D_{1}-D_{2}$. Since this differencing operation
effectively doubles the population variance we then estimate
$\sigma_{\mathrm{d}}^{2}$ by computing half of the sample variance of the
difference-frame

$\hat{\sigma}_{\mathrm{d}}^{2}=\frac{1}{2}\frac{1}{RC-1}\sum_{i=1}^{R}\sum_{j=1}^{C}(\Delta
D_{ij}-\overline{\Delta D})^{2},$

where $\Delta D_{ij}$ denotes the value of the difference-frame $\Delta D$ at
row $i$ and column $j$ and $\overline{\Delta
D}=\frac{1}{RC}\sum_{i=1}^{R}\sum_{j=1}^{C}\Delta D_{ij}$ denotes the sample
mean of the difference-frame. For a sensor with a very modest one megapixel
resolution $RC=10^{6}$ and so from just two images we can obtain an estimate
of $\sigma_{\mathrm{d}}^{2}$ based on an effective sample size of one-million
i.i.d. observations.

In actual sensors, the assumption of uniformity is rarely a good model and is
only useful for a small number of pixel architectures. Indeed, in Janesick’s
book _Photon Transfer_ he points out that nonuniformity is a key
characteristic of the _Complementary Metal-Oxide Semiconductor_ (cmos) sensor
architecture, which is by far the most prevalent architecture in modern
imaging systems [16, pg. 82]. When the uniformity assumption is violated we no
can longer provide global estimates of the necessary population parameters and
must turn to estimates for each individual pixel. Going back to the example of
estimating $\sigma_{\mathrm{d}}^{2}$, pixel-level estimation is performed by
capturing a sequence of images, i.e. an _image stack_ , and then computing the
pixelwise temporal variance of the stack. Supposing we capture a stack of $N$
dark images $(D_{1},\dots,D_{N})$ to form a $R\times C\times N$ dark image
stack, the estimate of $\sigma_{\mathrm{d}}^{2}$ for the pixel in row $i$ and
column $j$ becomes

$(\hat{\sigma}_{\mathrm{d}}^{2})_{ij}=\frac{1}{N-1}\sum_{k=1}^{N}(D_{ijk}-\overline{D}_{ij})^{2},$
(3.20)

where $D_{ijk}$ is the value of pixel $(i,j)$ in the $k$th image and
$\overline{D}_{ij}=\frac{1}{N}\sum_{k=1}^{N}D_{ijk}$ is the temporal sample
mean of that same pixel. From this, we can immediately see the challenge
introduced in pixelwise estimation, for in order to estimate
$\sigma_{\mathrm{d}}^{2}$ for each pixel with the same level of uncertainty as
in the uniform case a total of $N=10^{6}$ images must be captured! Setting
aside the shear ammount of data this entails we also run into the practical
problem of simply capturing the images in a short enough period of time so
that drift does not corrupt the population parameters we are attempting to
measure. So while the extension to pixelwise estimation is conceptually
straightforward, implementation can be quite challenging.

Given that uniformity is generally the exception and not the rule, a
generalized pixel-level approach to pt characterization is in order. Our goal
here is to take the theoretical results presented in this work and apply them
to develop a first attempt at the process of pixel-level conversion gain
estimation. Before this is done we must discuss the topic of _pixel grouping_.

####  3.6.1 Pixel grouping

When the assumption of uniformity is removed and pixel-level measurement is
performed, the characterization of each individual pixel effectively becomes
its own experiment with its own unique optimal sample sizes. However, since we
can only simultaneously sample all pixels by capturing full-frame images, the
conclusion of any pixel-level characterization procedure will be a single set
of sample sizes that must work for (nearly) all pixels in the sensor
array555We say nearly because typical sensors always contain some minority of
defect pixels which are either non-functioning or exhibit statistical behavior
too extreme to be properly characterized.. Knowing that estimates for each
pixel will ultimately be computed with the same number of samples,
individually estimating optimal sample sizes is nonsensical and thus we should
group pixels with similar characteristics to reduce the variance in estimates
of the optimal sample sizes. This point is made evident by first recalling the
asymptotic behavior of the optimal sample sizes in the low-illumination limit
$\zeta\searrow 1$

$(n_{1}^{\text{opt}},n_{2}^{\text{opt}})(\zeta)\sim(C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}(1-\zeta)^{-2},C_{\mathscr{T}_{\nu^{\smash{\dagger}}}}(1-\zeta)^{-2}).$

Letting $\hat{Y}=\hat{\sigma}_{\mathrm{d}}^{2}$ be an estimate of
$\sigma_{\mathrm{d}}^{2}$ according to (3.20) and
$\hat{X}=\hat{\sigma}_{\mathrm{p}+\mathrm{d}}^{2}$ be the analogous estimate
of $\sigma_{\mathrm{p}+\mathrm{d}}^{2}$ for a single pixel based on dark and
illuminated image stacks of $n_{2}$ and $n_{1}$ images, respectively, we can
estimate $\zeta$ with

$Z=\frac{n_{1}-3}{n_{1}-1}\frac{\hat{Y}}{\hat{X}}.$

Under the assumed normal model, $Z$ is the umvue for $\zeta$ with
$\mathsf{E}Z=\zeta$ and

$\mathsf{Var}Z=\left(\frac{(n_{1}-3)(n_{2}+1)}{(n_{1}-5)(n_{2}-1)}-1\right)\zeta^{2}.$

While $Z$ possesses good characteristics as an estimator for $\zeta$ its
density is nonzero in the neighborhood of $\zeta=1$ and so the moments of the
optimal sample size estimates
$(\hat{n}_{1}^{\text{opt}},\hat{n}_{2}^{\text{opt}})=(n_{1}^{\text{opt}},n_{2}^{\text{opt}})(Z)$
are strictly speaking undefined. However, if instead of measuring $Z$ for each
pixel individually we identify a group of $m$ pixels with similar population
values for $\sigma_{\mathrm{d}}^{2}$ and $\sigma_{\mathrm{p}+\mathrm{d}}^{2}$
we may estimate $\zeta$ for the entire group with

$Z=\frac{mn_{1}-3}{mn_{1}-1}\frac{\sum_{\text{group}}\hat{Y}_{\ell}}{\sum_{\text{group}}\hat{X}_{\ell}},$

where $\hat{Y}_{\ell}$ and $\hat{X}_{\ell}$ denotes the estimates for the
$\ell$th pixel in the group. Grouping in this manner does not alter the
unbiasedness of $Z$ and effectively increases the dark and illuminated sample
sizes by a factor of $m$. Hence, given sufficiently large $n_{1}$, $n_{2}$,
and group size $m$, one may drive down the asymptotic variance of the optimal
sample size estimates

$\mathsf{Var}\,\hat{n}_{i}^{\text{opt}}\sim(\partial_{\zeta}n_{i}^{\text{opt}}(\zeta))^{2}\mathsf{Var}Z,\quad
i=1,2$

to the point where they becomes useful.

Another important reason why grouping pixels in this manner is not only
important but practically necessary is in estimating the quantity
${\nu^{\smash{\dagger}}}(\zeta)=\log_{\zeta}(\mathsf{arb}_{0})+b$. In a
similar manner to the optimal sample sizes we see as $\zeta\searrow 1$

${\nu^{\smash{\dagger}}}(\zeta)\sim-\log\mathsf{arb}_{0}(1-\zeta)^{-1}.$

The pole of order one at $\zeta=1$ again implies that the estimate
${{\hat{\nu}}^{\smash{\dagger}}}={\nu^{\smash{\dagger}}}(Z)$ has undefined
moments although we may assign it an asymptotic variance

$\mathsf{Var}\,{{\hat{\nu}}^{\smash{\dagger}}}\sim\frac{\log^{2}_{\zeta}\mathsf{arb}_{0}}{\zeta^{2}\log^{2}\zeta}\mathsf{Var}Z.$

Since we ultimately want an estimate of ${\nu^{\smash{\dagger}}}$ with nearly
zero variance666Recall that the derivation of the mean and variance of the
estimator $\mathscr{G}_{\nu}$ assumed $\nu$ is a known constant., we need to
turn to pixel grouping to leverage the asymptotics of large sample sizes to
produce low-variance estimates of ${\nu^{\smash{\dagger}}}$ for groups of
similar pixels.

Identifying appropriate pixel groups for pixel-level characterization is
largely determined by the architecture of the sensor under test. In what
follows, we will demonstrate how to identify pixels groups through example as
we step through the process of pixel-level conversion gain estimation for a
real image sensor.

####  3.6.2 Experimental setup

For this experiment, the ON Semiconductor KAI-0407M monochrome interline
transfer _Charge-Coupled Deivce_ (ccd) sensor was chosen as an initial test
case. This particular sensor was chosen for two main reasons. First,
preliminary tests show that the noise it produces closely adheres to the
normal model assumed in the development of the estimator $\mathscr{G}_{\nu}$.
Secondly, interline transfer ccds typically exhibit only minor
nonuniformities. Since we expect significant nonuniformities to complicate the
experimental procedure, this particular sensor provides a gentle introduction
into pixel-level characterization.

The experimental setup began with a $650\,\mathrm{nm}$ superluminescent light
emitting diode (sled), which was collimated and attached to a _Variable
Optical Attenuator_ (voa) to facilitate control over total output power of the
light source. The intensity of the light exiting the voa could be adjusted
from $0-100\%$ power by rotating a half-wave plate inside the device. To
illuminate the sensor with a uniform field, the sled beam exiting the voa was
directed into a $12\,\mathrm{in}.$ diameter integrating sphere and the sensor
was placed at the output port of the sphere where uniformity is highest. The
sensor was operated at its full bit-depth of $14$-bits to minimize
quantization error and image data was read off the sensor using a single
readout register at its slowest setting ($40\,\mathrm{mhz}$) as not to
introduce additional nonuniformities or smear. For the sake of reducing the
total amount of data captured in the experiment, the source was turned on and
a live-stream of the sensor output was examined to select the most uniformly
illuminated $512\times 512\,\mathrm{px}$ subregion of the sensor array. All
subsequent images captured in the experiment were cropped to this subregion
prior to saving. Since image data needed to be captured under both dark and
illuminated conditions, a motorized mirror was placed next to the path of the
sled beam. Moving the mirror into the beam path effectively redirected the
light away from the integrating sphere and thus provided a dark environment
for the sensor.

####  3.6.3 Design of experiment

_Design of Experiment_ (doe) for pixel-level conversion gain measurement at
its core involves five major steps: (1) choosing values for
$\mathsf{arb}_{0}$, $b$ (bias profile), and $\mathsf{acv}_{0}$, (2) estimating
a global lower bound on the dark noise $\sigma_{\mathrm{d}}g$, (3) computing
$n_{1}^{\text{opt}}$, $n_{2}^{\text{opt}}$, and
$\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ on $\zeta\in[0,1)$, (4)
choosing an illumination level, and (5) identifying pixel groups.

Beginning with the first step it was decided to use $\mathsf{arb}_{0}=0.01$,
$b=0$ (constant bias profile), and $\mathsf{acv}_{0}=0.05$ as these choices
reflect typical values one might use. Then, to obtain a global lower bound
estimate of the dark noise, two $512\times 512\,\mathrm{px}$ dark images
$(D_{1},D_{2})$ were captured. The pixelwise difference $\Delta D=D_{1}-D_{2}$
of these two images was computed and $\sigma_{\mathrm{d}}^{2}$ was estimated
by half of the sample variance

$\hat{\sigma}_{\mathrm{d}}^{2}=\frac{1}{2(512^{2}-1)}\sum_{i=1}^{512}\sum_{j=1}^{512}(\Delta
D_{ij}-\overline{\Delta D})^{2}=39.94\,\mathrm{DN}^{2}.$

Since the value of the conversion gain is unknown, we make the safe assumption
$g\geq 1$ and then obtain a single lower bound estimate of the dark noise for
each pixel by

$\widehat{\sigma_{\mathrm{d}}g}=\hat{\sigma}_{\mathrm{d}}\,(\mathrm{DN})\times
1(e\mathchar 45\relax/\mathrm{DN})=\sqrt{39.94}\,e\mathchar 45\relax.$

For the next step we compute the optimal sample size curves according to
Algorithm 1. Figure 3.5 plots the computed curves for
$\mathrm{N}^{\text{opt}}=n_{1}^{\text{opt}}+n_{2}^{\text{opt}}$ and
$\hat{\mathscr{E}}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ given by

$\hat{\mathscr{E}}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}=\left(1+\frac{1+\mathsf{acv}_{0}^{-2}}{(\widehat{\sigma_{\mathrm{d}}g})^{2}}\frac{\zeta}{(1-\zeta)^{2}}\left(\frac{1}{n_{1}^{\text{opt}}}+\frac{\zeta}{n_{2}^{\text{opt}}}\right)\right)^{-1}$

along with the curves for $n_{1}^{\text{opt}}$ and $n_{2}^{\text{opt}}$ as
functions of $\zeta$.

Figure 3.5: $\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ and
$\mathrm{N}^{\text{opt}}$ versus $\zeta$ (top) with $n_{1}^{\text{opt}}$ and
$n_{2}^{\text{opt}}$ versus $\zeta$ (bottom) for $\mathsf{arb}_{0}=0.01$,
$\mathsf{acv}_{0}=0.05$, $b=0$, and
$\sigma_{\mathrm{d}}g=\sqrt{39.94}e\mathchar 45\relax$.

Due to the sufficiently large value of $\widehat{\sigma_{\mathrm{d}}g}$ we see
that $\hat{\mathscr{E}}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}$ is near unity
for virtually any illumination level so that the requirement
$\mathscr{E}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}\approx 1$ will not restrict
what illumination levels we can choose for the experiment. To select an
appropriate illumination level we first note that this sensor can record
images at $\approx 5\,\mathrm{fps}$ for the chosen readout rate of
$40\,\mathrm{mhz}$. Looking back at Figure 3.5 we observe that the
illumination level corresponding to $\zeta\approx 0.4$ is paired with an
optimal total sample size of $\mathrm{N}^{\text{opt}}\approx 4000$. At a
recording rate of $5\,\mathrm{fps}$ this number of images will take $\approx
13\,\mathrm{min.}$ to capture, which is short enough to avoid any significant
drift in the sensor or source.

Now equipped with a desired value for $\zeta$, the variable optical attenuator
was adjusted until the illumination level at the sensor plane resulted in a
mean pixel output of about $1\%$ of the sensors dynamic range. Using the same
process to estimate the dark noise, two illuminated frames $(I_{1},I_{2})$
were captured and their difference $\Delta I=I_{1}-I_{2}$ was used to obtain
the estimate

$\hat{\sigma}_{\mathrm{p}+\mathrm{d}}^{2}=\frac{1}{2(512^{2}-1)}\sum_{i=1}^{512}\sum_{j=1}^{512}(\Delta
I_{ij}-\overline{\Delta I})^{2}=112.89\,\mathrm{DN}^{2},$

which then gave an estimate $\hat{\zeta}=Z$ equal to

$Z=\frac{512^{2}-3}{512^{2}-1}\frac{\hat{\sigma}_{\mathrm{d}}^{2}}{\hat{\sigma}_{\mathrm{p}+\mathrm{d}}^{2}}=0.354.$

This illumination was sufficiently close to the target $\zeta=0.4$ and
corresponded to
$\hat{\mathscr{E}}_{\mathscr{T}_{\nu^{\smash{\dagger}}}}=0.9928$ and
$\mathrm{N}^{\text{opt}}=3237$ images, which needs only $\approx
11\,\mathrm{min}.$ to capture.

The last step in the doe process is to determine appropriate pixel groupings
for the data collection algorithm. As discussed in Section 3.6.1, identifying
these pixel groups is largely aided by knowledge of the sensor architecture
under test. In this example, we know we are working with an interline transfer
ccd architecture [20]. These devices work by transferring the packets of
charge collected by each pixel into columnwise vertical shift registers that
subsequently facilitate the transfer of the charge packets off the sensor to
downstream readout circuity. Since variations in the columnwise circuitry
supporting these vertical shift registers is expected, nonuniformities in
interline transfer ccds is typically observed between its columns. To see this
in action, a stack of dark images were captured with the KAI-0407M ccd and the
pixelwise temporal average of the stack was computed as seen in Figure 3.6.
From the figure we see clear columnwise variations in the estimated population
dark mean $\mu_{\mathrm{d}}$, which supports our conclusions about the sensor
nonuniformities. As such, we will group pixels according to their column
number for the data collection procedure.

Figure 3.6: Average of dark image stack generated from the KAI-0407M interline
transfer ccd showing columnwise nonuniformities in mean dark signal.

####  3.6.4 Data collection

The procedure for the data collection portion of this experiment centers
around successively capturing images under dark and illuminated conditions,
using these images to update sample statistics for each pixel, and determining
when the appropriate number of dark and illuminated images have been captured
based on these sample statistics. We will initialize the procedure by
capturing some predetermined number of dark images (i.e. a $Y$-image stack)
and illuminated images (i.e. a $X$-image stack) and then compute four _master
frames_ , denoted $\bar{Y}$, $\bar{X}$, $\hat{Y}$, and $\hat{X}$, from these
stacks. Using this nomenclature, the $\bar{Y}$ and $\bar{X}$ master frames
denote arrays computed by averaging the $Y$\- and $X$-image stacks in the
temporal dimension while the $\hat{Y}$ and $\hat{X}$ master frames denote
arrays computed by evaluating the variance of the $Y$\- and $X$-image stacks
in the temporal dimension. As an example, Figure 3.7 depicts the process of
computing a master $\hat{X}$-frame from a stack of $n$ $X$-images.

Figure 3.7: Master $\hat{X}$-frame computed from a stack of $n$ $X$-images.

With master frames computed from the two initial image stacks, we will then
capture individual $Y$\- and $X$-images and use these to continuously update
the master frames until a specified criteria is met to halt data capture. Due
to the large number of images that will be captured in this experiment,
updating the master frames by recomputing the sample statistics from the full
image stacks is too computationally expensive to work in practice. As such we
will use Welford’s algorithm to update the master frames when new image data
becomes available. Algorithm 2 presents the procedure UpdateStats(), which is
an adaptation of Welford’s algorithm to two-dimensional arrays [26]. This
procedure accepts as arguments a new two-dimensional array of data, newdata,
the current master sum of squares array, $M_{2,n}$, the current master average
frame $\bar{x}_{n}$, and the number of samples, $n$, the current master frames
are computed from. If no master frames have been initialized we may use the
syntax UpdateStats(newdata, $0$) to create them. For all subsequent $n\geq 1$,
the syntax UpdateStats(newdata, $n$, $M_{2,n}$, $\bar{x}_{n}$) takes in the
current master arrays and updates them with newdata according to Welford’s
algorithm. We note that all arithmetic operations in this procedure are
performed elementwise, where, in particular, the operator $\odot$ denotes the
Hadamard product (elementwise multiplication of matrices).

Algorithm 2 Welford’s online algorithm for two-dimensional arrays.

1:procedure UpdateStats(newdata, $n$, $M_{2,n}$, $\bar{x}_{n}$)

2: $n=n+1$;

3: if $n=1$ then

4: $\bar{x}_{n}=\text{{newdata}}$;

5: $M_{2,n}=\text{zeros(size({newdata}))}$;

6: $s_{n}^{2}=\text{zeros(size({newdata}))}$;

7: else

8: $\bar{x}_{n}=\bar{x}_{n-1}+\frac{1}{n}(\text{{newdata}}-\bar{x}_{n-1})$;

9: