# Estimation of the Diffusion Constant from Intermittent Trajectories with
Variable Position Uncertainties

Peter Relich, Mark Olah, Patrick Cutler, Keith Lidke

###### Abstract

The movement of a particle described by Brownian motion is quantified by a
single parameter, $D$, the diffusion constant. The estimation of $D$ from a
discrete sequence of noisy observations is a fundamental problem in biological
single particle tracking experiments since it can report on the environment
and/or the state of the particle itself via hydrodynamic radius. Here we
present a method to estimate $D$ that takes into account several effects that
occur in practice, that are important for correct estimation of $D$, and that
have hitherto not been combined together for estimation of $D$. These effects
are motion blur from finite integration time of the camera, intermittent
trajectories, and time-dependent localization uncertainty. Our estimation
procedure, a maximum likelihood estimation, follows directly from the
likelihood expression for a discretely observed Brownian trajectory that
explicitly includes these effects. The manuscript begins with the formulation
of the likelihood expression and then presents three methods to find the exact
solution. Each method has its own advantages in either computational
robustness, theoretical insight, or the estimation of hidden variables. We
then compare our estimator to previously published estimators using a squared
log loss function to demonstrate the benefit of including these effects.

##  1 Introduction

Single Particle Tracking (SPT) is a method to observe and classify the motion
of individual particles as trajectories: estimates of a particle’s position in
a sequence of discrete measurement times. In the field of biological
microscopy, SPT has been used for finding and analyzing protein motion in
heterogeneous environments like the cellular membrane (1, 2) and cytoplasm (3,
4). The SPT trajectory information can be used to resolve variations in the
individual motion of molecules that would otherwise be lost in ensemble
imaging techniques.

In the analysis of trajectories, the pure Brownian motion model is often the
first model used to describe a trajectory in the absence of prior information
about the movement. The behavior of a single particle dominated by Brownian
motion can be described by a normal distribution with the variance term
proportional to a single physical scale parameter $D$, the diffusion constant;
which makes Brownian motion the simplest model for describing stochastic
motion. More complicated behavior could potentially be modeled as Brownian
motion with discrete changes in the diffusion constant that could be
identified with change point analysis (5). Therefore, the estimation of the
diffusion constant of a particle from discrete, noisy, and possibly short
particle trajectories is a fundamental problem in single particle tracking.

In this manuscript, we focus on the likelihood distribution of $D$. We present
a maximum-likelihood-based approach for estimating the diffusion constant of a
particle given an SPT trajectory that includes the individual localization
error for each position in the trajectory, the time of the observation, and
the camera integration time. Our approach is based on a direct solution to the
likelihood equation for the observation of a particular trajectory. The need
for such an estimation procedure has evolved out of the rapid progress that
has been made in SPT analysis techniques over the last few years (6, 7, 8, 9,
10). In particular, some emitter localization techniques can not only
accurately resolve the location of an emitter to tens of nanometers, but can
also reliably estimate the localization error (11). Because the signal to
noise ratio of a particle can vary significantly from frame to frame in an
image sequence (e.g. from varying background, or photobleaching of the probe),
the localization error reported for each observation in a trajectory can also
vary significantly from frame to frame. We have therefore developed an
estimator that takes into account this information.

###  1.1 Background and Related Work

Historically, one of the primary techniques for estimating the diffusion
constant from trajectories relied on a linear regression of the mean-squared-
displacement (MSD) of the tracked particle coordinates as a function of time
lag (12). In the absence of measurement errors, the observed MSD for pure
Brownian motion scales linearly with time lag and intersects at the origin,
allowing the direct recovery of the diffusion constant from a linear
regression on the well sampled data points. It has been shown that a
regression of the MSD with an offset parameter can be interpreted to account
for the cumulative effects of static (13) and dynamic measurement errors (14).
If the MSD is built using the same data points for multiple time lags, the
correlation between MSD values must also be taken into account in the
regression (12, 15, 16). Although it seems theoretically possible to include
individual localization error into the MSD regression, to date this has not
been described.

A separate line of work has focused on maximum likelihood approaches to the
estimation procedure. A maximum likelihood estimator works by finding the
maximum of a likelihood function $\mathcal{L}(D)=P\\!\left({O}\,|\,{D}\right)$
that gives the probability of observing a particular trajectory $O$, given a
diffusion constant $D$. Ideally this probability should incorporate both the
variable localization errors of the trajectory and effect of motion-blur. The
motion-blur effect arises from the fact that each localization is performed on
data that is acquired over some non-zero exposure time. Typically camera
sensors integrate the signal over the exposure time resulting in a blurring of
the particle image. This blurring has important numerical effects on the
likelihood function (17). A specific solution to the likelihood function has
been accurately derived that incorporates the effects of motion-blur but with
the caveat that only a single global localization error estimate is used as an
input or estimated (18, 16). This estimator is a more robust alternative to
the MSD-based estimators because it can implement all trajectory information
without incurring systematic error when the data is not well conditioned for a
linear regression. Subsequent work has extended this approach to deal with
non-uniformly spaced or intermittent trajectories (19), however the particular
implementation in (19) did not account for motion blur. Maximum likelihood
estimators are not the only class of diffusion estimators that have evolved
recently; continued development on displacement-based estimators has resulted
in an estimator that incorporates the effects of covariances between
sequentially observed displacements (20).

In this work we provide a generalized solution to the likelihood function,
incorporating variable localization errors and variable displacement periods,
which results in an improvement in estimation accuracy for short trajectories,
trajectories with large variations in localization accuracy, and trajectories
with intermittently spaced measurements. In Sec. 2 we formulate the diffusion
likelihood function to directly incorporate the effects of motion-blur,
variable localization errors, and intermittent or non-uniformly spaced
observations in time. We present three independent solutions to this
likelihood function. The first derivation, the recursive method (Sec. 3), is a
sequential integration of the nuisance parameters and provides the fastest
numerical implementation. The second derivation, the Laplace method (Sec. 4),
utilizes a second order Taylor expansion to express the likelihood as a
multivariate Gaussian in the basis of integration. The Laplace method
additionally returns the maximum likelihood values of the true positions given
a sampled $D$. The third derivation, the Markov method (Sec. 5), calculates
the characteristic function in order to express the likelihood in the basis of
displacements. The Markov method allows us to verify that the generalized form
of the expression derived in (18) is the same distribution as the expressions
derived in this manuscript. The Markov method was also instrumental in
determining the coefficients necessary to reduce the computational complexity
of all the methods (Sec. 13.2). Each of these derivations leads to an
independent, numerically accurate computational algorithm for estimating the
likelihood of $D$ (Sec. 6), making full use of all the information contained
in a noisy trajectory. The resulting likelihood calculation allows for robust
computations in specific problems, such as a maximum likelihood estimator,
maximum a posteriori estimate, or change point analysis. We compare the
results of our maximum likelihood estimator (MLE) to the current state of the
art estimation software (16) with the squared log loss function and
demonstrate that the additional information provided from the localization
errors allows for better estimates of $D$ with trajectories parameterized by
any non-constant, but otherwise arbitrary distribution of localization
variances.

##  2 Theory

If a diffusing particle is accurately and exactly observed at a discrete
sequence of $N+1$ positions $\mathbf{X}=\\{\mathbf{x}_{i}\\}^{N+1}_{i=1}$ at
times $t_{i}$, then $P\\!\left({\mathbf{X}}\,|\,{D}\right)$, the probability
of sequence $\mathbf{X}$ given diffusion constant $D$, is

$P\\!\left({\mathbf{X}}\,|\,{D}\right)=\prod^{N}_{i=1}P\\!\left({\mathbf{x}_{i+1}}\,|\,{\mathbf{x}_{i}}\right).$
(1)

In Eq. 1,
$P\\!\left({\mathbf{x}_{i+1}}\,|\,{\mathbf{x}_{i}}\right)=P\\!\left({\mathbf{x}_{i+1}}\,|\,{\mathbf{x}_{i},D}\right)$
is the probability density of each discrete jump from
$\mathbf{x}_{i}\to\mathbf{x}_{i+1}$ over time step $\delta
t_{i}=t_{i+1}-t_{i}$, given diffusion constant $D$.

When measured experimentally, however, the true positions $\mathbf{X}$ are
never known exactly, but are related to $N$ observed positions $\mathbf{O}$ by
some distribution
$P\\!\left({\mathbf{o}_{i}}\,|\,{\mathbf{x}_{i},{\mathbf{x}_{i+1}}}\right)$,
where the dependence on both $\mathbf{x}_{i}$ and $\mathbf{x}_{i+1}$ arises
from the effects of exposure time integration by the observation apparatus
which will be dealt with in detail later. Under this experimental model,
$P\\!\left({\mathbf{O},\mathbf{X}}\,|\,{D}\right)$, the combined likelihood of
the observed positions $\mathbf{O}$ and the actual positions $\mathbf{X}$ is a
product of the observation probability densities
$P\\!\left({\mathbf{o}_{i}}\,|\,{\mathbf{x}_{i},\mathbf{x}_{i+1}}\right)$ and
the diffusion transition probability densities
$P\\!\left({\mathbf{x}_{i+1}}\,|\,{\mathbf{x}_{i}}\right)$ for each of the $N$
observed positions and displacements,

$P\\!\left({\mathbf{O},\mathbf{X}}\,|\,{D}\right)=\prod^{N}_{i=1}P\\!\left({\mathbf{o}_{i}}\,|\,{\mathbf{x}_{i},\mathbf{x}_{i+1}}\right)P\\!\left({\mathbf{x}_{i+1}}\,|\,{\mathbf{x}_{i}}\right).$
(2)

Since $\mathbf{X}$ is unknown for experimental data, we integrate Eq. 2 over
all possible $\mathbf{X}$ to marginalize out the dependence on $\mathbf{X}$,
and write the diffusion likelihood as an integral over the space of all
$\mathbf{X}$-values,

$P\\!\left({\mathbf{O}}\,|\,{D}\right)=\int\mathrm{d}{\mathbf{X}}\,P\\!\left({\mathbf{O},\mathbf{X}}\,|\,{D}\right)=\int\mathrm{d}{\mathbf{X}}\,\prod^{N}_{i=1}P\\!\left({\mathbf{o}_{i}}\,|\,{\mathbf{x}_{i},\mathbf{x}_{i+1}}\right)P\\!\left({\mathbf{x}_{i+1}}\,|\,{\mathbf{x}_{i}}\right).$

Experimental data typically involves trajectories with two or three spatial
dimensions. For diffusion in an isotropic medium and particle uncertainties
given as normal distributions with no covariance among the spatial dimensions,
the probability distribution of a particular displacement in each dimension is
separable. Thus, if $\Upsilon$ is the number of dimensions, then

$P\\!\left({\mathbf{O}}\,|\,{D}\right)=\prod_{n=1}^{\Upsilon}P\\!\left({O_{n}}\,|\,{D}\right).$
(3)

Hence, it is sufficient to only consider the estimation problem in the one-
dimensional (1D) case $O=\\{o_{i}\\}_{i=1}^{N}$, and

$P\\!\left({O}\,|\,{D}\right)=\int_{\mathbb{R}^{N+1}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right).$
(4)

###  2.1 Accounting for the effects of exposure time integration

Equation 40 is the fundamental description of the likelihood of diffusion
constant $D$ given observations $O$. Unfortunately, solving for this
expression explicitly is difficult because every $o_{i}$ term is dependent on
both $x_{i}$ and $x_{i+1}$. This is because the estimate of $o_{i}$’s position
is typically made from data collected over an exposure time
$0<t_{\epsilon}\leq t_{i+1}-t_{i}$. If the observational apparatus is a camera
sensor, the signal will be integrated over the frame, resulting in a motion-
blurred image of the moving particle, hence the observed location is
conditional upon the particle’s true position at the beginning ($x_{i}$) and
end ($x_{i+1}$) of the frame.

In the case where exposure time $t_{\epsilon}$ goes to 0, but $\delta
t_{i}=t_{i+1}-t_{i}$ remains constant, the motion-blur effect is no longer
present, so the observed location $o_{i}$ depends only on position $x_{i}$,

$P\\!\left({O}\,|\,{D}\right)=\int_{\mathbb{R}^{N}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}P\\!\left({o_{i}}\,|\,{x_{i}}\right)\prod^{N-1}_{j=1}P\\!\left({x_{j+1}}\,|\,{x_{j}}\right).$
(5)

Without the additional dependence on $x_{i+1}$, the methods required to solve
the integral in Eq. 5 are simpler. In order to use this simpler
representation, we will transform Eq. 40 into a form which resembles Eq. 5,
and seek functions $\mathcal{M}(o_{i},x_{i})$ and $\mathcal{T}(x_{j+1},x_{j})$
such that

$P\\!\left({O}\,|\,{D}\right)=\int_{\mathbb{R}^{N+1}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right)=\int_{\mathbb{R}^{N}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}\mathcal{M}(o_{i},x_{i})\prod^{N-1}_{j=1}\mathcal{T}(x_{j+1},x_{j}).$
(6)

The function $\mathcal{T}(x_{i+1},x_{i})$ stands for the transition
probability; it is simply the probability of a particle diffusing with
constant $D$ moving from $x_{i}$ to $x_{i+1}$, over time $\delta t_{i}$. The
function $\mathcal{M}(o_{i},x_{i})$ stands for the measurement probability and
it encapsulates the net effect of both the measurement localization error and
the motion-blur. The details of the representation equivalence of Eq. 6 are
important for correctness, but they also unnecessarily complicate the
exposition, and so can be found in Sec. 13. Other authors (14, 18) have
investigated the motion-blur effects of exposure time integration, and found
that the effect can be approximated by an effective decrease in variance of
the measurement localization error, dependent on diffusion constant $D$ and
exposure time $t_{\epsilon}$. Our derivations in Sec. 13 agrees with the
effective correction factor in (14, 18), and more importantly provides a form
for the diffusion likelihood that is directly amenable to the solution
techniques we employ in Secs. 3,4, and 5.

The result of the transformation of Eq. 6 is that the effective measurement
function $\mathcal{M}_{i}$ and the transition function $T_{i}$ take the form
of normalized Gaussians. We use the notation

$\,\mathcal{N}(a,a_{0},\eta)=\frac{1}{\sqrt{2\pi\eta}}\exp{\left[-\frac{(a-a_{0})^{2}}{2\eta}\right]}.$

to represent the normalized Gaussian function with variance $\eta=\sigma^{2}$
centered around mean $a_{0}$ considered as a function of $a,a_{0},$ and
$\eta$. Using this notation, we can succinctly represent the measurement and
transition functions as,

$\displaystyle\mathcal{T}_{i}=\mathcal{T}_{i}(x_{i+1},x_{i})=$
$\displaystyle\,\mathcal{N}(x_{i+1},x_{i},\omega_{i}(D)),\quad\textrm{for
}1\leq i\leq N-1,\textrm{ and }$ (7)
$\displaystyle\mathcal{M}_{i}=\mathcal{M}_{i}(o_{i},x_{i})=$
$\displaystyle\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i}(D)),\quad\textrm{for
}1\leq i\leq N.$ (8)

The transition functions (Eq. 46), are unaffected by the motion blur
transformation and their Gaussian representation follows directly from the
normally distributed displacements of diffusive processes, hence the variance
is

$\omega_{i}(D)=2D\delta t_{i}.$

For the measurement functions (Eq. 45), the variance $\varepsilon_{i}(D)$, is
the variance due to measurement error, $v_{i}$, combined with a correction for
the effect of motion-blur that is dependent on diffusion constant $D$ and
exposure time $t_{\epsilon}$,

$\varepsilon_{i}(D)=v_{i}-Dt_{\varepsilon}/3.$

where the factor of $1/3$ comes from the continuous limit integration of
photon emissions for averaged Brownian trajectories (Sec. 13.1). It is
important to note that the independence of $t_{\varepsilon}$ and $\delta
t_{i}$ allows for gaps in the trajectories, where $\delta t_{i}$ could span a
duration of multiple frames but $t_{\varepsilon}$ is the exposure time of a
single frame.

The result is that Eq. 6 allows us to express the likelihood function exactly
in a simple form that deals directly with variable localization error, motion-
blur effects, and missing or irregularly spaced trajectory localizations,

$P\\!\left({O}\,|\,{D}\right)=\int_{\mathbb{R}^{N}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}\mathcal{M}_{i}\prod^{N-1}_{j=1}\mathcal{T}_{j}.$
(9)

##  3 Recursive Method

The notation for the transition and measurement functions allows us to define
the likelihood function $\mathcal{L}(D)$, by writing Eq. 9 in a form that
emphasizes the dependencies on each marginalized position $x_{i}$,

$\mathcal{L}(D)=P\\!\left({O}\,|\,{D}\right)=\int\mathrm{d}{x_{N}}\,\mathcal{M}_{N}\int\mathrm{d}{x_{N-1}}\,\mathcal{M}_{N-1}\mathcal{T}_{N-1}\ldots\int\mathrm{d}{x_{2}}\,\mathcal{M}_{2}\mathcal{T}_{2}\int\mathrm{d}{x_{1}}\,\mathcal{M}_{1}\mathcal{T}_{1}.$
(10)

The form of Eq. 10 leads to a direct recursive solution, taking into account
the properties of integrals over products of normalized Gaussian functions.
Define $\mathcal{L}_{i}$ as the sub-integrand of $\mathcal{L}(D)$ considering
only the first $i$ observations,

$\displaystyle\mathcal{L}_{1}(D,x_{2})$
$\displaystyle=\int\\!\mathcal{M}_{1}\mathcal{T}_{1}\mathrm{d}{x_{1}}\,,$
$\displaystyle\mathcal{L}_{i}(D,x_{i+1})$
$\displaystyle=\int\\!\mathcal{M}_{i}\mathcal{T}_{i}\mathcal{L}_{i-1}\mathrm{d}{x_{i}}\,,\quad
2\leq i\leq N-1$ $\displaystyle\mathcal{L}(D)=\mathcal{L}_{N}(D)$
$\displaystyle=\int\\!\mathcal{M}_{N}\mathcal{L}_{N-1}\mathrm{d}{x_{N}}\,.$

Now, consider that the integral of a product of two normalized Gaussians
sharing a parameter $x$, with means and variances denoted by $c_{i}$ and
$\varphi_{i}$ respectively, is itself a normalized Gaussian (Sec. 11),

$\int\\!\mathrm{d}{x}\,\prod^{2}_{i=1}\,\mathcal{N}(x,c_{i},\varphi_{i})=\,\mathcal{N}(c_{1},c_{2},\varphi_{1}+\varphi_{2}).$
(11)

Hence, $\mathcal{L}_{1}$ is a normalized Gaussian in parameter $x_{2}$,

$\mathcal{L}_{1}(D,x_{2})=\int\mathrm{d}{x_{1}}\,\mathcal{M}_{1}\mathcal{T}_{1}=\,\mathcal{N}(o_{1},x_{2},\varepsilon_{1}+\omega_{1}).$

This implies that $\mathcal{L}_{2}(D,x_{3})$ which is an integral over
positions $x_{2}$, can now be written as an integral over three normalized
Gaussians, all of which share parameter $x_{2}$,

$\mathcal{L}_{2}(D,x_{3})=\int\mathrm{d}{x_{2}}\,\mathcal{M}_{2}\mathcal{T}_{2}\int\mathrm{d}{x_{1}}\,\mathcal{M}_{1}\mathcal{T}_{1}=\int\mathrm{d}{x_{2}}\,\mathcal{M}_{2}\mathcal{T}_{2}\mathcal{L}_{1}.$
(12)

Similarly, the integral of a product of three normalized Gaussians sharing the
integrated parameter is itself a product of two normalized Gaussians (Sec.
11),

$\int\\!\mathrm{d}{x}\,\prod^{3}_{i=1}\,\mathcal{N}(c_{i},x,\varphi_{i})=\,\mathcal{N}(c_{1},c_{2},\varphi_{1}+\varphi_{2})\,\mathcal{N}(c_{3},c^{\prime},\gamma),$
(13)

where,

$\displaystyle
c^{\prime}=\frac{c_{1}\varphi_{2}+c_{2}\varphi_{1}}{\varphi_{1}+\varphi_{2}},\quad\text{and}\quad\gamma=\frac{\varphi_{1}\varphi_{2}+\varphi_{1}\varphi_{3}+\varphi_{2}\varphi_{3}}{\varphi_{1}+\varphi_{2}}.$

Hence, applying Eq. 13 to Eq. 12, we find that

$\mathcal{L}_{2}=\int\mathrm{d}{x_{2}}\,\\\
M_{2}T_{2}\mathcal{L}_{1}=\,\mathcal{N}(o_{1},o_{2},(\varepsilon_{1}+\omega_{1})+\varepsilon_{2})\,\mathcal{N}(x_{3},\mu_{2},\eta_{2})$

is a product of two normalized Gaussians, one of which depends on $x_{3}$ and
the other is a constant with respect to $X$. The variables $\mu_{2}$ and
$\eta_{2}$ follow from Eq. 11 and Eq. 13, and are the second pair of values in
a recursive solution. Since all subsequent integrals, barring the last
integral, can be expressed as a product of three normalized Gaussians, we can
express the recursion variables as

$\mu_{1}=o_{1},\quad\ \eta_{1}=\varepsilon_{1}+\omega_{1},\quad\
\text{and}\quad\alpha_{1}=\eta_{1}+\varepsilon_{2},$ (14)

and for $2\leq i\leq N-1$,

$\mu_{i}=\frac{\mu_{i-1}\varepsilon_{i}+\eta_{i-1}o_{i}}{\alpha_{i-1}},\quad\eta_{i}=\frac{\eta_{i-1}\varepsilon_{i}}{\alpha_{i-1}}+\omega_{i},\quad\text{and}\quad\alpha_{i}=\eta_{i}+\varepsilon_{i+1}.$
(15)

Finally, this allows us to express our integrands $\mathcal{L}_{i}$ as

$\displaystyle\mathcal{L}_{1}$
$\displaystyle=\,\mathcal{N}(x_{2},\mu_{1},\eta_{1})$ (16)
$\displaystyle\mathcal{L}_{i}$
$\displaystyle=\,\mathcal{N}(x_{i+1},\mu_{i},\eta_{i})\prod_{k=1}^{i-1}\,\mathcal{N}(o_{k+1},\mu_{k},\alpha_{k}),\quad
2\leq i\leq N-1$ $\displaystyle\mathcal{L}(D)=\mathcal{L}_{N}$
$\displaystyle=\prod_{k=1}^{N-1}\,\mathcal{N}(o_{k+1},\mu_{k},\alpha_{k}).$

Equation 16 is the final form of the recursive solution for $\mathcal{L}(D)$
which is simply the product of $N-1$ normalized Gaussians each of which has
parameters which come from a recursive relationship on $o_{i}$,
$\varepsilon_{i}$, and $\omega_{i}$. The value of $D$ that maximizes
$\mathcal{L}(D)$ is the maximum likelihood estimate.

##  4 Laplace Method

An independent solution for Eq. 9 can be obtained using the Laplace method
which is based on integrating the second moment of the Taylor expansion of the
exponential component of a function. Given that the second moment of a Taylor
expansion is quadratic, this ensures that the function under the integral is
always a Gaussian function (21). Another caveat to the Laplace method is that
the Taylor expansion has to occur about the peak of the exponential, so that
the first moment of the Taylor expansion goes to 0. To perform the Laplace
method, we express our likelihood
$\mathcal{L}(D)=P\\!\left({O}\,|\,{D}\right)$ in terms of exponential and non-
exponential components

$\mathcal{L}(D)=\int\\!\mathrm{d}{X}\,f(X)=\int\\!\mathrm{d}{X}\,h(X)\exp{\left[-g(X)\right]},$

where $f(x)$ is simply the integrand of Eq. 9,

$f(X)=h(X)\exp{\left[-g(X)\right]}=\prod_{i=1}^{N}\mathcal{M}_{i}\prod_{j=1}^{N-1}\mathcal{T}_{j}.$
(17)

Thus, using equations 45 and 46, we see that $h=h(X)$ is independent of $X$
and $g(X)$ is quadratic in $X$,

$\displaystyle h=$
$\displaystyle\prod^{N}_{i=1}\frac{1}{\sqrt{2\pi\varepsilon_{i}}}\prod^{N-1}_{i=1}\frac{1}{\sqrt{2\pi\omega_{i}}},$
$\displaystyle g(X)=$
$\displaystyle\sum^{N}_{i=1}\frac{(o_{i}-x_{i})^{2}}{2\varepsilon_{i}}+\sum^{N-1}_{i=1}\frac{(x_{i+1}-x_{i})^{2}}{2\omega_{i}}.$

The maximum likelihood estimate $\widehat{X}$ of the actual positions $X$,
given $D$ and $O$ will be wherever the integrand is maximized, and since
$g(X)\geq 0$,

$\widehat{X}=\operatornamewithlimits{\text{argmax}}_{X}f(X)=\operatornamewithlimits{\text{argmin}}_{X}g(X).$

Now, given that $g(X)$ is quadratic, a second order Taylor expansion of $g$
about $\widehat{X}$ is exact and the Laplace method will provide an exact
solution for $\mathcal{L}(D)$ as the integral can be shown to take the form of
a standard Gaussian integral. To see this, we write out the second order
Taylor expansion

$\int\\!\mathrm{d}{X}\,f(X)=h\int\\!\mathrm{d}{X}\,\exp{\left[-g(\widehat{X})-\nabla
g(\widehat{X})(X-\widehat{X})-\frac{1}{2}(X-\widehat{X})^{\intercal}\nabla\nabla
g(\widehat{X})(X-\widehat{X})\right]}.$ (18)

Since $\widehat{X}$ is the minima of $g(X)$, the gradient $\nabla
g(\widehat{X})=0$, and we can rearrange Eq. 18 to extract all terms
independent of $X$,

$\int\\!\mathrm{d}{X}\,f(X)=f(\widehat{X})\int\\!\mathrm{d}{X}\,\exp{\left[-\frac{1}{2}(X-\widehat{X})^{\intercal}\nabla\nabla
g(\widehat{X})(X-\widehat{X})\right]}.$ (19)

Furthermore, since $h$ is independent of $X$, we know that
$-\nabla\nabla\ln{f(X)}=\nabla\nabla g(X)=M$, where $M$ can be thought of as
the inverse of the covariance matrix for the multivariate Gaussian, or
equivalently as the Hessian matrix of $-\ln{f(X)}$. Substituting $M$ for
$\nabla\nabla g(X)$ in Eq. 19 we are left with a Gaussian integral with the
solution,

$\mathcal{L}(D)=f(\widehat{X})\int\\!\mathrm{d}{X}\,\exp{\left[-\frac{1}{2}(X-\widehat{X})^{\intercal}M(X-\widehat{X})\right]}=f(\widehat{X})\sqrt{\frac{(2\pi)^{N}}{\det{M}}}.$
(20)

The Hessian matrix $M$ is independent of $X$ and is symmetric tri-diagonal
with non-zero elements,

$\displaystyle M_{1,1}$ $\displaystyle=-\frac{\partial^{2}\ln f}{\partial
x_{1}^{2}}=\frac{1}{\varepsilon_{1}}+\frac{1}{\omega_{1}}$ (21) $\displaystyle
M_{i,i}$ $\displaystyle=-\frac{\partial^{2}\ln f}{\partial
x_{i}^{2}}=\frac{1}{\varepsilon_{i}}+\frac{1}{\omega_{i}}+\frac{1}{\omega_{i-1}},\quad
2\leq i\leq N-1$ $\displaystyle M_{N,N}$ $\displaystyle=-\frac{\partial^{2}\ln
f}{\partial x_{N}^{2}}=\frac{1}{\varepsilon_{N}}+\frac{1}{\omega_{N-1}}$
$\displaystyle M_{i,i+1}=M_{i+1,i}$ $\displaystyle=-\frac{\partial^{2}\ln
f}{\partial x_{i}\partial x_{i+1}}=-\frac{1}{\omega_{i}},\quad 1\leq i\leq
N-1.$

We also require $\widehat{X}$ to compute Eq. 20, which can be solved for with
the relation $-\nabla\ln f(\widehat{X})=0$ (Sec. 14), giving

$\widehat{X}=M^{-1}\Theta,$ (22)

where $\Theta=\\{\theta_{i}=o_{i}/\varepsilon_{i}\\}_{i=1}^{N}$.

##  5 Markov Method

In this section we present a derivation of the likelihood function
$\mathcal{L}(D)$ utilizing a technique developed by Andrey Markov (22) and
generalized by Chandrasekhar (23). Markov’s method allows us to transform
$P\\!\left({O}\,|\,{D}\right)$ (Eq. 9) from a function of the $N$ observed
positions $O=\\{o_{i}\\}_{i=1}^{N}$ into a function of the $N-1$ discrete
steps (displacements) between subsequent observations

$S=\\{s_{i}=o_{i+1}-o_{i}\\}_{i=1}^{N-1}.$

This is possible because the spatial invariance of diffusion means
$\mathcal{L}(D)$ depends not on the absolute spatial positions $O$, but only
their relative displacements, $S$. Thus we should expect that
$P\\!\left({O}\,|\,{D}\right)$ can also be expressed as
$P\\!\left({S}\,|\,{D}\right)$, however Eq. 10, which defines
$P\\!\left({O}\,|\,{D}\right)$, cannot be directly transformed into a function
on $S$. This is where Markov’s method allows us to solve for a function
$P\\!\left({S}\,|\,{D}\right)=P\\!\left({O}\,|\,{D}\right)$ for a given $D$
value. For a particular fixed $S$ and $\mathrm{d}{S}\,$ of interest, the value
of $P\\!\left({S}\,|\,{D}\right)\mathrm{d}{S}\,$ gives the probability that
variable $S^{\prime}$ is within the bounds

$S-\frac{1}{2}\mathrm{d}{S}\,\leq S^{\prime}\leq
S+\frac{1}{2}\mathrm{d}{S}\,.$ (23)

More formally, $P\\!\left({S}\,|\,{D}\right)\mathrm{d}{S}\,$ is the integral
over the volume $\mathrm{d}{S}\,$ around the point of interest $S$, and we let
$S^{\prime}$ represent the variable integrated over,

$P\\!\left({S}\,|\,{D}\right)\mathrm{d}{S}\,=\int^{S+\frac{1}{2}\mathrm{d}{S}\,}_{S-\frac{1}{2}\mathrm{d}{S}\,}\mathrm{d}{S^{\prime}}\,P\\!\left({O}\,|\,{D}\right).$

The issue remains that $P\\!\left({O}\,|\,{D}\right)$ is expressed in a basis
of $O$ rather than of $S$, and integrating with respect to bounds in a
different basis is non-trivial. In order to circumvent this issue, Markov
utilized a product of Dirichlet integrals,
$\Xi(S^{\prime})=\prod^{N-1}_{k=1}\xi_{k}(s_{k}^{\prime})$, to expand the
limits of integration to all space

$P\\!\left({S}\,|\,{D}\right)\mathrm{d}{S}\,=\int\\!\mathrm{d}{S^{\prime}}\,\Xi(S^{\prime})P\\!\left({O}\,|\,{D}\right).$
(24)

The idea is that for each dimension of $S$, the Dirichlet integral
$\xi_{k}(s_{k}^{\prime})$ acts like a continuous indicator function
determining if $s_{k}^{\prime}$ is within the bounds of Eq. 23,

$\xi_{k}(s_{k}^{\prime})=\frac{1}{\pi}\int\\!\mathrm{d}{\rho_{k}}\,\frac{\sin(\frac{1}{2}\mathrm{d}{s_{k}}\,\rho_{k})}{\rho_{k}}\exp{\left[\imath\rho_{k}(s_{k}^{\prime}-s_{k})\right]},$

so that,

$\xi_{k}(s_{k}^{\prime})=\begin{cases}1&s_{k}-\frac{1}{2}\mathrm{d}{s_{k}}\,\leq
s_{k}^{\prime}\leq s_{k}+\frac{1}{2}\mathrm{d}{s_{k}}\\\
0&\text{otherwise}\end{cases}.$

Therefore $\Xi(S^{\prime})$ is the indicator function acting over the whole
space of $S^{\prime}$, and determining if $S^{\prime}$ is within the volume
$\mathrm{d}{S}\,$ around our point of interest $S$,

$\Xi(S^{\prime})=\prod_{k=1}^{N-1}\xi_{k}(s_{k}^{\prime})=\begin{cases}1&\displaystyle\bigwedge_{k=1}^{N-1}\textstyle
s^{\prime}_{k}\in\left[s_{k}-\frac{1}{2}\mathrm{d}{s_{k}}\,,s_{k}+\frac{1}{2}\mathrm{d}{s_{k}}\,\right]\\\
0&\text{otherwise}\end{cases}.$

This puts Eq. 24 in the form

$P\\!\left({S}\,|\,{D}\right)\mathrm{d}{S}\,=\int\mathrm{d}{S^{\prime}}\,\frac{1}{\pi^{N-1}}\int\mathrm{d}{\bm{\rho}}\,\left[\prod^{N-1}_{k=1}\frac{\sin(\frac{1}{2}\mathrm{d}{s_{k}}\,\rho_{k})}{\rho_{k}}\right]\exp{\left[\imath\bm{\rho}^{\intercal}(S^{\prime}-S)\right]}P\\!\left({O}\,|\,{D}\right),$
(25)

where $\bm{\rho}=\\{\rho_{k}\\}_{k=1}^{N-1}$ is a vector of conjugate
coordinates to $S^{\prime}$. We can then rearrange Eq. 25 to move the integral
over $\mathrm{d}{S^{\prime}}\,$ and all factors dependent on $S^{\prime}$ into
function $\Lambda({\bm{\rho}})$,

$P\\!\left({S}\,|\,{D}\right)\mathrm{d}{S}\,=\frac{1}{\pi^{N-1}}\int\mathrm{d}{\bm{\rho}}\,\left[\prod^{N-1}_{k=1}\frac{\sin(\frac{1}{2}\mathrm{d}{s_{k}}\,\rho_{k})}{\rho_{k}}\right]\exp{\left[-\imath\bm{\rho}^{\intercal}S\right]}\Lambda(\bm{\rho}).$
(26)

Now, we can interpret $\Lambda({\bm{\rho}})$ as the characteristic function of
$P\\!\left({O}\,|\,{D}\right)$ in the $S^{\prime}$ basis and it has the form

$\Lambda(\bm{\rho})=\int\mathrm{d}{S^{\prime}}\,\exp{\left[\imath\bm{\rho}^{\intercal}S^{\prime}\right]}P\\!\left({O}\,|\,{D}\right).$
(27)

The form of Eq. 27 implies that $\Lambda(\bm{\rho})$ is the inverse Fourier
transform of $P\\!\left({O}\,|\,{D}\right)$. Due to the properties of the
Fourier transform, we assert $\Lambda(\bm{\rho})$ is a bounded function with a
finite integral because $P\\!\left({O}\,|\,{D}\right)$ is expressed as a
finite product of Gaussians with non-zero variance, hence it is a continuous
function (24). Given that $\int\mathrm{d}{\bm{\rho}}\,\Lambda(\bm{\rho})$ is
bounded and $\mathrm{d}{S}\,$ is small we can approximate the product of sinc
functions in Eq. 26 as

$\prod_{k=1}^{N-1}\frac{\sin(\frac{1}{2}\mathrm{d}{s_{k}}\,\rho_{k})}{\rho_{k}}=\prod_{k=1}^{N-1}\frac{\mathrm{d}{s_{k}}\,}{2}=\frac{\mathrm{d}{S}\,}{2^{N-1}}.$

Thus Eq. 26 becomes the Fourier transform of $\Lambda(\bm{\rho})$

$P\\!\left({S}\,|\,{D}\right)\mathrm{d}{S}\,=\frac{\mathrm{d}{S}\,}{(2\pi)^{N-1}}\int\mathrm{d}{\bm{\rho}}\,\exp{\left[-\imath\bm{\rho}^{\intercal}S\right]}\Lambda(\bm{\rho}).$
(28)

We are now interested in evaluating $\Lambda(\bm{\rho})$ explicitly. To do so
we note that $\int\mathrm{d}{S^{\prime}}\,P\\!\left({O}\,|\,{D}\right)=1$
since $P\\!\left({O}\,|\,{D}\right)$ is a probability distribution that, as we
have argued, can be equivalently expressed in the $S$ basis as
$P\\!\left({S}\,|\,{D}\right)$ (Sec. 14). With this understanding we evaluate
$\Lambda(\bm{\rho})$, by expanding the exponential under integration in Eq. 27
as a Taylor series about the origin

$\exp{\left[\imath\bm{\rho}^{\intercal}S^{\prime}\right]}=1+\imath\left(\sum^{N-1}_{j=1}\rho_{j}s_{j}^{\prime}\right)-\frac{1}{2}\left(\left[\sum^{N-1}_{j=1}\rho_{j}^{2}s_{j}^{\prime\,2}\right]+\left[\sum^{N-2}_{j=1}\sum^{N-1}_{k=j+1}2\rho_{j}\rho_{k}s^{\prime}_{j}s^{\prime}_{k}\right]\right)+\mathcal{O}(\bm{\rho}^{3}).$

Solving the integral for $\Lambda(\bm{\rho})$ with the Taylor expanded
exponential term given that we know $P\\!\left({O}\,|\,{D}\right)$ is a
normalized probability density under $S^{\prime}$ allows us to write
$\Lambda(\bm{\rho})$ in terms of expected values for $s_{k}$ (using
$\left\langle{\cdot}\right\rangle$ to represent expectation),

$\Lambda(\bm{\rho})=1+\imath\left(\sum^{N-1}_{j=1}\rho_{j}\left\langle{s_{j}^{\prime}}\right\rangle\right)-\frac{1}{2}\left(\left[\sum^{N-1}_{j=1}\rho_{j}^{2}\left\langle{s_{j}^{\prime\,2}}\right\rangle\right]+\left[\sum^{N-2}_{j=1}\sum^{N-1}_{k=j+1}2\rho_{j}\rho_{k}\left\langle{s^{\prime}_{j}s^{\prime}_{k}}\right\rangle\right]\right)+\mathcal{O}(\bm{\rho}^{3}).$
(29)

We take the approach of solving for the expectation values on $S^{\prime}$ to
see if those results will simplify the expression in Eq. 29. If we set one of
the observations, $o_{i}$, as a constant, we find that we can marginalize all
of $P\\!\left({O}\,|\,{D}\right)$ except for the terms that are independently
represented by the basis we are interested in. In other words, we find that

$\displaystyle\left\langle{s_{i}^{\prime}}\right\rangle$
$\displaystyle=\int\mathrm{d}{s_{i}^{\prime}}\,s_{i}^{\prime}\,\mathcal{N}(s_{i}^{\prime},0,\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i})=0$
$\displaystyle\left\langle{s_{i}^{\prime\,2}}\right\rangle$
$\displaystyle=\int\mathrm{d}{s_{i}^{\prime}}\,s_{i}^{\prime\,2}\,\mathcal{N}(s_{i}^{\prime},0,\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i})=\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i}.$
$\displaystyle\left\langle{s_{i}^{\prime}s_{i+1}^{\prime}}\right\rangle$
$\displaystyle=\int\mathrm{d}{s_{i}^{\prime}}\,\mathrm{d}{s_{i+1}^{\prime}}\,s_{i}^{\prime}s_{i+1}^{\prime}\,\mathcal{N}(S_{b},0,\Sigma_{b})=-\varepsilon_{i+1}.$

Where the substitution of variables required for integrating the expectation
of $\left\langle{s_{i}^{\prime}s_{i+1}^{\prime}}\right\rangle$ induces a
bivariate Gaussian function with location parameters
$S_{b}=[s_{i},s_{i+1}]^{\intercal}$ and covariance matrix

$\Sigma_{b}=\begin{bmatrix}\omega_{i}+\epsilon_{i}+\epsilon_{i+1}&-\epsilon_{i+1}\\\
-\epsilon_{i+1}&\omega_{i+1}+\epsilon_{i+1}+\epsilon_{i+2}\end{bmatrix}.$

Furthermore, we find that the separability of two non-adjacent displacements
results in the relation
$\left\langle{s_{i}^{\prime}s_{i+k}^{\prime}}\right\rangle=\left\langle{s_{i}^{\prime}}\right\rangle\left\langle{s_{i+k}^{\prime}}\right\rangle=0$
for $k>1$. Since $P\\!\left({O}\,|\,{D}\right)$ is a multivariate Gaussian
with 0 mean, we can apply Isserlis’s theorem to solve for all of the moments
of the distribution (25) given knowledge of the second moments of the Fourier
transform on $P\\!\left({O}\,|\,{D}\right)$, which can be expressed as a
covariance matrix. This allows us to express Eq. 29 as

$\displaystyle\Lambda(\bm{\rho})=\exp{\left[-\frac{1}{2}\bm{\rho}^{\intercal}\Sigma\bm{\rho}\right]}.$
(30)

The covariance matrix $\Sigma$ is symmetric tri-diagonal, with non-zero
elements

$\displaystyle\Sigma_{i,i}=\omega_{i}+\varepsilon_{i}+\varepsilon_{i+1}$ (31)
$\displaystyle\Sigma_{i,i+1}=\Sigma_{i+1,i}=-\varepsilon_{i+1}.$

The expression in Eq. 30 is well known as the characteristic function of a
multivariate Gaussian. Substituting Eq. 30 into Eq. 28 and factoring out
$\mathrm{d}{S}\,$ gives

$\displaystyle\mathcal{L}(D)=P\\!\left({S}\,|\,{D}\right)=\frac{1}{\sqrt{(2\pi)^{N-1}\det{\Sigma}}}\exp{\left[-\frac{1}{2}S^{\intercal}\Sigma^{-1}S\right]}.$
(32)

##  6 Implementation

We have presented three independent solutions to the experimental diffusion
likelihood $\mathcal{L}(D)$ (Eq. 10): the recursive method (Sec. 3), the
Laplace method (Sec. 4), and the Markov method (Sec. 5). While each method
requires separate consideration, several features are common to all of the
implementations. The separability of the problem allows us to estimate
diffusion constants for any dimensional inputs inputs using the 1D algorithms
(Eq. 3). The inputs to the algorithms are: (1) the observed particle
locations, $\mathbf{O}=\\{\mathbf{o}_{i}\\}_{i=1}^{N}$; (2) the observation
times $T=\\{t_{i}\\}_{i=1}^{N}$; (3) the measurement variance for each
observation $\mathbf{V}=\\{\mathbf{v}_{i}\\}_{i=1}^{N}$; (4) the exposure of
each frame $t_{\epsilon}$; and (5) one or more diffusion constants $D$ at
which to evaluate the likelihood. The output for each $D$ value is
$\ln(\mathcal{L}(D))$. The logarithm of the likelihood makes the computation
of products and exponentials much faster, and avoids the problem of numerical
underflow for very small values of $\mathcal{L}(D)$. Additionally, because the
logarithm is a strictly monotonically increasing function,
$\operatornamewithlimits{\text{argmax}}_{D}{\mathcal{L}(D)}=\operatornamewithlimits{\text{argmax}}_{D}{\ln(\mathcal{L}(D))}$,
so the maximum likelihood estimate is identical for the log-likelihood.

###  6.1 Recursive Method

The recursive algorithm follows directly from the recursively defined
variables (Eqs. 14 and 15), and the expression of $\mathcal{L}(D)$ as a
product of Gaussians (Eq. 16). The recursive expressions for $\alpha_{i}$,
$\eta_{i}$, and $\mu_{i}$, are causal (the $i$-terms depend only on the
$(i-1)$-terms), enabling their computation in a simple for loop over $N$.
Noting that the logarithm of a normalized Gaussian is

$\ln\,\mathcal{N}(a,b,v)=-\frac{1}{2}\left[\ln(2\pi)+\ln(v)+\frac{(a-b)^{2}}{v}\right],$
(33)

we apply Eq. 33 directly to Eq. 16 to arrive at a computationally efficient
form for the recursive solution of the log-likelihood

$\ln\mathcal{L}(D)=\sum_{i=1}^{N-1}\ln\,\mathcal{N}(o_{i+1},\mu_{i},\alpha_{i})=-\frac{1}{2}\left[(N-1)\ln(2\pi)+\sum_{i=1}^{N-1}\ln(\alpha_{i})+\sum_{i=1}^{N-1}\frac{(o_{i+1}-\mu_{i})^{2}}{\alpha_{i}}\right].$

Of all the methods, the recursive method is the simplest to implement and the
most computationally efficient and numerically stable.

###  6.2 Laplace Method

The computational core of the Laplace method centers around the Hessian matrix
$M$ (Eq. 21). This matrix is symmetric tri-diagonal, which means all non-zero
elements are on the main diagonal and the diagonals immediately above and
below. Using $M$ we can solve the linear system $\widehat{X}=M^{-1}\Theta$
(Eq. 22) to obtain the maximum likelihood estimates $\widehat{X}$ for the true
particle locations. Typically, solving large linear systems is expensive but
since $M$ is tri-diagonal there are algorithms to solve this system in linear
time (26). We refer the reader to Sec. 15 for the details of tri-diagonal
matrix algorithms and our implementation.

Given a solution for $\widehat{X}$, we can use the definition of $f(X)$ in Eq.
17 along with Eq. 33 to compute

$\ln
f(\widehat{X})=-\frac{1}{2}\left[\sum_{i=1}^{N}\ln(2\pi\varepsilon_{i})+\sum_{i=1}^{N}\frac{(o_{i}-\widehat{x}_{i})^{2}}{\varepsilon_{i}}+\sum_{i=1}^{N-1}\ln(2\pi\omega_{i})+\sum_{i=1}^{N-1}\frac{(\widehat{x}_{i+1}-\widehat{x}_{i})^{2}}{\omega_{i}}\right].$
(34)

Finally we can compute the log-likelihood using the Laplace solution of Eq.
20, finding that

$\displaystyle\ln\mathcal{L}(D)$ $\displaystyle=\ln
f(\widehat{X})+\frac{N}{2}\ln(2\pi)-\frac{1}{2}\ln{\det{M}}$ (35)
$\displaystyle=-\frac{1}{2}\left[(N-1)\ln(2\pi)+\sum_{i=1}^{N-1}\ln(\omega_{i})+\sum_{i=1}^{N-1}\frac{(\widehat{x}_{i+1}-\widehat{x}_{i})^{2}}{\omega_{i}}+\sum_{i=1}^{N}\ln(\varepsilon_{i})+\sum_{i=1}^{N}\frac{(o_{i}-\widehat{x}_{i})^{2}}{\varepsilon_{i}}+\ln{\det{M}}\right]$

###  6.3 Markov Method

Finally, the Markov method computation, like the Laplace method, is centered
around matrix computations. In this case, the matrix of interest is the $N-1$
dimensional covariance matrix $\Sigma$ (Eq. 31), which also happens to be
symmetric tri-diagonal, so the same linear-time algorithms used in the Laplace
method are applicable (Sec. 15).

For the Markov method computation we first solve the linear system
$\Phi=\Sigma^{-1}S$, then apply this solution along with the tri-diagonal log-
determinant algorithm to compute the logarithm of the likelihood expression
from Eq. 32, giving

$\ln\mathcal{L}(D)=-\frac{1}{2}\left[(N-1)\ln(2\pi)+S^{\intercal}\Phi+\ln{\det{\Sigma}}\right].$

##  7 Results

To demonstrate the benefits of including individual localization errors in the
estimation, we opt for a loss function to evaluate the quality of the maximum
likelihood estimators (MLE) under consideration. The mean squared error is a
popular choice for evaluating estimators, but it is not a sufficient loss
function for the estimation of $D$, where the mean squared error shows over
penalization for higher estimations as it has a bounded loss at $\hat{D}=0$
and an unbounded loss at $\hat{D}=\infty$. Instead, we will use the squared
log loss function (27)

$\ell(D,\hat{D})=\left(\ln(D)-\ln(\hat{D})\right)^{2}.$

The squared log loss function is similar to the mean squared error, except
that the squared distance is between the logarithms of $D$ and $\hat{D}$. From
(28), we see that the variance term, $D$, scales with the observed data as a
logarithm, hence the choice of the squared log loss function represents a
metric between the expected data given by $D$ and $\hat{D}$.

For one set of SPT simulations, the true trajectory coordinates were first
generated with the pure diffusion model. Full frame motion blur was accounted
for and localization errors were independently and identically drawn from
either a Uniform or a Gamma distribution to test the effects of variable
localization errors without attributing success to a particular choice of
distribution. As a control, a set of simulations with a constant localization
error for all observations was generated to show that the likelihood
distribution represented by the class of diffusion estimators that only
recognize a Scalar Localization Error (SLE) was exactly the same as the
distribution represented by either of our three derived methods which
recognize a Vector of Localization Errors (VLE). We then used the Gamma and
Uniform distribution based localization errors to compare the new VLE based
estimators that account for individually measured localization errors over the
SLE based estimators, for which we had input the square root of the mean
localization variance as the best representation for the scalar error in all
trajectory observations. We performed 10,000 trajectory trials for each data
point to estimate the risk of using a particular estimator, where risk is
defined as $R(D,\hat{D})=\left\langle{\ell(D,\hat{D})}\right\rangle$. The
squared log risk is the variance of $\ln(\hat{D})$ for an unbiased estimator.

The simulations with the gamma distributed localization errors were generated
according to

$\sqrt{V}=\sigma\sim\mathrm{Gamma}(4,\langle{\sqrt{V}}\rangle/4),$

where the standard gamma distribution p.d.f. is
$\mathrm{Gamma}(k,\theta)=\theta^{-k}\sqrt{V}^{k-1}\exp(-\sqrt{V}/\theta)/\Gamma(k)$,
and $\langle{\sqrt{V}}\rangle$ represents the mean error. The simulations with
the uniform distributed localization errors were generated according to

$\sqrt{V}=\sigma\sim\mathrm{Uniform}(\frac{1}{2}\langle{\sqrt{V}}\rangle,\frac{3}{2}\langle{\sqrt{V}}\rangle),$

where the uniform distribution p.d.f. is $\mathrm{Uniform}(a,b)=1/(b-a)$ for
$\sqrt{V}\in[a,b]$ and $\mathrm{Uniform}(a,b)=0$ for all other values of
$\sqrt{V}$.

(A) | (B)
---|---
|
(C) | (D)
|
Figure 1: Comparison of estimation risk of the Vector of Localization Errors
(VLE) and Scalar Localization Error (SLE) MLEs with localization errors drawn
from Gamma or Uniform distributions. The MLE was found with the fminbnd
command from Matlab on the calculated likelihood distributions with a lower
bound of $10^{-8}$. Log-log plots (A) and (C) show the risk of using a
particular estimator on trajectories of various lengths with other constant
parameters set to typical experimental values with standard errors drawn from
gamma (A) or uniform (C) distributions. Log-log plots (B) and (D) show the
risk of using a particular estimator on trajectories of length 50 with various
standard errors drawn from gamma (B) or uniform distributions and all other
parameters were set to 1 to study the effects of relative localization error.
In both (A) and (C) a fiducial line shows that after 30 observations the risk
for the methods decreases approximately $\propto 1/N$ with trajectory length;
in this regime the SLE estimators perform worse by a constant factor (to
simulation precision) relative to the VLE estimators. In (B) and (D) the risk
for SLE estimators increase at a faster rate than the VLE estimators,
indicating that the localization errors provide valuable information for
improving estimator reliability.

Figures 1A and C are shown in log space to show that in the presence of
sufficient information the risk of the VLE estimators is less than the risk of
the SLE estimators by a constant proportionality factor when a scalar
localization error is used to parameterize a continuous distribution of
localization errors. In other words, for a given parameterized scalar
localization error, the amount of observations required to generate the same
quality $\hat{D}$ estimate is always less by a proportional factor for the VLE
estimators method compared to the SLE estimators. Figures 1B and D are shown
in log space to show how each estimator begins to fail in the presence of
increasing relative localization error given a fixed set of trajectory
observations (50). For these subplots, the VLE estimators show a noticeable
improvement in estimator reliability when the relative error is equal to or
greater than the true underlying $D$.

In an experimental trajectory, the distribution that parameterizes the
localization error is typically a function of several environmental variables
so that it can often appear arbitrary or specific to a particular experimental
trial. In this manuscript, we focus on two simple distributions to provide a
fair metric for validation over thousands of simulated trajectories. It is
worthwhile to further investigate the the precision increase of the VLE
estimators over their SLE counterparts with localization error distributions
of varying error variances. We do so by performing trials on trajectories with
localization errors parameterized by the Gamma and Uniform distribution, but
this time we vary the parameters that characterize the variance of these
distributions without altering the value of the mean localization variance. To
do so, we run simulations where the shape parameter, k, of the gamma
distribution is altered so that our expression looks like

$\sqrt{V}=\sigma\sim\mathrm{Gamma}(k,\langle{\sqrt{V}}\rangle/k),$

and the bounds of the uniform distribution are altered so the expression
becomes

$\sqrt{V}=\sigma\sim\mathrm{Uniform}([1-b]\langle{\sqrt{V}}\rangle,[1+b]\langle{\sqrt{V}}\rangle).$

(A) | (B)
---|---
|
Figure 2: Risk ratio of the Scalar Localization Error (SLE) and Vector of
Localization Errors (VLE) MLEs for various trajectories parameterized by
localization errors drawn from Gamma and Uniform distributions of a single
varying parameter but the same mean $\langle\sqrt{V}\rangle$. Log-log plot (A)
shows the risk ratio of SLE and VLE with 5 Gamma distributions of the same
mean localization error and different shape parameters, k. The increased value
of k reduces the variance of the gamma distribution; k = 1 is an exponential
distribution and k = 13 is approaching a gaussian distribution. All of the
gamma distributions in (A) have a variance of $\langle\sqrt{V}\rangle^{2}/k$.
Log-log plot (B) shows the risk ratio of SLE and VLE with 5 Uniform
distributions of the same localization variance and different sampling
boundaries, b. The reduced value of b reduces the variance of the uniform
distribution; a given b in plot (B) corresponds to a variance of
$\left(\langle\sqrt{V}\rangle b\right)^{2}/3$, hence b = 0.1 is nearly a
variance of 0. As long as the variance of $\langle\sqrt{V}\rangle$ is greater
than 0, the VLE estimators outperform the SLE estimators.

We see in Fig. 2 that the effect of increasing the variance relative to the
mean of the localization errors in these test distributions results in a
growing disparity between the two classes of estimators. From Fig. 1 the
effect of increasing $\langle\sqrt{V}\rangle$ sets a minimum trajectory length
where estimates become reasonable; e.g the linear decrease in risk for the
gamma distribution of Fig. 1A is seen at shorter trajectories than the uniform
distribution of Fig. 1C even though both the gamma and uniform distributions
have the same variance because the uniform distribution has a larger
$\langle\sqrt{V}\rangle$. In Fig. 2 the constant precision improvement of the
VLE estimators are recognized in the regime where both estimators start
reporting risk values that scale linearly with trajectory length. Prior to
that, the VLE estimators start performing significantly better at shorter
trajectories, which is represented by the peaks seen in Fig. 2. These
simulations, while overly simplified versions of real SPT trajectories,
highlight the importance of characterizing each localization error accordingly
to their associated localizations in a trajectory.

##  8 Discussion and Conclusion

Starting from the fundamental diffusion likelihood expression (Eq. 9), we
presented three independent solutions, each of which has different benefits,
and leads to a computational algorithm with different advantages. The
recursive method presents the simplest solution that is numerically more
stable than the other methods for estimating likelihoods when $D$ approaches
0. The Laplace method has the advantage that the expected true positions
$\hat{X}$ are computed along with the likelihood, which may be useful in some
applications. The Markov method was crucial for deriving the terms
$\epsilon_{i}$ in the components $\mathcal{M}_{i}$ given knowledge of the true
underlying probability distribution, hence the generality of the Markov method
is its main advantage. In terms of numerical accuracy and computational
efficiency the Markov method is better than the Laplace method especially for
very small $D$ values, but it remains computationally inferior to the
recursive method. For practical implementations, we recommend the recursive
method unless the MLE of the true positions is also desired.

The method described here naturally allows for trajectories with missing or
irregularly spaced localizations by decoupling the concept of observation
times $t_{i}$ from the exposure time $t_{\epsilon}$. This is important when
some localizations are missed because the gap between $t_{i}$ and $t_{i+1}$
becomes larger, but the exposure time remains the same. Thus to correctly
account for the motion-blur, the effective weighting of the dependence of
observed position $o_{i}$ on $x_{i}$ and $x_{i+1}$ changes, and our technique
directly incorporates this effect. Trajectory intermittency has been accounted
for in prior studies (19) and in those same studies extensions to dynamic
errors were suggested, but a convenient computational framework for an
estimator that seamlessly factors in both trajectory intermittencies and
dynamic error had not been explicitly worked out until now.

Numerical implementations of the likelihood forms resulting from the three
derivations were tested to justify the equivalence among the likelihood forms.
Since all three derivations began from the same set of first principles, the
three likelihood calculations are essentially equivalent. We note however that
our implementation remains unit agnostic, so that the trajectory with a very
small $D$ value can be easily scaled up to appropriate units so that the value
of $D$ is closer to 1, where the numerical calculations will be more robust.

Variable localization uncertainties could occur in practice from variable
background intensities or photobleaching of the fluorescent label. We compared
the performance of our VLE estimator to the current state-of-the-art SLE
estimator using the squared log loss function and found a clear performance
benefit when trajectories had variable localization uncertainties.

##  9 Author Contributions

PKR and MJO contributed equally to this work. KAL, PJC and PKR conceived the
project. PJC initiated the formulation of the estimation problem as a MLE
given a set of observations. PKR derived the recursive, Markov, and Laplace
methods. MJO derived the efficient algorithms for the three solution methods
and the C++ and MATLAB implementations of the methods, generated the
estimation accuracy results and helped to simplify the presentation. All
authors contributed to the writing and editing of the manuscript.

##  10 Acknowledgments

We wish to acknowledge Stan Steinberg and Michael Wester for reading our
manuscript and providing enlightening discussions and helpful comments.
Financial support for this work was provided primarily by the National Science
Foundation grant 0954836. Additional support was provided by The New Mexico
Spatiotemporal Modeling Center: NIH P50GM085273 (KAL), NIH grant 1R01GM100114
(KAL,MJO) and NIH grant 1R01NS071116 (KAL, MJO).

## References

  * (1) Michael J Saxton and Ken Jacobson.  SINGLE-PARTICLE TRACKING:Applications to Membrane Dynamics.  Annual Review of Biophysics and Biomolecular Structure, 26:373–399, 1997.
  * (2) Michael J Saxton.  Two-dimensional continuum percolation threshold for diffusing particles of nonzero radius.  Biophysical journal, 99(5):1490–1499, 2010.
  * (3) Arash Sanamrad, Fredrik Persson, Ebba G Lundius, David Fange, Arvid H Gynnå, and Johan Elf.  Single-particle tracking reveals that free ribosomal subunits are not excluded from the escherichia coli nucleoid.  Proceedings of the National Academy of Sciences, 111(31):11413–11418, 2014.
  * (4) Christopher P Calderon, Michael A Thompson, Jason M Casolari, Randy C Paffenroth, and WE Moerner.  Quantifying transient 3d dynamical phenomena of single mrna particles in live yeast cell measurements.  The Journal of Physical Chemistry B, 117(49):15701–15713, 2013\.
  * (5) Nilah Monnier, Syuan-Ming Guo, Masashi Mori, Jun He, Péter Lénárt, and Mark Bathe.  Bayesian approach to msd-based analysis of particle motion in live cells.  Biophysical journal, 103(3):616–626, 2012.
  * (6) Khuloud Jaqaman, Dinah Loerke, Marcel Mettlen, Hirotaka Kuwata, Sergio Grinstein, Sandra L Schmid, and Gaudenz Danuser.  Robust single-particle tracking in live-cell time-lapse sequences.  Nature methods, 5(8):695–702, 2008.
  * (7) Arnauld Sergé, Nicolas Bertaux, Hervé Rigneault, and Didier Marguet.  Dynamic multiple-target tracing to probe spatiotemporal cartography of cell membranes.  Nature Methods, 5(8):687–694, 2008.
  * (8) Nicolas Chenouard, Ihor Smal, Fabrice De Chaumont, Martin Maška, Ivo F Sbalzarini, Yuanhao Gong, Janick Cardinale, Craig Carthel, Stefano Coraluppi, Mark Winter, et al.  Objective comparison of particle tracking methods.  Nature methods, 2014.
  * (9) Alexander D Mont, Christopher P Calderon, and Aubrey B Poore.  A new computational method for ambiguity assessment of solutions to assignment problems arising in target tracking.  In SPIE Defense+ Security, pages 90920J–90920J. International Society for Optics and Photonics, 2014.
  * (10) Ji Won Yoon, Andreas Bruckbauer, William J Fitzgerald, and David Klenerman.  Bayesian inference for improved single molecule fluorescence tracking.  Biophysical journal, 94(12):4932–4947, 2008.
  * (11) Carlas S Smith, Nikolai Joseph, Bernd Rieger, and Keith a Lidke.  Fast, single-molecule localization that achieves theoretically minimum uncertainty.  Nature methods, 7(5):373–5, May 2010.
  * (12) H Qian, M P Sheetz, and E L Elson.  Single particle tracking. Analysis of diffusion and flow in two-dimensional systems.  Biophysical journal, 60(4):910–21, October 1991.
  * (13) Douglas S Martin, Martin B Forstner, and Josef A Käs.  Apparent subdiffusion inherent to single particle tracking.  Biophysical journal, 83(4):2109–2117, 2002.
  * (14) Thierry Savin and Patrick S Doyle.  Static and dynamic errors in particle tracking microrheology.  Biophysical journal, 88(1):623–38, January 2005.
  * (15) Xavier Michalet.  Mean square displacement analysis of single-particle trajectories with localization error: Brownian motion in an isotropic medium.  Phys Rev E Stat Nonlin Soft Matter Phys, 82(4):041914, October 2010\.
  * (16) Xavier Michalet and Andrew J Berglund.  Optimal diffusion coefficient estimation in single-particle tracking.  Physical review. E, Statistical, nonlinear, and soft matter physics, 85(6 Pt 1):061916, June 2012.
  * (17) D Montiel, H Cang, and H Yang.  Quantitative characterization of changes in dynamical behavior for single-particle tracking studies.  Journal of Physical Chemistry B, 110:19763–70, 2006.
  * (18) A J Berglund.  Statistics of camera-based single-particle tracking.  Phys Rev E Stat Nonlin Soft Matter Phys, 82:011917, 2010.
  * (19) Bo Shuang, Chad P Byers, Lydia Kisley, Lin-Yung Wang, Julia Zhao, Hiroyuki Morimura, Stephan Link, and Christy F Landes.  Improved analysis for determining diffusion coefficients from short, single-molecule trajectories with photoblinking.  Langmuir : the ACS journal of surfaces and colloids, 29(1):228–34, January 2013.
  * (20) Christian L Vestergaard, Paul C Blainey, and Henrik Flyvbjerg.  Optimal estimation of diffusion coefficients from single-particle trajectories.  Physical Review E, 89(2):022726, 2014.
  * (21) PS de Laplace.  Mémoire sur les suites récurro-récurrentes et sur leurs usages dans la théorie des hasards.  Mém. Acad. Roy. Sci. Paris, 6:353–371, 1774.
  * (22) A.A. Markov.  Wahrscheinlichkeitsrechnung.  BG Teubner, 1912.
  * (23) S. Chandrasekhar.  Stochastic problems in physics and astronomy.  Rev. Mod. Phys., 15:1–89, Jan 1943.
  * (24) Kevin Cahill.  Physical Mathematics.  Cambridge University Press, 2013.
  * (25) L. Isserlis.  On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables.  Biometrika, 12(1/2):pp. 134–139, 1918.
  * (26) Moawwad EA El-Mikkawy.  On the inverse of a general tridiagonal matrix.  Applied Mathematics and Computation, 150(3):669–679, 2004.
  * (27) L Brown.  Inadmissibility of the usual estimators of scale parameters in problems with unknown location and scale parameters.  The Annals of Mathematical Statistics, pages 29–48, 1968.
  * (28) Andrew Gelman, John B Carlin, Hal S Stern, and Donald B Rubin.  Bayesian data analysis, volume 2.  Taylor & Francis, 2014.
  * (29) Sheldon M Ross et al.  Stochastic processes, volume 2.  John Wiley & Sons New York, 1996.
  * (30) Conte de Boor.  Elementary numerical analysis.  McGraw-Hill, 1972.
  * (31) Moawwad E.A. El-Mikkawy.  On the inverse of a general tridiagonal matrix.  Applied Mathematics and Computation, 150(3):669 – 679, 2004.

SUPPORTING MATERIAL: Estimation of the Diffusion Constant from Intermittent
Trajectories with Variable Position Uncertainties

###### Contents

  1. 1 Introduction
    1. 1.1 Background and Related Work
  2. 2 Theory
    1. 2.1 Accounting for the effects of exposure time integration
  3. 3 Recursive Method
  4. 4 Laplace Method
  5. 5 Markov Method
  6. 6 Implementation
    1. 6.1 Recursive Method
    2. 6.2 Laplace Method
    3. 6.3 Markov Method
  7. 7 Results
  8. 8 Discussion and Conclusion
  9. 9 Author Contributions
  10. 10 Acknowledgments
  11. 11 Implemented Notation and Gaussian Functions
  12. 12 Problem Formulation
    1. 12.1 The Probability of a Single Observation and the True Start Positions
    2. 12.2 The Probability of a Trajectory of Observations
  13. 13 Explicit Derivations on the Problem Formulation
    1. 13.1 The Probability Density of a Time Averaged Position
    2. 13.2 Functional Form: From Products of Probability Components to Simpler Expressions
  14. 14 Method Component Derivations
    1. 14.1 Laplace Method: Maximum Likelihood of True Positions
    2. 14.2 Laplace Method: Direct Integration of the Probability Distribution Components
    3. 14.3 Markov Method: Marginalizing the Likelihood Function
    4. 14.4 Markov Method: Expectation Calculations
  15. 15 Implementation
    1. 15.1 Log–product Computation
    2. 15.2 Computation of Variance terms
    3. 15.3 The Recursive Method Algorithm
    4. 15.4 Tri-Diagonal matrix algorithms
    5. 15.5 DST Algorithm

##  11 Implemented Notation and Gaussian Functions

The derivations presented in the main text make heavy use of normalized
Gaussian functions, which we denote as a function of three arguments,

$\,\mathcal{N}(a,b,v)=\frac{1}{\sqrt{2\pi
v}}\exp{\left[-\frac{(a-b)^{2}}{2v}\right]}.$

The Gaussian function defined this way is symmetric with respect to the first
two position parameters, so that $\,\mathcal{N}(a,b,v)=\,\mathcal{N}(b,a,v)$,
and the variance of the function is given by the third parameter. Also, the
normalization factor ensures that the Gaussian integrated over all space with
respect to either of its position parameters is unity,

$\int^{\infty}_{-\infty}\mathrm{d}{a}\,\,\,\mathcal{N}(a,b,v)=\int^{\infty}_{-\infty}\mathrm{d}{b}\,\,\,\mathcal{N}(a,b,v)=1.$

As a corollary, a normalized Gaussian has a useful scaling identity, for
$q>0$,

$\,\mathcal{N}(a,b,v)=q\,\mathcal{N}(qa,qb,q^{2}v).$

Next, consider the case of the product of two normalized Gaussians sharing a
common position parameter, which can be rewritten as a product of two
normalized Gaussians where the common parameter only appears in one of the
two,

$\displaystyle\,\mathcal{N}(x,\mu_{1},\eta_{1})\,\mathcal{N}(x,\mu_{2},\eta_{2})$
$\displaystyle=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2})\,\mathcal{N}(x,\mu^{\prime},\eta^{\prime}),$
(36)

where,

$\mu^{\prime}=\frac{\mu_{1}\eta_{2}+\mu_{2}\eta_{1}}{\eta_{1}+\eta_{2}}\quad\mbox{and,}\quad\eta^{\prime}=\frac{\eta_{1}\eta_{2}}{\eta_{1}+\eta_{2}}.$

Using Eq. 36, the integral of the product of two normalized Gaussians over a
shared position parameter is itself a normalized Gaussian in the other two
(unintegrated) position parameters,

$\displaystyle\int\mathrm{d}{x}\,\,\mathcal{N}(x,\mu_{1},\eta_{1})\,\mathcal{N}(x,\mu_{2},\eta_{2})=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2})\int\mathrm{d}{x}\,\,\mathcal{N}(x,\mu^{\prime},\eta^{\prime})=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2}).$
(37)

With two successive applications of Eq. 36, the product of three Gaussians
which share a common position parameter becomes,

$\displaystyle\,\mathcal{N}(x,\mu_{1},\eta_{1})\,\mathcal{N}(x,\mu_{2},\eta_{2})\,\mathcal{N}(x,\mu_{3},\eta_{3})$
$\displaystyle=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2})\,\mathcal{N}(x,\mu^{\prime},\eta^{\prime})\,\mathcal{N}(x,\mu_{3},\eta_{3}),$
(38)
$\displaystyle=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2})\,\mathcal{N}(\mu_{3},\mu^{\prime},\eta^{\prime}+\eta_{3})\,\mathcal{N}(x,\mu^{\prime\prime},\eta^{\prime\prime}),$
$\displaystyle\mbox{where,}\quad\mu^{\prime\prime}$
$\displaystyle=\frac{\mu^{\prime}\eta_{3}+\mu_{3}\eta^{\prime}}{\eta^{\prime}+\eta_{3}},$
$\displaystyle\mbox{and,}\quad\eta^{\prime\prime}$
$\displaystyle=\frac{\eta^{\prime}\eta_{3}}{\eta^{\prime}+\eta_{3}}.$

Finally, Eq. 38 allows the integral of three normalized Gaussians over a
shared position parameter to reduce to

$\displaystyle\int\mathrm{d}{x}\,\,\mathcal{N}(x,\mu_{1},\eta_{1})\,\mathcal{N}(x,\mu_{2},\eta_{2})\,\mathcal{N}(x,\mu_{3},\eta_{3})$
$\displaystyle=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2})\,\mathcal{N}(\mu_{3},\mu^{\prime},\eta^{\prime}+\eta_{3})\int\mathrm{d}{x}\,\,\mathcal{N}(x,\mu^{\prime\prime},\eta^{\prime\prime})$
(39)
$\displaystyle=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2})\,\mathcal{N}(\mu_{3},\mu^{\prime},\eta^{\prime}+\eta_{3})$
$\displaystyle=\,\mathcal{N}(\mu_{1},\mu_{2},\eta_{1}+\eta_{2})\,\mathcal{N}(\mu_{3},\mu^{\prime},\gamma),$
$\displaystyle\mbox{where,}\quad\gamma$
$\displaystyle=\frac{\eta_{1}\eta_{2}+\eta_{1}\eta_{3}+\eta_{2}\eta_{3}}{\eta_{1}+\eta_{2}}.$

##  12 Problem Formulation

In the main manuscript, the likelihood distribution for a 1D random walk given
a set of observations, $O=\\{o_{i}\\}_{i=1}^{N}$, is described as an integral
that marginalizes over the $N+1$ unknown true positions
$X=\\{x_{i}\\}_{i=1}^{N+1}$,

$P\\!\left({O}\,|\,{D}\right)=\int_{\mathbb{R}^{N+1}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right).$
(40)

We seek to rigorously define our formalism to justify the analysis used in the
rest of the manuscript.

###  12.1 The Probability of a Single Observation and the True Start Positions

Under typical biological SPT experiments, the temporal resolution of a probe
is high enough that the effects of diffraction are much greater than the
effects of particle motion in generating the point spread function of the
image. Following the assumption that the variance due to diffraction is
significantly greater than the variance due to motion, the point spread
function can be defined as a stationary Gaussian function centered about the
average position of the particle during a single frame with some background
offset. With these considerations, we now set to define
$P\\!\left({o_{i},x_{i},x_{i+1}}\,|\,{D}\right)$: the probability of an
observation and the particle’s true start and end points in a frame given the
free diffusion model. Since the maximum likelihood estimator for a Gaussian
function returns the peak (the maximum likelihood) and variance (the error) of
a gaussian distribution, there is sufficient information from the estimator to
build a probability distribution relating the localization to the true
averaged position of a particle. The probability of obtaining a localized
position, $o_{i}$, given the true averaged position of the particle,
$\bar{y}_{i}$, and the estimator variance, $v_{i}$, is

$P\\!\left({o_{i}}\,|\,{\bar{y}_{i}}\right)=\,\mathcal{N}(o_{i},\bar{y}_{i},v_{i}).$
(41)

We incorporate more information on $\bar{y}_{i}$ by relating it to the start
positions at the considered frame, $i$, and the subsequent frame, $i+1$; where
it is assumed that a frame begins immediately after the prior frame ends.
Therefore the probability distribution of $o_{i}$ with frame start
coordinates, $x_{i}$ and $x_{i+1}$, is expressed as

$P\\!\left({o_{i},x_{i},x_{i+1}}\,|\,{D}\right)=\int\mathrm{d}{\bar{y}_{i}}\,\;P\\!\left({o_{i}}\,|\,{\bar{y}_{i}}\right)P\\!\left({\bar{y}_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right)P(x_{i}),$

where the diffusion constant, D, is implicitly included in the probability
distributions. The probabilities $P\\!\left({x_{i+1}}\,|\,{x_{i}}\right)$
represent a pure diffusion process.

$P\\!\left({x_{i+1}}\,|\,{x_{i}}\right)=\,\mathcal{N}(x_{i+1},x_{i},\omega_{i}(D)),$
(42)

where $\omega_{i}(D)=2D\delta t_{i}$ and $\delta t_{i}=t_{i+1}-t_{i}$.

In section 13, we will explicitly derive
$P\\!\left({\bar{y}_{i}}\,|\,{x_{i},x_{i+1}}\right)$, but here we will state
the result as

$P\\!\left({\bar{y}_{i}}\,|\,{x_{i},x_{i+1}}\right)=\,\mathcal{N}\left(\bar{y}_{i},\left(1-\frac{t_{\epsilon}}{2\delta
t_{i}}\right)x_{i}+\left(\frac{t_{\epsilon}}{2\delta
t_{i}}\right)x_{i+1},2Dt_{\epsilon}\left[\frac{1}{3}-\frac{t_{\epsilon}}{4\delta
t_{i}}\right]\right),$

where $t_{\epsilon}$ is the exposure time of a frame; in other words, the time
the last photon observed in one frame can be arbitrarily spaced from the first
photon in the next frame. Hence, trajectory intermittencies can be accounted
for by redefining $\delta t_{i}$ so that frame $i+1$ is the next frame that
observes a photon from the particle under consideration, omitting all frames
that do not provide measurement information. Given that
$P\\!\left({\bar{y}_{i}}\,|\,{x_{i},x_{i+1}}\right)$ is in the form of a
Gaussian function with $\bar{y}_{i}$ as one of the location parameters,
$\bar{y}_{i}$ can be effectively marginalized in
$P\\!\left({o_{i},x_{i},x_{i+1}}\,|\,{D}\right)$ so that the expression
reduces to

$P\\!\left({o_{i},x_{i},x_{i+1}}\,|\,{D}\right)=P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right)P(x_{i}).$

When obtaining the global probability of all $O$ given $D$, the probabilities
$P(x_{i})$ will be conditioned on prior observations, so it becomes necessary
to define the expression

$P\\!\left({o_{i},x_{i+1}}\,|\,{x_{i},D}\right)=P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right).$

###  12.2 The Probability of a Trajectory of Observations

We now wish to get the full expression for $P\\!\left({O}\,|\,{D}\right)$,
which has no dependency on the true positions, $X$. To do this, first define
the expression

$P\\!\left({O,X}\,|\,{D}\right)=P(x_{1})\prod^{N}_{i=1}P\\!\left({o_{i},x_{i+1}}\,|\,{x_{i},D}\right)=P(x_{1})\prod^{N}_{i=1}P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right).$

Since $X$ are nuisance parameters, integrate the $X$ variables over all
configuration space and set $P(x_{1})=1$ because it is assumed that $x_{1}$
must already be known, since the inference of a diffusion probability must
have an origin to relate all subsequent coordinates.

$P\\!\left({O}\,|\,{D}\right)=\int_{\mathbb{R}^{N+1}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i+1}}\,|\,{x_{i}}\right).$
(43)

We will see in section 13 that the following expression, which is the focus of
the main manuscript, is equivalent and easier to evaluate

$P\\!\left({O}\,|\,{D}\right)=\int_{\mathbb{R}^{N}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i}(D))\prod^{N-1}_{j=1}\,\mathcal{N}(x_{i+1},x_{i},\omega_{i}(D))=\int_{\mathbb{R}^{N}}\\!\\!\mathrm{d}{X}\,\prod^{N}_{i=1}\mathcal{M}_{i}\prod^{N-1}_{j=1}\mathcal{T}_{j}.$
(44)

where

$\displaystyle\mathcal{M}_{i}=\mathcal{M}_{i}(o_{i},x_{i})=$
$\displaystyle\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i}(D)),\quad\textrm{for
}1\leq i\leq N,\textrm{ and }$ (45)
$\displaystyle\mathcal{T}_{i}=\mathcal{T}_{i}(x_{i+1},x_{i})=$
$\displaystyle\,\mathcal{N}(x_{i+1},x_{i},\omega_{i}(D)),\quad\textrm{for
}1\leq i\leq N-1.$ (46)

##  13 Explicit Derivations on the Problem Formulation

There were a few results presented in section 12 that are clarified in this
section. The probability distribution of an averaged position given the start
and end points of a frame are discussed and analytically solved in the
continuous limit. Then the relation between representation of
$P\\!\left({O}\,|\,{D}\right)$ with the probability distribution and reduced
measurement functions is investigated.

###  13.1 The Probability Density of a Time Averaged Position

Given $D$, the probability density of a transition from point $a$ to to point
$b$ separated by a time $T$ is

$P\\!\left({b}\,|\,{a}\right)=\,\mathcal{N}(b,a,2DT).$

If an intermediate point, $y(t)$ sampled at a time $t<T$ is considered, the
joint probability density of a transition from $a$ to $y(t)$ and then from
$y(t)$ to $b$ is

$\displaystyle P(y(t),a,b)$
$\displaystyle=P(a)\,\mathcal{N}(y(t),a,2Dt)\,\mathcal{N}(b,y(t),2D(T-t))$
$\displaystyle=P(a)\,\mathcal{N}(b,a,2DT)\,\mathcal{N}\left(y(t),a\left(1-\frac{t}{T}\right)+b\frac{t}{T},2D\frac{t}{T}(T-t)\right).$

The probability density of the variable $y(t)$, preconditioned on the end
points, a and b, is then

$P\\!\left({y(t)}\,|\,{b,a}\right)=\,\mathcal{N}\left(y(t),a\left(1-\frac{t}{T}\right)+b\frac{t}{T},2D\frac{t}{T}(T-t)\right).$
(47)

Eq. 47 is the probability density for what is known as a Brownian Bridge (29),
or Brownian motion with preconditioned end points. It has a mean and
covariance defined as

$\displaystyle\left\langle{y(t)}\right\rangle$
$\displaystyle=a\left(1-\frac{t}{T}\right)+b\frac{t}{T}$
$\displaystyle\text{cov}[y(t),y(s)]$
$\displaystyle=2D\left(s-\frac{st}{T}\right)\quad\text{for}\quad s\leq t\leq
T$

It is now of interest to find the probability density for a quantity that
describes an integrated average of $y(t)$ such that

$\bar{y}=\frac{1}{t_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,y(t).$

Since each $y(t)$ is a normally distributed random variable, y(t) is a
gaussian process, then the time averaged integral of a Gaussian process,
$\bar{y}$, is also a normally distributed random variable. Therefore, from
Isserlis theorem (25) only the first two moments of $\bar{y}$ are needed to
determine its probability distribution. The first moment is

$\displaystyle\left\langle{\bar{y}}\right\rangle$
$\displaystyle=\left\langle{\frac{1}{t_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,y(t)}\right\rangle=\frac{1}{t_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\left\langle{y(t)}\right\rangle$
$\displaystyle=\frac{1}{t_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\left[a\left(1-\frac{t}{T}\right)+b\frac{t}{T}\right]=a\left(1-\frac{t_{\epsilon}}{2T}\right)+b\frac{t_{\epsilon}}{2T}$
$\displaystyle=(1-\alpha)a+\alpha b,$

where for notational convenience, we define $\alpha=t_{\epsilon}/2T$. It
follows that the second moment is

$\displaystyle\left\langle{\bar{y}^{2}}\right\rangle$
$\displaystyle=\left\langle{\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,y(t)\int^{t_{\varepsilon}}_{s=0}\mathrm{d}{s}\,y(s)}\right\rangle=\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\int^{t_{\varepsilon}}_{s=0}\mathrm{d}{s}\,\left\langle{y(t)y(s)}\right\rangle$
$\displaystyle=\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\int^{t_{\varepsilon}}_{s=0}\mathrm{d}{s}\,\left\langle{[y(t)-\left\langle{y(t)}\right\rangle+\left\langle{y(t)}\right\rangle][y(s)-\left\langle{y(s)}\right\rangle+\left\langle{y(s)}\right\rangle]}\right\rangle$
$\displaystyle=\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\int^{t_{\varepsilon}}_{s=0}\mathrm{d}{s}\,\text{cov}[y(t),y(s)]+\left\langle{y(t)}\right\rangle\left\langle{y(s)}\right\rangle$
$\displaystyle=\left\langle{\bar{y}}\right\rangle^{2}+\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\int^{t_{\varepsilon}}_{s=0}\mathrm{d}{s}\,\text{cov}[y(t),y(s)].$

Therefore, the variance for $\bar{y}$ is

$\displaystyle\left\langle{\bar{y}^{2}}\right\rangle-\left\langle{\bar{y}}\right\rangle^{2}$
$\displaystyle=\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\int^{t_{\varepsilon}}_{s=0}\mathrm{d}{s}\,\text{cov}[y(t),y(s)]$
$\displaystyle=\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,\left[\int^{t_{\varepsilon}}_{s=t}\mathrm{d}{s}\,2D\left(t-\frac{st}{T}\right)+\int^{t}_{s=0}\mathrm{d}{s}\,2D\left(s-\frac{st}{T}\right)\right]$
$\displaystyle=\frac{1}{t^{2}_{\varepsilon}}\int^{t_{\varepsilon}}_{t=0}\mathrm{d}{t}\,2D\left[t_{\epsilon}t-t^{2}-\frac{t^{2}_{\epsilon}t}{2T}+\frac{t^{2}}{2}\right]$
$\displaystyle=2D\left[\frac{t_{\epsilon}}{2}-\frac{t_{\epsilon}}{3}-\frac{t^{2}_{\epsilon}}{4T}+\frac{t_{\epsilon}}{6}\right]=2Dt_{\epsilon}\left[\frac{1}{3}-\frac{\alpha}{2}\right].$

Given that $\bar{y}$ must be a normally distributed and its first two moments
are known, then

$P\\!\left({\bar{y}}\,|\,{a,b}\right)=\,\mathcal{N}(\bar{y},\left\langle{\bar{y}}\right\rangle,\left\langle{\bar{y}^{2}}\right\rangle-\left\langle{\bar{y}}\right\rangle^{2})=\,\mathcal{N}\left(\bar{y},(1-\alpha)a+\alpha
b,2Dt_{\varepsilon}\left[\frac{1}{3}-\frac{\alpha}{2}\right]\right)$

Furthermore, it was established in Eq. 41 that a time averaged position was
related to a localized observation by a normal Gaussian function so that

$P\\!\left({o}\,|\,{a,b}\right)=\int\mathrm{d}{\bar{y}}\,P\\!\left({o}\,|\,{\bar{y}}\right)P\\!\left({\bar{y}}\,|\,{a,b}\right)=\,\mathcal{N}\left(o,(1-\alpha)a+\alpha
b,v+2Dt_{\varepsilon}\left[\frac{1}{3}-\frac{\alpha}{2}\right]\right)$ (48)

###  13.2 Functional Form: From Products of Probability Components to Simpler
Expressions

In the limit where the camera exposure time goes to 0, the probability of
$o_{i}$ is dependent on one coordinate, $x_{i}$. However, as the camera
exposure time becomes non-negligible with respect to the time spacing between
frames, the probability of $o_{i}$ becomes increasingly dependent on the
subsequent coordinate, $x_{i+1}$. The $o_{i}$ dependence on both $x_{i}$ and
$x_{i+1}$ make a direct approach to solving the integral computationally
difficult. In the main manuscript, Markov’s method approach showed the
following relationship for a multivariate gaussian

$\Sigma_{i,j}=\left\langle{s_{i}s_{j}}\right\rangle,$

where $\Sigma$ is the covariance matrix of a multivariate Gaussian function
describing the vector of random variables
$S=\\{s_{i}=o_{i+1}-o_{i}\\}_{i=1}^{N-1}$. Given that our probability
distribution is obtained by integrating several Gaussian functions, the result
of the distribution is a Gaussian function. Therefore, if the moments of $S$
are known, the parameter $\epsilon_{i}$ for the simpler expression can be
derived given that

$\left\langle{s_{i}s_{i+1}}\right\rangle=-\epsilon_{i+1}$

was previously shown to be true for the simpler expression in the main
manuscript. Starting from our derived probability expression in Eq. 48 for an
arbitrary $o_{i}$

$P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)=\,\mathcal{N}\left(o_{i},(1-\alpha_{i})x_{i}+\alpha_{i}x_{i+1},v_{i}+2Dt_{\epsilon}\left[\frac{1}{3}-\frac{\alpha_{i}}{2}\right]\right)$

where $v_{i}$ is defined as the localization variance and

$\alpha_{i}=\frac{t_{\epsilon}}{2\delta t_{i}}.$

From the properties of a Gaussian function with 0 mean

$\displaystyle\left\langle{o_{i}-(1-\alpha_{i})x_{i}-\alpha_{i}x_{i+1}}\right\rangle$
$\displaystyle=0,$
$\displaystyle\left\langle{(o_{i}-(1-\alpha_{i})x_{i}-\alpha_{i}x_{i+1})^{2}}\right\rangle$
$\displaystyle=v_{i}+2Dt_{\epsilon}\left[\frac{1}{3}-\frac{\alpha_{i}}{2}\right].$
(49)

Which implies

$\displaystyle\left\langle{o_{i}-x_{i}}\right\rangle$ $\displaystyle=0$
$\displaystyle\left\langle{o_{i}-x_{i}|X}\right\rangle$
$\displaystyle=\alpha_{i}(x_{i+1}-x_{i})$
$\displaystyle\left\langle{(o_{i}-x_{i})^{2}}\right\rangle$
$\displaystyle=v_{i}+2Dt_{\epsilon}\frac{1}{3}.$ (50)

Also

$\left\langle{s_{i}}\right\rangle=\left\langle{o_{i+1}-o_{i}}\right\rangle=\left\langle{o_{i+1}-x_{i+1}}\right\rangle-\left\langle{o_{i}-x_{i}}\right\rangle+\left\langle{x_{i+1}-x_{i}}\right\rangle=0.$
(51)

Given the relations in Eq. 13.2 and Eq. 51

$\displaystyle\left\langle{s_{i}s_{i+1}}\right\rangle$
$\displaystyle=\left\langle{(o_{i+2}-o_{i+1})(o_{i+1}-o_{i})}\right\rangle=-\epsilon_{i+1}$
$\displaystyle=\left\langle{(o_{i+2}-o_{i+1}+x_{i+2}-x_{i+2}+x_{i+1}-x_{i+1})(o_{i+1}-o_{i}+x_{i+1}-x_{i+1}+x_{i}-x_{i})}\right\rangle$
$\displaystyle=\left\langle{(o_{i+1}-x_{i})(x_{i+1}-x_{i})}\right\rangle-\left\langle{(o_{i+1}-x_{i+1})^{2}}\right\rangle+\cdots$
$\displaystyle=\alpha_{i+1}(2D\delta
t_{i+1})-2Dt_{\epsilon}\frac{1}{3}-v_{i+1}=2Dt_{\epsilon}\frac{1}{6}-v_{i+1}.$

Where all the other terms in the ellipsis ($\cdots$) go to 0. Additionally, it
follows that

$\left\langle{s_{i}^{2}}\right\rangle=w_{i}+\epsilon_{i}+\epsilon_{i+1}=2D\delta
t_{i}-4Dt_{\epsilon}\frac{1}{6}+v_{i}+v_{i+1}.$

Therefore

$\varepsilon_{i}(D)=v_{i}-2Dt_{\epsilon}\frac{1}{6}.$

Which is the variance correction discovered in earlier diffusion estimation
papers (14, 17, 18).

##  14 Method Component Derivations

The following sub-sections explain some of the relations that were explicitly
stated to complete the derivations in the main text.

###  14.1 Laplace Method: Maximum Likelihood of True Positions

Recalling the objective function in the Laplace method

$-\text{ln}(f(\textbf{X}))=\sum^{N}_{i=1}\left[\frac{1}{2}\text{ln}(2\pi\varepsilon_{i})+\frac{(o_{i}-x_{i})^{2}}{2\varepsilon_{i}}\right]+\sum^{N-1}_{i=1}\left[\frac{1}{2}\text{ln}(2\pi\omega_{i})+\frac{(x_{i+1}-x_{i})^{2}}{2\omega_{i}}\right],$
(52)

the gradient of Eq. 52 is

$\displaystyle-\frac{\partial\text{ln}f}{\partial x_{1}}$
$\displaystyle=\frac{(x_{1}-o_{1})}{\varepsilon_{1}}+\frac{(x_{1}-x_{2})}{\omega_{1}}$
$\displaystyle-\frac{\partial\text{ln}f}{\partial x_{i}}$
$\displaystyle=\frac{(x_{i}-o_{i})}{\varepsilon_{i}}+\frac{(x_{i}-x_{i-1})}{\omega_{i-1}}+\frac{(x_{i}-x_{i+1})}{\omega_{i}}$
$\displaystyle-\frac{\partial\text{ln}f}{\partial x_{N}}$
$\displaystyle=\frac{(x_{N}-o_{N})}{\varepsilon_{N}}+\frac{(x_{N}-x_{N-1})}{\omega_{N-1}},$
(53)

where $i\in 2:N-1$. The Hessian $-\text{ln}\nabla\nabla f(\widehat{X})=M$, of
Eq. 52 has the non-zero elements

$\displaystyle M_{1,1}$ $\displaystyle=-\frac{\partial^{2}\text{ln}f}{\partial
x_{1}\partial x_{1}}=\frac{1}{\varepsilon_{1}}+\frac{1}{\omega_{1}}$
$\displaystyle M_{i,i}$ $\displaystyle=-\frac{\partial^{2}\text{ln}f}{\partial
x_{i}\partial
x_{i}}=\frac{1}{\varepsilon_{i}}+\frac{1}{\omega_{i-1}}+\frac{1}{\omega_{i}}$
$\displaystyle M_{N,N}$ $\displaystyle=-\frac{\partial^{2}\text{ln}f}{\partial
x_{N}\partial x_{N}}=\frac{1}{\varepsilon_{N}}+\frac{1}{\omega_{N-1}}$
$\displaystyle M_{i,i+1}$
$\displaystyle=-\frac{\partial^{2}\text{ln}f}{\partial x_{i}\partial
x_{i+1}}=-\frac{1}{\omega_{i}}$ $\displaystyle M_{i,i-1}$
$\displaystyle=-\frac{\partial^{2}\text{ln}f}{\partial x_{i}\partial
x_{i-1}}=-\frac{1}{\omega_{i-1}}$

Setting the gradient in Eq. 14.1 equal to 0 and moving the constants to the
left hand side of the equation gives

$\displaystyle\frac{o_{1}}{\varepsilon_{1}}$
$\displaystyle=\frac{\widehat{x}_{1}}{\varepsilon_{1}}+\frac{(\widehat{x}_{1}-\widehat{x}_{2})}{\omega_{1}}$
$\displaystyle\frac{o_{i}}{\varepsilon_{i}}$
$\displaystyle=\frac{\widehat{x}_{i}}{\varepsilon_{i}}+\frac{(\widehat{x}_{i}-\widehat{x}_{i-1})}{\omega_{i-1}}+\frac{(\widehat{x}_{i}-\widehat{x}_{i+1})}{\omega_{i}}$
$\displaystyle\frac{o_{N}}{\varepsilon_{N}}$
$\displaystyle=\frac{\widehat{x}_{N}}{\varepsilon_{N}}+\frac{(\widehat{x}_{N}-\widehat{x}_{N-1})}{\omega_{N-1}}.$

With additional factoring, the expression looks like

$\displaystyle\frac{o_{1}}{\varepsilon_{1}}$
$\displaystyle=\widehat{x}_{1}\cdot\left(\frac{1}{\varepsilon_{1}}+\frac{1}{\omega_{1}}\right)+\widehat{x}_{2}\cdot\left(\frac{-1}{\omega_{1}}\right)$
$\displaystyle\frac{o_{i}}{\varepsilon_{i}}$
$\displaystyle=\widehat{x}_{i}\cdot\left(\frac{1}{\varepsilon_{i}}+\frac{1}{\omega_{i-1}}+\frac{1}{\omega_{i}}\right)+\widehat{x}_{i-1}\cdot\left(\frac{-1}{\omega_{i-1}}\right)+\widehat{x}_{i+1}\cdot\left(\frac{-1}{\omega_{i}}\right)$
$\displaystyle\frac{o_{N}}{\varepsilon_{N}}$
$\displaystyle=\widehat{x}_{N}\cdot\left(\frac{1}{\varepsilon_{N}}+\frac{1}{\omega_{N-1}}\right)+\widehat{x}_{N-1}\cdot\left(\frac{-1}{\omega_{N-1}}\right).$

The factored expression on the right can be expressed in terms of a vector
product of the Hessian matrix and the maximum likelihood of the true
positions, $M\cdot\widehat{X}$. We invert the Hessian matrix to bring it to
the other side of the equation so that the resulting expression for the
maximum likelihood looks like

$\widehat{X}=M^{-1}\Theta,$

where the components of $\Theta$ are

$\theta_{i}=o_{i}/\varepsilon_{i}.$

###  14.2 Laplace Method: Direct Integration of the Probability Distribution
Components

Starting from the probability component formalism

$f(X)=\prod^{N}_{i=1}P\\!\left({o_{i}}\,|\,{x_{i},x_{i+1}}\right)P\\!\left({x_{i}}\,|\,{x_{i+1}}\right)=\prod^{N}_{i=1}\,\mathcal{N}(o_{i},(1-\alpha_{i})x_{i}+\alpha_{i}x_{i+1},q_{i})\,\mathcal{N}(x_{i},x_{i+1},\omega_{i}),$

where $q_{i}$ is the variance due to the observation and
$\alpha_{i}=\frac{t_{\epsilon}}{2\delta t_{i}}$. We solve for the Hessian of
our objective function $M=-\nabla\nabla\text{ln}f(X)$

$\displaystyle M_{1,1}$ $\displaystyle=-\frac{\partial^{2}\ln f}{\partial
x_{1}^{2}}=\frac{(1-\alpha_{1})^{2}}{q_{1}}+\frac{1}{\omega_{1}}$
$\displaystyle M_{i,i}$ $\displaystyle=-\frac{\partial^{2}\ln f}{\partial
x_{i}^{2}}=\frac{(1-\alpha_{i})^{2}}{q_{i}}+\frac{(\alpha_{i-1})^{2}}{\varepsilon_{i-1}}+\frac{1}{\omega_{i}}+\frac{1}{\omega_{i-1}},\quad
2\leq i\leq N$ $\displaystyle M_{N+1,N+1}$
$\displaystyle=-\frac{\partial^{2}\ln f}{\partial
x_{N}^{2}}=\frac{\alpha_{N}^{2}}{q_{N}}+\frac{1}{\omega_{N}}$ $\displaystyle
M_{i,i+1}=M_{i+1,i}$ $\displaystyle=-\frac{\partial^{2}\ln f}{\partial
x_{i}\partial
x_{i+1}}=\frac{(1-\alpha_{i})(\alpha_{i})}{q_{i}}-\frac{1}{\omega_{i}},\quad
1\leq i\leq N.$

The form is similar to our functions, so we can solve for the maximum
likelihood in the same fashion

$\widehat{X}=M^{-1}\Theta$

Where the components of $\Theta$ are modified as

$\theta_{i}=\frac{(1-\alpha_{i})o_{i}}{q_{i}}+\frac{\alpha_{i-1}o_{i-1}}{q_{i-1}},$

and for completeness we set $\alpha_{0}=0$ and $\alpha_{N+1}=0$.

###  14.3 Markov Method: Marginalizing the Likelihood Function

To understand how the integration of $S^{\prime}$ on
$P\\!\left({O}\,|\,{D}\right)$ behaves, lets first consider the integration of
$O$ on $P\\!\left({O}\,|\,{D}\right)$ with one $o_{i}$ held constant, which is
equivalent to multiplying $P\\!\left({O}\,|\,{D}\right)$ with a delta
distribution $\delta(o_{i}^{\prime}-o_{i})$. The integral of
$P\\!\left({O}\,|\,{D}\right)$ and the delta distribution with respect to $O$
is of the form

$\displaystyle\int\mathrm{d}{O}\,\;\delta(o_{i}^{\prime}-o_{i})P(O|D)=\int\mathrm{d}{O}\,\mathrm{d}{X}\,\;\delta(o_{i}^{\prime}-o_{i})\prod^{N}_{i=1}\mathcal{M}_{i}\prod^{N-1}_{j=1}\mathcal{T}_{j}$
$\displaystyle=\int\mathrm{d}{X}\,\mathrm{d}{O}\,\;\delta(o_{i}^{\prime}-o_{i})\prod^{N}_{i=1}\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i})\prod^{N-1}_{j=1}\,\mathcal{N}(x_{j+1},x_{j},\omega_{j}).$

Shuffling the order of integration, so that the $O$ basis is integrated first
allows us to marginalize all $O$ except for $o_{i}$. There are then N terms of
$X$ which can be effectively marginalized

$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{O}\,\;\delta(o_{i}^{\prime}-o_{i})\prod^{N}_{i=1}\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i})\prod^{N-1}_{j=1}\,\mathcal{N}(x_{j+1},x_{j},\omega_{j})$
$\displaystyle=\int\mathrm{d}{X}\,\,\mathcal{N}(o_{i}^{\prime},x_{i},\varepsilon_{i})\prod^{N-1}_{j=1}\,\mathcal{N}(x_{j+1},x_{j},\omega_{j})=1.$

Now we wish to perform a similar integration, but this time on the basis of
$S^{\prime}$. The integration of $P\\!\left({O}\,|\,{D}\right)$ with respect
to $O$ is a completely different function than an integration with respect to
$S^{\prime}$, but we wish to show that the integration on $S^{\prime}$ yields
analogous to results to integration on $O$ with one $o_{i}$ held constant. If
$o_{i}=o_{i}^{\prime}$ is held constant, than we can directly express
$s^{\prime}_{i}=o_{i+1}-o_{i}^{\prime}$ in terms of one variable, $o_{i+1}$,
if $i<N$. We can also express $s^{\prime}_{i-1}$ in terms of $o_{i-1}$ if
$i>1$. Analogously, every $o_{i+k}$ can be expressed as

$o_{i+k}=\sum^{i+k-1}_{j=i}s^{\prime}_{j}+o_{i}^{\prime}=s^{\prime}_{i+k-1}+o_{i}^{\prime}+g(s^{\prime}_{i+k-2},i)$

and every $o_{i-l}$ can be expressed as

$o_{i-l}=\sum^{i-1}_{j=i-l}-s^{\prime}_{j}+o_{i}^{\prime}=-s^{\prime}_{i-l}+o_{i}^{\prime}+h(s^{\prime}_{i-l+1},i).$

We perform this substitution with an arbitrary $o_{i}$ held fixed to express
the integral over $S^{\prime}$ as

$\displaystyle\int\mathrm{d}{S^{\prime}}\,\;P(O|D)=\int\mathrm{d}{S^{\prime}}\,\mathrm{d}{X}\,\;\prod^{N}_{i=1}\mathcal{M}_{i}\prod^{N-1}_{j=1}\mathcal{T}_{j}$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{S^{\prime}}\,\;N(o_{i}^{\prime},x_{i},\varepsilon_{i})N(s_{i}^{\prime}+o_{i}^{\prime},x_{i+1},\varepsilon_{i+1})\prod^{N}_{j=i+2}\,\mathcal{N}(s_{j-1}^{\prime}+o_{i}^{\prime}+g(s_{j-2},i),x_{i},\varepsilon_{i})$
$\displaystyle\prod^{i-1}_{k=1}\,\mathcal{N}(-s^{\prime}_{1}+o_{i}^{\prime}+h(s^{\prime}_{2},i),x_{1},\varepsilon_{1})\prod^{N-1}_{l=1}\,\mathcal{N}(x_{l+1},x_{l},\omega_{l}).$

We can then shuffle the order of integration so that we can iteratively
integrate the components of the $S$ basis, essentially the components furthest
from $o_{i}^{\prime}$, that are expressed in only one of the univariate
gaussian functions that comprise $P\\!\left({O}\,|\,{D}\right)$, effectively
performing a sequential marginalization

$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{S^{\prime}}\,\;N(o_{i}^{\prime},x_{i},\varepsilon_{i})N(s_{i}^{\prime}+o_{i}^{\prime},x_{i+1},\varepsilon_{i+1})\prod^{N}_{j=i+2}\,\mathcal{N}(s_{j-1}^{\prime}+o_{i}^{\prime}+g(s_{j-2},i),x_{i},\varepsilon_{i})$
$\displaystyle\prod^{i-1}_{k=1}\,\mathcal{N}(-s^{\prime}_{1}+o_{i}^{\prime}+h(s^{\prime}_{2},i),x_{1},\varepsilon_{1})\prod^{N-1}_{l=1}\,\mathcal{N}(x_{l+1},x_{l},\omega_{l})$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{s_{i}}\,\mathrm{d}{s_{i-1}}\,\,\mathcal{N}(o_{i}^{\prime},x_{i},\varepsilon_{i})N(s_{i}^{\prime}+o_{i}^{\prime},x_{i+1},\varepsilon_{i+1})\,\mathcal{N}(-s_{i-1}^{\prime}+o_{i}^{\prime},x_{i},\varepsilon_{i})\prod^{N-1}_{l=1}\,\mathcal{N}(x_{l+1},x_{l},\omega_{l})$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\,\mathcal{N}(o_{i}^{\prime},x_{i},\varepsilon_{i})\prod^{N-1}_{l=1}\,\mathcal{N}(x_{l+1},x_{l},\omega_{l})=1.$
(54)

From this result, we see that holding a single $o_{i}$ fixed is quite
arbitrary, as the term is eventually marginalized by its associated $x_{i}$.
It is also apparent that fixing a particular $o_{i}$ allows complete isolation
of a particular $s_{i}$ basis if all other $s_{k}$ bases are marginalized.
However, if we wish to evaluate $s_{i}$ and $s_{i-1}$ components with $o_{i}$
fixed, we see in the integral expression that there will be some correlation
between adjacent displacements. Most importantly, we see that
$P\\!\left({O}\,|\,{D}\right)$ is a normalized probability density under
$S^{\prime}$.

###  14.4 Markov Method: Expectation Calculations

We shall solve for the expectation values on $S^{\prime}$ over
$P\\!\left({O}\,|\,{D}\right)$. In order to do so, we will use the trick in
Eq.14.3 where we hold one observation constant until the appropriate
substitutions are taken. Solving for $\langle s_{i}\rangle$ in this spirit
yields

$\displaystyle\langle s_{i}\rangle=$
$\displaystyle\int\mathrm{d}{S^{\prime}}\,\mathrm{d}{X}\,s_{i}\prod^{N}_{j=1}\mathcal{M}_{j}\prod^{N-1}_{k=1}\mathcal{T}_{k}$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{s^{\prime}_{i}}\,\;s^{\prime}_{i}N(s_{i}^{\prime}+o_{i},x_{i+1},\varepsilon_{i+1})\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i})\prod^{N-1}_{j=1}\,\mathcal{N}(x_{j+1},x_{j},\omega_{j})$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{s^{\prime}_{i}}\,\;s^{\prime}_{i}N(s_{i}^{\prime},0,\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i})\,\mathcal{N}(x_{i},\beta_{i},\gamma_{i})\,\mathcal{N}(x_{i+1},\beta_{i+1},\gamma_{i+1})\prod^{N-1}_{j=1,j\neq
i}\,\mathcal{N}(x_{j+1},x_{j},\omega_{j})$ $\displaystyle=$
$\displaystyle\int\mathrm{d}{s^{\prime}_{i}}\,\;s^{\prime}_{i}N(s_{i}^{\prime},0,\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i})=0,$
(55)

where the terms $\beta_{i}$ and $\gamma_{i}$ are intermediate variables that
are generated from rearranging Gaussian functions. These variables are
immediately marginalized by integration on X, so we shall omit their explicit
representation. Aside from the fact that the expectation value on the
displacement means are all 0, we also discover another striking property from
Eq. 14.4, that is, the marginalization of all other $s^{\prime}$ and all $X$
isolates the expectation calculation to components of $s^{\prime}_{i}$ that
are independent of all other marginalized variables. The same results from Eq.
14.4 are used so that the expectation of the variance is

$\displaystyle\langle s_{i}^{2}\rangle=$
$\displaystyle\int\mathrm{d}{S^{\prime}}\,\mathrm{d}{X}\,\;s^{{}^{\prime}2}_{i}\prod^{N}_{j=1}\mathcal{M}_{j}\prod^{N-1}_{k=1}\mathcal{T}_{k}$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{s^{\prime}_{i}}\,\;s^{{}^{\prime}2}_{i}N(s_{i}^{\prime}+o_{i},x_{i+1},\varepsilon_{i+1})\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i})\prod^{N-1}_{j=1}\,\mathcal{N}(x_{j+1},x_{j},\omega_{j})$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{s^{\prime}_{i}}\,\;s^{{}^{\prime}2}_{i}N(s_{i}^{\prime},0,\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i})\,\mathcal{N}(x_{i},\beta_{i},\gamma_{i})\,\mathcal{N}(x_{i+1},\beta_{i+1},\gamma_{i+1})\prod^{N-1}_{j=1,j\neq
i}\,\mathcal{N}(x_{j=1},x_{j},\omega_{j})$ $\displaystyle=$
$\displaystyle\int\mathrm{d}{s^{\prime}_{i}}\,\;s^{{}^{\prime}2}_{i}N(s_{i}^{\prime},0,\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i})=\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i}.$

Now on to the covariance terms; starting with adjacent displacements $\langle
s_{i}s_{i+1}\rangle$. However, if a given $o_{i}$ is fixed, there is a
dependence between adjacent displacements. With the techniques implemented in
Eq. 14.4 a bivariate gaussian function is isolated from the rest of the
integral expression.

$\displaystyle\langle s_{i}^{\prime}s_{i+1}^{\prime}\rangle=$
$\displaystyle\int\mathrm{d}{S^{\prime}}\,\mathrm{d}{X}\,\;s_{i}^{\prime}s_{i+1}^{\prime}\prod^{N}_{j=1}\mathcal{M}_{j}\prod^{N-1}_{k=1}\mathcal{T}_{k}$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{X}\,\mathrm{d}{s_{i}^{\prime}}\,\mathrm{d}{s_{i-1}^{\prime}}\,\;s_{i}^{\prime}s_{i+1}^{\prime}N(s_{i}^{\prime}+o_{i},x_{i+1},\varepsilon_{i+1})\,\mathcal{N}(o_{i},x_{i},\varepsilon_{i})\,\mathcal{N}(s_{i-1}^{\prime}+o_{i},x_{i},\varepsilon_{i})\prod^{N-1}_{j=1}\,\mathcal{N}(x_{j+1},x_{j},\omega_{j})$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{s^{\prime}_{i}}\,\mathrm{d}{s^{\prime}_{i+1}}\,\;s_{i}^{\prime}s_{i+1}^{\prime}\,\mathcal{N}(S_{b},0,\Sigma_{b})=-\varepsilon_{i+1}.$

Where the substitution of variables required for integrating the expectation
of $\left\langle{s_{i}^{\prime}s_{i+1}^{\prime}}\right\rangle$ induces a
bivariate Gaussian function with location parameters
$S_{b}=[s_{i},s_{i+1}]^{\intercal}$ and covariance matrix

$\Sigma_{b}=\begin{bmatrix}\omega_{i}+\epsilon_{i}+\epsilon_{i+1}&-\epsilon_{i+1}\\\
-\epsilon_{i+1}&\omega_{i+1}+\epsilon_{i+1}+\epsilon_{i+2}\end{bmatrix}.$

Conversely, when evaluating $\langle s_{i}s_{i+k}\rangle$ where $k>1$, it is
necessary to temporarily fix $o_{i}$ as well as $o_{i+k}$, this results in
isolating two univariate Gaussian functions which results in

$\displaystyle\langle s_{i}^{\prime}s_{i+k}^{\prime}\rangle=$
$\displaystyle\int\mathrm{d}{S^{\prime}}\,\mathrm{d}{X}\,\;s_{i}^{\prime}s_{i+k}^{\prime}\prod^{N}_{j=1}\mathcal{M}_{j}\prod^{N-1}_{k=1}\mathcal{T}_{k}$
$\displaystyle=$
$\displaystyle\int\mathrm{d}{s^{\prime}_{i}}\,\mathrm{d}{s^{\prime}_{i+k}}\,\;s_{i}^{\prime}s_{i+k}^{\prime}N(s_{i}^{\prime},0,\varepsilon_{i}+\varepsilon_{i+1}+\omega_{i})N(s_{i+k}^{\prime},0,\varepsilon_{i+k}+\varepsilon_{i+k+1}+\omega_{i+k})=0.$

##  15 Implementation

###  15.1 Log–product Computation

A common algorithmic problem shared by each of the three methods is the need
to compute logarithms of products. A straight forward method is to utilize the
identity

$\ln\left(\prod_{i=1}^{N}a_{i}\right)=\sum_{i=1}^{N}\ln a_{i},$ (56)

the left side of which requires one logarithm and $N-1$ multiplications, while
the right side requires $N$ logarithms and $N-1$ additions. Directly using the
left side of Eq. 56 in computations can lead to numerical over- or underflow,
while the $N$ logarithms required for the right side can dominate the
computational costs, taking up the majority of the computational cycles for
each of the three algorithms. Thus, to minimize the number of logarithm
computations, yet still maintain numerical accuracy, we utilize a hybrid log-
product implementation to evaluate forms like Eq. 56. The log–product method
builds up a product of $a_{i}$-values, only taking a logarithm when
multiplying by the next $a_{i}$ would lead to loss of precision, overflow, or
underflow.

###  15.2 Computation of Variance terms

Given the inputs, each of the method firsts compute the $N$ variance terms due
to measurement, $\varepsilon_{i}(D)=v_{i}-Dt_{\epsilon}/3$, and the $N-1$
variance terms due to diffusion, $\omega_{i}(D)=2D\delta t_{i}$. Any of these
variance terms can be arbitrarily close to zero, so to prevent numerical
instabilities we bound these terms away from zero by at least machine epsilon,
while preserving the sign which can be negative for $\varepsilon_{i}(D)$.

###  15.3 The Recursive Method Algorithm

Listing 1: The Recursive Algorithm for the log-likelihood calculation in C++.
This fundemental algorithm can be impemented with just the log function from
the C++ standard math library. We rely on the variances to be computed as
described in Sec. 15.2, and the log–product computation as described in Sec.
15.1

⬇

FloatT recusiveLLH(int N,

const FloatT Obs[],

const FloatT dT[],

const FloatT vD[],

const FloatT vM[])

{

FloatT alpha[N-1];

FloatT eta = vD[0]+vM[0];

FloatT mu = Obs[0];

FloatT LLH = 0;

for(int n=1;n<N-1;n++){

FloatT temp_alpha = vM[n]+eta;

alpha[n-1] = temp_alpha;

FloatT temp_diff = (Obs[n]-mu);

LLH += temp_diff*temp_diff/temp_alpha;

mu = (mu*vM[n]+Obs[n]*eta)/temp_alpha;

eta = vM[n]*eta/temp_alpha+vD[n];

}

alpha[N-2] = vM[N-1]+eta;

LLH += (N-1)*log2pi;

LLH += logprod(alpha);

FloatT temp_diff = (Obs[N-1]-mu);

LLH += temp_diff*temp_diff/alpha[N-2];

LLH *= -0.5;

return LLH;

}

###  15.4 Tri-Diagonal matrix algorithms

The Laplace and Markov method both require solving linear systems of the form
$Ax=b$, where $A$ is symmetric tri-diagonal (all non-zero terms are on the
main diagonal or those diagonals immediately above and below the main
diagonal). A naive solution based on inverting matrix $A$ has computational
complexity $\mathcal{O}\left({N^{3}}\right)$, where $N$ is the length of the
trajectory. However because $A$ is tri-diagonal there are established
algorithms (30) that use the Gaussian elimination strategy to solve the system
in time $\mathcal{O}\left({N}\right)$. Furthermore, the determinant of the
matrix can be shown to follow a recurrence relation (31) that leads to a
linear time computation. Combined with the log–product computation described
in Sec. 15.1, this leads to a fast algorithm for computing $\log(\det(A))$
with a minimum number of calls to the logarithm function, and time complexity
$\mathcal{O}\left({N}\right)$. We make use of these algorithms for the Laplace
and Markov method implementations.

###  15.5 DST Algorithm

The DST algorithm code is based on the Matlab code provided by the authors
(16). Our implementation is as faithful as possible to the original
implementation but there are a few caveats that should be mentioned. First
because we assume that the localization variances are known, we use the form
of the estimator that takes in the mean localization variance as an input
parameter. Because the DST can only incorporate a single localization variance
we use the mean of the given localization variances. For the $R$ parameter we
use $t_{\epsilon}/6$ as was derived for uniform exposure intervals (16). The
original implementation also leaves off the constant term
$-(1/2)(N-1)\log(2\pi)$ from the likelihood calculation. While for MLE
estimation this is not important, it does become important for other
downstream analysis using the likelihood, and our algorithms do include this
term, so we have added that in to the DST to make the plotted comparisons more
fair. Also, we ensure that the resulting values are always real as the
original code can return complex floating point numbers, but only the real
parts have meaning. Finally, because of the reliance on the discrete sine
transform function which is provided by Matlab, but not available natively in
the C++ standard math library, the speed results from the main paper only test
the provided Matlab implementation.