# On the optimal measurement of conversion gain in the presence of dark noise

Aaron Hendrickson  U.S. Navy, NAWCAD DAiTA Group, Atlantic Ranges & Targets,
Optical Systems, Electro-Optical Tracking Systems, 23013 Cedar Point Road,
Bldg. 2118, Patuxent River MD, 20670  Corresponding author<EMAIL_ADDRESS>David P. Haefner  U.S. Army Combat Capabilities Development Command (DEVCOM),
C5ISR Center, Research & Technology Integration (RTI), 10221 Burbeck Road,
Fort Belvoir, VA 22060  Bradley L. Preece  U.S. Army Combat Capabilities
Development Command (DEVCOM), C5ISR Center, Research & Technology Integration
(RTI), 10221 Burbeck Road, Fort Belvoir, VA 22060

###### Abstract

Working from a model of Gaussian pixel noise, we present and unify over
twenty-five years of developments in the statistical analysis of the photon
transfer conversion gain measurement. We then study a two-sample estimator of
the conversion gain that accounts for the general case of non-negligible dark
noise. The moments of this estimator are ill-defined (their integral
representations diverge) and so we propose a method for assigning
pseudomoments, which are shown to agree with actual sample moments under mild
conditions. A definition of optimal sample size pairs for this two-sample
estimator is proposed and used to find approximate optimal sample size pairs
that allow experimenters to achieve a predetermined measurement uncertainty
with as little data as possible. The conditions under which these
approximations hold are also discussed. Design and control of experiment
procedures are developed and used to optimally estimate a per-pixel conversion
gain map of a real image sensor. Experimental results show excellent agreement
with theoretical predictions and are backed up with Monte Carlo simulation.
The per-pixel conversion gain estimates are then applied in a demonstration of
per-pixel read noise estimation of the same image sensor. The results of this
work open the door to a comprehensive pixel-level adaptation of the photon
transfer method.

††journal: josaa

##  1 Introduction

Photon Transfer (pt) is a methodology developed in the 1970s to aid in the
design, characterization, and optimization of solid state image sensors [1].
Since its inception, pt has evolved to become the standard approach to image
sensor characterization for manufacturers and consumers alike, culminating in
its use as the basis for the European Machine Vision Association (emva) 1288
standard in 2005 [2, 3]. To fully characterize the performance of an image
sensor with pt, many performance parameters are measured including, but not
limited to, conversion gain, read noise, dynamic range, and quantum
efficiency.

Of all performance parameters prescribed by the pt method, the conversion
gain, $g\,(e\text{-}/\mathrm{DN})$, is the most critical as it is the unit
conversion constant needed to convert sensor measurements from device specific
units of Digital Numbers (dn) into physical units of electrons $(e\text{-})$.
Since units of $\mathrm{DN}$ are device specific, it is only after multiplying
by $g$ that $\mathrm{DN}$ measurements represent a physical quantity that can
be compared between different devices. For this reason, many of the
performance parameters measurable by the pt method, e.g. read noise, dynamic
range, and quantum efficiency, at some point require multiplying quantities in
units of $\mathrm{DN}$ by $g$.

Naturally, $g$ must be estimated through measurement; thus, the precision and
accuracy of its measurement fundamentally limits the precision and accuracy of
all measured pt parameters converted to units of electrons by $g$. To see why
this is, suppose $G$ is an estimator of $g$ and $T$ is an estimator of some
parameter $\tau$ with units of $\mathrm{DN}$. Then $\mathscr{T}=T\times G$ is
an estimator of $\tau$ in units of electrons and we have for the absolute
coefficient of variation
$\mathsf{ACV}\mathscr{T}=\sqrt{\mathsf{Var}\mathscr{T}}/|\mathsf{E}\mathscr{T}|$:

$\mathsf{ACV}^{2}\mathscr{T}=\mathsf{ACV}^{2}T+(\mathsf{ACV}^{2}T)(\mathsf{ACV}^{2}G)+\mathsf{ACV}^{2}G.$
(1)

It follows that $\mathsf{ACV}\mathscr{T}\geq\mathsf{ACV}G$ showing that the
relative uncertainty in the estimator $G$ represents a lower bound for the
relative uncertainty of all pt measurements converted to units of electrons.

Given the central role conversion gain plays in the pt method, much work has
been conducted into investigating various estimators for $g$, their
statistical properties, and procedures for performing $g$-estimation [4, 2, 3,
5, 6, 7, 1, 8, 9, 10, 11, 12, 13, 14]. From this body of research, only a
small subset have studied how sample size(s) relate to estimator uncertainty
for the purpose of properly designing and controlling conversion gain
measurement experiments [4, 1, 10, 12]. Perhaps of most influence to this
correspondence is the work of Beecken & Fossum (1996) who studied a shot-
noise-limited estimator under a model of Gaussian noise [4]. When the shot-
noise-limited assumption is made under the Gaussian model, statistical
analysis of the resulting estimator becomes much more tractable and parameters
such as the optimal (minimal) sample size needed to control the estimator’s
uncertainty can be derived. Of course, results derived from this assumption
only hold in the limiting case of shot-noise-limited estimation and is
generally not applicable for measuring $g$ in the presence of dark noise which
is the more general situation. Such a situation can occur when dealing with
sensors of lower-quality or when attempting to measure conversion gain under
low-illumination conditions.

An important consideration in how conversion gain is measured depends on if a
particular sensor under test exhibits uniform characteristics, whereby, each
pixel in the sensor array is assumed to exhibit identical values of conversion
gain, read noise, etc. When a sensor conforms to the uniform assumption, only
a single, global, value of each pt parameter needs to be reported to describe
the performance of the device. Since most sensors are comprised of many
pixels, only a few frames of data are needed to obtain large sample estimates
of these parameters. For example, Janesick showed that when working in the
shot-noise-limited regime, one can measure the conversion gain to a relative
uncertainty of $1\%$ by sampling $\mathcal{O}(10^{4})$ (approximately
$20\,000$) pixels from just a few frames of image data [1]. When the
assumption of sensor uniformity is not held, it is necessary to measure each
pt performance parameter on a per-pixel basis, requiring many frames of data.
Returning back to the example given by Janesick, the $\mathcal{O}(10^{4})$
pixels needed to obtain a global estimate of $g$ within $1\%$ uncertainty now
turns into $\mathcal{O}(10^{4})$ full-resolution image frames needed to
measure the conversion gain of each pixel to the same uncertainty.

The requirement of large datasets for per-pixel conversion gain estimation is
further compounded in the case of Complementary Metal-Oxide Semiconductor
(cmos) Active Pixel Sensors (aps), which not only exhibit per-pixel
nonuniformities, but also are generally nonlinear devices. To circumvent the
problem of nonlinearity, Janesick _et al_. proposed an extension of pt which
requires measuring $g$ at the low-illumination end of it’s dynamic range,
where the sensor exhibits linear characteristics (see Section $7.3$ of [7] and
[8]). As sensor dark noise dominates signal noise at low-illumination, the
number of samples needed to measure $g$ again dramatically increases. For
example, in this correspondence we show that if photon induced signal noise is
approximately equal to sensor dark noise, then to measure $g$, per-pixel, to
$1\%$ uncertainty one needs no less than $\mathcal{O}(10^{5})$ (approximately
$180\,000$) total image frames of data. The dramatic increase in the amount of
data needed for per-pixel estimation, especially when measuring $g$ outside
the shot-noise-limited regime, drives the need to find a means for estimating
$g$ with as little data as possible; this is the main goal of this
correspondence.

From an experimental perspective, capturing large datasets is not just time
consuming but also introduces sensitivity to drift in the sensor and/or light
source, which unchecked, will corrupt per-pixel estimates of $g$. Through
utilizing optimal sampling, one is able to measure $g$ with as few frames as
possible, reducing the time to capture data and mitigating drift. As such, our
goal here is to develop a general method for optimally estimating $g$ that
holds in the shot-noise-limit as well as the more general case where dark
noise is non-negligible.

We will organize this paper by first introducing a Gaussian model of sensor
noise along with the basic assumptions under which the model holds (Section
2). We follow with a condensed review of theory and statistical analysis for
estimators of the conversion gain in the shot-noise-limited case (Section 3)
and general case which accounts for the presence of dark noise (Section 4).
Section 5 then studies optimal sample size pairs for the general conversion
gain estimator, which is an estimator based on two independent samples.
Because analytical expressions for the exact optimal sample size pairs cannot
be derived, this section derives approximations and determines the conditions
under which these approximations are useful. From these analytical expressions
for the approximate optimal sample size pairs, Section 6 describes a method
for design and control of experiment for per-pixel conversion gain estimation
on a real image sensor. Lastly, Section 7 presents an application to per-pixel
conversion gain estimation by performing per-pixel estimation of read noise on
the same sensor.

##  2 Sensor noise model

Consider a sensor observing a constant irradiance, monochromatic light source.
The expected number of interacting photons, $\mu_{\gamma}$, received by each
pixel per fixed integration time is given by

$\mu_{\gamma}=\frac{AEt_{\mathrm{exp}}}{h\nu}QE_{\mathrm{int}},$ (2)

where $A\,(\mathrm{m}^{2})$ is the pixel area,
$E\,(\mathrm{W}/\mathrm{m}^{2})$ is the irradiance,
$t_{\mathrm{exp}}\,(\mathrm{s})$ is the integration time, $h\nu\,(\mathrm{J})$
is the quantization constant for photons of frequency $\nu\,(\mathrm{Hz})$,
and $QE_{\mathrm{int}}\,(-)$ is the interacting quantum efficiency
representing the probability an incident photon is detected [2].

The actual number of observed interacting photons for any given integration
time will vary randomly and is accurately modeled by the Poisson distribution
[1]. Assuming the transfer of interacting photons to photoelectrons is one-to-
one, variations in the number of photoelectrons generated in each pixel per
integration time can therefore be modeled as a Poisson random variable
$\mathscr{P}\sim\mathcal{P}(\mu_{e\text{-}})$. As the irradiance–which we
shall refer to as the illumination level–increases, $\mu_{e\text{-}}\to\infty$
leading to the asymptotic result

$\frac{\mathscr{P}-\mu_{e\text{-}}}{\sqrt{\mu_{e\text{-}}}}\overset{d}{\to}\mathcal{N}(0,1);$
(3)

however, for even relatively small values of $\mu_{e\text{-}}$, say
$\mu_{e\text{-}}>30$, the error in this approximation is generally acceptable
for applied purposes and so we can assume
$\mathscr{P}\sim\mathcal{N}(\mu_{e\text{-}},\mu_{e\text{-}})$.

Additionally, dark noise present in the sensor is comprised of dark current
shot noise and read noise with the read noise further comprised of several
other noise sources, e.g. source follower noise, reset noise, etc., so by the
central limit theorem we have the reasonable model
$\mathscr{D}\sim\mathcal{N}(\mu_{\mathscr{D}},\sigma_{\mathscr{D}}^{2})$,
where $\mu_{\mathscr{D}}$ and $\sigma_{\mathscr{D}}$ represent the sensor bias
and dark noise in units of electrons, respectively. We will further assume
$\mathscr{P}$ and $\mathscr{D}$ are independent. In considering the combined
photon induced and dark signal $\mathscr{P}+\mathscr{D}$, one might be tempted
to assume so long as $\mu_{e\text{-}}$ is large, that the approximate
normality of $\mathscr{P}$ implies $\mathscr{P}+\mathscr{D}$ must also be
approximately normal; however, this is not always the case. Once again using
$\mathscr{P}\sim\mathcal{P}(\mu_{e\text{-}})$ the density of
$\mathscr{P}+\mathscr{D}$ can be formally expressed by the convolution
integral
$f_{\mathscr{P}+\mathscr{D}}(x)=\int_{\mathbb{N}_{0}}f_{\mathscr{D}}(x-n)\,\mathrm{d}F_{\mathscr{P}}(n)$
leading to the explicit form (see [9] and Section 7.2 of [1]):

$f_{\mathscr{P}+\mathscr{D}}(x)=\sum_{n=0}^{\infty}\phi(x-n;\mu_{\mathscr{D}},\sigma_{\mathscr{D}})\frac{e^{-\mu_{e\text{-}}}\mu_{e\text{-}}^{n}}{n!},$
(4)

where $\phi(\cdot;\mu,\sigma)$ is the normal probability density with mean
$\mu$ and standard deviation $\sigma$. Due to the emerging importance of (4)
in the literature [9, 14] we also note that an alternative expression can be
obtained by writing $f_{\mathscr{P}+\mathscr{D}}$ in terms of the exponential
square series function [15]

$f_{\mathscr{P}+\mathscr{D}}(x)=\phi(x;\mu_{\mathscr{D}},\sigma_{\mathscr{D}})e^{-\mu_{e\text{-}}}E_{sq}(e^{-1/2\sigma_{\mathscr{D}}^{2}},\mu_{e\text{-}}e^{(x-\mu_{\mathscr{D}})/\sigma^{2}},1),$
(5)

where $E_{sq}(q,r,z)\coloneqq\sum_{n=0}^{\infty}q^{n^{2}}r^{n}z^{n}/n!$.
Proposition $5.2$ in Schmidt (2017) along with the identities $e^{iz}=\cos
z+i\sin z$ and $2\cos z=e^{iz}+e^{-iz}$ give

$E_{sq}(e^{-1/2\sigma_{\mathscr{D}}^{2}},\Omega,1)=\\\
2\int_{0}^{\infty}\phi(t;0,1)\exp(\Omega\cos(t/\sigma_{\mathscr{D}}))\cos(\Omega\sin(t/\sigma_{\mathscr{D}}))\,\mathrm{d}t,$
(6)

which subsequently provides a novel integral representation for (4).

Figure 1 plots $f_{\mathscr{P}+\mathscr{D}}$ for
$\mu_{e\text{-}}=30\,e\text{-}$, $\mu_{\mathscr{D}}=0\,e\text{-}$ and
$\sigma_{\mathscr{D}}=0.3,1.0\,e\text{-}$. Upon inspection, we see that
despite $\mu_{e\text{-}}$ being large, so that $\mathscr{P}$ is approximately
normal, only the density for $\mathscr{P}+\mathscr{D}$ corresponding to
$\sigma_{\mathscr{D}}=1.0\,e\text{-}$ can be accurately modeled as normal. As
such, on top of the restriction $\mu_{e\text{-}}>30\,e\text{-}$ we will also
assume $\sigma_{\mathscr{D}}>1.0\,e\text{-}$ so that we may use the model
$\mathscr{P}+\mathscr{D}\sim\mathcal{N}(\mu_{e\text{-}}+\mu_{\mathscr{D}},\mu_{e\text{-}}+\sigma_{\mathscr{D}}^{2})$.
This additional assumption excludes photon counting devices such as Deep-Sub-
Electron-Read Noise (dsern) image sensors. We do note that in the case where a
large number of sample are collected these assumptions can be loosened as the
distributions of the sample statistics, e.g. mean and variance, will agree
with our model even if $\mathscr{D}$ and $\mathscr{P}+\mathscr{D}$ deviate
from normality.

Figure 1: Probability density of $\mathscr{P}+\mathscr{D}$ for
$\mu_{e\text{-}}=30\,e\text{-}$, $\mu_{\mathscr{D}}=0.0\,e\text{-}$, and
$\sigma_{\mathscr{D}}=0.3\,e\text{-}$ (solid) and
$\sigma_{\mathscr{D}}=1.0\,e\text{-}$ (dash-dot).

In addition to being accurate for data from real image sensors, the normal
model

$\displaystyle\mathscr{D}\sim\mathcal{N}(\mu_{\mathscr{D}},\sigma_{\mathscr{D}}^{2})$
(7)
$\displaystyle\mathscr{P}+\mathscr{D}\sim\mathcal{N}(\mu_{e\text{-}}+\mu_{\mathscr{D}},\mu_{e\text{-}}+\sigma_{\mathscr{D}}^{2})$
(8)

is mathematically convenient due to the tractability of normal moments.
However, the continuous signals $\mathscr{D}$ and $\mathscr{P}+\mathscr{D}$
are never directly observed because they are quantized via an analog-to-
digital converter. This quantization step transforms the corresponding
continuous density functions into discrete probability mass functions, which
distorts the shape of the distributions and thus alters the moments. By
further imposing $g\leq\sigma_{\mathscr{D}}$, the effects of quantization are
negligible so that we can reasonably assume a normal model for the quantized
signal as well [5, 1].

##  3 Review of gain estimation theory: shot-noise-limited case

Our first goal is to present the major statistical results of conversion gain
estimation for the special case of a shot-noise-limited response. The results
in this section will serve as a starting point for the general case of gain
estimation in Section 4. While we focus on $g$-estimation for a single pixel
by repeated sampling of the pixel in time, the analysis is equally valid in
the case of spatially sampling an array of identical pixels exposed to a
uniform light source. To help keep things organized, Table 1 lists key symbols
pertaining to shot-noise-limited estimation and their associated formulae.
Note that all these symbols are built up from only three fundamental
quantities: $\mu_{e\text{-}}$, $g$, and $n$.

Table 1: List of symbols and corresponding formulae associated with shot-noise-limited estimation. symbol | formula | symbol | formula
---|---|---|---
$\mu_{P}$ | $\mu_{e\text{-}}/g$ | $\alpha$ | $(n-1)/2$
$\sigma_{P}^{2}$ | $\mu_{e\text{-}}/g^{2}$ | $\beta$ | $\alpha/\sigma_{P}^{2}$
$\mu_{\bar{P}}$ | $\mu_{P}$ | $g$ | $\mu_{P}/\sigma_{P}^{2}$
$\sigma_{\bar{P}}^{2}$ | $\sigma_{P}^{2}/n$ | $-$ | $-$

###  3.1 Estimator derivation

A pixel can be modeled as a transfer function
$\mathcal{T}:e\text{-}\to\mathrm{DN}$ mapping photoelectrons to a digital
number output. In general, each pixel comprising the active sensor array is
assumed to have it’s own unique transfer function. Suppose
$\mathcal{T}(e\text{-})=e\text{-}/g$ and our pixel has zero bias and dark
noise so that the dark signal can be formally represented by the degenerate
variable $\mathscr{D}\sim\delta(0)$. Because the only noise in the pixel
output will come from photon shot noise we say the pixel exhibits a shot-
noise-limited response. If $P=\mathcal{T}(\mathscr{P})$ is the random variable
representing the photon induced output signal in $\mathrm{DN}$, then it’s easy
to see

$\mu_{P}\coloneqq\mathsf{E}P=\mathsf{E}(\mathscr{P}/g)=\mu_{e\text{-}}/g$ (9)

and

$\sigma_{P}^{2}\coloneqq\mathsf{Var}P=\mathsf{Var}(\mathscr{P}/g)=\mu_{e\text{-}}/g^{2}.$
(10)

Combining these two results give us the fundamental photon transfer relation

$g=\mu_{P}/\sigma_{P}^{2}.$ (11)

This fundamental relation implies a natural estimator for $g$. Let
$\\{P_{1},\dots,P_{n}\\}$ be a sample of $n$ i.i.d. observations of our pixel
exposed to some constant level of incident illumination for a fixed, nonzero
integration time. Then we can estimate $g$ with

$G=\bar{P}/\hat{P},$ (12)

where $\bar{P}=\frac{1}{n}\sum_{k=1}^{n}P_{k}$ and
$\hat{P}=\frac{1}{n-1}\sum_{k=1}^{n}(P_{k}-\bar{P})^{2}$ are the sample mean
and sample variance, respectively. Under the normal model
$P_{k}\sim\mathcal{N}(\mu_{P},\sigma_{P}^{2})$, $(\bar{P},\hat{P})$ is a
complete sufficient statistic of the unknown parameter
$(\mu_{P},\sigma_{P}^{2})$ so that $G$ also happens to be the Uniformly
Minimum-Variance Unbiased Estimator (umvue) of its expected value [16].

###  3.2 Historical developments

Statistical analysis of the estimator (12) has been previously conducted by
Beecken & Fossum (1996) as well as Janesick (2001) [4, 5]. In both works the
moments of $G$ were approximated with the moments of it’s first-order Taylor
polynomial about
$(\mathsf{E}\bar{P},\mathsf{E}\hat{P})=(\mu_{P},\sigma_{P}^{2})$

$G\approx
g+\frac{g}{\mu_{P}}(\bar{P}-\mu_{P})-\frac{g}{\sigma_{P}^{2}}(\hat{P}-\sigma_{P}^{2}).$
(13)

Using these approximate moments Beecken & Fossum were able to show under the
normal model of sensor noise (c.f. Eq. 20 in [4] using
$S_{g}/g\mapsto\mathsf{ACV}G$, $g\mapsto 1/g$, $\bar{x}\mapsto\mu_{P}$,
$N\mapsto n$, and $\sigma/S\mapsto 1$)

$\mathsf{ACV}^{2}G\approx\frac{2}{n-1}+\frac{1}{n}\frac{1}{\mu_{e\text{-}}}.$
(14)

For clarity we note that the paper by Beecken & Fossum actually studied the
estimator $G^{-1}=\hat{P}/\bar{P}$, which is the conversion gain in units of
$\mathrm{DN}/e\text{-}$; however, applying their noise model and statistical
analysis to $G$ as given by (12) gives the result in (14). Furthermore, as $n$
becomes large we may replace $2/(n-1)$ with $2/n$, which is the same estimate
given by Janesick (c.f. Eq. 2.18 in [5] using
$\sigma_{K}^{2}\mapsto\mathsf{Var}G$, $N_{pix}\mapsto n$,
$S(\mathrm{DN})\mapsto\mu_{P}$, and $K\mapsto g$).

In both works it was noted that at typical illumination levels where $g$ is
measured, (14) is very well approximated by its first term, which happens to
be the first-order Taylor approximation of $\mathsf{ACV}^{2}\hat{P}^{-1}$. In
other words, for sufficiently large illumination, we have the approximate
relation $\mathsf{ACV}G\approx\mathsf{ACV}\hat{P}^{-1}$. Such an approximation
is useful because it tells us that the number of samples needed to measure $g$
to a given uncertainty can be approximated by the number of samples needed to
measure $1/\sigma_{P}^{2}$ to the same uncertainty. This is a key insight that
will appear several more times throughout this correspondence. Using the high
illumination approximation $\mathsf{ACV}G\approx\sqrt{2/n}$ [1, c.f. Eq.
$6.12$,] we subsequently obtain Janesick’s approximation for the optimal
(minimal) number of samples needed to estimate $g$ to a desired relative
uncertainty $\mathsf{acv}_{0}$:

$n^{\mathrm{opt}}\approx\frac{2}{\mathsf{acv}_{0}^{2}}.$ (15)

###  3.3 Further developments

Statistical analysis of the shot-noise-limited estimator in (12) yields
tractable results for the density function, moments, and optimal sample size
without the need to invoke approximate methods. For a normal model, we have
for the distributions of the photon induced signal and it’s sample statistics
: $P_{k}\sim\mathcal{N}(\mu_{P},\sigma_{P}^{2})$,
$\bar{P}\sim\mathcal{N}(\mu_{\bar{P}},\sigma_{\bar{P}}^{2})$, and
$\hat{P}\sim\mathcal{G}(\alpha,\beta)$, with the latter being a gamma variable
parameterized in terms of shape $\alpha$ and rate $\beta$ (see Table 1 for
parameter formulae). Since the $P_{k}$ are normal, $\bar{P}$ and $\hat{P}$ are
also independent. For the density function we use change of variables to write
$f_{G}(g)=\int_{0}^{\infty}t\phi(gt;\mu_{\bar{P}},\sigma_{\bar{P}})f_{\hat{P}}(t)\,\mathrm{d}t$,
which after substituting $u=|g|t/\sigma$ gives

$f_{G}(g)=\frac{\alpha}{|g|}\left(\frac{\beta\sigma_{\bar{P}}}{|g|}\right)^{\alpha}\frac{e^{z^{2}(g)/4-\mu_{\bar{P}}^{2}/(2\sigma_{\bar{P}}^{2})}}{\sqrt{2\pi}}D_{-\alpha-1}(z(g)).$
(16)

Here,
$z(g)=\beta\sigma/|g|-\mu_{\bar{P}}\operatorname{sign}(g)/\sigma_{\bar{P}}$
and
$D_{\nu}(z)\coloneqq\frac{e^{-z^{2}/4}}{\Gamma(-\nu)}\int_{0}^{\infty}t^{-\nu-1}e^{-t^{2}/2-zt}\,\mathrm{d}t$,
which is the parabolic cylinder function.

As for the moments of $G$ we have by the independence of $\bar{P}$ and
$\hat{P}$

$\mathsf{E}G^{k}=(\mathsf{E}\bar{P}^{k})\mathsf{E}\hat{P}^{-k}.$ (17)

The moments of $\bar{P}$ are easily found by comparing the generating function
$e^{2zt-t^{2}}=\sum_{n=0}^{\infty}H_{n}(z)t^{n}/n!$ to the moment generating
function of $\bar{P}$ yielding

$\mathsf{E}\bar{P}^{k}=(i\sigma_{\bar{P}}/\sqrt{2})^{k}H_{k}\left(-i\frac{\mu_{\bar{P}}}{\sqrt{2}\sigma_{\bar{P}}}\right),$
(18)

where $H_{k}(z)\coloneqq(2z-\partial_{z})^{k}\cdot 1$ denotes the $k$th degree
Hermite polynomial and $i$ the imaginary unit. Likewise, we have for the
moments of $\hat{P}^{-1}$

$\mathsf{E}\hat{P}^{-k}=\beta^{k}(\alpha)_{-k}$ (19)

with $(s)_{n}\coloneqq\Gamma(s+n)/\Gamma(s)$ denoting the Pochhammer symbol.

From here, the approximate results of Beecken, Fossum, and Janesick can be
derived rigorously from an exact expression for $\mathsf{ACV}G$. The following
lemma will aid us in this goal and also be used extensively throughout the
rest of this work. All proofs can be found in Section 9.

######  Lemma 1.

Let $T=XY$. If $X$ and $Y$ are independent then

$\mathsf{ACV}^{2}T=\mathsf{ACV}^{2}X+(\mathsf{ACV}^{2}X)(\mathsf{ACV}^{2}Y)+\mathsf{ACV}^{2}Y,$
(20)

with $\mathsf{ACV}T\coloneqq\sqrt{\mathsf{Var}T}/|\mathsf{E}T|$.

With the help of Lemma 1 and the moment expressions given above we deduce the
exact expression (c.f. (14))

$\mathsf{ACV}^{2}G=\frac{2}{n-5}+\frac{2}{n(n-5)}\frac{1}{\mu_{e\text{-}}}+\frac{1}{n}\frac{1}{\mu_{e\text{-}}}.$
(21)

Setting $\mathsf{ACV}^{2}G=\mathsf{acv}_{0}^{2}$ yields a quadratic equation
in $n$, which upon solving gives the optimal sample size needed to measure $g$
to a desired relative uncertainty $\mathsf{acv}_{0}$:

$n^{\mathrm{opt}}=\frac{2+5\mathsf{acv}_{0}^{2}+\frac{1}{\mu_{e\text{-}}}+\left((2+5\mathsf{acv}_{0}^{2}+\frac{1}{\mu_{e\text{-}}})^{2}-12\frac{\mathsf{acv}_{0}^{2}}{\mu_{e\text{-}}}\right)^{1/2}}{2\mathsf{acv}_{0}^{2}}.$
(22)

Unfortunately, this exact expression is not of great use in practice because
it depends on the unknown quantity $\mu_{e\text{-}}$, which also cannot be
directly measured without _a priori_ knowledge of $g$. To obtain an
approximation that is independent of $\mu_{e\text{-}}$, we first consider the
following result showing that $\mathsf{ACV}G$ is dominated by
$\mathsf{ACV}\hat{P}^{-1}$ at high-illumination.

######  Theorem 1.

Let $G$ be as given in (12). As illumination increases,
$\mu_{e\text{-}}\to\infty$ and

$\mathsf{ACV}G=\mathsf{ACV}\hat{P}^{-1}\left(1+\frac{n-3}{4n}\frac{1}{\mu_{e\text{-}}}+\mathcal{O}(\mu_{e\text{-}}^{-2})\right),$
(23)

with $\mathsf{ACV}\hat{P}^{-1}=\sqrt{2/(n-5)}$.

Theorem 1 confirms the observations of [4, 5] in that

$\mathsf{ACV}G\sim\mathsf{ACV}\hat{P}^{-1}$ (24)

at high-illumination. Because of this finding, the optimal sample size for $G$
can be approximated by the optimal sample size for $\hat{P}^{-1}$ when
restricted to high-illumination conditions. Setting
$\mathsf{ACV}^{2}\hat{P}^{-1}=\mathsf{acv}_{0}^{2}$ yields a linear equation
in $n$, which upon solving for $n$ subsequently gives us the high-
illumination, asymptotic approximation for the optimal sample size of $G$

$n^{\mathrm{opt}}\sim\frac{2}{\mathsf{acv}_{0}^{2}}+5,\quad\mu_{e\text{-}}\to\infty.$
(25)

For example, choosing a desired relative uncertainty of $1\%$ we have
$\mathsf{acv}_{0}=0.01$ and $n^{\mathrm{opt}}\approx 20\,005$, which agrees
with Janesick’s approximation of $n^{\mathrm{opt}}\approx 20\,000$ as given by
(15). Unfortunately, all of the results in this section break down when sensor
dark noise is non-negligible and thus have limited applications. We are now
ready to move onto the more general case.

##  4 Review of gain estimation theory: general case

As was the case in the previous section, general gain estimation requires the
use of many symbols that can be combined and manipulated. To stay organized,
Table 2 lists many of the key symbols used along with their corresponding
formulae. Note that all of these symbols are constructed from six fundamental
quantities: $\mu_{e\text{-}}$, $\mu_{\mathscr{D}}$,
$\sigma_{\mathscr{D}}^{2}$, $g$, $n_{1}$, and $n_{2}$.

Table 2: List of symbols and corresponding formulae associated with general estimation. symbol | formula | symbol | formula
---|---|---|---
$\mu_{D}$ | $\mu_{\mathscr{D}}/g$ | $\mu_{P}^{2}$ | $\mu_{e\text{-}}/g$
$\sigma_{D}^{2}$ | $\sigma_{\mathscr{D}}^{2}/g^{2}$ | $\sigma_{P}^{2}$ | $\mu_{e\text{-}}/g^{2}$
$\mu_{P+D}$ | $\mu_{P}+\mu_{D}$ | $\alpha_{1}$ | $(n_{1}-1)/2$
$\sigma_{P+D}^{2}$ | $\sigma_{P}^{2}+\sigma_{D}^{2}$ | $\alpha_{2}$ | $(n_{2}-1)/2$
$\mu_{\bar{P}}$ | $\mu_{P}$ | $\beta_{1}$ | $\alpha_{1}/\sigma_{P+D}^{2}$
$\sigma_{\bar{P}}^{2}$ | $\sigma_{P+D}^{2}/n_{1}+\sigma_{D}^{2}/n_{2}$ | $\beta_{2}$ | $\alpha_{2}/\sigma_{D}^{2}$
$\mu_{\hat{P}}$ | $\sigma_{P}^{2}$ | $g$ | $\frac{\mu_{P+D}-\mu_{D}}{\sigma_{P+D}^{2}-\sigma_{D}^{2}}$
$\sigma_{\hat{P}}^{2}$ | $\alpha_{1}/\beta_{1}^{2}+\alpha_{2}/\beta_{2}^{2}$ | $-$ | $-$

###  4.1 Estimator derivation

We now consider the more general case where the pixel exhibits both a bias and
non-negligible dark noise:
$\mathscr{D}\sim\mathcal{N}(\mu_{\mathscr{D}},\sigma_{\mathscr{D}}^{2})$ with
$\mu_{\mathscr{D}},\sigma_{\mathscr{D}}\neq 0$. We again note that in general,
each pixel in the active sensor array will exhibit unique values of
$\mu_{\mathscr{D}}$ and $\sigma_{\mathscr{D}}$. By linearity of the transfer
function, the digital output of a pixel in the absence of illumination for
some fixed, nonzero integration time is $D=\mathcal{T}(\mathscr{D})$ with

$\mu_{D}\coloneqq\mathsf{E}D=\mu_{\mathscr{D}}/g$ (26)

and

$\sigma_{D}^{2}\coloneqq\mathsf{Var}D=\sigma_{\mathscr{D}}^{2}/g^{2}.$ (27)

Likewise, the digital output of the same pixel exposed to incident
illumination for the same fixed, nonzero integration time gives the combined
dark and photon induced signal
$P+D=\mathcal{T}(\mathscr{P}+\mathscr{D})=\mathcal{T}(\mathscr{P})+\mathcal{T}(\mathscr{D})$,
which by the assumed independence of $\mathscr{P}$ and $\mathscr{D}$ further
gives

$\mu_{P+D}\coloneqq\mathsf{E}(P+D)=\mu_{P}+\mu_{D}$ (28)

and

$\sigma_{P+D}^{2}\coloneqq\mathsf{Var}(P+D)=\sigma_{P}^{2}+\sigma_{D}^{2}.$
(29)

Noting that $\mu_{P}=\mu_{P+D}-\mu_{D}$ and
$\sigma_{P}^{2}=\sigma_{P+D}^{2}-\sigma_{D}^{2}$ leads to the modified gain
relation

$g=(\mu_{P+D}-\mu_{D})/(\sigma_{P+D}^{2}-\sigma_{D}^{2}).$ (30)

As before, this relation suggests a natural estimator for $g$. Let
$\\{X_{1},\dots,X_{n_{1}}\\}$ be a sample of $n_{1}$ i.i.d. observations of a
pixel exposed to a constant level of illumination for a fixed, nonzero
integration time and $\\{Y_{1},\dots,Y_{n_{2}}\\}$ be an independent sample of
$n_{2}$ i.i.d. observations of the same pixel in the absence of illumination
for the same fixed, nonzero integration time. Then we can estimate $g$ for our
pixel with

$G=(\bar{X}-\bar{Y})/(\hat{X}-\hat{Y}),$ (31)

with bar $(\bar{\cdot})$ and hat $(\hat{\cdot})$ accents again denoting sample
means and sample variances, respectively. In this context we see

$\bar{P}\coloneqq\bar{X}-\bar{Y}$ (32)

estimates the photon induced signal $\mu_{P}$, while

$\hat{P}\coloneqq\hat{X}-\hat{Y}$ (33)

estimates the photon induced signal variance $\sigma_{P}^{2}$, in a similar
manner to that of the estimators $\bar{P}$ and $\hat{P}$ in (12). Note that
now $\bar{P}$ no longer represents a single sample mean but rather the
difference of sample means and likewise for $\hat{P}$.

###  4.2 Historical developments

Janesick (2007) was the first work to include the contributions of dark noise
in conversion gain estimator uncertainty. He did this by studying a variant of
(31) where the numerator is substituted with the random variable $\bar{P}$
($\mathsf{E}\bar{P}=\mu_{P}$, $\mathsf{Var}\bar{P}=\sigma_{P}^{2}/n_{1})$. The
relative uncertainty of the resulting estimator was then approximated via the
moments of its first-order Taylor polynomial [1]. Further imposing large and
equal numbers of samples for both the bright and dark collection,
$n_{1}=n_{2}=n$, Janesick obtained the approximation (c.f. Eq. 6.7-6.8 in [1]
using $N_{\mathrm{PIX}}\mapsto n$, $K\mapsto g$, $S\mapsto\mu_{P}$,
$N\mapsto\sigma_{P+D}$, and $R\mapsto\sigma_{D}$)

$\mathsf{ACV}^{2}G\approx\frac{2}{n}+\frac{1}{n}\frac{1}{\mu_{e\text{-}}}+\frac{4}{n}\frac{\sigma_{\mathscr{D}}^{2}}{\mu_{e\text{-}}}+\frac{4}{n}\left(\frac{\sigma_{\mathscr{D}}^{2}}{\mu_{e\text{-}}}\right)^{2}.$
(34)

It was noted that when photon shot noise is large ($\mu_{e\text{-}}\to\infty)$
and dominates sensor dark noise
($\mu_{e\text{-}}\gg\sigma_{\mathscr{D}}^{2}$), this estimate reduces to the
high illumination, shot-noise-limited estimate
$\mathsf{ACV}G\approx\sqrt{2/n}$ discussed in Section 3.

Following Janesick’s work, Hendrickson (2017) was the first attempt to draw
exact statistical conclusions about the full estimator (31) by applying the
normal model $X_{k}\sim\mathcal{N}(\mu_{P+D},\sigma_{P+D}^{2})$ and
$Y_{k}\sim\mathcal{N}(\mu_{D},\sigma_{D}^{2})$ to derive the density of $G$ in
the form of the Centralized Inverse-Fano (cif) distribution

$f_{G}(g)=\int_{-\infty}^{\infty}|t|\phi(gt;\mu_{\bar{P}},\sigma_{\bar{P}})f_{\hat{P}}(t)\,\mathrm{d}t,$
(35)

where

$f_{\hat{P}}(t)=C\times\begin{cases}\frac{e^{\beta_{2}t}}{\Gamma(\alpha_{2})}U\left({1-\alpha_{2}\atop
2-\alpha_{1}-\alpha_{2}};-(\beta_{1}+\beta_{2})t\right)&t<0\\\
\frac{e^{-\beta_{1}t}}{\Gamma(\alpha_{1})}U\left({1-\alpha_{1}\atop
2-\alpha_{1}-\alpha_{2}};(\beta_{1}+\beta_{2})t\right)&t\geq 0,\end{cases}$
(36)

with
$C=\beta_{1}^{\alpha_{1}}\beta_{2}^{\alpha_{2}}(\beta_{1}+\beta_{2})^{1-\alpha_{1}-\alpha_{2}}$
is the gamma-difference distribution [17, 18, 19, 20], and $U(a,b,z)$ is
Kummer’s confluent hypergeometric function of the 2nd-kind [10]. Under this
model it was shown that $G$ has ill-defined moments which is due to the fact
that the tails of the probability density $f_{G}$ decay like those of the
Cauchy density [21]. For this reason, Hendrickson (2019) [11] extended the
notion of statistical moments in the same manner as Peng [22, 23] by deriving
the first moment of $G$ in the sense of the Cauchy principal value
$\mathsf{E}_{\mathcal{P}}G=\lim_{R\to\infty}\int_{-R}^{R}tf_{G}(t)\,\mathrm{d}t$

$\mathsf{E}_{\mathcal{P}}G=g\frac{(\frac{\alpha_{1}}{\beta_{1}}-\frac{\alpha_{2}}{\beta_{2}})\beta_{1}^{\alpha_{1}}\beta_{2}^{\alpha_{2}}(\beta_{1}+\beta_{2})^{1-\alpha_{1}-\alpha_{2}}}{(\alpha_{1}+\alpha_{2}-1)\operatorname{B}(\alpha_{1},\alpha_{2})}\biggl{(}\\\
\psi(\alpha_{1})-\log\beta_{1}+\frac{(\alpha_{1}-1)\beta_{2}}{\alpha_{2}\beta_{1}}{{}_{3}F_{2}}\left({2-\alpha_{1},1,1\atop
1+\alpha_{2},2};-\frac{\beta_{2}}{\beta_{1}}\right)\\\
-\psi(\alpha_{2})+\log\beta_{2}-\frac{(\alpha_{2}-1)\beta_{1}}{\alpha_{1}\beta_{2}}{{}_{3}F_{2}}\left({2-\alpha_{2},1,1\atop
1+\alpha_{1},2};-\frac{\beta_{1}}{\beta_{2}}\right)\biggr{)}.$ (37)

Here, $\log z$ is the natural logarithm, $\psi(z)$ is the digamma function,
$\operatorname{B}(\alpha,\beta)$ is the beta function, and
${{}_{p}F_{q}}(\mathbf{a};\mathbf{b};z)$ is the generalized hypergeometric
function. It was shown that $\mathsf{E}_{\mathcal{P}}G$ agrees with actual
sample means of conversion gain data when $\mathsf{P}(\hat{P}\leq 0)\approx
0$. Additionally, Hendrickson (2021) showed that no unbiased, finite variance
estimator of the modified gain relation (30) exists for all possible
parameters under the normal model of noise despite
$(\bar{Y},\hat{Y},\bar{X},\hat{X})$ constituting a complete sufficient
statistic for the parameter
$(\mu_{D},\sigma_{D}^{2},\mu_{P+D},\sigma_{P+D}^{2})$ [12, Theorem 14,].

###  4.3 Further developments

The nonexistence of $G$’s moments lies in the fact that $\hat{P}$ has positive
and continuous probability density at zero, which manifests as non-integrable
singularities in the integral representations of these moments. To assign some
notion of higher-order moments to $G$ we may take advantage of the
independence of $\bar{P}$ and $\hat{P}$ and the concept of regularization to
define pseudomoments as

$\mathsf{E}_{\mathcal{P}}G^{k}\coloneqq(\mathsf{E}\bar{P}^{k})\mathsf{E}_{\mathcal{P}}\hat{P}^{-k},\quad
k\in\mathbb{N},$ (38)

where $\mathcal{P}$ denotes the principal-value regularization

$\mathsf{E}_{\mathcal{P}}\hat{P}^{-k}\coloneqq\lim_{\epsilon\to
0^{+}}\int_{\mathbb{R}\setminus(-\epsilon,\epsilon)}\frac{f_{\hat{P}}(t)}{t^{k}}\,\mathrm{d}t-h_{k}(\epsilon),$
(39)

with $h_{1}(\epsilon)=0$ and

$h_{k}(\epsilon)=\sum_{\ell=0}^{k-2}\frac{f_{\hat{P}}^{(\ell)}(0)}{\ell!}\left(\frac{1-(-1)^{k-\ell-1}}{(k-\ell-1)\epsilon^{k-\ell-1}}\right)$
(40)

for $k\geq 2$ [24, 25]. Alternative definitions and methods also exist for
evaluating the pseudomoments. For example, we may express them via a moment
generating function as

$\mathsf{E}_{\mathcal{P}}\hat{P}^{-k}=\frac{1}{(k-1)!}\partial_{\omega}^{k-1}\mathcal{H}[f_{\hat{P}}](\omega)\Big{|}_{\omega=0}$
(41)

with

$\mathcal{H}[f_{\hat{P}}](\omega)=\lim_{\epsilon\to
0^{+}}\int_{\mathbb{R}\setminus(\omega-\epsilon,\omega+\epsilon)}\frac{f_{\hat{P}}(t)}{t-\omega}\,\mathrm{d}t,$
(42)

denoting the Hilbert transform of the density $f_{\hat{P}}$ [26]. Regardless
of the method used to evaluate them, the subscript $\mathcal{P}$ is there to
remind the reader that these moments are regularized and do not exist in the
traditional sense because the integrals representing them diverge.

Evaluating the pseudomoments proves to be quite challenging. As pointed out
earlier, the work in [11] was able to obtain an expression for the special
case $k=1$, as given in (37), through the use of complex methods involving
contour integration. For $k\geq 2$, one could implement (39) numerically;
however, when $n_{1}$ and $n_{2}$ become even moderately large the evaluation
of $f_{\hat{P}}$ becomes unstable causing numerical computation to fail. In
such a case we may approximate the distribution of $\hat{P}$ by a normal one
with mean $\mu_{\hat{P}}$ and variance $\sigma_{\hat{P}}^{2}$ so that we may
use the normal approximation given by Quenouille [27]

$\mathcal{H}[f_{\hat{P}}](\omega)\sim\frac{\sqrt{2}}{\sigma_{\hat{P}}}\mathcal{D}\left(\frac{\mu_{\hat{P}}-\omega}{\sqrt{2}\sigma_{\hat{P}}}\right).$
(43)

Here, $\mathcal{D}(z)\coloneqq e^{-z^{2}}\int_{0}^{z}e^{t^{2}}\,\mathrm{d}t$
denotes the Dawson integral. Higher-order derivatives for the Dawson integral
are given by Barakat (1971), which allow us to deduce a closed-form,
asymptotic approximation for the pseudomoments of $\hat{P}^{-1}$ given large
$n_{1}$ and $n_{2}$ [28]:

$\mathsf{E}_{\mathcal{P}}\hat{P}^{-k}\sim\frac{2}{(k-1)!}\frac{1}{(\sqrt{2}\sigma_{\hat{P}})^{k}}(H_{k-1}(z_{\hat{P}})\mathcal{D}(z_{\hat{P}})-P_{k-2}(z_{\hat{P}})),$
(44)

where $z_{\hat{P}}=\mu_{\hat{P}}/(\sqrt{2}\sigma_{\hat{P}})$, and $P_{n}$ is a
polynomial satisfying $P_{n}(t)=2tP_{n-1}(t)-2nP_{n-2}(t)$ with $P_{-1}(t)=0$
and $P_{0}(t)=1$. Combining this result with the moments of $\bar{P}$
subsequently give us the large sample size asymptotic approximation of the
pseudomoments for $G$

$\mathsf{E}_{\mathcal{P}}G^{k}\sim\frac{2}{(k-1)!}\left(\frac{i\sigma_{\bar{P}}}{2\sigma_{\hat{P}}}\right)^{k}H_{k}(-iz_{\bar{P}})\\\
\cdots\times(H_{k-1}(z_{\hat{P}})\mathcal{D}(z_{\hat{P}})-P_{k-2}(z_{\hat{P}})),$
(45)

where $z_{\bar{P}}=\mu_{\bar{P}}/(\sqrt{2}\sigma_{\bar{P}})$. For example, in
the case $k=1$ we obtain the large sample size (large $\alpha$) asymptotic
approximation of (37)

$\mathsf{E}_{\mathcal{P}}G\sim
g\frac{\sqrt{2}(\frac{\alpha_{1}}{\beta_{1}}-\frac{\alpha_{2}}{\beta_{2}})}{\sqrt{\frac{\alpha_{1}}{\beta_{1}^{2}}+\frac{\alpha_{2}}{\beta_{2}^{2}}}}\mathcal{D}\left(\frac{\frac{\alpha_{1}}{\beta_{1}}-\frac{\alpha_{2}}{\beta_{2}}}{\sqrt{2}\sqrt{\frac{\alpha_{1}}{\beta_{1}^{2}}+\frac{\alpha_{2}}{\beta_{2}^{2}}}}\right).$
(46)

As a verification of the accuracy of these approximations, using the
parameters $\mu_{P}=9\,\mathrm{DN}$, $\sigma_{P+D}^{2}=10\,\mathrm{DN}^{2}$,
$\sigma_{D}^{2}=1\,\mathrm{DN}^{2}$ (so that $g=1$), $n_{1}=101$, and
$n_{2}=51$ we calculated $\mathsf{E}_{\mathcal{P}}G$ using the exact
expression (37) as well as the normal approximation (46) yielding
$\mathsf{E}_{\mathcal{P}}G=1.02604$ and $\mathsf{E}_{\mathcal{P}}G\approx
1.02738$, respectively. This resulted in only a $0.13\%$ approximation error
showing that the normal approximation in (45) will approximate
$\mathsf{E}_{\mathcal{P}}G$ as well as the higher-order pseudomoments for the
chosen sample sizes. We typically deal with much larger sample sizes in pt
conversion gain estimation and therefore expect (46) to be a good
approximation to the exact pseudomoments of $G$ in most scenarios.

Through the use of pseudomoments we can subsequently derive other quantities
of interest for $G$ like the (pseudo) absolute coefficient of variation

$\mathsf{ACV}_{\mathcal{P}}^{2}G=\mathsf{ACV}_{\mathcal{P}}^{2}\hat{P}^{-1}+(\mathsf{ACV}_{\mathcal{P}}^{2}\hat{P}^{-1})(\mathsf{ACV}^{2}\bar{P})+\mathsf{ACV}^{2}\bar{P}$
(47)

where

$\mathsf{ACV}^{2}_{\mathcal{P}}\hat{P}^{-1}=\frac{\mathsf{E}_{\mathcal{P}}\hat{P}^{-2}-(\mathsf{E}_{\mathcal{P}}\hat{P}^{-1})^{2}}{(\mathsf{E}_{\mathcal{P}}\hat{P}^{-1})^{2}},$
(48)

and $\mathsf{ACV}^{2}\bar{P}=\mathsf{Var}\bar{P}/(\mathsf{E}\bar{P})^{2}$ is
defined in the traditional sense. Likewise, we have the (pseudo) absolute
relative bias

$\mathsf{ARB}_{\mathcal{P}}G=\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}$ (49)

with

$\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}=\left|\frac{\mathsf{E}_{\mathcal{P}}\hat{P}^{-1}-(\mathsf{E}\hat{P})^{-1}}{(\mathsf{E}\hat{P})^{-1}}\right|.$
(50)

What remains is to address why these pseudomoments are useful for describing
moments of actual data. The motivation for introducing pseudomoments was to
assign analytical expressions to the moments $\mathsf{E}\hat{P}^{-k}$, which
diverge on sets of the form $|\hat{P}|<\epsilon$. This is clear from writing

$\mathsf{E}\hat{P}^{-k}=\int_{\mathbb{R}\setminus(-\epsilon,\epsilon)}\frac{f_{\hat{P}}(t)}{t^{k}}\,\mathrm{d}t+\int_{-\epsilon}^{\epsilon}\frac{f_{\hat{P}}(t)}{t^{k}}\,\mathrm{d}t,$
(51)

where the first integral in this decomposition always converges, while the
second integral always diverges for any choice of $k\in\mathbb{N}$ and
$\epsilon>0$. As such, the principal value regularization provides a means of
discarding the divergent terms arising out of the second integral to provide a
finite expression for $\mathsf{E}\hat{P}^{-k}$. In practice we don’t ever
observe $\hat{P}$ out in the extreme tails of its assumed distribution due to
the small probabilities of such events, and even more so, inaccuracies in our
assumed noise model. In particular, recall that we assumed a normal
distribution for the digital signal $X_{k}$ and $Y_{k}$ but in practice these
quantities can only take on values between $0$ and $2^{N_{\mathrm{bits}}}-1$
with $N_{\mathrm{bits}}$ denoting the bit-depth of the analog-to-digital
converter. Because the normal model assigns positive density outside this
interval, we see the model inherently overestimates the tails of the actual
data. Therefore, if $\mathsf{ACV}\hat{P}$ is small we won’t observe
$|\hat{P}|<\epsilon$ and the sample moments will agree with the pseudomoments.
Again using the parameters following (46), $10^{9}$ pseudorandom observations
of $\bar{P}$ and $\hat{P}$ were generated, which were then used to compute a
sample of $10^{9}$ observations of $G$. Due to a sufficiently small value of
$\mathsf{ACV}\hat{P}$, all observed values of $\hat{P}$ were strictly
positive, so we should expect the sample statistics to agree with the
theoretical pseudomoments. Computing the sample mean yielded
$\bar{G}=1.02604\dots$, which agreed with the exact value of
$\mathsf{E}_{\mathcal{P}}G$ to six significant digits. Likewise, we would
expect the higher-order sample moments to agree with their corresponding
higher-order pseudomoments for these parameters.

##  5 Optimal sample size pairs for conversion gain estimation

Now that we have a thorough understanding of the statistical characteristics
for the general conversion gain estimator, we may begin to tackle the problem
of optimal measurement. Recall from Section 3 that the optimal sample size for
the one-sample, shot-noise-limited estimator $G$ satisfied
$\mathsf{ACV}G(n^{\mathrm{opt}})=\mathsf{acv}_{0}$ for any choice of
$\mathsf{acv}_{0}>0$. In the present problem, we are now working with an
estimator of two samples and must first define what is meant by optimal sample
sizes in this two-sample case.

######  Definition 1 (Optimal sample size pairs).

Let $T=f(\mathbf{X},\mathbf{Y})$ be a statistic of a sample $\mathbf{X}$ of
size $n_{1}$ and another sample $\mathbf{Y}$ of size $n_{2}$. Furthermore, let
$\mathsf{ACV}T(n_{1},n_{2})$ denote the absolute coefficient of variation for
$T$ as a function of the sample sizes. Then the optimal sample size pairs for
$T$ shall be defined as the ordered pair
$(n_{1}^{\mathrm{opt}},n_{2}^{\mathrm{opt}})$, which satisfies the system of
equations

$\displaystyle\inf_{n_{2}}\,\mathsf{ACV}T(N-n_{2},n_{2})\Big{|}_{N=n_{1}^{\mathrm{opt}}+n_{2}^{\mathrm{opt}},\,n_{2}=n_{2}^{\mathrm{opt}}}$
(52)
$\displaystyle\mathsf{ACV}T(n_{1}^{\mathrm{opt}},n_{2}^{\mathrm{opt}})=\mathsf{acv}_{0}.$
(53)

In this system of equations we see (52) fixes the total number of samples to
$N=n_{1}+n_{2}$ and solves for the $n_{2}$ that minimizes $\mathsf{ACV}T$. In
the special case where $\mathsf{ACV}T(N-n_{2},n_{2})$ is strictly convex in
$n_{2}$, this minimization can be solved via equating the derivative with
zero: $\partial_{n_{2}}\mathsf{ACV}T(N-n_{2},n_{2})=0$. Solving (52) and then
substituting $N\mapsto n_{1}^{\mathrm{opt}}+n_{2}^{\mathrm{opt}}$ and
$n_{2}\mapsto n_{2}^{\mathrm{opt}}$ implicitly defines the optimal sample
sizes as a function of each other, which we will call the _optimality
relation_. For example, in the proof of Lemma 3 we derive the optimality
relation for the estimator $\bar{P}$ in the form
$n_{2}^{\mathrm{opt}}=\frac{\sigma_{D}}{\sigma_{P+D}}n_{1}^{\mathrm{opt}}$,
which expresses the relationship between the optimal sample sizes for
$\bar{P}$. Substituting the optimality relation into (53) then scales the
optimal sample sizes so that they not only satisfy the optimality relation but
also achieve a prescribed final absolute coefficient of variation equal to
$\mathsf{acv}_{0}$. In this way, $n_{1}^{\mathrm{opt}}$ and
$n_{2}^{\mathrm{opt}}$ represent the sample sizes whose sum is the minimal
possible number of total samples needed to force $\mathsf{ACV}T$ equal to
$\mathsf{acv}_{0}$ and therefore serves as a good generalization of optimal
sample size to the two-sample case.

In Section 4 we were unable to derive an exact expression for
$\mathsf{ACV}_{\mathcal{P}}G$ and thus are forced to make some form of
approximation to get a handle on this problem. Even if we did have an exact
expression, substituting it in Definition 1 would almost certainly yield an
intractable system of equations. To overcome these barriers, we will take the
approach of linearizing $\hat{P}^{-1}$ by replacing it with it’s first-order
Taylor polynomial about $\mathsf{E}\hat{P}=\sigma_{P}^{2}$ and instead focus
on

$G_{\delta}=\bar{P}\times\hat{P}^{-1}_{\delta},$ (54)

with

$\hat{P}^{-1}_{\delta}=\frac{1}{\sigma_{P}^{2}}-\frac{\hat{P}-\sigma_{P}^{2}}{\sigma_{P}^{4}}.$
(55)

The advantage we gain from this linearization is that $\hat{P}^{-1}_{\delta}$
and $G_{\delta}$ have simple and well-defined moments that can be used to make
concrete conclusions about their statistical properties including their
optimal sample sizes. Furthermore, so long as $\mathsf{ACV}\hat{P}$ is small
we have $\hat{P}^{-1}_{\delta}\overset{d}{\approx}\hat{P}^{-1}$, which implies
$G_{\delta}\overset{d}{\approx}G$ so that any conclusions we make about the
random variables $\hat{P}^{-1}_{\delta}$ and $G_{\delta}$ apply to
$\hat{P}^{-1}$ and $G$ when $\mathsf{ACV}\hat{P}$ is small. In particular,
note that $\mathsf{ACV}\hat{P}=\mathsf{ACV}\hat{P}^{-1}_{\delta}$ so we can
conclude that the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ should be
good approximations for those of $G$ when $\mathsf{acv}_{0}$ is chosen to be
small and the optimal sample sizes of $\hat{P}^{-1}_{\delta}$ are
approximately equal to those of $G_{\delta}$.

Our angle of attack from here will be as follows. We will first use Definition
1 to derive properties pertaining to the optimal sample sizes of $\bar{P}$
(Subsection 5.1), which will be used later on to derive similar properties of
the optimal sample sizes for $G_{\delta}$. Following analysis of $\bar{P}$
will be analogous study of $\hat{P}^{-1}_{\delta}$ (Subsection 5.2), where we
derive the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ and their
properties. Subsection 5.3 then studies the estimator $G_{\delta}$ with the
goal of understanding when the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ are a good approximation for those of $G_{\delta}$. We
first show that $\mathsf{ACV}G_{\delta}$ is dominated by
$\mathsf{ACV}\hat{P}^{-1}_{\delta}$ at high-illumination (Theorem 3) so that
the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ are asymptotically equal
to those of $G_{\delta}$ at high-illumination. At the other end of the
illumination range, Theorem 4 shows that the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ are asymptotically proportional to those of
$G_{\delta}$ with the constant of proportionality approaching one with
increasing dark noise $\sigma_{\mathscr{D}}$. What this demonstrates is that
the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ should be excellent
approximations to those of $G_{\delta}$ at any illumination level given
sufficiently large dark noise. This observation leads us to construct a
metric, $\mathcal{E}_{\mathrm{opt}}$, which indicates how good the optimal
sample sizes of $\hat{P}^{-1}_{\delta}$ are as approximations for those of
$G_{\delta}$ as a function of illumination level. By studying
$\mathcal{E}_{\mathrm{opt}}$ at low-illumination we are able to put a rule-of-
thumb on how large the dark noise must be to obtain a good approximation. In
particular, results will show that for sensors with moderate dark noise of
$\sigma_{\mathscr{D}}\geq 5\,e\text{-}$, the optimal sample sizes of
$\hat{P}^{-1}_{\delta}$ are excellent approximations for those of $G_{\delta}$
at any illumination level. Theorem 5 then confirms these findings by showing
for small $\mathsf{acv}_{0}$

$\mathsf{ACV}_{\mathcal{P}}G(n_{1}^{\mathrm{opt}},n_{2}^{\mathrm{opt}})=c\cdot\mathsf{acv}_{0}+\mathcal{O}(\mathsf{acv}_{0}^{3}),$
(56)

with $n_{i}^{\mathrm{opt}}$ again representing the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ and $c\to 1$ with increasing illumination and/or
increasing dark noise. Consequently, we show that the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ serve as excellent approximations for those of $G$ at
any illumination level when $\mathsf{acv}_{0}$ is small and
$\sigma_{\mathscr{D}}\geq 5\,e\text{-}$.

Before proceeding, it will also be convenient to introduce the dimensionless
quantity

$\zeta=\sigma_{D}^{2}/\sigma_{P+D}^{2}.$ (57)

We will reparameterize all subsequent analysis in terms of this quantity via
the following lemma.

######  Lemma 2.

Under the assumed model we have $\mu_{P}=\sigma_{D}^{2}\zeta^{-1}(1-\zeta)g$
and $\sigma_{P}^{2}=\sigma_{D}^{2}\zeta^{-1}(1-\zeta)$.

Because $0<\sigma_{D}<\sigma_{P+D}$, it follows that $\zeta\in(0,1)$. As
illumination decreases to zero, we find $\zeta\to 1^{-}$. Likewise, as
illumination increases without bound $\zeta\to 0^{+}$, with $\zeta=0$ denoting
a mathematical definition of the shot-noise-limit. Real sensors always contain
some dark noise and are limited in well capacity; therefore, one can only
achieve the shot-noise-limit in theory and never in real experiments. Lastly,
the reader should take note that $\zeta$ is increasing as illumination
decreases.

###  5.1 Statistical analysis of bar-$P$

We begin our analysis by studying the random variable $\bar{P}$ with the
primary purpose of understanding its optimal sample sizes for later use in
studying $G_{\delta}$. Using Lemma 2, the dark noise relation
$\sigma_{\mathscr{D}}=\sigma_{D}\times g$, and the distributional result
$\bar{P}\sim\mathcal{N}(\mu_{P},\sigma_{P+D}^{2}/n_{1}+\sigma_{D}^{2}/n_{2})$
we are able to write the squared absolute coefficient of variation as

$\mathsf{ACV}^{2}\bar{P}=\frac{1}{\sigma_{\mathscr{D}}^{2}}\frac{\zeta}{(1-\zeta)^{2}}\left(\frac{1}{n_{1}}+\frac{\zeta}{n_{2}}\right).$
(58)

With this expression in hand, we present some key properties of the optimal
sample sizes for $\bar{P}$.

######  Lemma 3.

Let $n_{1}^{\mathrm{opt}}$ and $n_{2}^{\mathrm{opt}}$ denote the optimal
sample sizes for $\bar{P}$. As the illumination decreases, $\zeta\to 1^{-}$,
$n_{2}^{\mathrm{opt}}/n_{1}^{\mathrm{opt}}\to 1^{-}$, and
$n_{i}^{\mathrm{opt}}\sim C_{\bar{P}}(1-\zeta)^{-2}$ for $i=1,2$ with
$C_{\bar{P}}=2/(\sigma_{\mathscr{D}}^{2}\mathsf{acv}_{0}^{2})$.

At low-illumination, from Lemma 3, the optimal sample sizes for $\bar{P}$ are
asymptotically proportional to $(1-\zeta)^{-2}$, demonstrating that as
expected, infinite sample sizes are needed in the zero illumination limit.

###  5.2 Statistical analysis of hat-$P^{-1}_{\delta}$

We now carry out a similar but more detailed analysis for
$\hat{P}^{-1}_{\delta}$. Under the normal model of Section 2 we have the
distributional results $\hat{X}\sim\mathcal{G}(\alpha_{1},\beta_{1})$ and
$\hat{Y}\sim\mathcal{G}(\alpha_{2},\beta_{2})$ (see Table 2). As such, the
difference $\hat{P}=\hat{X}-\hat{Y}$ is distributed as
$\hat{P}\sim\mathcal{GD}(\alpha_{1},\alpha_{2},\beta_{1},\beta_{2})$, which is
a gamma-difference variable. Working with the properties of the gamma-
difference distribution as well as Lemma 2 we can write

$\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta}=\frac{2}{(1-\zeta)^{2}}\left(\frac{1}{n_{1}-1}+\frac{\zeta^{2}}{n_{2}-1}\right),$
(59)

which leads to the following expressions for optimal sample sizes of
$\hat{P}^{-1}_{\delta}$.

######  Lemma 4.

The optimal sample sizes for $\hat{P}^{-1}_{\delta}$ are given by,

$(n_{1}^{\mathrm{opt}},n_{2}^{\mathrm{opt}})=\left(\frac{2(1+\zeta)}{\mathsf{acv}_{0}^{2}(1-\zeta)^{2}}+1,\frac{2\zeta(1+\zeta)}{\mathsf{acv}_{0}^{2}(1-\zeta)^{2}}+1\right).$
(60)

######  Remark 1.

Our goal for deriving the optimal sample sizes of $\hat{P}^{-1}_{\delta}$ was
for the purpose of optimally estimating $\hat{P}^{-1}$. It turns out that we
can alter the sample sizes of Lemma 4 to make them exact in the shot-noise-
limit. See Remark 2 in the appendix for more information.

######  Proposition 1.

Let $n_{1}^{\mathrm{opt}}$ and $n_{2}^{\mathrm{opt}}$ denote the optimal
sample sizes for $\hat{P}^{-1}_{\delta}$, then:

  1. (i)

$n_{1}^{\mathrm{opt}}$ and $n_{2}^{\mathrm{opt}}$ are strictly increasing with
decreasing illumination ($\zeta$ increasing),

  2. (ii)

$n_{2}^{\mathrm{opt}}<n_{1}^{\mathrm{opt}}$ for all $\zeta\in[0,1)$,

  3. (iii)

As illumination decreases, $\zeta\to 1^{-}$,
$n_{2}^{\mathrm{opt}}/n_{1}^{\mathrm{opt}}\to 1^{-}$, and
$n_{i}^{\mathrm{opt}}\sim C_{\hat{P}^{-1}_{\delta}}(1-\zeta)^{-2}$ for $i=1,2$
where $C_{\hat{P}^{-1}_{\delta}}=4/\mathsf{acv}_{0}^{2}$.

######  Corollary 1.

Let $n_{1}^{\mathrm{opt}}(\zeta)$ and $n_{2}^{\mathrm{opt}}(\zeta)$ denote the
optimal sample sizes for $\hat{P}^{-1}_{\delta}$ as a function of $\zeta$. To
estimate $\hat{P}^{-1}_{\delta}$ to a relative uncertainty of
$\mathsf{acv}_{0}$, one needs at minimum a total of
$N^{\mathrm{opt}}(0)=n_{1}^{\mathrm{opt}}(0)+n_{2}^{\mathrm{opt}}(0)=2/\mathsf{acv}_{0}+2$
observations.

Proposition 1 demonstrates that the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ have similar characteristics for those of $\bar{P}$.
In particular, we see at low-illumination that the optimal sample sizes are
asymptotically equal and asymptotically proportional to $(1-\zeta)^{-2}$
showing that infinite sample sizes are again needed in the zero illumination
limit. Furthermore, Corollary 1 takes advantage of the monotonicity of the
optimal sample sizes in $\zeta$ in order to derive a lower bound to the
minimum total number of samples needed to achieve an absolute coefficient of
variation for $\hat{P}^{-1}_{\delta}$ equal to $\mathsf{acv}_{0}$. Note that
this lower bound closely matches the asymptotic result for shot-noise-limited
estimator of $g$ given in (25). Choosing $\mathsf{acv}_{0}=0.01$ gives the
lower bound $N^{\mathrm{opt}}(0)=20\,002$, which serves to show the large
sample sizes needed to obtain typically desired measurement uncertainties.

Now that we have explicit expressions for the optimal sample sizes of
$\hat{P}^{-1}_{\delta}$, and know their properties, we can assess the accuracy
of the approximation obtained by using the optimal sample sizes of
$\hat{P}^{-1}_{\delta}$ in place of those for $\hat{P}^{-1}$.

######  Theorem 2.

Let $\hat{P}^{-1}_{\mathrm{opt}}$ denote the estimator $\hat{P}^{-1}$ as a
function of the optimal sample sizes for $\hat{P}^{-1}_{\delta}$. Then as
$\mathsf{acv}_{0}\to 0^{+}$

$\displaystyle\mathsf{ACV}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}$
$\displaystyle=\mathsf{acv}_{0}+3\mathsf{acv}_{0}^{3}+\mathcal{O}(\mathsf{acv}_{0}^{5}),$
(61) $\displaystyle\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}$
$\displaystyle=\mathsf{acv}_{0}^{2}+3\mathsf{acv}_{0}^{4}+\mathcal{O}(\mathsf{acv}_{0}^{6}).$

As Theorem 2 shows, when $\mathsf{acv}_{0}$ is small, the optimal samples
sizes for $\hat{P}^{-1}_{\delta}$ force
$\mathsf{ACV}_{\mathcal{P}}\hat{P}^{-1}$ to be nearly $\mathsf{acv}_{0}$,
which indicates a good approximation. Since the secondary term in the
expansion for $\mathsf{ACV}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}$ is
cubic, we expect the approximation to be quite good for any small choice of
$\mathsf{acv}_{0}$, say $\mathsf{acv}_{0}<0.1$. This restriction is not
problematic as one generally does not seek to measure $1/\sigma_{P}^{2}$ with
greater than $10\%$ uncertainty. Furthermore, Theorem 2 also quantifies the
bias of $\hat{P}^{-1}_{\mathrm{opt}}$ showing that the main term is quadratic
in $\mathsf{acv}_{0}$. As such, we expect the estimator $\hat{P}^{-1}$ to be
nearly unbiased for $1/\sigma_{P}^{2}$ (in the sense that
$\mathsf{E}_{\mathcal{P}}\hat{P}^{-1}\approx 1/\sigma_{P}^{2}$) when subject
to the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ for a small choice of
$\mathsf{acv}_{0}$.

###  5.3 Statistical analysis of $G_{\delta}$

Equipped with the results of Subsections 5.1-5.2, we are now ready to perform
a statistical analysis on $G_{\delta}$ with the goal of determining when the
optimal sample sizes for $\hat{P}^{-1}_{\delta}$ serve as good approximations
for those of $G_{\delta}$. Using the formula in Lemma 1 we first write the
squared absolute coefficient of variation for $G_{\delta}$ as

$\mathsf{ACV}^{2}G_{\delta}=\mathsf{ACV}^{2}\bar{P}+(\mathsf{ACV}^{2}\bar{P})(\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta})+\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta},$
(62)

with $\mathsf{ACV}^{2}\bar{P}$ and $\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta}$
given in (58) and (59), respectively. Our first goal here is to show the
estimator $G_{\delta}$ behaves like previous estimators for $g$ in the sense
that its uncertainty is dominated by that of $\hat{P}^{-1}_{\delta}$ at high-
illumination (c.f. Theorem 1).

######  Theorem 3.

Let $G_{\delta}$ and $\hat{P}^{-1}_{\delta}$ be as given in (54) and (55),
respectively. As illumination increases, $\zeta\to 0^{+}$ and

$\mathsf{ACV}G_{\delta}=\mathsf{ACV}\hat{P}^{-1}_{\delta}\left(1+\frac{1}{\sigma_{\mathscr{D}}^{2}}\frac{n_{1}+1}{4n_{1}}\zeta+\mathcal{O}(\zeta^{2})\right).$
(63)

Since $\mathsf{ACV}G_{\delta}\sim\mathsf{ACV}\hat{P}^{-1}_{\delta}$ we have by
the same reasoning in Section 3 that the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ are asymptotically equal to those of $G_{\delta}$ at
high-illumination with equality in the shot-noise-limit $(\zeta=0)$. However,
if our goal is optimal sampling for the general case, this approximation must
also hold in the low-illumination regime $(\zeta\to 1^{-})$. To determine how
good this approximation is at low-illumination, we may consider the following
theorem.

######  Theorem 4.

Let $n_{1}^{\mathrm{opt}}$ and $n_{2}^{\mathrm{opt}}$ denote the optimal
sample sizes for $G_{\delta}$. Then, as illumination decreases, $\zeta\to
1^{-}$, $n_{2}^{\mathrm{opt}}/n_{1}^{\mathrm{opt}}\to 1^{-}$, and
$n_{i}^{\mathrm{opt}}\sim C_{G_{\delta}}(1-\zeta)^{-2}$ for $i=1,2$ with

$C_{G_{\delta}}=\frac{2}{\mathsf{acv}_{0}^{2}}\left(1+\frac{1}{2\sigma_{\mathscr{D}}^{2}}+\left(\left(1+\frac{1}{2\sigma_{\mathscr{D}}^{2}}\right)^{2}+\frac{2}{\sigma_{\mathscr{D}}^{2}}\mathsf{acv}_{0}^{2}\right)^{1/2}\right).$
(64)

Comparing Theorem 4 with Proposition 1 $(iii)$ shows that at low-illumination
the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ are asymptotically
proportional to the optimal sample sizes for $G_{\delta}$. The constant of
proportionality, $C_{G_{\delta}}/C_{\hat{P}^{-1}_{\delta}}$, depends on the
dark noise $\sigma_{\mathscr{D}}$ and as the dark noise increases

$\frac{C_{G_{\delta}}}{C_{\hat{P}^{-1}_{\delta}}}\sim
1+\frac{1+\mathsf{acv}_{0}^{2}}{2\sigma_{\mathscr{D}}^{2}}+\mathcal{O}(\sigma_{\mathscr{D}}^{-4}),$
(65)

which further shows that this constant of proportionality approaches one with
increasing dark noise. Bringing the observations following Theorems 3-4
together, we expect the optimal sample sizes of $\hat{P}^{-1}_{\delta}$ to
serve as good approximations to those of $G_{\delta}$ at any illumination
level given enough dark noise.

To get a better grasp on these observations we will define the metric

$\mathcal{E}=\frac{\mathsf{ACV}\hat{P}^{-1}_{\delta}}{\mathsf{ACV}G_{\delta}}.$
(66)

and consider this metric as a function of the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$, that is,

$\mathcal{E}_{\mathrm{opt}}=\frac{\mathsf{acv}_{0}}{\mathsf{ACV}G_{\delta,\mathrm{opt}}}=\left(1+(\mathsf{ACV}^{2}\bar{P}_{\mathrm{opt}})(1+\mathsf{acv}_{0}^{-2})\right)^{-1/2}$
(67)

with

$\mathsf{ACV}^{2}\bar{P}_{\mathrm{opt}}=\frac{1}{\sigma_{\mathscr{D}}^{2}}\frac{\zeta}{(1-\zeta)^{2}}\left(\frac{1}{n_{1}^{\mathrm{opt}}}+\frac{\zeta}{n_{2}^{\mathrm{opt}}}\right).$
(68)

Both $\mathcal{E}$ and it’s counterpart $\mathcal{E}_{\mathrm{opt}}$ are
normalized in the sense that $0\leq\mathcal{E}_{\mathrm{opt}}\leq 1$. In
particular, whenever $\mathcal{E}_{\mathrm{opt}}\approx 1$ we can conclude
that the optimal sample sizes of $\hat{P}^{-1}_{\delta}$ are good
approximations for those of $G_{\delta}$ with equality achieved when
$\mathcal{E}_{\mathrm{opt}}=1$.

It is straightforward to show $\lim_{\zeta\to
0^{+}}\mathcal{E}_{\mathrm{opt}}=1$, which again indicates that the optimal
sample sizes for $\hat{P}^{-1}_{\delta}$ are exactly equal to those for
$G_{\delta}$ in the shot-noise-limit. Taking the limit now in the opposite
direction, we see
$\mathsf{ACV}^{2}\bar{P}_{\mathrm{opt}}\to\mathsf{acv}_{0}^{2}/(2\sigma_{\mathscr{D}}^{2})$
as $\zeta\to 1^{-}$ and so

$\overline{\mathcal{E}_{\mathrm{opt}}}\coloneqq\lim_{\zeta\to
1^{-}}\mathcal{E}_{\mathrm{opt}}=\left(1+\frac{1+\mathsf{acv}_{0}^{2}}{2\sigma_{\mathscr{D}}^{2}}\right)^{-1/2},$
(69)

which is nonzero. To get a sense of how close
$\overline{\mathcal{E}_{\mathrm{opt}}}$ is to one recall from the discussion
following Theorem 2, that we generally want to impose the restriction
$\mathsf{acv}_{0}\in(0,0.1)$. Since $\overline{\mathcal{E}_{\mathrm{opt}}}$ is
decreasing in $\mathsf{acv}_{0}$ we set $\mathsf{acv}_{0}=0.1$ giving the
lower bound

$\overline{\mathcal{E}_{\mathrm{opt}}}>\left(1+\frac{101}{200}\frac{1}{\sigma_{\mathscr{D}}^{2}}\right)^{-1/2}.$
(70)

Plotting this lower bound as a function of $\sigma_{\mathscr{D}}$ we find
$\overline{\mathcal{E}_{\mathrm{opt}}}>0.99$ for $\sigma_{\mathscr{D}}\geq
5\,e\text{-}$ and so for sensors with dark noise greater than $5\,e\text{-}$,
we can expect the optimal sample sizes of $\hat{P}^{-1}_{\delta}$ to give
excellent approximations to those of $G_{\delta}$ at any level of
illumination. In conclusion, if $\mathsf{acv}_{0}$ is small
$(\mathsf{acv}_{0}<0.1)$ and the dark noise is sufficiently large
$(\sigma_{\mathscr{D}}\geq 5\,e\text{-})$, the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ serve as excellent approximations for those of the
general conversion gain estimator $G$ at any illumination level.

We close this section with an analogous result to that of Theorem 2, which
confirms our findings.

######  Theorem 5.

Let $G_{\mathrm{opt}}$ denote the estimator $G$ as a function of the optimal
sample sizes for $\hat{P}^{-1}_{\delta}$. Then as $\mathsf{acv}_{0}\to 0^{+}$

$\displaystyle\mathsf{ACV}_{\mathcal{P}}G_{\mathrm{opt}}$
$\displaystyle=\sqrt{1+a}\,\mathsf{acv}_{0}+\frac{6+a+b}{2\sqrt{1+a}}\mathsf{acv}_{0}^{3}+\mathcal{O}(\mathsf{acv}_{0}^{5}),$
(71) $\displaystyle\mathsf{ARB}_{\mathcal{P}}G_{\mathrm{opt}}$
$\displaystyle=\mathsf{acv}_{0}^{2}+3\mathsf{acv}_{0}^{4}+\mathcal{O}(\mathsf{acv}_{0}^{6}),$

where $a=\frac{1}{\sigma_{\mathscr{D}}^{2}}\frac{\zeta}{1+\zeta}$ and
$b=-\frac{1}{\sigma_{\mathscr{D}}^{2}}\frac{(1-\zeta)^{2}}{4(1+\zeta)}$.

Notice that in particular, if $\sigma_{\mathscr{D}}$ is large, and/or $\zeta$
is small, the main term in the expansion for
$\mathsf{ACV}_{\mathcal{P}}G_{\mathrm{opt}}$ is close to $\mathsf{acv}_{0}$
showing that the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ are good
approximations for those of $G$ under these conditions.

##  6 Design and control of experiment for per-pixel conversion gain
estimation

In this section, we demonstrate how the derived expressions can be used in the
Design & Control of Experiment (doe & coe) for per-pixel conversion gain
estimation. For this example we will use the ON Semiconductor KAI-04070
monochrome interline transfer ccd sensor [29]. As we will demonstrate, the
entire experimental process for per-pixel conversion gain estimation centers
around the optimal sample sizes for $\hat{P}^{-1}_{\delta}$ in Lemma 4 and the
metric $\mathcal{E}_{\mathrm{opt}}$ given in (67).

A full Monte Carlo implementation of the doe and coe algorithms presented in
this section can be found in Code 1 [30] and Code 2 [31]. Code 1
(DOE_and_COE.m) is the primary script and requires matlab’s _Statistics and
Machine Learning_ toolbox to run. Since Code 1 is meant to guide the reader in
understanding the steps involved in doe and coe it is advised to run it one
section at a time and observe the results of each section.

###  6.1 Experimental Setup

Figure 2 shows a schematic diagram of the experimental setup used for
performing per-pixel conversion gain estimation. The experimental setup
consisted of a $650\,\mathrm{nm}$ Superluminescent Light Emitting Diode (sled)
passed into an integrating sphere with the Sensor Under Test (sut) placed in
the plane of the spheres output port ($f/0$ geometry) where uniformity is
highest [32]. To facilitate control of the illumination level, a Variable
Optical Attenuator (voa) was introduced between the sled and integrating
sphere. The sensor was configured at its full bit-depth of $14$-bits to
minimize quantization error and image data was read off the sensor at a
$512\times 512\,\mathrm{px}$ resolution using a single readout register
operating at $40\,\mathrm{mhz}$. By operating the sensor in a single-tap mode
like this, the gain of each pixel will be the same and so any excursions in
the measured gain of each pixel will be, in theory, completely due to sampling
error. This high level of uniformity will allow us to see if the proposed
algorithms are able to successfully measure the gain to the desired relative
uncertainty, verify our noise model, and compare experimental results to
theoretical predictions.

In order to capture imagery under both dark and illuminated conditions, a
Motorized Mirror (mm) was placed next to the path of the sled beam. Moving the
mirror into the beam path redirected the beam away from the integrating sphere
and into a Beam Dump (bd); thus, providing a dark environment for the sensor.

Figure 2: Schematic diagram of experimental setup for per-pixel conversion
gain estimation.

###  6.2 Design of experiment

doe for per-pixel conversion gain estimation begins with choosing a suitable
value for the desired relative uncertainty $\mathsf{acv}_{0}$. As a rough rule
of thumb, if we are operating under the condition
$\mathcal{E}_{\mathrm{opt}}\approx 1$, then for small $\mathsf{acv}_{0}$ we
may approximate
$G_{\mathrm{opt}}\sim\mathcal{N}(g,(g\times\mathsf{acv}_{0})^{2})$ so that
$\mathsf{P}(G_{\mathrm{opt}}\in g(1\pm\mathsf{acv}_{0}))=0.683$. This amounts
to the optimal estimate of $g$ being within $\mathsf{acv}_{0}\times 100\%$ of
its exact value $68.3\%$ of the time. For this particular experiment it was
decided to use $\mathsf{acv}_{0}=0.05$ as this represents a typical value one
might choose.

Next, we want to determine an appropriate illumination level for measuring the
conversion gain by selecting a value of $\zeta$ where: (1) the total number of
required observations is attainable and (2) our approximate optimal sample
sizes are valid. To determine when both conditions are satisfied, we will
create an $\mathcal{E}$-$N$ plot, which consists of plotting the functions
$\mathcal{E}_{\mathrm{opt}}$ and
$N^{\mathrm{opt}}=n_{1}^{\mathrm{opt}}+n_{2}^{\mathrm{opt}}$ as a function of
$\zeta$. To determine an illumination level, we select a value of $\zeta$
where $N^{\mathrm{opt}}$ is small enough and
$\mathcal{E}_{\mathrm{opt}}\approx 1$. However, because
$\mathcal{E}_{\mathrm{opt}}$ contains the unknown value of the dark noise, we
will need to provide an estimate of $\sigma_{\mathscr{D}}$ for this plot,
which can be accomplished through a preliminary measurement, a vendor
specification sheet, or an educated guess/lower bound.

In this example we will take the route of estimating a global lower bound for
the dark noise. For this particular sensor we expect $g\geq 1$ and since
$\sigma_{\mathscr{D}}=\sigma_{D}\times g$ it follows that $\sigma_{D}$
provides a lower bound on $\sigma_{\mathscr{D}}$. Since
$\mathcal{E}_{\mathrm{opt}}$ is an increasing function of
$\sigma_{\mathscr{D}}$, it follows that evaluating it at $\sigma_{D}$ provides
a lower bound for it’s exact value. To obtain a global estimate of
$\sigma_{D}^{2}$ we capture two dark frames $\mathbf{Y}_{1}$ and
$\mathbf{Y}_{2}$ and then compute half the sample variance of the difference-
frame $\Delta\mathbf{Y}=\mathbf{Y}_{1}-\mathbf{Y}_{2}$. Note the use of bold
symbols to denote $\mathbf{Y}_{k}$ as a two-dimensional array so that
$\Delta\mathbf{Y}$ is the pixel-wise difference of the two dark frames
$\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$. For the sensor under test we found

$\hat{\sigma}_{D}^{2}=\frac{1}{2}\mathrm{var}(\Delta\mathbf{Y})=39.94\,\mathrm{DN}^{2},$
(72)

where

$\mathrm{var}(\Delta\mathbf{Y})=\frac{1}{512^{2}-1}\sum_{i=1}^{512}\sum_{j=1}^{512}(\Delta\mathbf{Y}_{ij}-\overline{\Delta\mathbf{Y}})^{2}.$
(73)

Figure 3 plots

$\hat{\mathcal{E}}_{\mathrm{opt}}(\zeta)=\\\
\left(1+\frac{1+\mathsf{acv}_{0}^{-2}}{\hat{\sigma}_{D}^{2}}\frac{\zeta}{(1-\zeta)^{2}}\left(\frac{1}{n_{1}^{\mathrm{opt}}(\zeta)}+\frac{\zeta}{n_{2}^{\mathrm{opt}}(\zeta)}\right)\right)^{-1}$
(74)

along with $N^{\mathrm{opt}}(\zeta)$ for our choice of $\mathsf{acv}_{0}=0.05$
and estimate $\hat{\sigma}_{D}$.

Figure 3: $\mathcal{E}$-$N$ plot (top) with optimal sample sizes (bottom) for
$\mathsf{acv}_{0}=0.05$ versus $\zeta$.

Due to the sufficiently large value of $\hat{\sigma}_{D}$ we see that
$\hat{\mathcal{E}}_{\mathrm{opt}}$ is near unity for virtually any
illumination level so that the requirement $\mathcal{E}_{\mathrm{opt}}\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 we observe that
the illumination level corresponding to $\zeta=0.4$ is paired with an optimal
total sample size of $N^{\mathrm{opt}}\approx 4300$ and approximation quality
metric $\mathcal{E}_{\mathrm{opt}}\approx 0.993$. At a recording rate of
$5\,\mathrm{fps}$ this number of images will take $\approx 15\,\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$, we guessed the required
illumination level by adjusting the variable optical attenuator until the
illumination level at the sensor plane resulted in a mean pixel output of
about $1\%$ of the sensor’s dynamic range. Using the same process to estimate
the dark noise, two illuminated frames $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$
were captured and their difference
$\Delta\mathbf{X}=\mathbf{X}_{1}-\mathbf{X}_{2}$ was used to obtain the
estimate

$\hat{\sigma}_{P+D}^{2}=\frac{1}{2}\mathrm{var}(\Delta\mathbf{X})=112.89\,\mathrm{DN}^{2},$
(75)

which lead to a global estimate of $\zeta$ equal to

$\hat{\zeta}=\frac{\hat{\sigma}_{D}^{2}}{\hat{\sigma}_{P+D}^{2}}=0.354.$ (76)

For our purposes, this illumination level was sufficiently close to the target
$\zeta=0.4$ and corresponded to $\hat{\mathcal{E}}_{\mathrm{opt}}=0.996$ and
$N^{\mathrm{opt}}=3521$ images, which needs only $\approx 12\,\mathrm{min}.$
to capture. If this $\zeta$-value was not appropriate we could simply re-
adjust the illumination level, capture another two illuminated frames, and re-
estimate $\zeta$. We can repeat this procedure until we have found an
illumination level corresponding to a suitable $\zeta$-value.

Before moving on to data capture, its important to understand what (76) is an
estimate of. Recall that
$\mathbf{Y}_{ijk}\sim\mathcal{N}(\mu_{D},\sigma_{D}^{2})$ so that we have for
the difference frame
$\mathbf{Y}_{ij1}-\mathbf{Y}_{ij2}=\Delta\mathbf{Y}_{ij}\sim\mathcal{N}(0,2\sigma_{D}^{2})$.
In general $\sigma_{D}^{2}$ will vary from pixel-to-pixel (the $ij$-dimension)
so that we may treat it as a random variable and model it according to some
probability distribution $\sigma_{D}^{2}\sim F_{\sigma_{D}^{2}}$. By the law
of total variance one finds

$\mathsf{Var}(\Delta\mathbf{Y}_{ij})=\mathsf{Var}(\mathsf{E}(\Delta\mathbf{Y}_{ij}|\sigma_{D}^{2}))+\mathsf{E}(\mathsf{Var}(\Delta\mathbf{Y}_{ij}|\sigma_{D}^{2}))=2\,\mathsf{E}(\sigma_{D}^{2}),$
(77)

which shows that the estimator (72) gives an unbiased estimate of the average
value of $\sigma_{D}^{2}$ across the sensor array and likewise for the
estimator of $\sigma_{P+D}^{2}$ in (75). Thus, $\hat{\zeta}$ is a ratio of
unbiased estimates for the average noise values across the sensor array;
providing a useful global estimate of $\zeta$.

###  6.3 Control of experiment

With the illumination level set, the experiment was ready to commence. The
algorithm for data capture is presented in Algorithm 1. This algorithm
utilizes Welford’s online algorithm in UpdateStats() to iteratively update the
master frames $\bar{\mathbf{X}}$, $\bar{\mathbf{Y}}$, $\hat{\mathbf{X}}$, and
$\hat{\mathbf{Y}}$, which contain sample means and variances for each pixel in
both dark and illuminated conditions. In each iteration of the algorithm a
batch of illuminated frames and another batch of dark frames are captured and
used to update the master frames before recalculating a global estimate of
$\zeta$, which is then in turn used to update the next batch sizes. The
initial batch sizes are $\mathtt{batch1}=\lceil
n_{1}^{\mathrm{opt}}(0,\mathsf{acv}_{0})\rceil$ and $\mathtt{batch2}=2$, which
are the minimal possible number of each frame type needed. The batch sizes are
updated to capture $m\times 100\%$ ($m\in(0,1]$) of the remaining difference
between the current sample sizes and their estimates. The parameter $m$
controls how aggressive the algorithm is and ultimately how many times the
light source needs to be turned on and off. Large values of $m$ mean the
algorithm will iterate less times (putting more confidence in the estimated
sample sizes) while smaller values of $m$ result in more iterations (less
confidence in the estimates). In the limit $m\to 0$, the algorithm iterates
every time a dark and illuminated frame are captured. The algorithm terminates
when both batch sizes are nonpositive indicating that the current sample sizes
meet or exceed their respective estimates. This is then followed by the per-
pixel calculation of $g$, denoted $\mathbf{G}$, which we shall call the
$g$-map.

While the form of $\hat{\zeta}$ in Algorithm 1 has a much different form than
that of (76), we can show it still is estimating the same quantity. To see why
note that $\hat{\mathbf{Y}}_{ij}\sim\mathcal{G}(\alpha_{2},\beta_{2})$ with
$\sigma_{D}^{2}\sim F_{\sigma_{D}^{2}}$. By the law of total expectation

$\mathsf{E}(\hat{\mathbf{Y}}_{ij})=\mathsf{E}(\mathsf{E}(\hat{\mathbf{Y}}_{ij}|\sigma_{D}^{2}))=\mathsf{E}(\sigma_{D}^{2}),$
(78)

which shows $\sum\hat{\mathbf{Y}}_{ij}$ is an unbiased estimator for $I\cdot
J\cdot\mathsf{E}(\sigma_{D}^{2})$ with $I$ and $J$ being the vertical and
horizontal resolution of the sensor in units of pixels, respectively.
Likewise, $\sum\hat{\mathbf{X}}_{ij}$ is an unbiased estimator for $I\cdot
J\cdot\mathsf{E}(\sigma_{P+D}^{2})$ so that
$\hat{\zeta}=\sum\hat{\mathbf{Y}}_{ij}/\sum\hat{\mathbf{X}}_{ij}$ is
equivalent to a ratio of unbiased estimates for $\mathsf{E}(\sigma_{D}^{2})$
and $\mathsf{E}(\sigma_{P+D}^{2})$ just as (76) is.

Algorithm 1 coe algorithm.

1:procedure OptimalCOE($\mathsf{acv}_{0}$,$m$)

2: $n_{1}=0$; $n_{2}=0$; $\hat{\zeta}=0;$

3: $\mathtt{batch1}=\lceil
n_{1}^{\mathrm{opt}}(\hat{\zeta},\mathsf{acv}_{0})\rceil$;

4: $\mathtt{batch2}=2$;

5:

6: while $\mathtt{batch1}>0\lor\mathtt{batch2}>0$ do

7: if $\mathtt{batch1}>0$ then

8: Turn on light source.

9: $\mathtt{NewData}=$ Capture $\mathtt{batch1}$ $\mathbf{X}$-frame(s);

10: $n_{1}=n_{1}+\mathtt{batch1}$;

11:
$[\bar{\mathbf{X}},\hat{\mathbf{X}}]=\text{\sc{UpdateStats}}(\bar{\mathbf{X}},\hat{\mathbf{X}},\mathtt{NewData})$;

12: if $\mathtt{batch2}>0$ then

13: Turn off light source.

14: $\mathtt{NewData}=$ Capture $\mathtt{batch2}$ $\mathbf{Y}$-frame(s);

15: $n_{2}=n_{2}+\mathtt{batch2}$;

16:
$[\bar{\mathbf{Y}},\hat{\mathbf{Y}}]=\text{\sc{UpdateStats}}(\bar{\mathbf{Y}},\hat{\mathbf{Y}},\mathtt{NewData})$;

17: $\hat{\zeta}=\sum\hat{\mathbf{Y}}_{ij}/\sum\hat{\mathbf{X}}_{ij}$;

18: $\mathtt{batch1}=\lceil
m\times(n_{1}^{\mathrm{opt}}(\hat{\zeta},\mathsf{acv}_{0})-n_{1})\rceil$;

19: $\mathtt{batch2}=\lceil
m\times(n_{2}^{\mathrm{opt}}(\hat{\zeta},\mathsf{acv}_{0})-n_{2})\rceil$;

20:

21:
$\mathbf{G}=(\bar{\mathbf{X}}-\bar{\mathbf{Y}})/(\hat{\mathbf{X}}-\hat{\mathbf{Y}})$;
$\triangleright$ Calculation is per-pixel.

22: return $\mathbf{G}$

###  6.4 Experimental Data

Algorithm 1 was executed on the chosen ccd for the specified uncertainty
$\mathsf{acv}_{0}=0.05$ and multiplier $m=0.8$, which halted after capturing
$n_{1}=2602$ illuminated frames and $n_{2}=921$ dark frames. Due to the large
value of $m$, the mirror only needed to be moved eleven times throughout the
entire coe procedure. Figure 4 presents the final $g$-map created from
Algorithm 1 along with a histogram of the $g$-map values. Upon inspection of
the $g$-map, there is no evidence of structure or nonuniformities as was
expected. Because the excursions in the $g$-map are almost entirely due to
random sampling error, as opposed to nonuniformity in the sensor or source,
and the sensor noise obeys the assumed noise model, the estimated value of $g$
for each pixel represents an i.i.d. observation from the cif distribution
[10]. Given the large values of $n_{1}$ and $n_{2}$, the fitted cif
distribution for this example was approximated by Hinkley’s normal ratio
distribution [33]

$f_{G}(g)\sim\frac{e^{-c/2}}{\sigma_{\bar{P}}\sigma_{\hat{P}}a^{2}(g)}\left(\frac{b(g)e^{b^{2}(g)/2}}{\sqrt{2\pi}}(2\Phi(b(g))-1)+\frac{1}{\pi}\right)$
(79)

where $a(g)=(g^{2}/\sigma_{\bar{P}}^{2}+1/\sigma_{\hat{P}}^{2})^{1/2}$,
$b(g)=(\mu_{\bar{P}}g/\sigma_{\bar{P}}^{2}+\mu_{\hat{P}}/\sigma_{\hat{P}}^{2})/a(g)$,
and
$c=\mu_{\bar{P}}^{2}/\sigma_{\bar{P}}^{2}+\mu_{\hat{P}}^{2}/\sigma_{\hat{P}}^{2}$.
The parameters $\mu_{\bar{P}}$, $\sigma_{\bar{P}}$, $\mu_{\hat{P}}$, and
$\sigma_{\hat{P}}$ are easily estimated from sample means and sample standard
deviations of the $\bar{\mathbf{P}}=\bar{\mathbf{X}}-\bar{\mathbf{Y}}$ and
$\hat{\mathbf{P}}=\hat{\mathbf{X}}-\hat{\mathbf{Y}}$ master frames.

Given that the fitted distribution is that of a normal ratio, the
pseudomoments of the fitted distribution are exactly described by (45). Table
3 presents the exact values of the pseudomoments for the fitted distribution
and compares them to sample moments of the $g$-map. Upon inspection, we see
the theoretical pseudomoments show a high level of agreement with the sample
moments and thus demonstrate that the pseudomoments serve as a useful way to
characterize moments of actual sensor data.

Table 3: Comparison of $g$-map sample moments to fitted distribution pseudomoments. | fit | sample | error ($\%$)
---|---|---|---
$\mathsf{E}_{\mathcal{P}}G$ | $2.1971$ | $2.1972$ | $4.3508\times 10^{-3}$
$\sqrt{\mathsf{Var}_{\mathcal{P}}G}$ | $0.11069$ | $0.11122$ | $0.47739$
Figure 4: Histogram of $g$-map values fit with the cif distribution (top) and
$g$-map (bottom).

Recall that the coe algorithm was run with a design parameter of
$\mathsf{acv}_{0}=0.05$; comparing this to the measurement, we find the sample
absolute coefficient of variation of the $g$-map to be

$\mathrm{acv}(\mathbf{G})=0.05059.$ (80)

Here,
$\mathrm{acv}(\mathbf{G})=\sqrt{\mathrm{var}(\mathbf{G})}/|\mathrm{mean}(\mathbf{G})|$
with $\mathrm{mean}(\mathbf{G})$ and $\mathrm{var}(\mathbf{G})$ denoting the
sample mean and variance of the two-dimensional $g$-map, respectively. We note
that $\mathrm{acv}(\mathbf{G})>\mathsf{acv}_{0}$. The discrepancy between
$\mathrm{acv}(\mathbf{G})$ and $\mathsf{acv}_{0}$ is due in part to variance
in the coe procedure and other factors like having to estimate the optimal
sample sizes. But even if the algorithm was deterministic and we knew the
optimal sample sizes exactly, Theorem 5 tells us that the absolute coefficient
of variation for the $g$-map should have a positive bias away from the desired
value of $\mathsf{acv}_{0}$.

To study this effect more, a Monte Carlo experiment was setup to replicate the
sensor data presented above. The only parameters we need to estimate for the
experiment are $g$, $\mu_{\mathscr{D}}$, $\sigma_{\mathscr{D}}$, and
$\mu_{e\text{-}}$. To obtain a good estimate for $g$ we use the definition of
relative bias $\mathsf{RB}_{\mathcal{P}}G=(\mathsf{E}_{\mathcal{P}}G-g)/g$ to
write the first moment of $G$ as

$\mathsf{E}_{\mathcal{P}}G=g(1+\mathsf{RB}_{\mathcal{P}}G).$ (81)

In the context of conversion gain estimation we expect the relative bias to be
positive, i.e. $\mathsf{RB}_{\mathcal{P}}G=\mathsf{ARB}_{\mathcal{P}}G$ and
since the $g$-map, $\mathbf{G}$, was computed using the optimal sample sizes
we have according to Theorem 5:
$\mathsf{ARB}_{\mathcal{P}}G_{\mathrm{opt}}\sim\mathsf{acv}_{0}^{2}+3\mathsf{acv}_{0}^{4}$.
This leads to the bias corrected estimate

$\hat{g}=\mathrm{mean}(\mathbf{G})/(1+\mathsf{acv}_{0}^{2}+3\mathsf{acv}_{0}^{4})$
(82)

with $\mathsf{acv}_{0}=0.05$. Table 4 presents the estimate for $g$ along with
the other estimates needed for the experiment. We note that these parameter
estimates are used as the simulation parameters in Code 1 [30].

Table 4: Parameter estimates for Monte Carlo experiment. parameter | estimate | value
---|---|---
$\hat{g}$ | $\mathrm{mean}(\mathbf{G})/(1+\mathsf{acv}_{0}^{2}+3\mathsf{acv}_{0}^{4})$ | $2.1917$
$\hat{\mu}_{\mathscr{D}}$ | $\mathrm{mean}(\bar{\mathbf{Y}})\times\hat{g}$ | $92.858$
$\hat{\sigma}_{\mathscr{D}}$ | $\sqrt{\mathrm{mean}(\hat{\mathbf{Y}})}\times\hat{g}$ | $13.853$
$\hat{\mu}_{e\text{-}}$ | $\mathrm{mean}(\bar{\mathbf{X}}-\bar{\mathbf{Y}})\times\hat{g}$ | $350.03$

With the parameter estimates in Table 4, Algorithm 1 was run a total of ten
times using the design parameters $\mathsf{acv}_{0}=0.05$ and $m=0.8$. To
simulate the sensor data, $512\times 512\mathrm{px}$ dark images,
$\mathbf{Y}$, and illuminated images, $\mathbf{X}$, were generated according
to the model in Section 2 so that they contained elements of the form

$\displaystyle\mathbf{Y}_{ij}=\lceil\mathscr{D}/\hat{g}\rfloor$ (83)
$\displaystyle\mathbf{X}_{ij}=\lceil(\mathscr{P}+\mathscr{D})/\hat{g}\rfloor,$
(84)

with $\mathscr{P}\sim\mathcal{P}(\hat{\mu}_{e\text{-}})$ and
$\mathscr{D}\sim\mathcal{N}(\hat{\mu}_{\mathscr{D}},\hat{\sigma}_{\mathscr{D}}^{2})$.
Table 5 presents the results of the Monte Carlo experiment. We see that the
sample sizes resulting from the simulated coe algorithm were very close to
those from our actual experiment and have very little variance. This indicates
that: (1) our assumed noise model is effective and (2) the difference between
$\mathrm{acv}(\mathbf{G})$ and $\mathsf{acv}_{0}$ in our original experiment
is unlikely to be due to uncertainty in the optimal sample size estimates.
Furthermore, we see that the absolute coefficient of variation for the $g$-map
created by the algorithm shows a distinct positive bias away the desired value
$\mathsf{acv}_{0}$, also with small variance. Such results are indicative that
the difference between $\mathrm{acv}(\mathbf{G})$ and $\mathsf{acv}_{0}$ in
our original experiment is a result of approximating optimal sample sizes for
$G$ with those of $\hat{P}^{-1}_{\delta}$.

Table 5: Results of Monte Carlo experiment for ten runs. | $\mathrm{mean}(\cdot)$ | $\mathrm{var}(\cdot)$ | target
---|---|---|---
$n_{1}$ | $2606.9$ | $0.1$ | $2602$
$n_{2}$ | $923.9$ | $0.1$ | $921$
$\mathrm{acv}(\mathbf{G})$ | $0.050279$ | $5.0328\times 10^{-9}$ |  $0.050409$ (Thm. 5)

##  7 Application: Per-Pixel Read Noise Estimation

With a method for estimating per-pixel conversion gain, we can estimate
downstream parameters such as the read noise $\sigma_{\mathscr{R}}$ on a per-
pixel basis. Recall that the random variable $\mathscr{D}$ represents the dark
noise of the sensor for some nonzero integration time. As such, the noise
represented by $\mathscr{D}$ contains the combined noise from read noise and
dark current shot noise. To isolate the sensor read noise component of
$\mathscr{D}$, given by
$\mathscr{R}\sim\mathcal{N}(\mu_{\mathscr{R}},\sigma_{\mathscr{R}}^{2})$, we
simply set the sensor integration time to $t_{\mathrm{exp}}=0\,\mathrm{sec}$
(or the shortest allowable integration time), which eliminates dark current.
Again by linearity of the transfer function we find for the sensor output at
zero illumination for a zero second integration time to be
$R=\mathcal{T}(\mathscr{R})$ with

$\mu_{R}\coloneqq\mathsf{E}R=\mu_{\mathscr{R}}/g$ (85)

and

$\sigma_{R}^{2}\coloneqq\mathsf{Var}R=\sigma_{\mathscr{R}}^{2}/g^{2}.$ (86)

Upon inspection of the expression for $\sigma_{R}^{2}$ we find an equation for
the read noise as

$\sigma_{\mathscr{R}}=\sigma_{R}\times g.$ (87)

We already have an estimator for $g$ so all that is needed to estimate
$\sigma_{\mathscr{R}}$ is an estimator for $\sigma_{R}$.

To estimate $\sigma_{R}$ we let $\\{Z_{1},\dots,Z_{n_{3}}\\}$ denote a
sequence of $n_{3}$ i.i.d. observations of a pixel captured in the dark over a
zero second integration time. Under the assumed normal model,
$Z_{k}\sim\mathcal{N}(\mu_{R},\sigma_{R}^{2})$; thus we can obtain an unbiased
estimate of $\sigma_{R}^{2}$ via the sample variance

$\hat{Z}=\frac{1}{n_{3}-1}\sum_{k=1}^{n_{3}}(Z_{k}-\bar{Z})^{2}.$ (88)

It follows that we can estimate $\sigma_{R}$ via the sample standard deviation
$\breve{Z}=\sqrt{\hat{Z}}$ although this estimate will no longer be unbiased.
It’s easy to show that $\hat{Z}\sim\mathcal{G}(\alpha_{3},\beta_{3})$ with
$\alpha_{3}=(n_{3}-1)/2$ and $\beta_{3}=\alpha_{3}/\sigma_{R}^{2}$ so that
using the transformation $\breve{Z}=\sqrt{\hat{Z}}$ we find for the density of
$\breve{Z}$

$f_{\breve{Z}}(z)=\frac{2\beta^{\alpha}}{\Gamma(\alpha)}z^{2\alpha-1}e^{-\beta
z^{2}},$ (89)

which is the Nakagami distribution
$\breve{Z}\sim{\mathcal{N}a}(\alpha_{3},\beta_{3})$. One finds for the moments
of $\breve{Z}$

$\mathsf{E}\breve{Z}^{n}=\frac{(\alpha_{3})_{n/2}}{\alpha_{3}^{n/2}}\sigma_{R}^{n},$
(90)

with $(s)_{n}\coloneqq\Gamma(s+n)/\Gamma(s)$ again denoting the Pochhammer
symbol. Thus, an unbiased estimator for $\sigma_{R}$ can be given by

$\breve{Z}^{\ast}=\frac{\sqrt{\alpha_{3}}}{(\alpha_{3})_{1/2}}\sqrt{\hat{Z}}.$
(91)

The read noise $\sigma_{\mathscr{R}}$ of our pixel can therefore be estimated
via

$\breve{\mathscr{R}}=\breve{Z}^{\ast}\times G,$ (92)

with the estimator $G$ given by (31). Since $\breve{Z}^{\ast}$ is computed
from a the sample $\mathbf{Z}$, which is independent from the samples
$\mathbf{X}$ and $\mathbf{Y}$ used to calculate $G$, we further find for the
density of $\breve{R}$

$f_{\breve{R}}(r)=\int_{0}^{\infty}f_{\breve{Z}^{\ast}}(t)f_{G}(r/t)\frac{\mathrm{d}t}{t},$
(93)

with $f_{G}$ again being the cif distribution. We shall simply refer to
$f_{\breve{\mathscr{R}}}$ as the Read Noise (rn) distribution. Likewise, we
have for the pseudomoments of $\breve{\mathscr{R}}$

$\mathsf{E}_{\mathcal{P}}\breve{\mathscr{R}}^{n}=\frac{(\alpha_{3})_{n/2}}{(\alpha_{3})_{1/2}^{n}}\sigma_{R}^{n}\,\mathsf{E}_{\mathcal{P}}G^{n}.$
(94)

To experimentally measure the per-pixel read noise of our sensor, we again
employed an iterative algorithm. In each iteration of this algorithm a master
$\breve{\mathbf{Z}}^{\ast}$-frame was updated with a new $\mathbf{Z}$-frame
using Welford’s algorithm, which was then multiplied, per-pixel, by the
$g$-map estimated in the previous section to produce the
$\sigma_{\mathscr{R}}$-map. The algorithm was stopped when the sample absolute
coefficient of variation for the $\sigma_{\mathscr{R}}$-map satisfied
$\mathrm{acv}(\breve{\mathscr{R}}_{\mathrm{map}})\leq
1.05\,\mathrm{acv}(\mathbf{G})$. This algorithm halted after $n_{3}=2291$
$\mathbf{Z}$-frames were captured.

Figure 5 presents the $\sigma_{\mathscr{R}}$-map generated from this procedure
along with its histogram fit the with rn-distribution of (93). Unlike the
$g$-map, the $\sigma_{\mathscr{R}}$-map does show some column-wise
nonuniformities, which are linked to the interline transfer CCD architecture
of the sensor. However, these nonuniformities are not severe as indicated by
how well the rn-distribution fits the data. What this demonstrates is that
per-pixel maps allow further insight into how sensor architecture affects the
uniformity of key performance parameters across the sensor array.

Figure 5: Histogram of $\sigma_{\mathscr{R}}$-map values fit with the rn
distribution (top) and $\breve{R}$-map (bottom).

##  8 Conclusions

In this work we have presented a general method for sample size determination
of the photon transfer conversion gain measurement given a desired uncertainty
requirement. So long as this uncertainty requirement is small and sensor dark
noise is greater than $5\,e\text{-}$, this method of determining optimal
sample sizes works across the full dynamic range of a sensor. Additionally, we
have developed analytical expressions for the moments of the conversion gain
sampling distribution (the cif distribution) through the use of pseudomoments,
showing that these pseudomoments accurately describe the sampling moments of
conversion gain data under the proposed sensor noise model. With these
theoretical results, we were able to construct simple design and control of
experiment procedures that guides the number of samples required for both dark
and illuminated conditions based on iterative statistics and predicted
convergence. These experimental procedures were executed on a real image
sensor; the results of which agreed with our theoretical predictions and were
further confirmed through Monte Carlo simulation.

The ability to optimally measure per-pixel conversion gain is a key
development in a more comprehensive approach to per-pixel photon transfer
characterization. We have already shown how per-pixel conversion gain maps
enable per-pixel read noise estimation and we plan to extend this approach to
measure other important pt parameters such as well-capacity and dynamic range
on a per-pixel basis. The additional information that comes from these per-
pixel maps provides a richer characterization of the sensor and opens up the
idea of assigning quality metrics to a sensor based on the uniformity of its
per-pixel maps.

##  9 Proofs

###### Proof of Lemma 1.

By the law of total variance

$\mathsf{Var}T=(\mathsf{E}X^{2})\mathsf{Var}Y+(\mathsf{E}Y)^{2}\mathsf{Var}X.$
(95)

Substituting $\mathsf{E}X^{2}=\mathsf{Var}X+(\mathsf{E}X)^{2}$, expanding, and
dividing both sides by $(\mathsf{E}T)^{2}=(\mathsf{E}X)^{2}(\mathsf{E}Y)^{2}$
gives the desired result. ∎

###### Proof of Theorem 1.

Letting $T=G$, $X=\bar{P}$, and $Y=\hat{P}^{-1}$ we have after dividing both
sides of the relation in Lemma 1 by $\mathsf{ACV}^{2}\hat{P}^{-1}$ and
combining with

$\mathsf{ACV}^{2}\bar{P}=\frac{1}{n}\frac{\sigma_{P}^{2}}{\mu_{P}^{2}}=\frac{1}{n}\frac{1}{\mu_{P}g}=\frac{1}{n}\frac{1}{\mu_{e\text{-}}}$
(96)

and $\mathsf{ACV}^{2}\hat{P}^{-1}=2/(n-5)$:

$\frac{\mathsf{ACV}G}{\mathsf{ACV}\hat{P}^{-1}}=\left(1+\frac{n-3}{2n}\frac{1}{\mu_{e\text{-}}}\right)^{1/2}.$
(97)

The result then follows from noting that
$\sqrt{1+a/x}=1+a/(2x)+\mathcal{O}(x^{-2})$ as $x\to\infty$. ∎

###### Proof of Lemma 2.

The relation for $\sigma_{P}^{2}$ is derived by combining
$\sigma_{P}^{2}=\sigma_{P+D}^{2}-\sigma_{D}^{2}$ and
$\sigma_{P+D}^{2}=\sigma_{D}^{2}/\zeta$. To derive the relationship for
$\mu_{P}$ we use the fundamental gain relation in (30) to write
$\mu_{P}=\sigma_{P}^{2}g$. ∎

###### Proof of Lemma 3.

From (58) we observe $\mathsf{ACV}^{2}\bar{P}(N-n_{2},n_{2})$ is strictly
convex on $n_{2}\in(0,N)$ for all $\zeta\in(0,1)$; thus, the optimal sample
sizes for $\bar{P}$ satisfy

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

The first equation gives the optimality relation
$n_{2}^{\mathrm{opt}}/n_{1}^{\mathrm{opt}}=\sqrt{\zeta}$, which proves
$n_{2}^{\mathrm{opt}}/n_{1}^{\mathrm{opt}}\to 1$ as $\zeta\to 1^{-}$. Solving
the systems of equations then gives

$(n_{1}^{\mathrm{opt}},n_{2}^{\mathrm{opt}})=\left(\frac{1}{\sigma_{\mathscr{D}}^{2}\mathsf{acv}_{0}^{2}}\frac{\zeta(1+\sqrt{\zeta})}{(1-\zeta)^{2}},\frac{1}{\sigma_{\mathscr{D}}^{2}\mathsf{acv}_{0}^{2}}\frac{\zeta(\sqrt{\zeta}+\zeta)}{(1-\zeta)^{2}}\right).$
(99)

We now see that $n_{2}^{\mathrm{opt}}<n_{1}^{\mathrm{opt}}$, which allows us
to deduce the stronger result $n_{2}^{\mathrm{opt}}/n_{1}^{\mathrm{opt}}\to
1^{-}$ as $\zeta\to 1^{-}$. Then working directly with the expressions for the
optimal sample sizes it is straightforward to show $n_{i}^{\mathrm{opt}}\sim
2/(\sigma_{\mathscr{D}}^{2}\mathsf{acv}_{0}^{2})(1-\zeta)^{-2}$, which
completes the proof. ∎

###### Proof of Lemma 4.

Upon inspection of (59) we see
$\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta}(N-n_{2},n_{2})$ is strictly convex on
$n_{2}\in(1,N-1)$ for all $\zeta\in(0,1)$; thus, its optimal sample sizes
satisfy

$\displaystyle\partial_{n_{2}}\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta}(N-n_{2},n_{2})\Big{|}_{N=n_{1}^{\mathrm{opt}}+n_{2}^{\mathrm{opt}},\,n_{2}=n_{2}^{\mathrm{opt}}}$
$\displaystyle=0,$ (100)
$\displaystyle\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta}(n_{1}^{\mathrm{opt}},n_{2}^{\mathrm{opt}})$
$\displaystyle=\mathsf{acv}_{0}^{2}.$

Equating the derivative with zero we find for optimality relation
$1/(n_{1}^{\mathrm{opt}}-1)=\zeta/(n_{2}^{\mathrm{opt}}-1)$, which upon
substituting into the second equation gives us $n_{1}^{\mathrm{opt}}$.
Substituting the solution for $n_{1}^{\mathrm{opt}}$ back into the optimality
relation then gives $n_{2}^{\mathrm{opt}}$. The proof is now complete. ∎

######  Remark 2.

Suppose we define $\hat{P}^{-1}$ in terms of the optimal sample sizes for
$\hat{P}^{-1}_{\delta}$ and consider what happens as we pass to the shot-
noise-limit. Denoting

$\hat{P}^{-1}_{\mathrm{opt}}\coloneqq\frac{1}{\hat{X}(n_{1}^{\mathrm{opt}})-\hat{Y}(n_{2}^{\mathrm{opt}})},$
(101)

we see as $\zeta\to 0^{+}$ two simultaneous events occurring: (1)
$n_{2}^{\mathrm{opt}}\to 1$, which implies
$\hat{Y}(n_{2}^{\mathrm{opt}})\overset{d}{\to}\delta(0)$ and (2)
$\sigma_{P+D}\to\sigma_{P}$ so that
$\hat{X}\overset{d}{\to}\mathcal{G}(\alpha_{1},\alpha_{1}/\sigma_{P}^{2})$
with $\alpha_{1}=(n_{1}^{\mathrm{opt}}(0)-1)/2$. In conclusion,

$\lim_{\zeta\to
0^{+}}\hat{P}^{-1}_{\mathrm{opt}}\overset{d}{=}\frac{1}{\hat{X}(n_{1}^{\mathrm{opt}}(0))}$
(102)

with the r.h.s. being the shot-noise-limited estimator of $1/\sigma_{P}^{2}$
in Section 3 computed from $n_{1}^{\mathrm{opt}}(0)$ observations. Comparing
the shot-noise-limited optimal sample sizes for this estimator in (25) with
$n_{1}^{\mathrm{opt}}(0)=2/\mathsf{acv}_{0}+1$ we see that we can force our
optimal sample sizes to be exact in the shot-noise-limit by instead using

$(n_{1}^{\mathrm{opt}},n_{2}^{\mathrm{opt}})=\left(\frac{2(1+\zeta)}{\mathsf{acv}_{0}^{2}(1-\zeta)^{2}}+5,\frac{2\zeta(1+\zeta)}{\mathsf{acv}_{0}^{2}(1-\zeta)^{2}}+1\right).$
(103)

######  Lemma 5.

As $|z|\to\infty$

$\sqrt{2}\,\mathcal{D}(z/\sqrt{2})\sim\frac{1}{z}\sum_{k=0}^{\infty}(2k-1)!!\frac{1}{z^{2k}}.$
(104)

###### Proof of Lemma 5.

The proof follows from combining the relation

$\mathcal{D}(z)=\frac{\sqrt{\pi}}{2}e^{-z^{2}}\operatorname{erfi}(z)$ (105)

with the asymptotic expansion for $|z|\to\infty$

$\operatorname{erfi}(z)\sim\operatorname{sgn}(\Im
z)i+\frac{1}{\sqrt{\pi}\,z}e^{z^{2}}\sum_{k=0}^{\infty}\frac{(1/2)_{k}}{z^{2k}},$
(106)

and the relation between the Pochhammer symbol and double factorial
$(2k-1)!!=2^{k}(1/2)_{k}$. ∎

###### Proof of Theorem 2.

As $\mathsf{acv}_{0}\to 0^{+}$, $n_{1}^{\mathrm{opt}}\to\infty$ and
$n_{2}^{\mathrm{opt}}\to\infty$ so that we have by the central limit theorem
$\hat{P}_{\mathrm{opt}}\overset{d}{\to}\mathcal{N}(\mu_{\hat{P}_{\mathrm{opt}}},\sigma_{\hat{P}_{\mathrm{opt}}}^{2})$.
Noting that
$\mathsf{ACV}\hat{P}_{\mathrm{opt}}=\mathsf{ACV}\hat{P}^{-1}_{\delta,\mathrm{opt}}=\mathsf{acv}_{0}$
we have after combing (50) with (43)

$\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}\to\left|\frac{\sqrt{2}}{\mathsf{acv}_{0}}\mathcal{D}\left(\frac{1}{\sqrt{2}\,\mathsf{acv}_{0}}\right)-1\right|.$
(107)

Since $\mathsf{acv}_{0}$ is small, we subsequently have according to Lemma 5

$\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}\sim\mathsf{acv}_{0}^{2}\sum_{k=0}^{\infty}(2k+1)!!\,\mathsf{acv}_{0}^{2k},$
(108)

which leads to the desired result for
$\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}$.

By the same line of reasoning we combine (48) and (43)

$\mathsf{ACV}^{2}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}\to\frac{\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}}{(\sqrt{2}\,\mathcal{D}(\mathsf{acv}_{0}^{-1}/\sqrt{2}))^{2}}-1.$
(109)

Using standard results for the products and multiplicative inverses of power
series we write with the help of Lemma 5

$(\sqrt{2}\,\mathcal{D}(\mathsf{acv}_{0}^{-1}/\sqrt{2}))^{-2}\sim\frac{1}{\mathsf{acv}_{0}^{2}}\sum_{k=0}^{\infty}b_{k}\,\mathsf{acv}_{0}^{2k},$
(110)

with $b_{0}=1$, $b_{k}=-\sum_{i=1}^{k}a_{i}b_{k-i}$ and

$a_{i}=\sum_{j=0}^{i}(2(i-j)-1)!!\,(2j-1)!!.$ (111)

It follows

$\mathsf{ACV}^{2}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}\sim\sum_{k=1}^{\infty}\left(\sum_{\ell=0}^{k}(2(k-\ell)+1)!!\,b_{\ell}\right)\mathsf{acv}_{0}^{2k},$
(112)

which upon combining with
$\sqrt{1+ax^{2}+\mathcal{O}(x^{4})}=1+ax^{2}/2+\mathcal{O}(x^{4})$ as $x\to 0$
yields the desired asymptotic result for
$\mathsf{ACV}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}$. ∎

###### Proof of Theorem 3.

The proof follows much in the same way as that for Theorem 1. We write

$\frac{\mathsf{ACV}G_{\delta}}{\mathsf{ACV}\hat{P}^{-1}_{\delta}}=\left(1+\mathsf{ACV}^{2}\bar{P}+(\mathsf{ACV}^{2}\bar{P})(\mathsf{ACV}^{-2}\hat{P}^{-1}_{\delta})\right)^{1/2}.$
(113)

Combining this with the asymptotic approximations

$\mathsf{ACV}^{2}\bar{P}=\frac{1}{n_{1}\sigma_{\mathscr{D}}^{2}}\zeta+\mathcal{O}(\zeta^{2}),$
(114)
$\mathsf{ACV}^{-2}\hat{P}^{-1}_{\delta}=\frac{n_{1}-1}{2}+(n_{1}-1)\zeta+\mathcal{O}(\zeta^{2}),$
(115)

and $\sqrt{1+ax+\mathcal{O}(x^{2})}=1+ax/2+\mathcal{O}(x^{2})$ as $x\to 0$
gives the desired result. ∎

###### Proof of Theorem 4.

The results of Lemma 3 and Proposition 1 show that as $\zeta\to 1^{-}$, the
optimal sample sizes for $\hat{P}^{-1}_{\delta}$ and $\bar{P}$ are
asymptotically equal and of the form $n_{i}^{\mathrm{opt}}\sim
C(1-\zeta)^{-2}$. As such, the optimal sample sizes for $G_{\delta}$ must also
be asymptotically equal and of the form $n_{i}^{\mathrm{opt}}\sim
C_{G_{\delta}}(1-\zeta)^{-2}$ since

$\mathsf{ACV}^{2}G_{\delta}=\mathsf{ACV}^{2}\bar{P}+(\mathsf{ACV}^{2}\bar{P})(\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta})+\mathsf{ACV}^{2}\hat{P}^{-1}_{\delta}.$
(116)

To determine $C_{G_{\delta}}$ we substitute
$n_{1}=n_{2}=C_{G_{\delta}}(1-\zeta)^{-2}$ into the expression for
$\mathsf{ACV}^{2}G_{\delta}$, pass to the limit, and equate with
$\mathsf{acv}_{0}$ yielding

$\frac{2}{\sigma_{\mathscr{D}}^{2}}\frac{1}{C_{G_{\delta}}}+\frac{8}{\sigma_{\mathscr{D}}^{2}}\frac{1}{C_{G_{\delta}}^{2}}+4\frac{1}{C_{G_{\delta}}}=\mathsf{acv}_{0}.$
(117)

The resulting expression is a quadratic equation in $C_{G_{\delta}}$, the
roots of which are real and differing in sign. Taking the positive root then
yields the expression for $C_{G_{\delta}}$. ∎

###### Proof of Theorem 5.

Since $\mathsf{ARB}_{\mathcal{P}}G=\mathsf{ARB}_{\mathcal{P}}\hat{P}^{-1}$,
the result for $\mathsf{ARB}_{\mathcal{P}}G_{\mathrm{opt}}$ immediately
follows from Theorem 2. To obtain the result for
$\mathsf{ACV}_{\mathcal{P}}G_{\mathrm{opt}}$, combine

$\mathsf{ACV}_{\mathcal{P}}^{2}G=\mathsf{ACV}_{\mathcal{P}}^{2}\hat{P}^{-1}+(\mathsf{ACV}_{\mathcal{P}}^{2}\hat{P}^{-1})(\mathsf{ACV}^{2}\bar{P})+\mathsf{ACV}^{2}\bar{P}$
(118)

with

$\displaystyle\mathsf{ACV}^{2}_{\mathcal{P}}\hat{P}^{-1}_{\mathrm{opt}}$
$\displaystyle=\mathsf{acv}_{0}+6\mathsf{acv}_{0}^{3}+\mathcal{O}(\mathsf{acv}_{0}^{5})$
(119) $\displaystyle\mathsf{ACV}^{2}_{\mathcal{P}}\bar{P}_{\mathrm{opt}}$
$\displaystyle=a\mathsf{acv}_{0}+b\mathsf{acv}_{0}^{3}+\mathcal{O}(\mathsf{acv}_{0}^{5})$

and $\sqrt{1+ax^{2}+\mathcal{O}(x^{4})}=1+ax^{2}/2+\mathcal{O}(x^{4})$ as
$x\to 0$. ∎

##  10 Backmatter

Funding The authors report no funding for this work.

Acknowledgments The authors would like thank Nico Schlömer for his matlab2tikz
function used to create the figures throughout this work [34]. The authors
would also like to thank Paul Enta for pointing us to reference [15] on square
series generating functions.

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

Supplemental document Code for the doe and coe procedures is available on
figshare (Ref. [30, 31]).

## References

  * [1] J. R. Janesick, _Photon Transfer: $DN\to\lambda$_ (SPIE, 2007).
  * [2] EMVA 1288 working group, “EMVA Standard 1288: Standard for characterization and presentation of specification data for image sensors and cameras, release A1.00,” (2005).
  * [3] EMVA 1288 working group, “EMVA Standard 1288: Standard for characterization of image sensors and cameras, release 4.0 linear,” (2021).
  * [4] B. P. Beecken and E. R. Fossum, “Determination of the conversion gain and the accuracy of its measurement for detector elements and arrays,” Appl. Opt. 35, 3471–3477 (1996).
  * [5] J. R. Janesick, _Scientific Charge-Coupled Devices_ , Press Monographs (SPIE, 2001).
  * [6] B. Pain and B. R. Hancock, “Accurate estimation of conversion gain and quantum efficiency in CMOS imagers,” in _Sensors and Camera Systems for Scientific, Industrial, and Digital Photography Applications IV,_ (2003), pp. 94–103.
  * [7] J. Janesick, J. T. Andrews, and T. Elliott, “Fundamental performance differences between CMOS and CCD imagers: Part 1,” in _High Energy, Optical, and Infrared Detectors for Astronomy II,_ vol. 6276 D. A. Dorn and A. D. Holland, eds., International Society for Optics and Photonics (SPIE, 2006), pp. 208–226.
  * [8] S. E. Bohndiek, A. Blue, A. T. Clark, M. L. Prydderch, R. Turchetta, G. J. Royle, and R. D. Speller, “Comparison of methods for estimating the conversion gain of CMOS active pixel sensors,” IEEE Sensors Journal 8, 1734–1744 (2008).
  * [9] D. A. Starkey and E. R. Fossum, “Determining conversion gain and read noise using a photon-counting histogram method for deep sub-electron read noise image sensors,” IEEE Journal of the Electron Devices Society 4, 129–135 (2016).
  * [10] A. J. Hendrickson, “Centralized inverse-Fano distribution for controlling conversion gain measurement accuracy of detector elements,” J. Opt. Soc. Am. A 34, 1411–1423 (2017).
  * [11] A. Hendrickson, “The inverse gamma-difference distribution and its first moment in the Cauchy principal value sense,” Statistics and Its Interface 12, 467–478 (2019).
  * [12] A. Hendrickson, “A novel approach to photon transfer conversion gain estimation,” Cornell University arXiv 2106.14958 (2021).
  * [13] B. L. Preece and D. P. Haefner, “3d noise photon transfer curve,” Appl. Opt. 61, 6202–6212 (2022).
  * [14] K. Nakamoto and H. Hotaka, “Efficient and accurate conversion-gain estimation of a photon-counting image sensor based on the maximum likelihood estimation,” Opt. Express 30, 37493–37506 (2022).
  * [15] M. D. Schmidt, “Square series generating function transformations,” Journal of Inequalities and Special Functions 8, 125–156 (2017).
  * [16] G. Casella and R. Berger, _Statistical Inference_ , Duxbury advanced series in statistics and decision sciences (Thomson Learning, 2002), 2nd ed.
  * [17] A. Mathai, “On noncentral generalized Laplacianness of quadratic forms in normal variables,” Journal of Multivariate Analysis 45, 239–246 (1993).
  * [18] E. Krishna and K. Jose, “Marshall-Olkin generalized asymmetric Laplace distributions and processes,” Statistica 71, 453–467 (2011).
  * [19] B. Klar, “A note on gamma difference distributions,” Journal of Statistical Computation and Simulation 85, 3708–3715 (2015).
  * [20] M. Hančová, A. Gajdoš, and J. Hanč, “A practical, effective calculation of gamma difference distributions with open data science tools,” Journal of Statistical Computation and Simulation 92, 2205–2232 (2022).
  * [21] E. L. Lehmann and J. P. Shaffer, “Inverted distributions,” The American Statistician 42, 191–194 (1988).
  * [22] C. Y. Peng, “The first negative moment in the sense of the Cauchy principal value,” Statistics & Probability Letters 78, 1765–1774 (2008).
  * [23] C. Y. Peng, “The first negative moment of skew-t and generalized Student’s t-distributions in the principal value sense,” Journal of Applied Mathematics 2013 (2013).
  * [24] C. Fox, “A generalization of the Cauchy principal value,” Canadian Journal of Mathematics 9, 110–117 (1957).
  * [25] E. A. Galapon, “The Cauchy principal value and the Hadamard finite part integral as values of absolutely convergent integrals,” Journal of Mathematical Physics 57, 033502 (2016).
  * [26] G. Criscuolo, “A new algorithm for Cauchy principal value and Hadamard finite-part integrals,” Journal of Computational and Applied Mathematics 78, 255–275 (1997).
  * [27] M. H. Quenouille, “Notes on bias in estimation,” Biometrika 43, 353–360 (1956).
  * [28] R. Barakat, “The derivatives of Dawson’s function,” Journal of Quantitative Spectroscopy and Radiative Transfer 11, 1729–1730 (1971).
  * [29] ON Semiconductor, _KAI-04070: 2048 (H) x 2048 (V) Interline CCD Image Sensor_ (2013). Rev. 3.
  * [30] A. Hendrickson, D. P. Haefner, and B. L. Preece, “Supplementary material 1: On the optimal measurement of conversion gain in the presence of dark noise.” figshare (2022). https://opticapublishing.figshare.com/s/e7ce99e8995a5581cdb9.
  * [31] A. Hendrickson, D. P. Haefner, and B. L. Preece, “Supplementary material 2: On the optimal measurement of conversion gain in the presence of dark noise.” figshare (2022). https://opticapublishing.figshare.com/s/96ea4554dca4c769a04d.
  * [32] Labsphere, “Integrating sphere theory and applications,” (2017).
  * [33] D. V. Hinkley, “On the ratio of two correlated normal random variables,” Biometrika 56, 635–639 (1969).
  * [34] N. Schlömer, “matlab2tikz: A script to convert MATLAB/Octave into TikZ figures for easy and consistent inclusion into LaTeX.” GitHub. URL: https://github.com/matlab2tikz/matlab2tikz (retrieved May 8, 2021).