%!TEX root = ../dissertation.tex
%\begin{savequote}[75mm]
%Nulla facilisi. In vel sem. Morbi id urna in diam dignissim feugiat. Proin molestie tortor eu velit. Aliquam erat volutpat. Nullam ultrices, diam tempus vulputate egestas, eros pede varius leo.
%\qauthor{Quoteauthor Lastname}
%\end{savequote}

\chapter{A Bayesian approach to analyzing holograms of colloidal particles}
\label{chap:bayesian}

\newthought{Digital in-line holography} has emerged as a powerful tool for both
tracking and characterizing colloidal particles. Such particles are a
few hundred nanometers to a few micrometers in diameter and dispersed in
a fluid. In a typical tracking experiment, one measures the position of
the particle in three dimensions as a function of time by analyzing a
time-series of holograms; the particle might be subject to Brownian
motion, interactions, or external forces. This type of measurement is
important to a number of different experiments in soft-matter physics
and fluid mechanics, including
microrheology~\cite{cheong_holographic_2009}, the study of phase
transitions~\cite{perry_real-space_2012}, measuring interactions in
colloidal systems~\cite{fung_measuring_2011, kaz_physical_2012}, and
particle image velocimetry~\cite{ovryn_three-dimensional_2000,
  cheong_flow_2009}. In a characterization experiment, one measures the
size, refractive index, and shape of individual colloidal particles.
This type of experiment is important for quality control of
formulations~\cite{wang_holographic_2016} and studies of polymerization
kinetics~\cite{wang_celebrating_2015}, among other applications.

Many of these experiments have been enabled by a method of analyzing
holograms based on scattering theory. A hologram is a two-dimensional
interference pattern that encodes information about the
three-dimensional (3D) positions and optical properties of a particle or
set of particles. An inline hologram results from the interference of
the scattered fields from the sample and the undiffracted beam, which
acts as a reference wave (Figure~\ref{fig:setup}). Traditionally, the 3D
information is recovered through reconstruction---effectively, shining
light back through the hologram to generate a real 3D image. However,
when the particles are comparable to the wavelength, reconstruction
produces a number of imaging artifacts \cite{pu_intrinsic_2003}. An
alternative to reconstruction was first demonstrated by Ovryn and
Izen~\cite{ovryn_imaging_2000, ovryn_three-dimensional_2000}, who showed
that inline holograms could be directly predicted from Lorenz--Mie
scattering theory and a model for the propagation and interference of
the fields in the optical train of the microscope. Later work by Lee and
coworkers~\cite{lee_characterizing_2007}, using a simpler model of the
optical train, showed that one could fit a scattering model based on
Lorenz--Mie theory to both track and characterize spherical colloidal
particles. Further work from our research group has shown that other
models, based on different scattering theories, can be fit to holograms
of non-spherical particles\cite{wang_using_2014} and clusters of
spherical particles~\cite{perry_real-space_2012, fung_imaging_2012}.

\begin{figure}
  \centering
  \includegraphics[width=5in]{figures/bayesian-holography/diagram}
  \caption{A typical optical setup used to record inline digital
    holograms. We image with a collimated laser beam and record the
    hologram on a CMOS camera. The coordinate $z$ is defined as the
    distance from the particle to the focal plane of the objective
    ($z>0$ if the particle is above (upstream of) the focal plane). $x$
    and $y$, not shown, represent the position of the particle in the
    plane normal to the optical axis. On the right, a measured,
    background-divided hologram of a 1-\si{\micro\meter}-diameter silica particle
    in water, taken using a 660 nm laser (hologram courtesy of Anna
    Wang).}
  \label{fig:setup}
\end{figure}

With these techniques, the optical field is never reconstructed;
instead, a forward model---a model that can calculate a hologram based
on scattering theory and propagation along the optical train--- is fit
to the data. In this method, which we will refer to as ``least-squares
fitting,'' an iterative minimization algorithm is used to find the
parameters in the model that minimize the sum of the squared
differences, pixel-by-pixel, between the measured hologram and the
model. Although this process is computationally more expensive than
reconstruction, it yields much more accurate and precise measurements:
for tracking experiments, analysis of the mean-square displacement shows
that the uncertainty in the position can be on the order of a few
nanometers or less \cite{kaz_physical_2012}.

However, there are two problems with the least-squares technique. First,
it is difficult to obtain accurate uncertainties on parameters. In a
characterization experiment, for example, one would like to know the
uncertainty in the measured size or refractive index. The uncertainties
on the fit parameters are critical not only for characterization
experiments, but for any experiment in which the data will later be
compared to a physical theory---for example an equation of motion for
the particles---since the agreement with theory can only be judged when
uncertainties are known. Second, the iterative minimization schemes used
in least-squares fitting require good initial guesses for and
constraints on all the parameters~\cite{fung_measuring_2011,
  perry_real-space_2012, wang_using_2014}. Determining the initial
guesses and the constraints often requires a great deal of effort,
particularly for more complex scattering models.

Here we demonstrate a Bayesian analysis method that overcomes these
problems. Bayesian approaches to parameter estimation differ from
classical frequentist approaches such as least-squares fitting in that
we are allowed to speak of the probability distribution \emph{of a
  parameter}. As in the least-squares technique, we start with a
forward model, but instead of maximizing the likelihood, which is the
probability of the data given parameter values, we maximize the
posterior probability, which is the probability of the parameter values
given the data.  There are other Bayesian approaches to analyzing
holograms \cite{denis_inline_2009}, but none that we are aware of use
forward models based on scattering solutions.

There are two principle advantages of the Bayesian approach over
least-squares fitting: first, any prior information about the
particle, such as its position, size, or composition, can be
formalized as a prior probability distribution, rather than implicitly
encoded into initial guesses or constraints. Second, because
probabilities are assigned to parameters, we can marginalize, or
integrate out, the uncertainties on some of the parameters to obtain
the uncertainties on the others. For example, in the characterization
problem, we are interested in the uncertainty on the size or
refractive index after taking into account the uncertainties on the
positions. This type of analysis is not easily done using classical
frequentist methods, and hence the parameters that are not of interest
(but required by the model) are often called ``nuisance parameters.''

As an example of the utility of this technique, we show the results from
a Bayesian analysis of a hologram of a 1-\si{\micro\meter}-diameter silica particle
in water in Figure~\ref{fig:comp_mcmc}, compared to the results obtained
from least-squares fitting. The Bayesian approach yields nearly the same
``best fit'' values for the parameters, but it also gives us much more
information about the parameter uncertainties. From the Bayesian
results, we see that there are strong correlations between the estimated
particle refractive index, size, and axial position. The covariances
between these parameters are difficult to extract in a least-squares
approach. Because we have access to these covariances in the Bayesian
approach, we can integrate them out, or ``marginalize,'' to estimate the
uncertainties on any one parameter, as shown by the distributions along
the diagonal.

\begin{figure}
  \centering
  \includegraphics[width=5in]{figures/bayesian-holography/mle-comparison}
  \caption{Comparison of results obtained from least-squares fitting
    using the Levenberg--Marquardt method, as described in
    Section~\ref{sec:bh:background}, to those obtained from the Bayesian
    inference method described in this paper. Units of $x$, $y$, $z$ and
    $r$ are micrometers. The results from fitting are shown by the
    dashed orange lines. Fully marginalized distributions obtained from
    Markov-chain Monte Carlo (MCMC) sampling are shown along the
    diagonal, while the off-diagonal contour plots show the joint
    distributions of each pair of parameters, represented as Gaussian
    kernel density estimates. We explain the Bayesian approach and the
    MCMC technique later in the text.}
  \label{fig:comp_mcmc}
\end{figure}

The rest of this paper explains how we obtain these results and how we
interpret them. In practice, the Bayesian approach involves three steps:
specifying a forward model, specifying priors on all the parameters, and
calculating the posterior probability distribution of all the
parameters. The last step, calculating the posterior, is most
effectively done using Markov-chain Monte Carlo (MCMC) techniques. In
what follows, we limit our discussion to the measurement of a single
microsphere, and we use a simple model for the propagation and
interference along the optical train, equivalent to that of Lee and
coworkers\cite{lee_characterizing_2007}. Nonetheless, the approach can
easily be extended to more complex models of scattering or propagation.
After describing the model and the Bayesian approach to calculating the
posterior probability distribution, we apply our method to both the
particle tracking problem and the particle characterization problem. We
also show how to speed up the approach---and obviate the need for
initial starting guesses for our MCMC chains---by a type of tempering
based on choosing random subsets of pixels.

\section{Background}
\label{sec:bh:background}

Here we describe the forward model for hologram formation from a single
spherical colloidal particle and the least-squares method used to fit
this model to the data. The data are obtained from an in-line
holographic setup as shown in Figure \ref{fig:setup}. We follow the work
of Lee and coworkers\cite{lee_characterizing_2007} closely. In their
model, a hologram $\mathbf{H}$ is computed as a function of six
parameters:
%
\begin{equation}
  \mathbf{H}(x, y, z, n, r, \alpha) =
  \left\lvert \alpha E_\mathrm{scat}(x, y, z, r, n) + E_\mathrm{ref}\right\rvert^2
\end{equation}
%
where $x$, $y$, and $z$ are the coordinates of the center of the
particle, $n$ is its (possibly complex) refractive index, $r$ is its
radius, and $\alpha$ is an auxiliary parameter that rescales the
scattered field $E_\mathrm{scat}$. $E_\mathrm{ref}$ is the reference
field. We use boldface type for the hologram $\mathbf{H}$ to
indicate that it is a matrix consisting of measurements of the intensity
at each pixel $(i,j)$ on our detector. The parameter $\alpha$ is needed
to compensate for a number of different effects along the optical train
and to produce satisfactory agreement between the model and measured
holograms. Although a single scalar parameter is likely insufficient to
account for all of these effects, the Bayesian analysis framework we
develop below can easily be extended to more complex models, including
ones that model the propagation along the optical train.

We divide the measured hologram $\mathbf{H}$ by a measured background,
usually from a hologram captured without any particles in the field of
view. The measured background should be very close to the intensity of
the reference wave, $\left\lvert E_\mathrm{ref}\right\rvert^2$,
differing only by the light scattered by the particle of interest. We
therefore approximate the background as the reference intensity to
obtain an expression for the normalized hologram $\mathbf{h}$:
%
\begin{equation}
  \mathbf{h} = \frac{\mathbf{H}(x, y, z, r, n, \alpha)}
                    {\left\lvert E_\mathrm{ref} \right\rvert^2} =
  \left\lvert\frac{\alpha E_\mathrm{scat}(x, y, z, r, n)}
                  {E_\mathrm{ref}} + 1\right\rvert^2.
  \label{eq:forward_model}
\end{equation}

To infer the parameters $x$, $y$, $z$, $r$, $n$, and $\alpha$, we need a
model for $E_\mathrm{scat}$ and an inference method. For a single
sphere, the Lorenz--Mie theory can be used to calculate
$E_\mathrm{scat}$; other scattering theories\cite{fung_measuring_2011,
  fung_imaging_2012, wang_using_2014} can be used for different types of
particles or clusters of particles. In the inference method used by Lee
and coworkers\cite{lee_characterizing_2007} and found in many subsequent
papers, the parameters are inferred by minimizing the sum (over pixels)
of the squared differences between model and data. Letting $\theta$ be
the set of parameters $\left\{x,y,z,n,r,\alpha\right\}$, we have
%
\begin{equation}
  \label{eq:chisquared_min}
  \theta = \argmin_\theta\chi^2 =
  \argmin_\theta \sum_{i, j} \left[h_{ij} - h^M_{ij}(\theta)\right]^2
\end{equation}
%
where $h_{ij}$ is the intensity measured at pixel $(i,j)$ in the
normalized recorded hologram and $h^M_{ij}$ is the same for the hologram
calculated from the forward model $M$. Using this least-squares
formulation is equivalent to maximum-likelihood estimation when the
noise in the measurements has a Gaussian distribution.

Because evaluating the Mie fields involves many computationally
expensive special-function evaluations over a large number of pixels,
several alternatives and variations on this method have been developed.
Soulez and coworkers \cite{soulez_inverse-problem_2007} used a local
refinement approach to find particle positions more efficiently. Yevick
and coworkers~\cite{yevick_machine-learning_2014} used a
machine-learning method employing a trained support-vector machine,
instead of the full Lorenz--Mie solution, to infer the parameters. It is
also possible to obtain particle positions from reconstructed volumes
using deconvolution~\cite{latychevskaia_depth-resolved_2010} or
successive resampling~\cite{seifi_three-dimensional_2012}. In a recent
paper our research group showed that good fits could be obtained by
modeling only a small, randomly chosen subset of the
pixels\cite{dimiduk_random-subset_2014}. We later discuss how the
random-subset technique can speed up the Bayesian approach as well.

Least-squares fitting is generally done using an iterative non-linear
algorithm such as Levenberg--Marquardt~\cite{marquardt_algorithm_1963},
which requires the user to specify an initial guess for the parameters.
Good initial guesses are required for the algorithm to converge. To keep
the algorithm from wandering off into unphysical regions of parameter
space, one can fix certain parameters or allow them to vary only over a
certain range. One can also add arbitrary functions to the $\chi^2$
function in Eq.~(\ref{eq:chisquared_min}) to penalize certain
regions of parameter space. However, fixing parameters might not
accurately reflect our knowledge of the system. We might have some prior
information about the particle size (for example, from the manufacturer)
but not trust it enough to set the size to a certain value. A penalty
function can more accurately capture our state of knowledge, but as we
shall see, Bayesian priors offer an elegant alternative to these various
methods of limiting the parameter space.

\section{Bayesian Approach}

The problem of estimating parameters from a measured hologram and a
forward model is well suited to Bayesian inference, which allows us to
replace guesses, penalties, and constraints with a single prior
probability function. Here we describe our approach using standard
Bayesian terminology; for a more detailed discussion of Bayesian
inference and incorporation of priors, see chapter 3 of Gregory's book on Bayesian data analysis~\cite{gregory_bayesian_2005}.

We are interested in the posterior probability $p(\theta| \mathbf{h}, M,
I)$ of a set of parameters $\theta$ given the recorded hologram
$\mathbf{h}$, the forward model $M$, and known prior information $I$.
$I$ might include information about the apparatus, such as the
magnification, wavelength, and observation volume, and the fact that
we know we are looking at a single sphere. Applying Bayes' rule, we can
express the posterior in terms of a prior probability density
$p(\theta|M, I)$, the likelihood function $p(\mathbf{h}|\theta, M, I)$,
and a normalization factor $p(\mathbf{h}|M, I)$:
\begin{equation}
  p(\theta|\mathbf{h}, M, I) =
  \frac{p(\theta|M, I) p(\mathbf{h}|\theta, M, I)}{p(\mathbf{h}|M, I)}.
\end{equation}

Because the normalization factor $p(\mathbf{h}|M,I)$ has no effect on
the parameters we measure or their uncertainties, we work with the
unnormalized posterior probability distribution:
%
\begin{eqnarray}
  p'(\theta|M, \mathbf{h}, I) &=& p(\theta|\mathbf{h}, M, I) p(\mathbf{h}|M, I) \nonumber\\
  &=& p(\theta|M, I) p(\mathbf{h}|\theta, M, I). \label{eq:unnormalized}
\end{eqnarray}
%
As shown by Eq.~(\ref{eq:unnormalized}), we need to specify the prior
$p(\theta|M, I)$ and the likelihood function $p(\mathbf{h}|\theta, M,
I)$ to proceed.

\subsection{Priors}
\label{sec:priors}

The prior $p(\theta|M, I)$ is simply a mathematical statement of what we
know about the parameters before doing the inference calculation. For
example, priors for $r$ and $n$ might come from manufacturer's data
sheets for commercially sourced particles. If we know the material
composition of the particle, we can obtain a prior for $n$ from a
refractive index table and some estimated uncertainty. Similarly, we
might obtain a prior probability distribution for $r$ from the size
distribution of the particles, as measured by the manufacturer or by
some other technique, such as light scattering. In general, the prior
distribution reflects our state of knowledge about a parameter, with
wider distributions indicating greater uncertainty. We discuss specific
choices of priors in Section~\ref{sec:computation}.

\subsection{Likelihood function}
\label{sec:likelihood}

To evaluate the likelihood $p(\mathbf{h}| \theta, M, I)$, we must model
the noise in the measured holograms. There are a variety of noise
sources, including photon shot noise, electronic readout noise at the
camera, and stray fringes from other scatterers in the optical train. We
consider only the shot noise and readout noise in what follows. We then
write the measured hologram as the sum of the model hologram plus noise:
\begin{equation}
  h_{ij} = h^M_{ij} + u_{ij}
\end{equation}
where the model hologram $\mathbf{h}^M$ (the elements of which are
$h^M_{ij}$) is calculated from the forward model
(Eq.~(\ref{eq:forward_model})) and the Lorenz--Mie theory. The noise
in the intensity at pixel $(i,j)$ is represented by $u_{ij}$.

In writing this equation, we are assuming that the noise can be
represented by a single additive term that characterizes pixel-by-pixel
variations in the intensity. We now further assume that the noise
$\mathbf{u}$ is delta correlated in space and time and that each
$u_{ij}$ is Gaussian distributed with a variance $\sigma^2_{ij}$ at each
pixel. With these assumptions, we can write the likelihood as
%
\begin{equation}
  p(\mathbf{h}|\theta, M, I) = \prod_{i, j} p(u_{ij}|\theta, M, I)
  \label{eq:likelihood1}
\end{equation}
%
where
\begin{equation}
  p(u_{ij}|\theta, M, I) = \frac{1}{\sqrt{2\pi}\sigma_{ij}} \exp\left\{\frac{-\left[h_{ij}-h^M_{ij}(\theta)\right]^2}{2\sigma^2_{ij}}\right\}.
  \label{eq:likelihood2}
\end{equation}
%
Combining Eq.~(\ref{eq:likelihood1}) with
Eq.~(\ref{eq:likelihood2}), we obtain the full likelihood
%
\begin{equation}
  \label{eq:likelihood}
  p(\mathbf{h}|\theta, M, I) = \frac{1}{(2\pi)^{N/2}\prod_{i, j} \sigma_{ij}}\exp\left\{-\sum_{i, j}\frac{\left[h_{ij}-h^M_{ij}(\theta)\right]^2}{2\sigma^2_{ij}}\right\}
\end{equation}
%
where $N$ is the total number of pixels in the hologram. We can estimate
$\sigma^2_{ij}$ from a series of background holograms where no particles
are present.

To justify this form of the likelihood, we first note that at the
intensities we measure, the number of photons and readout electrons are
large enough that the both shot noise and readout noise are
approximately Gaussian, and hence their sum is also Gaussian. Also, the
non-linear least squares fitting methods used in previous work
implicitly assume that the noise is Gaussian, so our approach is
consistent with these methods. Other noise models that include the
effect of correlations (that is, models for speckle) are a topic for
future work.

\subsection{Putting it all together}
\label{sec:posterior}

We recognize the term in the exponent of the likelihood in
Eq.~(\ref{eq:likelihood}) as $\chi^2/2$, where
\begin{equation}
  \chi^2(\theta) =
  \sum_{i, j}  \frac{\left[h_{ij}-h^M_{ij}(\theta)\right]^2}{\sigma^2_{ij}}.
\end{equation}
%
In this work we assume that the standard deviation of the noise at each
pixel is identical and equal to a constant $\sigma$, so
that
%
\begin{equation}
  \prod_{i, j} \sigma_{ij} = \sigma^{N}.
\end{equation}
%
Substituting the above into Eq.~(\ref{eq:unnormalized}) yields an
equation for the unnormalized posterior:
%
\begin{equation}
  \label{eq:posterior}
  p'(\theta|, \mathbf{h}, M, I) =
  \frac{p(\theta|M, I)}{(2\pi)^{N/2}\sigma^{N}}
  \exp\left[-\frac{\chi^2(\theta)}{2}\right].
\end{equation}
%
Because the $\chi^2$ term depends on the parameters $\theta$ through the
forward model, it is difficult if not impossible to calculate the
posterior directly from Eq.~(\ref{eq:posterior}). We therefore use
a Markov-chain Monte Carlo technique as described in
Section~\ref{sec:computation}.

For the specific problems of particle tracking and particle
characterization, we are interested only in subsets of the parameters.
In particle tracking, we want to infer the position as a function of
time. The position is described by one or more of the parameters $x$,
$y$, and $z$ along with their uncertainties, which can be obtained by
marginalizing over $n$, $r$, and $\alpha$ (and over the position
variables that are not of interest). In particle characterization, we
want to infer the index of refraction $n$ and radius $r$ of the particle
along with their associated uncertainties, which can be obtained by
marginalizing over $x$, $y$, $z$, and $\alpha$. Formally,
marginalization involves integration of the posterior distribution over
the parameters to be marginalized. The marginalized distribution for the
particle tracking problem is
\begin{equation}
  p'(x, y, z|\mathbf{h}, M, I) = \int dr\, dn\, d\alpha\, p'(x, y, z, r, n, \alpha|\mathbf{h}, M, I)
\end{equation}
and that for the particle characterization problem is
\begin{equation}
  p'(n, r| \mathbf{h}, M, I) = \int dx\, dy\, dz\, d\alpha\, p'(x, y, z, n, r,
  \alpha| \mathbf{h}, M, I).
\end{equation}
In practice, we do not have to evaluate these integrals. As we describe
in Section~\ref{sec:computation}, we can sample the posterior using MCMC
techniques and then marginalize by binning samples or by kernel density
estimation.

\subsection{Time-Series}

In both particle tracking and particle characterization we work with
time-series of holograms. There are two different categories of
parameters within time-series: those that change over time and those
that do not. In most cases, the particle radius $r$ and refractive index
$n$ do not change as a function of time, whereas $x$, $y$, and $z$ do.
Therefore when calculating the marginalized posterior for the
characterization problem, $p'(n, r| \mathbf{h}, M, I)$, we can use
information from the entire time-series to minimize the uncertainty on
the parameters. When computing $p'(x, y, z | \mathbf{h}, M, I)$ we
analyze each frame individually.

There are at least two procedures to calculate the posterior $p'(n, r|
\mathbf{h}, M, I)$ using information from the entire time-series. One is
to evaluate it for a single frame and then use the result as the prior
for the next frame. This procedure is similar to using a Kalman filter.
It should yield the best estimate of $n$ and $r$ when the last frame is
analyzed. The other procedure is to do a global calculation over the
entire time-series. This procedure would involve expanding the model to
include parameters for $x$, $y$, and $z$ in each frame, then
marginalizing all of these position parameters. Here we use the first
approach because it is simpler to implement and computationally less
demanding.

\section{Results and discussion}

\subsection{Markov-chain Monte Carlo technique}
\label{sec:computation}

Because Eq.~(\ref{eq:posterior}) does not admit an analytic
solution, we must calculate the posterior probability distribution
numerically. Brute force evaluation of the posterior is prohibitively
expensive owing to the six parameters. Thus we use a sampling technique
based on a Markov-chain Monte Carlo algorithm.

Specifically, we use an affine-invariant ensemble
sampler\cite{goodman_ensemble_2010} as implemented in the Python library
Emcee \cite{foreman-mackey_emcee:_2013}. This sampler uses an ensemble of
``walkers'' to explore the posterior probability distribution. After a
sufficient number of steps, the ensemble will converge to a steady-state
distribution that is, by construction, equal to the posterior
probability density. Thus, in the long-time limit, the method yields a
set of samples directly from the posterior. Each sample consists of
values for all six parameters in our forward model. We can plot these
samples as in Figure~\ref{fig:comp_mcmc} to visualize the joint
distributions between the parameters, or we can marginalize simply by
binning the samples as a function of only the parameters we wish to
infer.

We demonstrate these techniques on holograms of 1-\si{\micro\meter}-diameter
silica spheres in water (index 1.333), illuminated with a 660 nm
laser. We assign priors as follows. For the refractive index $n$, we
choose a Gaussian prior with a mean of 1.5 and a standard deviation of
0.1, chosen based on the typical variation of refractive index from
particle to particle measured in holographic microscopy. Based on
manufacturer's specifications, we choose a Gaussian prior for the
particle radius with mean of 0.5~\si{\micro\meter} and standard deviation of 0.05
\si{\micro\meter}. We estimate the in-plane particle position $(x,y)$ using a
Hough-transform based algorithm \cite{cheong_flow_2009}, and we assign
Gaussian priors with standard deviations of 0.1~\si{\micro\meter}, chosen based
on our prior experience with this algorithm. For $\alpha$, we choose a
Gaussian prior centered at 0.7 with width 0.05, based on prior
experience from least-squares fitting. We have also used a uniform
prior ranging from 0 to 1, with no significant effect on the results.

Choosing a prior distribution for $z$ is the most difficult task; our
knowledge of $z$ is limited before we do the full inference calculation
(indeed, one of the main goals of the model-based approach to holography
is to make it possible to extract precise information about $z$). We
discuss various choices for the prior on $z$ below. In the end, we are
able to use a uniform prior from 0 to 100~\si{\micro\meter} from the focal plane.
We chose this width based on the sample height, which is on the order of
100~\si{\micro\meter}. The prior can be widened or narrowed with little effect on
the final results.

We compute the noise ($\sigma$ in Eq.~\ref{eq:posterior}) from a
background image $\mathbf{b}$ as
%
\begin{equation}
  \sigma = \frac{\mathrm{std}(\mathbf{b})}{\mathrm{mean}(\mathbf{b})}.
\end{equation}
Typically, $\sigma\approx0.1$ for holograms with no other particles
present. For the results presented here, $\sigma=0.119$.

\subsubsection{Equilibration of the MCMC chain}
\label{sec:burnin}

With these choices for priors, we run our MCMC algorithm and examine
the period of time required to reach steady state, called the
``burn-in'' period. Evaluating the burn-in period is important because
the samples obtained within it are not representative of the posterior
and therefore must be discarded.

We find that the number of samples needed to reach steady state is
strongly affected by the initial positions of the ensemble of walkers.
Figure \ref{fig:trace} shows that if the initial positions are not too
far from the maximum \textit{a posteriori} (MAP) estimate of $z$, the
burn-in is quick. However, if the initial positions are slightly further
from the MAP estimate, the MCMC walkers can get stuck in local minima of
the posterior distribution and take a long time to equilibrate
(Figure~\ref{fig:stuck}).

\begin{figure}
  \centering
  \begin{subfigure}[b]{2.5in}
    \includegraphics{figures/bayesian-holography/trace_5um_uniform}
    \caption{Walkers initially distributed uniformly between 0 and 5~\si{\micro\meter}}
    \label{fig:trace}
  \end{subfigure}
  \begin{subfigure}[b]{2.5in}
    \includegraphics{figures/bayesian-holography/trace_8um_uniform}
    \caption{Walkers initially distributed uniformly between 0 and 8~\si{\micro\meter}}
    \label{fig:stuck}
  \end{subfigure}
  \caption{Plots of walker position as a function of MCMC time step. As can
    be seen from the plot on the right, walkers that are further away
    from the MAP position get stuck and do not reach steady state.}
\end{figure}

\subsubsection{Random-subset tempering}
\label{sec:tempering}

The slow equilibration poses a chicken-and-egg problem for the MCMC
method: without a good guess for the position of the particle, the
method does not yield a reliable MAP estimate in any reasonable time,
but it is difficult to obtain a sufficiently precise guess (within a
few micrometers of the MAP position) without knowing the MAP position
ahead of time. One way to address this challenge is parallel
tempering\cite{earl_parallel_2005}, in which several different MCMC
chains are run in parallel at different ``temperatures'' (a parameter
used to vary the sharpness of the probability distributions).
Higher-temperature chains can more easily explore the entire probability
distribution.

Inspired by this approach, we use a tempering scheme, but one that is
simpler to implement and computationally less expensive. Our technique
is based on previous work \cite{dimiduk_random-subset_2014} showing that
only a small, randomly selected subset of the pixels in a hologram are
needed to extract information about the parameters in the model. We
therefore start by choosing a small fraction $f=0.0005$ of pixels (10
pixels out of 20,000 total in the hologram). We then initialize 500
walkers with positions sampled from the prior and run the MCMC algorithm
for a short period (30 samples after burn-in, which takes 15 seconds or
less on a modern CPU at these low fractions). We then calculate the
marginalized probability distribution of $z$ based on this MCMC chain,
select a Gaussian distribution around the most likely value of $z$, and
use this distribution to select the walker positions for the next
iteration at higher fraction. We then repeat this procedure at
successively higher fractions. At each stage, we change only the initial
distribution of walkers, and not the priors. For holograms of single
spheres, we find that three preliminary stages (at $f=0.00025$,
$f=0.00116$, and $f=0.00539$) are sufficient to obtain a narrow
distribution for $z$, which we then use for a final run at $f=0.025$. We
also find that the posterior distribution does not change significantly
at fractions above $f=0.025$ for the holograms we analyze in this paper.

To illustrate this approach, we analyze a hologram of a
1-\si{\micro\meter}-diameter polystyrene sphere in water (index 1.33). Figure
\ref{fig:subsets} shows that there is a peak in the posterior
probability of $z$ even at the lowest pixel fraction. The probability
distributions show that the fraction of pixels is analogous to inverse
temperature in a tempering scheme; as the fraction increases, the peak
in the posterior distribution becomes sharper, and the other modes
subside. The prior and posterior distributions at each stage for a
1-\si{\micro\meter}-diameter silica sphere in water are shown in
Figure~\ref{fig:tempering}. The plots show that one can start with a
very wide prior in $z$ (a range of 0 to 100~\si{\micro\meter} above the focal
plane, representing nearly complete ignorance about the position of the
particle within the sample volume), and refine the uncertainty in the
position down to a fraction of the particle size.

\begin{figure}
  \centering
  \includegraphics[width=1.7in]{figures/bayesian-holography/masked_f0_0005}
  \includegraphics[width=1.7in]{figures/bayesian-holography/masked_f0_01}
  \includegraphics[width=1.7in]{figures/bayesian-holography/masked_f0_1}

  \includegraphics{figures/bayesian-holography/prob_vs_z_f0_0005}
  \includegraphics{figures/bayesian-holography/prob_vs_z_f0_01}
  \includegraphics{figures/bayesian-holography/prob_vs_z_f0_1}
  \caption{Comparison of posterior probability distributions of the
    parameter $z$, as obtained through MCMC calculations at three subset
    fractions. Top row shows a hologram of a 1-\si{\micro\meter}-diameter
    polystyrene sphere in water (same hologram in all three columns).
    The red pixels are the randomly chosen pixels used in the
    calculations. The bottom row shows the posterior probabilities at the
    different random-subset fractions $f$. The peak in the probability
    gets sharper as the fraction increases, but it is present even at
    the lowest fraction, which represents only 10 pixels in the
    hologram.}
  \label{fig:subsets}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=4.9in]{figures/bayesian-holography/tempering}
  \caption{The random-subset tempering procedure involves refining
    position estimates at successively higher fractions, starting from
    coarse guesses. Top row shows the prior probability distributions of
    $z$ for each fraction and bottom row the posterior, as obtained from
    MCMC ensemble sampling.  The prior for the second and third stage are
    chosen from the peak in the posterior for the previous stage. The
    final posterior distribution is much narrower than the prior used
    for the first stage; note the change in scale in the horizontal
    axis, which shows the $z$ position in micrometers.}
  \label{fig:tempering}
\end{figure}

This tempering approach allows us to start with very little information
and---in under a minute---obtain an excellent starting point for a final
MCMC calculation that requires only a short burn-in. We speed up the
computation by doing this final run at a fraction $f=0.025$, for which
the calculation is about 40 times faster than that for the full
hologram. We use 500 walkers in the final run. Under these conditions,
it takes about 20 minutes on a modern desktop CPU to obtain 25000
samples, roughly 5000 of which are independent. Although the method
takes longer than least-squares fitting, it does not require a good
initial guess for $z$ (or $\alpha$). Thus there is less manual
intervention required to accurately estimate the parameters, and the
procedure can more easily be parallelized.

\subsection{Particle tracking}

With the tempering scheme in hand, we proceed to the particle-tracking
problem. We present results from the analysis of a trajectory of a
1-\si{\micro\meter} silica particle, as described in section \ref{sec:computation}.
First, we note that the MAP estimates for the position, as obtained from
a single frame of the trajectory, are close to those obtained from
least-squares fitting, but with some key differences as shown in Figure
\ref{fig:position_pairplots}, in which we have marginalized over $n$,
$r$, and $\alpha$. The MAP estimates, with boundaries of the 68\%
credibility interval shown in subscript and superscript, are
$x=17.48^{+0.024}_{-0.026}$~\si{\micro\meter}, $y=17.64^{+0.027}_{-0.028}$~\si{\micro\meter},
and $z=1.56^{+0.077}_{-0.072}$~\si{\micro\meter}. There is no appreciable
covariance between $x$ and $y$ and only slight covariance between the
in-plane coordinate $(x, y)$ and the height $z$. By comparison, a
least-squares fit yields $x=17.483\pm 0.0016$~\si{\micro\meter}, $17.640\pm
0.0026$~\si{\micro\meter}, and $z=1.569\pm 0.0019$~\si{\micro\meter}, where the uncertainties
represent one-sigma confidence intervals, as calculated from the
parameter covariance matrix.

The least-squares fits yield parameters with smaller uncertainties, but
this is because they do not account for the uncertainties in the other
parameters. For example, the uncertainty in $z$ is larger in the
Bayesian estimate because $z$ is correlated with $n$, $r$, and $\alpha$;
in the least-squares fit the uncertainty in $z$ is calculated based on
fixed (best-fit) values for these other parameters. The differences
between the MAP estimates and the best-fit values likely also arise
because the Bayesian approach accounts for the correlations. We
therefore argue that the marginalized distributions offer more realistic
estimates of the parameters and their uncertainties.

\begin{figure}
  \centering
  \begin{center}
    \includegraphics{figures/bayesian-holography/tracking_pairplot}
  \end{center}
  \caption{Results for particle position estimation from a single frame
    of a time-series. The marginalized posteriors for $x$, $y$, and $z$
    obtained from our MCMC scheme are shown along the diagonal, while
    the off-diagonal contour plots show the covariances (Gaussian kernel
    density estimates). Orange lines are results from least-squares
    fitting.}
  \label{fig:position_pairplots}
\end{figure}

With the posterior distributions obtained for each frame in a
time-series, we construct the trajectory of the particle along with the
uncertainty by sequential MCMC calculations, in which we use the
posterior from one frame in the time-series to inform priors for the
subsequent frame. Specifically, we approximate the posterior from one
frame as a Gaussian and use it as a prior for $n$, $r$, and $\alpha$ for
the next frame. For $x$, $y$, $z$, we use a Gaussian prior centered on
the MAP position of the previous frame and assign it a width of
$0.1$~\si{\micro\meter} to allow for diffusion. The results of our procedure are
shown in Figure \ref{fig:trajectory}.

\begin{figure}
  \centering
  \includegraphics{figures/bayesian-holography/trajectory_x}
  \includegraphics{figures/bayesian-holography/trajectory_y}
  \includegraphics{figures/bayesian-holography/trajectory_z}
  \caption{Particle trajectory inferred from analysis of a 300-frame
    time-series (3 seconds total). The blue curve shows the MAP position
    estimate at each frame, and the shaded magenta region indicates the
    uncertainty (68\% credibility region).}
  \label{fig:trajectory}
\end{figure}

The uncertainty in the trajectory is important for hypothesis testing or
model selection. In this data set the particle is subject to Brownian
motion and forces from an optical trap and gravity. In other, similar
experiments, one might wish to compare the actual trajectory to an
equation of motion derived from a model of the forces. To quantitatively
determine whether a particular equation of motion accurately models the
data, one needs to know the uncertainty on the trajectory.

\subsection{Particle characterization}

We apply the same tempered MCMC scheme to the particle characterization
problem, using the same data set as for the tracking analysis. We show
the results from a single frame and those from 10 frames in Figure
\ref{fig:characterization}. For the 10-frame analysis, we use the
posterior for one frame as the prior for the subsequent frame in the
time-series, then we combine the MCMC samples from the 10 individual
frames.  The resulting posterior distribution is significantly
narrower (by a factor of about 10 for $n$ and 26 for $r$), thus yielding
much better estimates for the refractive index and radius than could be
obtained from a single frame.

\begin{figure}
  \centering
  \includegraphics{figures/bayesian-holography/characterization}
  \caption{Posterior distributions of particle properties (refractive
    index $n$ and radius $r$) inferred from analysis of one frame (top
    row) and ten successive frames (bottom row) of a time-series of
    holograms. Plots are Gaussian kernel density estimates of MCMC
    samples.}
  \label{fig:characterization}
\end{figure}

The refinement of the particle parameters is important in applications
such as particle sizing. In holographic particle sizing, the size
distribution is measured by combining measurements for individual
particles. It is therefore important that the uncertainty on the each
measured particle be as small as possible, so that narrow size
distributions and multimodal distributions can be resolved. In the
10-frame measurement we show in Figure~\ref{fig:characterization}, the
uncertainty relative to the MAP radius is 0.016\%, which (as desired)
is much smaller than the width of even a monodisperse colloidal size
distribution.

\section{Conclusions}
\label{sec:conclusions}

The Bayesian approach we have demonstrated is well-suited to the problem
of extracting information from holograms of colloidal particles. One
reason is that we understand well how particles scatter light, and there
are accurate ways to calculate the scattering. Another is that the
technique captures essential details of the scattered field, such
as the phase information that is encoded in the interference fringes.
All of this means that accurate forward models for hologram formation can be
developed from first principles.  With a good forward model, MCMC
schemes like the one we demonstrate here can be brought to bear on
parameter estimation.

The main disadvantage of the Bayesian approach compared to least-squares
fitting is that the computations take much longer. This disadvantage is
outweighed, in our opinion, by the many advantages of the approach.
First, prior information can be specified explicitly at the outset, and
need not be spread over initial guesses, constraints, and penalty
functions. Second, the tempered MCMC scheme we have developed does not
require accurate initial guesses, thus eliminating much of the manual
effort required to tune the guesses. Third, the analysis yields the
entire joint probability distribution of the parameters, so that the
uncertainty can be quantified accurately and marginalized to account for
uncertainties in parameters that are not of interest. Thus the method
can be applied, as we have shown, to extract realistic uncertainties on
the trajectory of a particle or on its properties, taking into account
data over an entire time-series.

But we have not demonstrated what could be the most important advantage
of the Bayesian approach: its adaptability to other, more
complex forward models. In the future, we aim to apply this technique
to more complex scattering models such as the superposition solution for
multiple spheres, which allows us to infer the structure and dynamics of
colloidal clusters \cite{perry_real-space_2012}. We would also like to
model correlated noise (speckle) in the holograms as well as other
optical effects arising from propagation along the optical train of the
microscope. A more detailed description of noise and imaging effects
should allow us to model the hologram more accurately than the
current formulation, which lumps these effects into a single parameter,
$\alpha$.

All of these extensions to the forward model will require more
parameters, which is problematic for least-squares analysis but much
less so for the Bayesian approach. In least-squares techniques, most of
the computational time is spent evaluating the Hessian, which scales as
the number of parameters squared. Furthermore, adding parameters
complicates the chi-squared surface, producing local minima
corresponding to unphysical solutions. It becomes more and more
difficult to tune the least-squares algorithm as the number of
parameters increases, even when there is no risk of overfitting. In
contrast, the MCMC techniques used in a Bayesian approach do not rely on
gradients and are designed to handle models with hundreds or even
thousands of parameters. Also, tempering methods allow one to explore
the entire probability distribution, including local maxima. Thus we
believe that tempered MCMC approaches like the one we describe in this
paper will scale better to more complex models.

Finally, although the approach we describe is by no means simple, it is
easy to use because it requires very little tuning. We therefore hope
that it will be useful to other researchers considering the use of
holographic microscopy. The subset tempering algorithm makes it
straightforward for a scientist who is not experienced with holography
to make measurements. To facilitate the adoption of these techniques, we
have implemented them in our open-source hologram analysis library
HoloPy \cite{dimiduk_holopy_2016}.

\section*{Funding}

We acknowledge support from the National Science Foundation
through grant number DMR-1306410.

\section*{Acknowledgments}

We thank Anna Wang for providing the raw holograms that we analyze in
the paper.
