# Kernel Density Estimation with Linked Boundary Conditions

Matthew J. Colbrook Department of Applied Mathematics and Mathematical
Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK.
Email<EMAIL_ADDRESS>Zdravko I. Botev School of Mathematics and
Statistics, The University of New South Wales, Sydney, NSW 2052, Australia.
Karsten Kuritz Institute for Systems Theory and Automatic Control, University
of Stuttgart, 70569 Stuttgart, Germany.  Shev MacNamara ARC Centre of
Excellence for Mathematical and Statistical Frontiers, School of Mathematical
and Physical Sciences, University of Technology Sydney, NSW 2007, Australia.

###### Abstract

Kernel density estimation on a finite interval poses an outstanding challenge
because of the well-recognized bias at the boundaries of the interval.
Motivated by an application in cancer research, we consider a boundary
constraint linking the values of the unknown target density function at the
boundaries. We provide a kernel density estimator (KDE) that successfully
incorporates this linked boundary condition, leading to a non-self-adjoint
diffusion process and expansions in non-separable generalized eigenfunctions.
The solution is rigorously analyzed through an integral representation given
by the unified transform (or Fokas method). The new KDE possesses many
desirable properties, such as consistency, asymptotically negligible bias at
the boundaries, and an increased rate of approximation, as measured by the
AMISE. We apply our method to the motivating example in biology and provide
numerical experiments with synthetic data, including comparisons with state-
of-the-art KDEs (which currently cannot handle linked boundary constraints).
Results suggest that the new method is fast and accurate. Furthermore, we
demonstrate how to build statistical estimators of the boundary conditions
satisfied by the target function without apriori knowledge. Our analysis can
also be extended to more general boundary conditions that may be encountered
in applications.

Keywords: density estimation, diffusion, unified transform, linked boundary
conditions, boundary bias, biological cell cycle.

##  1 Introduction and Background

Suppose we are given an independent and identically distributed sample
$X_{1},\ldots,X_{n}$ from some unknown density function $f_{X}$. Throughout,
we will use a subscript $X$ in $f_{X}$ to indicate that $f_{X}$ is the
probability density function of the random variable $X$. We will also denote
expectation and variance with respect to $f_{X}$ by $\mathbb{E}_{f_{X}}$ and
$\mathrm{Var}_{f_{X}}$ respectively. Estimating the density $f_{X}$ is one of
the most common problems for discovering patterns in statistical data [7, 62,
63]. When the support of $f_{X}$ is the whole real line, a simple and popular
non-parametric method for estimating $f_{X}$ is the kernel density estimator
(KDE)

$\widehat{f}(x;t)=\frac{1}{n\sqrt{t}}\sum_{k=1}^{n}\varphi\left(\frac{x-X_{k}}{\sqrt{t}}\right),$
(1)

with a kernel $\varphi(x)$. A common choice is a Gaussian kernel
$\varphi(x)=\exp(-x^{2}/2)/\sqrt{2\pi}$. Here $\sqrt{t}$ is the so-called
bandwidth parameter that controls the smoothness of the estimator (see, for
example, [68, 67, 62, 63] and references therein). Another viewpoint is to
connect kernel density estimation to a diffusion equation, an approach
pioneered by the second author in [5]. Our goal in this article is to extend
this analysis to linked boundary conditions. A key tool in our analysis is the
unified transform (also known as the Fokas method), a novel transform for
analyzing boundary value problems for linear (and integrable non-linear)
partial differential equations [25, 26, 28, 27, 66, 17, 16, 12, 11, 10, 13, 9,
61]. An excellent pedagogical review of this method can be found in the paper
of Deconinck, Trogdon & Vasan [18].

It is well-known that $\widehat{f}(x;t)$ is not an appropriate kernel
estimator when $f_{X}$ has compact support [29], which (without loss of
generality) we assume to be the unit interval $[0,1]$. The main reason for
this is that $\widehat{f}(x;t)$ exhibits significant boundary bias at the end-
points of the interval. For example, with a Gaussian kernel, no matter how
small the bandwidth parameter, $\widehat{f}(x;t)$ will have non-zero
probability mass outside the interval $[0,1]$. Various solutions have been
offered to cope with this boundary bias issue, which may be classified into
three main types:

  1. (a)

Using special (non-Gaussian) kernels with support on $[0,1]$ or on
$[0,\infty)$, as in [6, 42, 56];

  2. (b)

Adding bias-correction terms to $\widehat{f}(x;t)$ as in [14, 37];

  3. (c)

Employing domain transformations [29, 44], which work by mapping the data to
$(-\infty,\infty)$, constructing a KDE on the whole real line, and finally
mapping the estimate back to $[0,1]$.

Additionally, sometimes we not only know that $f_{X}$ has support on $[0,1]$,
but also have extra information about the values of $f_{X}$ at the boundaries.
One example of this situation is what we will refer to as a _linked boundary
condition_ , where we know apriori that

$f_{X}(0)=rf_{X}(1)$

for some known given parameter $r\geq 0$. Most of our analysis also carries
over to complex $r$, as long as $r\neq-1$ (the PDE (2) is degenerate irregular
and the problem ill-posed when $r=-1$), but we focus on $r\geq 0$ since in
statistics $f_{X}\geq 0$. An example that motivated the current article arises
in the field of biology [39, 38], in particular cell cycle studies in cancer
research. The cell cycle itself is one of the fundamentals of biology and
knowledge about its regulation is crucial in the treatment of various
diseases, most prominently cancer. Cancer is characterized by an uncontrolled
cell growth and commonly treated with cytotoxic drugs. These drugs interfere
with the cell cycle and in this way cause cancer cells to die. By studying the
effect of chemicals on the cell cycle one can discover new drugs, identify
potential resistance mechanisms or evaluate combinatorial therapy. These kind
of studies have benefited from continued improvement in cell population
analysis methods like fluorescence microscopy, flow cytometry, CyTOF or
single-cell omics, where the abundance of up to thousands of cellular
components for every individual cell in a population is measured. In such an
experiment, cells in an unsynchronized cell population are spread over all
stages of the cell cycle. Trajectory inference algorithms then reduce the
dimensionality to a pseudotime scale by ordering cells in the population based
on their similarity in the dataset [53]. Subsequently, mathematical methods
based on ergodic principles infer molecular kinetics in the cell cycle from
the distribution of cells in pseudotime. The value at the left boundary of
this distribution must, because of cell division, be double the value at the
right boundary. In other words, we have linked boundary conditions with the
constant $r=2$, but otherwise, we do not know the value of the density at the
boundaries of the domain. The problem is described in more detail in Section
5.2, where we also demonstrate the estimator with linked boundary condition on
a real dataset. In particular, for this example, respecting the linked
boundary condition is crucial for generating the correct kinetics due to a
certain mapping between pseudotime and real time. See also [39, 38], for
example, for further motivation and discussion. In other applications, even if
we do not know the value of $r$, one can approximate the true value of $r$
which, together with the methods proposed in this article, leads to an
increase in the rate of approximation of $f_{X}$ as the sample size $n$
becomes large (we do this for an example in Section 5.1, see also §3 for some
results in this direction).

Unfortunately, to the best of our knowledge, all of the currently existing
kernel density estimation methods, bias-correcting or not, cannot
satisfactorily handle the linked boundary condition. Figure 1 shows a typical
example of what can go wrong when a standard density estimator is applied to
real biological data. The result is a smooth density with two unacceptable
features:

  * •

The domain $x\in[0,1]$ is not respected, and instead the solution has positive
density for negative values of $x$, and also for $x>1$, which are physically
unreasonable. This problem can be addressed using existing bias-correction
methods and is not the challenge that we had to overcome in this article.

  * •

The density does not respect the important biological constraint of the linked
boundary condition (that the left value should be double the right, in this
particular application), and instead the density decays to zero as $|x|$
becomes large. Existing bias-correction methods do not address this problem.

Figure 1: A typical example of output from a KDE (ksdensity from MATLAB)
applied to our real biological data. This does not respect the domain, and it
also does not respect the important linked boundary conditions. The methods
that we propose in this article address those issues simultaneously, with
results for this data set shown in Figure 7.

The purpose of this article is to describe a new KDE that can handle the more
general problem of linked boundary conditions with an arbitrary value of $r$;
the situation of interest in the biological application where $r=2$ is then
solved as an important special case. Figure 7 (C) shows a successful
application of our proposed method. The MAPiT toolbox for single-cell data
analysis [38] applies our new KDE with linked boundary conditions to analyze
cell cycle dependent molecular kinetics.

Our proposed estimator is of type (a), that is, we construct a special kernel
with support on $[0,1]$, and such that the linked boundary condition is
incorporated into the resulting estimator. Our kernel is inspired by the
solution of a diffusion-type PDE [1, 5, 45, 55, 69]. In particular, we modify
the diffusion model in [5] so that it satisfies the linked boundary
conditions. Unlike the case in [5], the non-self-adjoint initial-boundary
problem that arises cannot be diagonalized, meaning the solution cannot be
expressed as a series solution of eigenfunctions of the spatial differential
operator in the usual sense. Instead, we use the unified transform, which
provides an algorithmic recipe for solving these types of problems via an
integral solution. This was the way we first found the solution formula to our
diffusion model, and the integral representation simplifies many of the proofs
in our analysis. So far, the only case of our problem considered in the
literature on this method has been $r=1$ [66] (periodic). For the heat
equation with oblique Robin boundary conditions/non-local boundary conditions
we refer the reader to [43, 47, 51] and for interface problems we refer the
reader to [59, 60, 58]. Recently linked boundary conditions have been
considered for the Schrödinger equation in [50] (however, in [50], the
parameters were chosen such that the characteristic values were simple, in
other words the eigenvalues were simple, making the analysis easier and
leading to a series solution in terms of bona fide eigenfunctions).

We then construct a series expansion in non-separable generalized
eigenfunctions of the spatial derivative operator by deforming the contours in
the integral representation and applying Cauchy’s residue theorem. This formal
solution is then rigorously verified and studied via a non-symmetric heat
kernel. Each of these representations (integral and series) is beneficial for
different analysis. For instance, the integral representation is much easier
to construct and makes it easier to study regularity properties, as well as
some parts of the behavior as $t\downarrow 0$, whereas the kernel
representation is useful for proving conservation of mass (the solution
generates a true probability measure) and studying the asymptotic mean
integrated squared error (AMISE). Although it is not the goal of the present
article, we envisage that the method that we demonstrate here can also be
generalized to the multivariate case and to scenarios where other types of
boundary conditions (such as linked derivatives or on-local boundary
conditions) arise or can be estimated. In these situations, we recommend using
the unified transform to find the solution of the resulting PDE. For numerical
implementation of the unified transform, we refer the reader to [15].

We also consider the discrete counterpart of the continuous model for two
reasons. First, it is a numerical approximation to the continuous model and a
useful way to compute the solution of the PDE. Second, the discrete model is
relevant when we deal with data which is already pre-binned.

The rest of the article is organized as follows. In the next section, we
describe the continuous model for the application at hand. Our results provide
the necessary assurances that the PDE model is a valid and accurate density
estimator. We then discuss the issue of choosing an optimal bandwidth
(stopping time for the PDE model), including pointwise bias, asymptotic
properties and the AMISE. We briefly discuss numerical methods for calculating
the estimator and, in particular, a discretized version of the continuous PDE,
which we prove converges to the unique continuous solution. Finally, the new
method is applied to a real dataset from a biological application in Section
5.2, and we also provide an illustrative set of examples with synthetic
datasets. We compare our new estimator to several well-known methods and these
results suggest that our new method is typically more accurate and faster, and
that it does not suffer from boundary bias. We finish with a short conclusion.

All technical analysis and proofs are moved to the appendices to ensure that
the presentation flows more easily. Freely available code for the new kernel
methods is also provided at https://github.com/MColbrook/Kernel-Density-
Estimation-with-Linked-BCs.

##  2 The Continuous Linked–Boundaries Model

In this section, we present the continuous diffusion model that satisfies the
linked boundary condition and discuss the analytical properties of its
solution. Our proposed diffusion model for a linked-boundary KDE is the
solution of the formal PDE system:

$\begin{split}\frac{\partial f}{\partial
t}&=\frac{1}{2}\frac{\partial^{2}f}{\partial x^{2}},\qquad
x\in[0,1],\;\;\;t>0,\\\ \mathrm{IC:}\quad\lim_{t\downarrow
0}f(\cdot,t)&=f_{0},\\\ \mathrm{BCs:}\quad f(0,t)&=rf(1,t),\quad\frac{\partial
f}{\partial x}(0,t)=\frac{\partial f}{\partial x}(1,t),\quad\forall
t>0.\end{split}$ (2)

The boundary condition $\frac{\partial f}{\partial x}(0,t)=\frac{\partial
f}{\partial x}(1,t)$ is enforced so that the solution at any time $t\geq 0$
gives a probability measure (see Theorem 4). When considering the setup
described in the introduction, the initial condition is given by

$\textstyle f_{0}=\frac{1}{n}\sum_{k=1}^{n}\delta_{X_{k}},$ (3)

the empirical measure of the given sample $X_{1},\ldots,X_{n}$. In other
words, $f_{0}$ is a sum of Dirac delta distributions. However, in our analysis
we also consider more general initial conditions. Many of the existence and
uniqueness theorems carry over from the well-known $r=1$ (periodic) case. In
particular, the boundary conditions and PDE make sense when the initial data
is given by a finite Borel measure, which we also denote by $f_{0}$. Sometimes
we will also refer to a function $g$ as a measure through the formula
$g(U)=\int_{U}g(x)dx$ for Borel sets $U$. Therefore, since the initial data is
a distribution, we need to be precise by what we mean when writing
$\lim_{t\downarrow 0}f(\cdot,t)=f_{0}$.

######  Definition 1.

Denote the class of finite Borel measures on $[0,1]$ by $M([0,1])$ and equip
this space with the vague topology (i.e. weak∗ topology). We let
$C^{w}(0,T;M([0,1]))$ denote the space of all continuous maps

$\displaystyle\mu:[0,T)\rightarrow M([0,1]),$ $\displaystyle\mu(t)=\mu_{t}.$

In other words, $\mu_{t}$ is continuous as a function of $t$ in the vague
topology, meaning that for any given function $g$ that is continuous on the
interval $[0,1]$, the integral $\int_{0}^{1}g(x)d\mu_{t}(x)$ is continuous as
a function of $t\in[0,T)$.

We look for weak solutions of (2). In terms of notation, we will denote the
derivative with respect to $x$ by $g_{x}$ and use $\mu(g)$ to denote the
integration of a function $g$ against a measure $\mu$. The following adjoint
boundary conditions are exactly those that arise from formal integration by
parts.

######  Definition 2.

Let $\mathcal{F}(r)$ denote all $g\in C^{\infty}([0,1])$ satisfying the
adjoint linked boundary conditions

$g(1)=g(0),\quad g_{x}(1)=rg_{x}(0).$ (4)

######  Definition 3 (Weak Solution).

Let $f_{0}\in M([0,1])$ such that $f_{0}(\\{0\\})=rf_{0}(\\{1\\})$. We say
that $\mu\in C^{w}(0,T;M([0,1]))$ is a weak solution to (2) for $t\in[0,T)$ if
$\mu_{0}=f_{0}$ and for all $g\in\mathcal{F}(r)$, $\mu_{t}(g)$ is
differentiable for $t>0$ with

$\frac{d}{dt}\mu_{t}(g)=\frac{1}{2}\mu_{t}(g_{xx}).$ (5)

We can now precisely state the well-posedness of (2).

######  Theorem 1 (Well-Posedness).

Assume our initial condition $f_{0}$ lies in $M([0,1])$ and satisfies
$f_{0}(\\{0\\})=rf_{0}(\\{1\\})$. Then there exists a unique weak solution to
(2) for $t\in[0,T)$ for any $T\in(0,\infty]$, which we denote by $f(\cdot,t)$.
For $t>0$ this weak solution is a function that is smooth in $t$ and real
analytic as a function of $x$. Furthermore, the solution has the following
properties which generalize the classical periodic case of $r=1$:

  1. 1.

If $f_{0}\in C([0,1])$ (the space of continuous functions on $[0,1]$), then
for any $x\in(0,1)$, $f(x,t)$ converges to $f_{0}(x)$ as $t\downarrow 0$. If
$f_{0}(0)=rf_{0}(1)$ then $f(\cdot,t)$ converges to $f_{0}$ as $t\downarrow 0$
uniformly over the whole closed interval $[0,1]$.

  2. 2.

If $1\leq p<\infty$ and $f_{0}\in L^{p}([0,1])$, then $f$ is the unique weak
solution in $C(0,T;L^{p}([0,1]))$ and $f(\cdot,t)$ converges to $f_{0}$ as
$t\downarrow 0$ in $L^{p}([0,1])$.

###### Proof.

See Appendix A.2. ∎

The system (2) is a natural candidate for density estimation with such a
linked boundary condition. Whilst Theorem 1 is expected and analogous to the
$r=1$ case, due to the non-self-adjoint boundary conditions, it is not
immediately obvious what properties solutions of (2) have. For example, one
question is whether or not the solution is a probability density for $t>0$,
and what its asymptotic properties are. Moreover, we would like to be able to
write down an explicit solution formula (and ultimately use this to
numerically compute the solution), a formal derivation of which is given in
Appendix A.1 using the unified transform.

###  2.1 Solution formula and consistency of KDE at boundaries

If we ignore the constant $r$ in the boundary conditions of (2) (and replace
it by the special case $r=1$), then we would have the simple diffusion
equation with periodic boundary conditions. One can successfully apply Fourier
methods, separation-of-variables or Sturm–Liouville theory to solve the
periodic version of this PDE [24, 30]. However, when $r\neq 1$, making the
ansatz that a solution is of the ‘rank one’, separable form $f(x,t)=g(x)h(t)$
leads to a non-complete set of functions and separation of variables fails.
The differential operator associated with the evolution equation in (2) is
regular in the sense of Birkhoff [3] but not self-adjoint when $r\neq 1$, due
to the boundary conditions. Nevertheless, it is possible to generalize the
notion of eigenfunctions of the differential operator [8] and these
generalized eigenfunctions form a complete system in $L^{2}([0,1])$ [49, 40]
(and in fact form a Riesz basis). This allows us to obtain a series expansion
of the solution. The easiest way to derive this is through the unified
transform, which also generates a useful integral representation.

######  Theorem 2 (Integral and Series Representations of Diffusion
Estimator).

Suppose that the conditions of Theorem 1 hold. Then the the unique solution of
(2) has the following representations for $t>0$.

Integral representation:

$\begin{split}&2\pi
f(x,t)=\int_{-\infty}^{\infty}{\exp(ikx-k^{2}t/2)}\hat{f}_{0}(k)dk\\\
&-\textstyle\int_{\partial
D^{+}}\frac{\exp(ikx-k^{2}t/2)}{{\Upsilon(k)}}\left\\{\hat{f}_{0}(k)[(1+r)\exp(ik)-2r]+\hat{f}_{0}(-k)(1-r)\exp(-ik)\right\\}dk\\\
&-\textstyle\int_{\partial
D^{-}}\frac{\exp(ik(x-1)-k^{2}t/2)}{\Upsilon(k)}\left\\{\hat{f}_{0}(k)[2\exp(ik)-(1+r)]+\hat{f}_{0}(-k)(1-r)\right\\}dk.\end{split}$
(6)

Here the contours $\partial D^{\pm}$ are shown in Figure 8 and are
deformations of the boundaries of
$D^{\pm}=\\{k\in\mathbb{C}^{\pm}:\mathrm{Re}(k^{2})<0\\}$. The determinant
function is given by $\Upsilon(k)=2(1+r)(\cos(k)-1)$ and
$\hat{f}_{0}(k):=\int_{0}^{1}\exp(-ikx)f_{0}(x)dx.$

Series representation:

$\begin{split}f(x,t)=&\frac{2}{(1+r)}\hat{c}_{0}(0)\phi_{0}(x)\\\
&+\sum_{n\in\mathbb{N}}\frac{4\exp(-k_{n}^{2}t/2)}{(1+r)}\big{\\{}\hat{c}_{0}(k_{n})\phi_{n}(x)-k_{n}t(1-r)\hat{c}_{0}(k_{n})\sin(k_{n}x)\\\
&\quad\quad\quad\quad\quad+[\hat{s}_{0}(k_{n})-(1-r)\hat{s}_{1}(k_{n})]\sin(k_{n}x)\big{\\}},\end{split}$
(7)

where $k_{n}=2n\pi$ and

$\displaystyle\phi_{n}(x)=\left(r+(1-r)x\right)\cos(k_{n}x),$
$\displaystyle\hat{s}_{0}(k)=\int_{0}^{1}\sin(kx)f_{0}(x)dx,$
$\displaystyle\hat{c}_{0}(k)=\int_{0}^{1}\cos(kx)f_{0}(x)dx,$
$\displaystyle\hat{s}_{1}(k)=\int_{0}^{1}\sin(kx)xf_{0}(x)dx.$

###### Proof.

See Appendix A.2. ∎

In the case where $r\neq 1$, in addition to the usual separable solutions
$\exp(ik_{n}x-k_{n}^{2}t/2)$, the series expansion also includes the non-
separable solutions $\exp(ik_{n}x-k_{n}^{2}t/2)(x+ik_{n}t)$. We can understand
these as being generalized eigenfunctions in the following sense (see the
early papers [41, 65]). Define the operator

$\mathbb{A}=-\frac{d^{2}}{dx^{2}},\quad\mathcal{D}(\mathbb{A})=\\{u\in
H^{2}([0,1]):u(0)=ru(1),u_{x}(0)=u_{x}(1)\\},$ (8)

where $\mathcal{D}(\mathbb{A})$ denotes the domain of $\mathbb{A}$. We use
$\mathcal{N}$ to denote the null space, which is sometimes often termed the
kernel, of an operator, i.e. $\mathcal{N}(S)$ is the space of all vectors $v$
with $S(v)=0$. It is then easily checked that
$\phi_{n}\in\mathcal{N}((\mathbb{A}-k_{n}^{2}I)^{2})$. In particular, both
$\phi_{n}$ and $(\mathbb{A}-k_{n}^{2}I)\phi_{n}$ satisfy the required boundary
conditions. These functions block diagonalize the operator in an analogous
form to the Jordan normal form for finite matrices. If we consider any
generalized eigenspace $\mathcal{N}((\mathbb{A}-k_{n}^{2}I)^{2})$
corresponding to $k_{n}^{2}=4\pi^{2}n^{2}$ with $n>0$ and choose the basis
$\\{\sin(k_{n}x),\phi_{n}(x)/(2k_{n})\\}$, the operator acts on this subspace
as the matrix

$\left(\begin{tabular}[]{cc}$k_{n}^{2}$&$1-r$\\\
$0$&$k_{n}^{2}$\end{tabular}\right),$

which cannot be diagonalized for $r\neq 1$.

For our purposes of kernel density estimation, we define an integral kernel
$K$ so that we can write the solution as

$f(x,t)=\int_{0}^{1}K(r;x,y,t)f_{0}(y)dy.$

After some residue calculus (see (34) in the Appendix), this is given by the
somewhat complicated expression:

$\begin{split}K(r;x,y,t)&=\sum_{n\in\mathbb{Z}}{\exp(ik_{n}x-k_{n}^{2}t/2)}\Big{[}\exp(-ik_{n}y)+\frac{1-r}{1+r}(x+ik_{n}t)\exp(-ik_{n}y)\\\
&+\frac{1-r}{1+r}(x+ik_{n}t-1)\exp(ik_{n}y)+\frac{1-r}{1+r}y(\exp(ik_{n}y)-\exp(-ik_{n}y))\Big{]},\end{split}$
(9)

which can be re-expressed in terms of the more common $r=1$ kernel and its
derivative, as in (40). For the initial data (3) this gives the density
estimate

$f(x,t)=\frac{1}{n}\sum_{k=1}^{n}K(r;x,X_{k},t),$

a generalization of (1). A key consequence of the solution from Theorem 2 is
that the pointwise bias of the corresponding diffusion estimator vanishes if
$f_{X}$ is continuous with $f_{X}(0)=rf_{X}(1)$. Namely, we have the
following.

######  Theorem 3 (Consistency of Estimator at Boundaries).

Suppose that the initial data is given by (3) and that $f_{X}\in C([0,1])$
with $f_{X}(0)=rf_{X}(1)$. Then the solution of the PDE (2) satisfies

$\lim_{t\downarrow 0}\mathbb{E}_{f_{X}}(f(x,t))=f_{X}(x),$ (10)

uniformly in $x$. Further, if in addition $f_{X}\in C^{1}([0,1])$ and
$x_{t}=x+\mathcal{O}(\sqrt{t})$, then our estimator satisfies

$\left|\mathbb{E}_{f_{X}}(f(x_{t},t))-f_{X}(x)\right|\leq C(f_{X})\sqrt{t},$

with $C(f_{X})$ a constant independent of $x\in[0,1]$, but dependent on the
true $f_{X}$.

###### Proof.

See Appendix A.3. ∎

######  Remark 1.

For consistency in $L^{p}$ spaces, we refer the reader to Proposition 3 in
Appendix A.3.

###  2.2 Conservation of probability and non-negativity

In addition to establishing that the behavior of the PDE solution near the
boundaries is satisfactory, we also want the PDE solution to be a proper bona
fide probability density — a non-negative function integrating to unity. The
main tool in the proof of this is the Maximum Principle [24, 30] for parabolic
PDEs. The Maximum Principle states that a solution of the diffusion equation
attains a maximum on the ‘boundary’ of the two-dimensional region in space
$x\in[0,1]$ and time $t\geq 0$. If our initial condition is given by a
continuous function, then the maximum principle gives the following.

######  Proposition 1 (Bounds on Diffusion Estimator).

Suppose that the conditions of Theorem 1 hold and that $f_{0}$ is a continuous
function with $f_{0}(0)=rf_{0}(1)$ and non-negative with $0\leq a\leq
f_{0}(x)\leq b$ for all $x\in[0,1]$. Then for any $t>0$ and $x\in[0,1]$ we
have

$\textstyle\min\left\\{\frac{2r}{1+r},\frac{2}{1+r}\right\\}a\leq
f(x,t)\leq\max\left\\{\frac{2r}{1+r},\frac{2}{1+r}\right\\}b.$ (11)

In particular, $f$ remains bounded away from $0$ if $a>0$ and $r>0$.

###### Proof.

See Appendix A.4. ∎

However, we also want this to hold when $f_{0}$ is given by (3). Furthermore,
if we start with a probability measure as our initial condition, then we want
the solution to be the density function of a probability distribution for any
$t>0$. In the context of density estimation, this essential property
corresponds to conservation of probability. This is made precise in the
following theorem, which does not require continuous initial data.

######  Theorem 4 (A Bona Fide Kernel Density Estimator).

Suppose that the conditions of Theorem 1 hold and that the initial condition
$f_{0}$ is a probability measure. Then,

  1. 1.

$\int_{0}^{1}f(x,t)dx=1$, for $t>0$,

  2. 2.

$f(x,t)\geq 0$ for $t>0$ and $x\in[0,1]$.

###### Proof.

See Appendix A.5. ∎

From the solution formula (7), we can also characterize the behavior of the
solution for large bandwidths (large $t$), that is, when the estimator
oversmooths the data. An example of this is given in Figure 2.

######  Corollary 1 (Oversmoothing Behavior with Large Bandwidth).

Suppose that the conditions of Theorem 1 hold, then as $t\rightarrow\infty$,
$f$ converges uniformly on $[0,1]$ to the linear function

$\textstyle f_{\infty}(x):=\frac{2}{(1+r)}\hat{c}_{0}(0)\phi_{0}(x).$ (12)

This linear function is the unique stationary function that obeys the boundary
conditions and has the same integral over $[0,1]$ as $f_{0}$.

Figure 2: An example of the solution of the continuous PDE (2) at three time
points, with $f_{0}(x)=\frac{6}{11}(-2x^{2}+x+2)$. The values at the
boundaries change with time, but the ratio remains a constant with
$f(0,t)=2f(1,t)$. As $t\rightarrow\infty$, the solution converges to a
straight line.

##  3 Asymptotic Properties and Bandwidth Choice

An important issue in kernel density estimation is how to choose the bandwidth
parameter or, equivalently, the final or stopping time $T$ at which we compute
the solution of the PDE. This issue has already received extensive attention
in the literature [34, 57, 19, 35]. We now give a brief summary of that issue,
and we also make two suggestions for known methods already available. After
that, we address the issue specifically in the context of our linked
boundaries model.

At one extreme, if we choose $T=0$, then we recover the initial condition,
which is precisely the raw data, with an estimator with zero bias, but
infinite variance. At the other extreme, if we let $T\rightarrow\infty$, then
we obtain a stationary density that is a straight line (see Corollary 1),
which contains no information whatsoever about the raw data (other than the
empirical mean), giving an estimator of zero variance, but significant bias.
In between, $0<T<\infty$, we have some smoothing effect while also retaining
some information from the original data — an optimal balance between the
variance and the bias of the estimator.

One would also like a consistent estimator — as more and more data are
included, it must converge to the true density (for instance, in the mean
squared sense). Various proposals for the stopping times and their properties
are available. One of the most common choices is ‘Silverman’s rule of thumb’
[62], which works very well when the data is close to being normally
distributed. We expect that this choice is fine for the simpler datasets and
examples that we consider in this article. Another possible approach is to use
the output from the freely available software of one of the authors:
https://au.mathworks.com/matlabcentral/fileexchange/14034-kernel-density-
estimator. This is expected to be a better choice than Silverman’s rule in
situations where there are many widely separated peaks in the data. In
particular, [5] introduced a non-parametric selection method that avoids the
so-called _normal reference rules_ that may adversely affect plug-in
estimators of the bandwidth.

We now give a more precise treatment of the choice of smoothing bandwidth for
the linked boundaries model, as well as discussing the Mean Integrated Squared
Error (MISE) defined by

$\displaystyle\mathrm{MISE}\\{f\\}(t)$
$\displaystyle=\mathbb{E}_{f_{X}}\left\\{\int_{0}^{1}[f(x,t)-f_{X}(x,t)]^{2}dx\right\\}$
(13)
$\displaystyle=\int_{0}^{1}\mathbb{\\{}\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)\\}^{2}dx+\int_{0}^{1}\mathrm{Var}_{f_{X}}[f(x,t)]dx.$
(14)

Often one is interested in the asymptotic approximation to the MISE, denoted
AMISE, under the requirements that $t=t_{n}\downarrow 0$ and
$n\sqrt{t_{n}}\rightarrow\infty$, which ensure consistency of the estimator.
The asymptotically optimal bandwidth is then the minimizer of the AMISE. For
our continuous model of kernel density estimation we have the following result
(proven in Appendix B) which gives the same $\mathcal{O}(n^{-4/5})$ rate of
convergence as the Gaussian KDE on the whole real line.

######  Theorem 5 (Asymptotic Bias and Variance of Diffusion Estimator).

Let $t_{n}$ be such that $\lim_{n\rightarrow\infty}t_{n}=0$ and
$\lim_{n\rightarrow\infty}n\sqrt{t_{n}}=\infty$ and suppose that $f_{X}\in
C^{2}([0,1])$ (twice continuously differentiable) with $f_{X}(0)=rf_{X}(1)$.
Then the following hold as $n\rightarrow\infty$:

  1. 1.

The integrated variance has the asymptotic behavior

$\int_{0}^{1}\mathrm{Var}_{f_{X}}[f(x,t_{n})]dx\sim\frac{1}{2n\sqrt{\pi
t_{n}}}.$ (15)

  2. 2.

If $f_{X}^{\prime}(0)=f_{X}^{\prime}(1)$ then the integrated squared bias is

$\int_{0}^{1}\left\\{\mathbb{E}_{f_{X}}[f(x,t_{n})]-f_{X}(x)\right\\}^{2}dx\sim
t_{n}^{2}\int_{0}^{1}\frac{1}{4}\left[f_{X}^{{}^{\prime\prime}}(x)\right]^{2}dx.$
(16)

  3. 3.

If $f_{X}^{\prime}(0)\neq f_{X}^{\prime}(1)$ then the integrated squared bias
is

$\int_{0}^{1}\\{\mathbb{E}_{f_{X}}[f(x,t_{n})]-f_{X}(x)\\}^{2}dx\sim
t_{n}^{3/2}\frac{4-2\sqrt{2}}{3\sqrt{\pi}}\frac{r^{2}+1}{(1+r)^{2}}[f_{X}^{\prime}(1)-f_{X}^{\prime}(0)]^{2}.$
(17)

###### Proof.

See Appendix B. ∎

A direct consequence of this result is that we can select the stopping time
$t$ or bandwidth to minimize the AMISE.

######  Corollary 2 (Asymptotically Optimal Bandwidth Choices).

Combining the leading order bias and variance terms gives the asymptotic
approximation to the MISE:

  1. 1.

If $f_{X}^{\prime}(1)=f_{X}^{\prime}(0)$ then

$\mathrm{AMISE}\\{f\\}(t_{n})=\frac{1}{2n\sqrt{\pi
t_{n}}}+t_{n}^{2}\int_{0}^{1}\frac{1}{4}\left[f_{X}^{{}^{\prime\prime}}(x)\right]^{2}dx.$
(18)

Hence, the square of the asymptotically optimal bandwidth is

$t^{*}=(2n\sqrt{\pi}\|f_{X}^{{}^{\prime\prime}}\|_{L^{2}}^{2})^{-2/5}$

with the minimum value

$\min_{t}\mathrm{AMISE}\\{f\\}(t)=\frac{5\|f_{X}^{{}^{\prime\prime}}\|_{L^{2}}^{2/5}}{2^{14/5}\pi^{2/5}}n^{-4/5}.$

  2. 2.

If $f_{X}^{\prime}(1)\neq f_{X}^{\prime}(0)$ then

$\begin{split}\mathrm{AMISE}\\{f\\}(t_{n})&=\frac{1}{2n\sqrt{\pi
t_{n}}}+t_{n}^{3/2}\frac{4-2\sqrt{2}}{3\sqrt{\pi}}\frac{r^{2}+1}{(1+r)^{2}}[f_{X}^{\prime}(1)-f_{X}^{\prime}(0)]^{2}\\\
&=\frac{1}{2n\sqrt{\pi
t}}+t^{3/2}\frac{A(r)}{3}\left[f_{X}^{\prime}(1)-f_{X}^{\prime}(0)\right]^{2}.\end{split}$
(19)

Hence, the square of the asymptotically optimal bandwidth is

$t^{*}=(2n\sqrt{\pi}A(r))^{-1/2}\left|f_{X}^{\prime}(1)-f_{X}^{\prime}(0)\right|^{-1}$

with the minimum value

$\min_{t}\mathrm{AMISE}\\{f\\}(t)=\frac{2^{5/4}\sqrt{\left|f_{X}^{\prime}(1)-f_{X}^{\prime}(0)\right|}}{3\pi^{3/8}}A(r)^{1/4}n^{-3/4}.$

A few remarks are in order. First, it is interesting to note that in the case
of $f_{X}^{\prime}(1)=f_{X}^{\prime}(0)$, the optimum choice $t^{*}$ and the
minimum AMISE do not depend on $r$, and are the same as the more familiar
‘whole line’ situation — in other words, we can confidently use existing
methods in the literature (such as recommended above) to choose a stopping
time. Second, it seems plausible that we could estimate
$f_{X}^{\prime}(1)-f_{X}^{\prime}(0)$ (or the value of $r$) adaptively and
change the boundary conditions in the model (2) accordingly. A full discussion
of solving the heat equation with linked boundary conditions for the first
spatial derivative is beyond the scope of this paper but can be done using the
same methods we present here. Future work will aim to incorporate an adaptive
estimate of the true boundary conditions (both for the density function and
its first derivative - we do this for the density function in §5.1) and
resulting adaptive boundary conditions. We mention a result in this direction
which will appear when we compare our model to that of [5], whose code is
based around the discrete cosine transform, the continuous version of which
solves the heat equation subject to the boundary conditions
$f^{\prime}_{c}(0)=f_{c}^{\prime}(1)=0.$ We have used the subscript $c$ to
avoid confusion with our solution $f$ to (2). The analogous result to Theorem
5 is the following theorem which can be proven using the same techniques and
hence we have omitted the proof. Similarly, one can then derive the optimum
choice of $t$ and the minimum AMISE $\mathcal{O}(n^{-3/4})$ (slower rate)
under the condition that $(f_{X}^{\prime}(1),f_{X}^{\prime}(0))\not=(0,0)$.

######  Theorem 6 (Boundary Effects on Asymptotic Bias).

Let $t_{n}$ be such that $\lim_{n\rightarrow\infty}t_{n}=0$ and also
$\lim_{n\rightarrow\infty}n\sqrt{t_{n}}=\infty$. Suppose that $f_{X}\in
C^{2}([0,1])$. Then the following hold as $n\rightarrow\infty$:

  1. 1.

The integrated variance has the asymptotic behavior

$\int_{0}^{1}\mathrm{Var}_{f_{X}}[f_{c}(x,t_{n})]dx\sim\frac{1}{2n\sqrt{\pi
t_{n}}}.$ (20)

  2. 2.

If $f_{X}^{\prime}(0)=f_{X}^{\prime}(1)=0$ then

$\int_{0}^{1}\\{\mathbb{E}_{f_{X}}[f_{c}(x,t_{n})]-f_{X}(x)\\}^{2}dx\sim
t_{n}^{2}\int_{0}^{1}\frac{1}{4}\left[f_{X}^{{}^{\prime\prime}}(x)\right]^{2}dx.$
(21)

  3. 3.

If $f_{X}^{\prime}(0)\neq 0$ or $f_{X}^{\prime}(1)\neq 0$ then

$\int_{0}^{1}\\{\mathbb{E}_{f_{X}}[f_{c}(x,t_{n})]-f_{X}(x)\\}^{2}dx\sim
t_{n}^{3/2}\frac{4-2\sqrt{2}}{3\sqrt{\pi}}[f_{X}^{\prime}(1)^{2}+f_{X}^{\prime}(0)^{2}].$
(22)

##  4 Numerical Approximations of the PDE Estimator

Before giving numerical examples with the new estimator, we consider practical
methods for solving the PDE (2), in order to evaluate the KDE
$f(x,t)=\frac{1}{n}\sum_{k=1}^{n}K(r;x,X_{k},t),$ on a regular grid. There are
two different practical computational methods to compute the density estimator
based on the PDE (2):

1\. Series Expansion:

Essentially solving the continuous model (2) via the series or contour
integral representation in Theorem 2.

2\. Backward Euler method:

Solving a discretized or binned version of (2), as explained in the rest of
this section. In Theorem 7, we show that this binned estimator converges to
the continuous PDE estimator.

The two methods have relative advantages and disadvantages. The backward Euler
method is a first order finite difference method (however, this is not a
problem in practice as argued below), but it is simple and easy to use,
especially if the initial data is already discretely binned. The backward
Euler method also maintains the key property of positivity and satisfies the
same maximum principle properties as the continuous solution (see Appendix C
and Lemma 3). The reason for not using second order methods such as
Crank–Nicolson is that for large time steps this would not preserve non-
negativity of the solution. In other words, the discrete solution can no
longer be interpreted as a probability distribution (a well-known result says
that any general linear method that is unconditionally positivity preserving
for all positive ODEs must have order $\leq 1$ [4]). However, methods such as
Crank–Nicolson can also easily be used for the discrete model if desired, but
for brevity we do not discuss such methods further. The series expansion of
the continuous PDE model is typically highly accurate for $t>0$, but less easy
to implement. We provide MATLAB codes for both methods:
https://github.com/MColbrook/Kernel-Density-Estimation-with-Linked-Boundary-
Conditions.

To derive the appropriate time-stepping method, we do the following:

  1. 1.

We approximate the exact solution $f$ by a vector $\bm{u}$. That is,
$u(x_{i};\cdot)\approx f(x_{i},\cdot)$. Here $x_{i}=ih$ is the $i$th grid
point on the grid of $m+2$ equally spaced points in the domain $[0,1]$, for
$i=0,1,\ldots,m,m+1$. The spacing between two consecutive grid points is
$h=\frac{1}{m+1}.$ Note here that $m$ is typically smaller than $n$, the
number of samples that form the empirical measure.

  2. 2.

The two boundary conditions in (2) give two equations involving values at the
two boundary nodes, i.e. at node $0$ and at node $m+1$. That is,

$\displaystyle\quad u_{0}$ $\displaystyle=$ $\displaystyle ru_{m+1},$ (23)
$\displaystyle u_{1}-u_{0}$ $\displaystyle=$ $\displaystyle u_{m+1}-u_{m}.$
(24)

This motivates us to make the following definitions for the boundary nodes:

$u_{0}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{r}{r+1}(u_{1}+u_{m}),\qquad
u_{m+1}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{r+1}(u_{1}+u_{m}).$
(25)

We are left with a set of $m$ equations involving $m$ unknown values
$u_{1},\ldots,u_{m}$, at the $m$ interior nodes $1,\ldots,m$, where we use a
standard second-order finite difference approximation of the (spatial) second
derivative.

  3. 3.

We consider the corresponding $m\times m$ _four-corners matrix_ with the
following structure:

$\mathbf{A}=\left(\begin{tabular}[]{ccccc}$2-\frac{r}{r+1}$&$-1$&&&$-\frac{r}{r+1}$\\\
$-1$&$2$&$-1$\\\ &$\ddots$&$\ddots$&$\ddots$\\\ &&$-1$&$2$&$-1$\\\
$-\frac{1}{r+1}$&&&$-1$&$-\frac{1}{r+1}$\end{tabular}\right).$ (26)

Given a time $T$ at which we wish to evaluate the solution, we consider a time
step $\Delta t=2h^{2}$. For ease of the analysis, we assume that $T$ is a
multiple of $\Delta t$, though his can be avoided by making the last time step
smaller if needed. We use a superscript $k$ to denote the solution at time
$k\Delta t$ (i.e. the $k$th step), then the backwards Euler method can be
written as

$\bm{u}^{k+1}=\left(\mathbf{I}+\mathbf{A}\right)^{-1}\bm{u}^{k},\quad
k=0,...,T/\Delta t-1,$ (27)

where $\mathbf{I}$ denotes the $m\times m$ identity matrix. The matrix inverse
can be applied in $\mathcal{O}(m)$ operations using the fact that $\mathbf{A}$
is a rank one perturbation of a tridiagonal matrix. Even though we take small
time steps, the total time $T=\mathcal{O}(n^{-2/5})$ is small. It follows that
the total complexity is $\mathcal{O}(m^{3}n^{-2/5})$, giving an error (in the
interior) of order $\mathcal{O}(h^{2})=\mathcal{O}(m^{-2})$. The error of the
continuous model scales as $\mathcal{O}(n^{-2/5})$. If there is freedom in
selecting the number of bins $m+2$, this suggests choosing
$m=\mathcal{O}(n^{1/5})$ which leads to a modest
$\mathcal{O}(n^{1/5})=\mathcal{O}(m)$ complexity. A key property of the matrix
(26) is that it has zero column sum, off-diagonals are negative or zero, and
the main diagonal entries are positive. This allows the interpretation of (27)
as a discrete-time Markov process. In Appendix C, we prove the following
theorem for completeness (using explicit formulae for the eigenvalues and
eigenvectors of $\mathbf{A}$).

######  Theorem 7 (Convergence of Binned to Diffusion Estimator).

The solution of the binned estimator (27) with the four corner matrix in (26)
converges to the solution of the continuous problem (2) as
$m\rightarrow\infty$:

$\sup_{\epsilon\leq t\leq T}\sup_{0\leq k\leq
m+1}|u(k/(m+1);t)-f(k/(m+1);t)|\rightarrow 0,\qquad n\rightarrow\infty.$

###### Proof.

See Appendix C. ∎

Further interesting properties of the discrete system are discussed in
Appendix C. In Theorem 7, we have restricted $t\geq\epsilon>0$ to include the
possibility that the initial condition may not be a proper function, but an
empirical measure. We finally remark that sometimes the solution is needed at
later times (e.g. $\mathcal{O}(1)$), for example when querying the solution at
various times $t$ as part of minimizing least squares cross validation to
determine a good choice of $T$. In that case, we recommend computing the
matrix exponential

$\textstyle\bm{u}(t)=\exp\left(-\frac{t}{2h^{2}}\mathbf{A}\right)\bm{u}(0).$

There are many possible methods to compute the matrix exponential [48], such
as MATLAB’s expm code based on [31, 2].

##  5 Numerical Experiments

###  5.1 Numerical examples with synthetic data

First, we test the estimator on examples where the true density $f_{X}$ is
known. We begin with the trimodal distribution shown in Figure 3. We will
demonstrate two versions of the method. First, when the exact value of $r$ is
known (labelled “Linked 1”), and second where we estimate the value of $r$ by

$r_{\mathrm{est}}=\frac{\sum_{j=1}^{n}\chi_{<n^{-1/2}}(X_{j})}{\sum_{j=1}^{n}\chi_{>1-n^{-1/2}}(X_{j})}$

(labelled “Linked 2”). We expect both to perform similarly for sufficiently
large $n$. For stopping times, we have used the software that adaptively
chooses the bandwidth, discussed in Section 3. In other words, we do not give
our algorithms any information other than the given sample. We compare with
three other methods. The first is the density estimation proposed in [5] based
on the discrete cosine transform (labelled “Cosine”). The second is the well-
known and arguably state-of-the-art beta kernel method of [6], which we label
“Beta” in the plots. This method is free from boundary bias, at the cost of an
increased boundary variance. Finally, we also compare with a method which uses
copula kernels [33] and which has been found to be competitive with the beta
kernel approach of [6]. This method has an automatic bandwidth selector which
we shall use, and we label it “Copula” in the plots. The latter two methods
are freely available in the R package evmix [32] which can be found at
https://CRAN.R-project.org/package=evmix.

We estimate the error using the $L^{2}$ and $L^{\infty}$ norms at the points
$l\times 10^{-3}$ for $l=0,...,10^{3}$. The only change is when considering
the copula method, where we take $l=1,...,10^{3}-1$ instead since we found
this method to be unstable near the boundaries. Figure 3 shows a typical
approximation of the distribution function using our proposed method and the
other methods for a sample size of $n=10^{4}$. Our proposed method is more
accurate near the boundaries of the domain (see magnified section of plots)
and behaves similarly in the middle of the domain. We found that using the
estimate $r_{\mathrm{est}}$ instead of the exact value of $r$ did not have a
great effect on the error. In other words, we can apply our model without
needing to know the value of $r$.

Figure 4 (left) shows the $L^{2}$ measure of error averaged over $100$
independent samples for each $n$. The $L^{2}$ errors for both “Linked” methods
and the “Cosine” method agreed almost perfectly with the minimum AMISE and the
analysis in Section 3 for large $n$. Using our model with an estimate of $r$
increases the convergence rate from $\mathcal{O}(n^{-3/4})$ to
$\mathcal{O}(n^{-4/5})$. Both “Linked” methods and the “Cosine” method are
found to be more accurate than the “Beta” and “Copula” methods. The tailing-
off convergence for the “Copula” method was due to a need to implement a lower
bound for the bandwidth. Below this limit, we found the “Copula” method to be
unstable. Figure 4 (right) shows the same plot but now for the $L^{\infty}$
measure of error. Here we see a more pronounced difference between the
methods, with both “Linked” methods producing much smaller errors than the
other methods. We found the same behavior in these plots for a range of other
tested distributions. Finally, we comment on the CPU times for each method,
shown in Figure 5 (averaged over the 100 samples for each $n$). In order to
produce a fair comparison, we have included the CPU time taken for automatic
bandwidth selection when using the “Linked” methods. All methods appear to
have CPU times that grow linearly with $n$. The “Cosine” method in fact scales
like $\mathcal{O}(n\log(n))$ due to the use of the discrete cosine transform.
The linked estimator is faster by about an order of magnitude than the other
methods. This is due to the exponential decay of the series for $t>0$ \- only
a small number of terms need to be summed in order to get very accurate
results.

Figure 3: Example of different methods for a sample size $n=10^{4}$. The
proposed diffusion model (“Linked”) is much more accurate near the boundaries
than the cosine model (“Cosine”) as highlighted by the magnified sections. The
method “Copula” is found to be unstable near the boundaries.

Figure 4: Left: $L^{2}$ errors of methods averaged over $100$ samples for each
$n$. Right: $L^{\infty}$ errors of methods averaged over $100$ samples for
each $n$. The $L^{2}$ errors agree well with the minimum AMISE from Section 3,
whereas the increased accuracy gained near the boundary by using the linked
boundary model is highlighted by the $L^{\infty}$ errors. Figure 5: CPU times
for each method averaged over $100$ samples for each $n$. Experiments were
performed on a basic four year old laptop. Each method appears to grow almost
linearly (up to logarithmic factors), with the linked boundary estimator an
order of magnitude faster than the other methods.

$a$ | $1.1$ | $1.2$ | $1.3$ | $1.4$ | $1.5$
---|---|---|---|---|---
Linked | $2.98\text{\times}{10}^{-3}$ | $1.33\text{\times}{10}^{-3}$ | $6.82\text{\times}{10}^{-4}$ | $3.22\text{\times}{10}^{-4}$ | $2.38\text{\times}{10}^{-4}$
LC | $1.05\text{\times}{10}^{-3}$ | $1.14\text{\times}{10}^{-3}$ | $1.26\text{\times}{10}^{-3}$ | $1.38\text{\times}{10}^{-3}$ | $1.52\text{\times}{10}^{-3}$
LCS | $1.23\text{\times}{10}^{-3}$ | $1.03\text{\times}{10}^{-3}$ | $9.42\text{\times}{10}^{-4}$ | $1.04\text{\times}{10}^{-3}$ | $1.19\text{\times}{10}^{-3}$

$a$ | $1.6$ | $1.7$ | $1.8$ | $1.9$ | $2$
---|---|---|---|---|---
Linked | $1.58\text{\times}{10}^{-4}$ | $1.13\text{\times}{10}^{-4}$ | $8.01\text{\times}{10}^{-5}$ | $5.96\text{\times}{10}^{-5}$ | $5.05\text{\times}{10}^{-5}$
LC | $1.65\text{\times}{10}^{-3}$ | $1.80\text{\times}{10}^{-3}$ | $1.94\text{\times}{10}^{-3}$ | $2.09\text{\times}{10}^{-3}$ | $2.27\text{\times}{10}^{-3}$
LCS | $1.30\text{\times}{10}^{-3}$ | $1.40\text{\times}{10}^{-3}$ | $1.66\text{\times}{10}^{-3}$ | $1.74\text{\times}{10}^{-3}$ | $2.16\text{\times}{10}^{-3}$

Table 1: Mean $L^{2}$ squared error over 10 simulations for different $a$.

$a$ | $1.1$ | $1.2$ | $1.3$ | $1.4$ | $1.5$
---|---|---|---|---|---
Linked | $7.32\text{\times}{10}^{-2}$ | $4.19\text{\times}{10}^{-2}$ | $2.52\text{\times}{10}^{-2}$ | $1.31\text{\times}{10}^{-2}$ | $7.97\text{\times}{10}^{-3}$
LC | $5.34\text{\times}{10}^{-1}$ | $6.40\text{\times}{10}^{-1}$ | $7.51\text{\times}{10}^{-1}$ | $8.71\text{\times}{10}^{-1}$ | $1.00$
LCS | $1.84\text{\times}{10}^{-1}$ | $1.26\text{\times}{10}^{-1}$ | $1.42\text{\times}{10}^{-1}$ | $1.72\text{\times}{10}^{-1}$ | $1.95\text{\times}{10}^{-1}$

$a$ | $1.6$ | $1.7$ | $1.8$ | $1.9$ | $2$
---|---|---|---|---|---
Linked | $4.42\text{\times}{10}^{-3}$ | $2.85\text{\times}{10}^{-3}$ | $1.18\text{\times}{10}^{-3}$ | $4.78\text{\times}{10}^{-4}$ | $2.39\text{\times}{10}^{-4}$
LC | $1.14$ | $1.28$ | $1.44$ | $1.60$ | $1.78$
LCS | $2.27\text{\times}{10}^{-1}$ | $2.47\text{\times}{10}^{-1}$ | $2.82\text{\times}{10}^{-1}$ | $3.22\text{\times}{10}^{-1}$ | $3.70\text{\times}{10}^{-1}$

Table 2: Mean $L^{\infty}$ squared error over 10 simulations for different
$a$.

Figure 6: Typical estimates for $n=10^{4}$ and $a=1.1$, $a=2$. We used the R
package logcondens for the log-concave projection method.

Next, we consider the case when $f_{X}$ is log-concave and not necessarily
smooth. Denoting the PDF of the beta distribution with parameters
$(\alpha,\beta)$ by $b(\alpha,\beta;x)$, we let

$f_{X}(x)=\frac{b(1,2;x)+2b(a,1;x)}{3}.$

The parameter $a$ controls the smoothness of $f_{X}$ near $x=0$. We have
compared our method to a method that computes log-concave maximum likelihood
estimators [20, 21]. This seeks to compute the log-concave projection of the
empirical distribution through an active set approach. Code is freely
available in logcondens [22] which can be found at
https://CRAN.R-project.org/package=logcondens. Details on such methods can be
found in [54], with a study of the more involved case of censored data in
[23]. Tables 1 and 2 show the mean squared $L^{2}$ and $L^{\infty}$ errors
respectively over $10$ simulations for $n=10^{5}$, as we vary $a$ for the
linked boundary diffusion estimator and the log-concave projection method
(abbreviated to LC), as well as its smoothed version (LCS). In each case, we
have shaded the most accurate estimator. The linked boundary diffusion
estimator performs much better when measured in the uniform norm but is
slightly worse in the $L^{2}$ sense when the distribution function becomes
less smooth. This is demonstrated in Figure 6 for a typical estimation using
$n=10^{4}$. To produce the tables, the linked boundary diffusion estimator
took about 0.5s on average per simulation, the log-concave projection took
about 5s, but its smoothed version was much slower, taking about 73s.

###  5.2 Numerical example with cell data

Figure 7: (A) Schematic cell cycle with geminin expression starting at the end
of G1. (B) DNA and geminin signal from individual cells can be used to obtain
a pseudo-temporal ordering of the population. An average cell follows the
indicated path (red) through the dataset. (C) Pseudotime values (gray), binned
data (blue) and kernel density estimate (red). The kernel density estimate was
obtained by solving our continuous PDE (2) by our discrete numerical method
with the ‘four corners matrix’ in (26). The stopping time, $t=0.00074$, came
from the stopping time software of one of the authors:
https://au.mathworks.com/matlabcentral/fileexchange/14034-kernel-density-
estimator.

This section demonstrates the application of the methods that we propose to a
problem in biology with the data taken from [39]. As mentioned in the
introduction, Figure 1 shows an example of what goes wrong when current
methods are applied. Figure 7 C demonstrates our proposed method, which
successfully incorporates the desired linked boundary condition.

This example originates from the study of biological processes, in particular,
cell cycle studies in cancer research (Figure 7 A). A recently developed
theory [36, 39] which relies on the distribution of cells along the cell cycle
enables the study of entire cell cycle progression kinetics. The method
utilizes data from single cell experiments like flow cytometry or single cell
RNA sequencing, where the abundance of up to thousands of cellular components
for every individual cell in a population is measured. Cells in a
unsynchronized cell population are spread over all stages of the cell cycle,
which can be seen in the exemplary dataset where levels of DNA and geminin in
single cells were measured by flow cytometry (Figure 7 B). The red curve in
Figure 7 B indicates the path that the average cell takes when it goes through
the cell cycle. Pseudotime algorithms perform a dimensionality reduction by
assigning a pseudotime value to each cell, which can be interpreted as its
position on the average curve. In this example, the pseudotime is a
quantitative value of the progression through the cell cycle. However, it is
in general not equal to real time. As the number of cells in a particular
stage is related to the average transit time through that stage, one can
derive a mapping from pseudotime to real time based on ergodic principles [36,
39, 38]. This mapping relies on the distribution of cells on the pseudotime
scale. As mentioned in the introduction, the distribution at the beginning and
the end of the cell cycle are linked due to cell division by

$f(0,t)=2\,f(1,t)\;.$ (28)

Ignoring this fact when estimating the density on the pseudotime scale results
in an erroneous transformation and thus inaccurate kinetics. The KDE with
linked boundary condition ($r=2$) produces a distribution that satisfies the
conditions (28) on the density due to cell division (Figure 7 C). The MAPiT
toolbox for single-cell data analysis [38] applies our new KDE with linked
boundary conditions to analyze cell cycle dependent molecular kinetics.

##  6 Conclusion

Our study was motivated by a dataset from a biological application. This
biological application required a method of density estimation that can handle
the situation of linked boundaries, which are crucial for gaining correct
kinetics. More broadly, boundary bias issues are known to be a difficult
problem in the context of kernel density estimation. To our knowledge, the
linked boundary conditions that we handle here have not been previously
addressed. We have proposed a new diffusion KDE that can successfully handle
the linked boundary conditions. By using the unified transform, we obtained an
explicit solution. In particular, we proved that this diffusion estimator is a
bona fide probability density, which is also a consistent estimator at the
linked boundaries, and derived its asymptotic integrated squared bias and
variance (which shows an increase in the rate of convergence with sample
size).

We also proposed two numerical methods to compute the estimator — one is based
on its series or integral representation and the other on the backward Euler
method. We proved that the discrete/binned estimator converges to the
continuous estimator. We found the new method competes well with other
existing methods, including state-of-the-art methods designed to cope with
boundary bias, both in terms of speed and accuracy. In particular, the new
method is more accurate close to the boundary. Our new KDE with linked
boundary conditions is now used in the MAPiT toolbox for single-cell data
analysis [38] to analyze cell cycle dependent molecular kinetics.

There remain some open questions regarding the proposed models. First, it is
possible to adapt the methods in this paper to multivariate distributions.
Second, it is possible to adapt these methods to other types of boundary
conditions such as constraints on the moments of the distribution (and other
non-local conditions). In this regard, we expect that the flexibility of the
unified transform in PDE theory will be useful in designing smoothing kernel
functions with the desired statistical properties.

#### Acknowledgments & Contributions:

MJC was supported by EPSRC grant EP/L016516/1. ZIB was supported by ARC grant
DE140100993. KK was supported by DFG grant AL316/14-1 and by the Cluster of
Excellence in Simulation Technology (EXC 310/2) at the University of
Stuttgart. SM was supported by the ARC Centre of Excellence for Mathematical
and Statistical Frontiers (ACEMS). MJC performed the theoretical
PDE/statistical analysis of both the continuous and discrete models, and the
numerical tests. SM developed and tested the binned version of the estimator.
ZIB proposed the PDE model and assisted MJC and SM in writing the paper. KK
provided the cell data and assisted in the writing of the numerical section.
MJC is grateful to Richard Samworth, Tom Trogdon and David Smith for comments,
and to Arieh Iserles for introducing him to the problem. The authors are
grateful to the referees for comments that improved the manuscript.

##  Appendix A Proofs of Results in Section 2

###  A.1 Formal derivation of solution formula

We begin with a formal description of how to obtain the solution formulae in
Theorem 2. The most straightforward way to construct the solution is via the
unified transform, and the following steps provide a formal solution which we
must then rigorously prove is indeed a solution.

The first step is to write the PDE in divergence form:

$[\exp(-ikx+k^{2}t/2)f]_{t}-\frac{1}{2}[\exp(-ikx+k^{2}t/2)(f_{x}+ikf)]_{x}=0,\quad
k\in\mathbb{C}.$

We will employ Green’s theorem,

$\textstyle\iint_{\Omega}\Big{(}\frac{\partial F}{\partial x}-\frac{\partial
G}{\partial y}\Big{)}dxdy=\int_{\partial\Omega}\big{(}Gdx+Fdy\big{)},$ (29)

over the domain $(0,1)\times(0,t)$. Here one must assume apriori estimates on
the smoothness of the solution $f$ which will be verified later using the
candidate solution. Define the transforms:

$\displaystyle\textstyle\hat{f}_{0}(k):=\int_{0}^{1}\exp(-ikx)f_{0}(x)dx,\quad$
$\displaystyle\textstyle\hat{f}(k,t):=\int_{0}^{1}\exp(-ikx)f(x,t)dx,$
$\displaystyle\textstyle\tilde{g}(k,t):=\int_{0}^{t}\exp(k\tau)f(1,\tau)d\tau,\quad$
$\displaystyle\textstyle\tilde{h}(k,t):=\int_{0}^{t}\exp(k\tau)f_{x}(1,\tau)d\tau,$

where again we assume these are well defined. Green’s theorem and the boundary
conditions imply (after some small amount of algebra) the so called ‘global
relation’, coupling the the transforms of the solution and initial data:

$\begin{split}\hat{f}(k,t)\exp(k^{2}t/2)=&\textstyle\hat{f}_{0}(k)-\frac{1}{2}[\tilde{h}(k^{2}/2,t)+ikr\tilde{g}(k^{2}/2,t)]\\\
&+\textstyle\frac{\exp(-ik)}{2}[\tilde{h}(k^{2}/2,t)+ik\tilde{g}(k^{2}/2,t)],\quad
k\in\mathbb{C}.\end{split}$ (30)

The next step is to invert via the inverse Fourier transform, yielding

$\begin{split}f(x,t)=\textstyle\frac{1}{2\pi}\int_{-\infty}^{\infty}&\exp(ikx-k^{2}t/2)\big{\\{}\hat{f}_{0}(k)-\frac{1}{2}[\tilde{h}(k^{2}/2,t)+ikr\tilde{g}(k^{2}/2,t)]\\\
&+\textstyle\frac{\exp(-ik)}{2}[\tilde{h}(k^{2}/2,t)+ik\tilde{g}(k^{2}/2,t)]\big{\\}}dk.\end{split}$
(31)

However, this expression contains the unknown functions $\tilde{g}$ and
$\tilde{h}$. To get rid of these we use some complex analysis and symmetries
of the global relation (30). Define the domains

$D^{+}=\\{k\in\mathbb{C}^{+}:\mathrm{Re}(k^{2})<0\\},\quad
D^{-}=\\{k\in\mathbb{C}^{-}:\mathrm{Re}(k^{2})<0\\},\quad D=D^{+}\cup D^{-}.$
(32)

These are shown in Figure 8. A quick application of Cauchy’s theorem and
Jordan’s lemma means we can re-write our solution as

$\begin{split}f(x,t)=&\textstyle\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp(ikx-k^{2}t/2)\hat{f}_{0}(k)dk\\\
&\textstyle-\frac{1}{2\pi}\int_{\partial
D^{+}}\frac{\exp(ikx-k^{2}t/2)}{2}[\tilde{h}(k^{2}/2,t)+ikr\tilde{g}(k^{2}/2,t)]dk\\\
&\textstyle-\frac{1}{2\pi}\int_{\partial
D^{-}}\frac{\exp(ik(x-1)-k^{2}t/2)}{2}[\tilde{h}(k^{2}/2,t)+ik\tilde{g}(k^{2}/2,t)]dk.\end{split}$
(33)

We now use the symmetry under $k\rightarrow-k$ of the global relation (30) and
the fact that the argument in each of $\tilde{g}$ and $\tilde{h}$ is $k^{2}/2$
to set up the linear system:

$\begin{split}\frac{1}{2}\begin{pmatrix}[\exp(-ik)-1]&ik[\exp(-ik)-r]\\\
[\exp(ik)-1]&-ik[\exp(ik)-r]\end{pmatrix}\begin{pmatrix}\tilde{h}(\frac{k^{2}}{2},t)\\\
\tilde{g}(\frac{k^{2}}{2},t)\end{pmatrix}=\begin{pmatrix}\hat{f}(k,t)\exp(\frac{tk^{2}}{2})-\hat{f}_{0}(k)\\\
\hat{f}(-k,t)\exp(\frac{tk^{2}}{2})-\hat{f}_{0}(-k)\end{pmatrix}.\end{split}$

Defining the determinant function $\Upsilon(k)=2(1+r)(\cos(k)-1),$ solving the
linear system leads to the relations:

$\displaystyle\textstyle\frac{\tilde{h}(k^{2},t)+ikr\tilde{g}(k^{2},t)}{2}$
$\displaystyle=\textstyle\frac{1}{\Upsilon(k)}\Big{\\{}\hat{f}_{0}(k)[(1+r)\exp(ik)-2r]$
$\displaystyle\quad\quad\quad\quad+\textstyle\hat{f}_{0}(-k)(1-r)\exp(-ik)$
$\displaystyle\quad\quad\quad\quad\quad-\textstyle\exp(k^{2}t/2)\hat{f}(k,t)[(1+r)\exp(ik)-2r]$
$\displaystyle\quad\quad\quad\quad\quad\quad-\textstyle\exp(k^{2}t/2)\hat{f}(-k,t)(1-r)\exp(-ik)\Big{\\}},$
$\displaystyle\textstyle\frac{\tilde{h}(k^{2}/2,t)+ik\tilde{g}(k^{2}/2,t)}{2}$
$\displaystyle=\frac{1}{\Upsilon(k)}\Big{\\{}\hat{f}_{0}(k)[2\exp(ik)-(1+r)]+\hat{f}_{0}(-k)(1-r)$
$\displaystyle\textstyle\quad\quad\quad\quad-\exp(k^{2}t/2)\hat{f}(k,t)[2\exp(ik1)-(1+r)]$
$\displaystyle\textstyle\quad\quad\quad\quad\quad-\exp(k^{2}t/2)\hat{f}(-k,t)(1-r)\Big{\\}}.$

Since $\Upsilon(k)$ is zero whenever $\cos(k)=1$, before we substitute these
relations into our integral solution we deform the contours $\partial D^{+}$
and $\partial D^{-}$ as shown in Figure 8 to avoid the poles of
$\Upsilon(k)^{-1}$ along the real line.

Figure 8: Left: The domains $D^{\pm}$ as well as the orientation of the
boundaries $\partial D^{\pm}$. Right: The deformed contours to avoid the
singularity at $k=0$. The bold arrow shows a path on which both the $x$ and
$t$ exponential parts of the integrand are exponentially decaying which can be
used for efficient numerical evaluation.

Upon substitution, we are still left with unknown contributions proportional
to

$\displaystyle\textstyle I_{1}(x,t):=\int_{\partial
D^{+}}\frac{\exp(ikx)}{\Upsilon(k)}\big{\\{}\hat{f}(k,t)[(1+r)\exp(ik)-2r]\\!+\hat{f}(-k,t)(1-r)\exp(-ik)\big{\\}}dk$
$\displaystyle\textstyle I_{2}(x,t):=\int_{\partial
D^{-}}\frac{\exp(ik(x-1))}{\Upsilon(k)}\big{\\{}\hat{f}(k,t)[2\exp(ik)-(1+r)]+\hat{f}(-k,t)(1-r)\big{\\}}dk.$

We will argue that the integral $I_{1}(x,t)$ along $\partial D^{+}$ vanishes
and the argument for $I_{2}(x,t)$ follows the same reasoning. First observe
that as $k\rightarrow\infty$ in $\mathbb{C}^{+}$,
$\Upsilon(k)^{-1}\sim\exp(ik)/(1+r)$. Also, we must have that

$\displaystyle\textstyle\exp(ik)\hat{f}(k,t)=\int_{0}^{1}\exp(ik(1-x))f(x,t)dx$

is bounded in $\mathbb{C}^{+}$. $\hat{f}(-k,t)$ is also bounded in
$\mathbb{C}^{+}$ and hence the function

$\textstyle\frac{\hat{f}(k,t)[(1+r)\exp(ik)-2r]+\hat{f}(-k,t)(1-r)\exp(-ik)}{\Upsilon(k)}$

is bounded in $\mathbb{C}^{+}$. It follows that we can close the contour in
the upper half plane and use Jordan’s lemma to see that $I_{1}(x,t)$ vanishes.
We then obtain the integral form of the solution in Theorem 2.

To obtain the series form we can write
$2\exp(ik)-(1+r)=-\exp(ik)\Upsilon(k)+\exp(ik)[(1+r)\exp(ik)-2r],$ which
implies

$\displaystyle\textstyle\int_{\partial
D^{-}}\frac{\exp(ik(x-1)-k^{2}t/2)}{\Upsilon(k)}\hat{f}_{0}(k)[2\exp(ik)-(1+r)]dk=$
$\displaystyle\textstyle\int_{\partial
D^{-}}\exp(ikx-k^{2}t/2)\hat{f}_{0}(k)dk-\\!\\!\int_{\partial
D^{-}}\\!\\!\\!\frac{\exp(ikx-k^{2}t/2)}{{\Upsilon(k)}}\hat{f}_{0}(k)[(1+r)\exp(ik)-2r]dk.$

Taking into account the orientation of $\partial D^{-}$, upon deforming the
first of these integrals back to the real line, we see that it cancels the
first integral in (6). Hence we have

$\textstyle 2\pi f(x,t)=-\int_{\partial
D}\\!\\!\\!\frac{e^{ikx-k^{2}t/2}}{{\Upsilon(k)}}\big{\\{}\hat{f}_{0}(k)[(1+r)e^{ik}-2r]+\hat{f}_{0}(-k)(1-r)e^{-ik}\big{\\}}dk.$
(34)

Define the function

$\displaystyle\textstyle
F(x,t;k):=\frac{e^{ikx-k^{2}t/2}\big{\\{}\hat{f}_{0}(k)[(1+r)e^{ik}-2r]+\hat{f}_{0}(-k)(1-r)e^{-ik}\big{\\}}}{2(1+r)}.$

The integrand in (34) has a double pole at $k_{n}=2n\pi$ so we deform the
contour $\partial D$ to $\partial\tilde{D}$ shown in Figure 9. Cauchy’s
residue theorem then implies that

$\textstyle f(x,t)=-\frac{1}{2\pi}\int_{\partial
D}\frac{F(x,t;k)}{\cos(k)-1}dk=\sum_{n\in\mathbb{Z}}-2iF^{\prime}(k_{n}).$
(35)

It is then straightforward to check the equality of (35) and (7).

Figure 9: Deformation of the contour to circle the poles. The contributions
along the real line between these circles cancel.

###  A.2 Proof of Theorems 1 and 2

###### Proof of Theorems 1 and 2.

For $t>0$, it is clear that the function $f$ given by (6) is smooth in $x,t$
and real analytic in $x$, as well as solving the heat equation. This follows
from being able to differentiate under the integral sign due to the
$\exp(-k^{2}t/2)$ factor and the fact that extending $x$ to a complex argument
yields an analytic function. Note also that the argument in Section A.1 does
rigorously show equivalence between the series and integral forms of $f$. It
is easy to check via the series (7) that the function $f$ satisfies the
required boundary conditions and hence (5) also holds by simple integration by
parts. Regarding the convergence properties as $t\downarrow 0$ when extra
regularity of the initial condition is assumed, Proposition 2 deals with the
case of continuous $f_{0}$, whilst Proposition 3 deals with $f_{0}\in
L^{p}([0,1])$ for $1\leq p<\infty$.

Hence there are two things left to prove; the fact that
$\mu_{t}:=f(\cdot,t)dx$ lies in $C^{w}(0,T;M([0,1]))$ as well as uniqueness in
$C^{w}(0,T;M([0,1]))$ (and $C(0,T;L^{p}([0,1]))$ for $1\leq p<\infty$).

To prove that $\mu_{t}\in C^{w}(0,T;M([0,1]))$, let $g\in C([0,1])$ and
consider the integral kernel defined by (9). By Fubini’s theorem we have

$\textstyle\int_{0}^{1}f(x,t)g(x)dx=\int_{0}^{1}\int_{0}^{1}K(r;x,y,t)g(x)dxdf_{0}(y).$

By Proposition 2 the integral

$\textstyle\int_{0}^{1}K(r;x,y,t)g(x)dx$

converges for all $x$ and is uniformly bounded as $t\downarrow 0$. We will use
the explicit calculation of the endpoints limits at $x=0,1$. By the dominated
convergence theorem, we have

$\displaystyle\textstyle\lim_{t\downarrow 0}$
$\displaystyle\textstyle\int_{0}^{1}f(x,t)g(x)dx=\big{(}\frac{r}{1+r}g(0)+\frac{1}{1+r}g(1)\big{)}f_{0}(\\{0\\})$
$\displaystyle+\textstyle\big{(}\frac{r}{1+r}g(0)+\frac{1}{1+r}g(1)\big{)}f_{0}(\\{1\\})+\int_{x\in(0,1)}g(x)df_{0}(x)$
$\displaystyle=\textstyle
f_{0}(g)+\frac{g(1)-g(0)}{1+r}[f_{0}(\\{0\\})-rf_{0}(\\{1\\})]=f_{0}(g),$

which proves the required weak continuity.

To prove uniqueness, suppose that there exists $\mu_{t},\tau_{t}\in M([0,1])$
which are both weak solutions with $\mu_{0}=\tau_{0}=f_{0}$. Set
$m_{t}=\mu_{t}-\tau_{t}$. We will consider expansions of functions in the
generalized eigenfunctions of the adjoint problem. It is straightforward to
check that the adjoint problem (with the boundary conditions in (4)) is
Birkhoff regular and hence the generalized eigenfunctions are complete in
$L^{2}([0,1])$. In fact we can show that any continuous function $g\in
C([0,1])$ of bounded variation with $g(0)=g(1)$ can be approximated uniformly
by linear combinations of these functions. This follows by either arguing as
we did in the proof of Proposition 2 (the case of non-matching derivatives
holds but is more involved) or follows from Theorem 7.4.4 of [46]. Now suppose
that $\lambda$ lies in the spectrum of the adjoint $\mathbb{A}^{*}$ defined by

$\textstyle\mathbb{A}^{*}=-\frac{d^{2}}{dx^{2}},\quad\mathcal{D}(\mathbb{A}^{*})=\\{u\in
H^{2}([0,1]):u_{x}(1)=ru_{x}(0),u(0)=u(1)\\}.$ (36)

In our case, the generalized eigenfunctions associated with $\lambda$
correspond to a basis of $\mathcal{N}((\mathbb{A}-\lambda I)^{l})$ where $l=1$
or $2$. If $l=2$, and the nullity of $(\mathbb{A}-\lambda I)^{2}$ is greater
than $\mathbb{A}-\lambda I$, we can choose a basis $\\{g_{1},g_{2}\\}$ such
that $(\mathbb{A}-\lambda I)g_{2}=g_{1}$. For the general case and chains of
generalized eigenfunctions, we refer the reader to [46]. Now suppose that
$g\in\mathcal{N}(\mathbb{A}-\lambda I)$, then $g$ must be smooth on $[0,1]$.
It follows that for $t>0$

$\textstyle 2\frac{d}{dt}m_{t}(g)=-\lambda m_{t}(g).$

Note that $m_{0}(g)=0$ and hence we must have that $m_{t}(g)=0$ for all $t\geq
0$. Similarly, suppose that
$\\{g_{1},g_{2}\\}\subset\mathcal{N}((\mathbb{A}-\lambda I)^{2})$ with
$(\mathbb{A}-\lambda I)g_{2}=g_{1}$. Then by the above reasoning we have
$m_{t}(g_{1})=0$ for all $t\geq 0$ and hence

$\textstyle 2\frac{d}{dt}m_{t}(g_{2})=-\lambda
m_{t}(g_{2})-m_{t}(g_{1})=-\lambda m_{t}(g_{2}).$

Again we see that $m_{t}(g_{2})=0$ for all $t\geq 0$. Though we don’t have to
consider it in our case, it is clear that the same argument would work for
chains of longer lengths. The expansion theorem discussed above together with
the dominated convergence theorem shows that if $g\in C([0,1])$ of bounded
variation with $g(0)=g(1)=0$, then $m_{t}(g)=0$ for all $t\geq 0$. This
implies that if $U\subset(0,1)$ is open then $m_{t}(U)=0$ for all $t\geq 0$.
In particular, we must have

$m_{t}=a(t)\delta_{0}+b(t)\delta_{1}$

with $a,b$ continuous. In fact, for any $f\in\mathcal{F}(r)$ we have

$\textstyle\frac{d}{dt}[a(t)+b(t)]f(1)=\frac{a(t)}{2}f_{xx}(0)+\frac{b(t)}{2}f_{xx}(1),$

from which we easily see that $a=b=0$ and hence uniqueness follows. This also
shows uniqueness in the space $C(0,\mathbb{A};L^{p}([0,1]))$, where no
argument at the endpoints is needed. ∎

###  A.3 Proof of Theorem 3

The proof requires that we study the solution of the PDE as $t\downarrow 0$.
We break down the proof into a number of smaller results, which allows us to
use them elsewhere. Recall the definition in (9). We shall also need the
function

$\textstyle K_{1}(x,t):=\sum_{n\in\mathbb{Z}}{\exp(ik_{n}x-k_{n}^{2}t/2)},$
(37)

defined for $t>0$. Using the Poisson summation formula, we can write $K_{1}$
as

$\textstyle K_{1}(x,t)=\frac{1}{\sqrt{2\pi
t}}\sum_{n\in\mathbb{Z}}\exp\Big{(}-\frac{(x-n)^{2}}{2t}\Big{)},$

a periodic summation of the heat kernel. The following lemma is well-known and
hence stated without proof.

######  Lemma 1.

Let $w\in C([0,1])$ then

$\textstyle\int_{0}^{1}K_{1}(x-y,t)w(y)dy$ (38)

is bounded by $\|w\|_{\infty}$ and converges pointwise to $w(x)$ for any
$x\in(0,1)$ and to $(w(0)+w(1))/2$ for $x=0,1$ as $t\downarrow 0$. If
$w(0)=w(1)$ then (38) converges to $w(x)$ uniformly over the interval $[0,1]$.

We will also need the following.

######  Lemma 2.

Let $f_{0}\in L^{1}([0,1])$, then

$\textstyle
t\sum_{n\in\mathbb{N}}\left|\exp(-k_{n}^{2}t/2)k_{n}\hat{f}_{0}(k_{n})\right|\rightarrow
0\quad\text{ as }\quad t\downarrow 0.$

###### Proof.

By the Riemann–Lebesgue lemma, we have that
$\lim_{n\rightarrow\infty}\hat{f}_{0}(k_{n})=0$. So given $\epsilon>0$, let
$N$ be large such that if $n\geq N$ then
$\left|\hat{f}_{0}(k_{n})\right|\leq\epsilon$. Then

$\textstyle
t\sum_{n>N}\left|\exp(-k_{n}^{2}t/2)k_{n}\hat{f}_{0}(k_{n})\right|\leq\frac{t\epsilon}{2\pi}\sum_{n>N}\exp(-2n^{2}\pi^{2}t)4n\pi^{2}.$

Let $h=2\pi$. The sum is an approximation of the integral
$\int_{h(N+1)}^{\infty}\exp(-y^{2}t/2)ydy$ and we have

$t\sum_{n>N}\left|\exp(-k_{n}^{2}t/2)k_{n}\hat{f}_{0}(k_{n})\right|\leq\frac{\epsilon}{2\pi}\int_{0}^{\infty}\exp(-y^{2}t/2){t(y+2h)}dy<\tilde{C}\epsilon,$

for some constant $\tilde{C}$. It follows that

$\limsup_{t\downarrow
0}t\sum_{n\in\mathbb{N}}\left|\exp(-k_{n}^{2}t/2)k_{n}\hat{f}_{0}(k_{n})\right|\leq\tilde{C}\epsilon.$

Since $\epsilon>0$ was arbitrary, the lemma follows. ∎

The following Proposition then describes the limit properties of our
constructed solution as $t\downarrow 0$ in the case of continuous initial
data.

######  Proposition 2.

Let $f_{0}\in C([0,1])$ and $K$ be given by (9). For $t\in(0,1]$ define

$\textstyle f(x,t):=\int_{0}^{1}K(r;x,y,t)f_{0}(y)dy,\quad
q(x,t):=\int_{0}^{1}K(r;y,x,t)f_{0}(y)dy,$

(note the interchange of $x,y$ as arguments of $K$ for the definition of $q$).
Then there exists a constant $C$ (dependent on $r$) such that

$\sup_{x\in[0,1],t\in(0,1]}\max\\{\left|f(x,t)\right|,\left|q(x,t)\right|\\}\leq
C\|f_{0}\|_{\infty}.$ (39)

Furthermore,

$\displaystyle\lim_{t\downarrow 0}f(x,t)$
$\displaystyle=\begin{cases}f_{0}(x),\quad x\in(0,1)\\\
\frac{r}{1+r}[f_{0}(0)+f_{0}(1)],\quad x=0\\\
\frac{1}{1+r}[f_{0}(0)+f_{0}(1)],\quad x=1\end{cases}$
$\displaystyle\lim_{t\downarrow 0}q(x,t)$
$\displaystyle=\begin{cases}f_{0}(x),\quad x\in(0,1)\\\
\frac{1}{1+r}[rf_{0}(0)+f_{0}(1)],\quad x=0,1\end{cases}.$

Finally, in the case that $f_{0}(0)=rf_{0}(1)$, $f(x,t)$ converges to
$f_{0}(x)$ uniformly over $x\in[0,1]$ as $t\downarrow 0$.

###### Proof.

We can write

$\begin{split}\textstyle K(r;x,y,t)&=\textstyle
K_{1}(x-y,t)\big{[}1+(x-y)\frac{1-r}{1+r}\big{]}+K_{1}(x+y,t)(x+y-1)\frac{1-r}{1+r}\\\
&\quad\quad\quad\quad+\frac{t(1-r)}{1+r}\big{[}K_{1}^{\prime}(x+y,t)+K_{1}^{\prime}(x-y,t)\big{]}.\end{split}$
(40)

Here ′ means the derivative with respect to the spatial variable.

To study the limit as $t\downarrow 0$, we note that we can ignore the terms
with a factor of $t$ using Lemma 2. By a change of variable we have

$\textstyle\int_{0}^{1}K_{1}(x+y,t)(x+y-1)f_{0}(y)dy=\int_{0}^{1}K_{1}(x-y,t)(x-y)f_{0}(1-y)dy.$

The bound (39) now follows from Lemma 1, as do the pointwise limits from a
straightforward somewhat tedious calculation.

Now suppose that $f_{0}(0)=rf_{0}(1)$ and split the initial data as follows:

$f_{0}(x)=f_{0}(0)+x(1-r)f_{0}(1)+p_{0}(x).$ (41)

Then $p_{0}\in C([0,1])$ with the crucial property that $p_{0}(0)=p_{0}(1)=0$.
Arguing as above and using Lemma 1, we see that the following limit holds
uniformly

$\textstyle\lim_{t\downarrow 0}\int_{0}^{1}K(r;x,y,t)p_{0}(y)dy=p_{0}(x).$

So it only remains to show that

$\textstyle\int_{0}^{1}K(r;x,y,t)[f_{0}(0)+y(1-r)f_{0}(1)]dy=f_{0}(0)+x(1-r)f_{0}(1).$
(42)

Let $l(x)=f_{0}(0)+x(1-r)f_{0}(1)$ and set $a=f_{0}(1)$. An explicit
calculation yields

$\textstyle\hat{l}(k)=\frac{i}{k}(\exp(-ik)-r)a+\frac{1}{k^{2}}(\exp(-ik)-1)a(1-r).$

We then have

$\displaystyle\textstyle\hat{l}_{0}(k)[(1+r)\exp(ik)-2r]+\hat{l}_{0}(-k)(1-r)\exp(-ik)=$
$\displaystyle\textstyle-a\big{[}\frac{ri}{k}+\frac{1-r}{k^{2}}\big{]}\Upsilon(k).$

We can then apply the residue theorem to the representation (34) to obtain
(42). ∎

In the case where the true density is not continuous but belongs to
$L^{p}([0,1])$ for $p\geq 1$, we have the following.

######  Proposition 3.

Let $1\leq p<\infty$, $f_{0}\in L^{p}([0,1])$ and $K$ be given by (9). For
$t\in(0,1]$ define

$\displaystyle f(x,t):=$
$\displaystyle\textstyle\int_{0}^{1}K(r;x,y,t)f_{0}(y)dy.$

Then $f(\cdot,t)$ converges to $f_{0}$ in $L^{p}([0,1])$ as ${t\downarrow 0}$.

###### Proof.

Note that the case $r=1$ is well-known. The fact that $f_{0}\in L^{1}([0,1])$
by Hölder’s inequality together with Lemma 2 show that we can ignore the parts
multiplied by $t$ in the kernel representations (40). The fact that
$yf_{0}(y)\in L^{p}([0,1])$ implies the convergence by simply summing the
parts in (40) and using the $r=1$ case with a change of variable for the
$K_{1}(x+y,t)$ part. ∎

###### Proof of Theorem 3.

We have that

$\mathbb{E}_{f_{X}}(f(x,t))=\frac{1}{n}\sum_{k=1}^{n}\int_{0}^{1}K(r;x,y,t)dy=\int_{0}^{1}K(r;x,y,t)dy.$

The first part of the theorem therefore follows from Proposition 2. For the
second part, assume that $f_{X}\in C^{1}([0,1])$ and
$x_{t}=x+\mathcal{O}(\sqrt{t})$. The relation (49) implies the result since we
have that

$\int_{0}^{t}\left|K_{1}(x,s)\right|ds\leq C\sqrt{t}$

for some $C$ independent of $x$ and
$\left|f_{X}(x)-f_{X}(x_{t})\right|=\mathcal{O}(\sqrt{t})$. ∎

###  A.4 Proof of Proposition 1

###### Proof of Proposition 1.

We first show that in this case the solution is continuous on
$[0,1]\times[0,T)$ for any $T\in(0,\infty]$. The case of continuity at points
$t>0$ has already been discussed so suppose that
$(x_{n},t_{n})\rightarrow(x,0)$ then

$\left|f(x_{n},t_{n})-f_{0}(x)\right|\leq\left|f_{0}(x_{n})-f_{0}(x)\right|+\left|f(x_{n},t_{n})-f_{0}(x_{n})\right|.$

The first term converges to zero by continuity of $f_{0}$ whilst the second
term converges to zero by the proven uniform convergence as $t\downarrow 0$.
Using the limit given by Proposition 1, we will take $T=\infty$ without loss
of generality.

Since the solution is regular in the interior and continuous on the closure,
this immediately means that we can apply the maximum principle to deduce that

$\sup_{(x,t)\in\overline{\Omega}}f(x,t)=\sup_{(x,t)\in\partial{}\Omega}f(x,t),$

where $\Omega=(0,1)\times(0,T)$. A similar result holds for the infinum.
Evaluating (35) at $x=0$ leads to

$\textstyle
f(0,t)=\frac{2r}{1+r}\sum_{n\in\mathbb{Z}}{\exp(-k_{n}^{2}t/2)}\hat{f}_{0}(k_{n})=\frac{r\sqrt{2}}{(1+r)\sqrt{\pi
t}}\int_{-\infty}^{\infty}\exp(-x^{2}/(2t))f_{0}(y)dy,$

where we have used the function $K_{1}$ defined by (37) and extended $f_{0}$
periodically (values at the endpoints contributed nothing). Hence

$\textstyle\frac{2ar}{1+r}\leq f(0,t)\leq\frac{2br}{1+r}.$

Similar calculations yield

$\textstyle\frac{2a}{1+r}\leq f(1,t)\leq\frac{2b}{1+r}.$

The fact $\max\\{2r/(1+r),2/(1+r)\\}\geq 1$ (recall $r\geq 0$) finishes the
proof. ∎

###  A.5 Proof of Theorem 4

###### Proof of Theorem 4.

We have that

$\int_{0}^{1}f(x,t)dx=\int_{0}^{1}\int_{0}^{1}K(r;x,y,t)df_{0}(y)dx=\int_{0}^{1}\int_{0}^{1}K(r;x,y,t)dxdf_{0}(y)$

by Fubini’s theorem. Using the series representation (40) and integrating term
by term (justified due to the exponential decaying factors) we have

$\int_{0}^{1}K(r;x,y,t)dx=1+\frac{1-r}{1+r}\sum_{n\in\mathbb{Z}}\exp(-k_{n}^{2}t/2)\int_{0}^{1}\big{[}x\exp(ik_{n}(x-y))+(x-1)\exp(ik_{n}(x+y))\big{]}dx.$

All other terms vanish since the integral of $\exp(ik_{n}x)$ is $0$ unless
$n=0$. We can change variables for the second term in the integrand to see
that the above is equal to

$\displaystyle
1+\frac{1-r}{1+r}\sum_{n\in\mathbb{Z}}\exp(-k_{n}^{2}t/2)\int_{0}^{1}\big{[}x\exp(ik_{n}(x-y))-x\exp(-ik_{n}(x-y))\big{]}dx$
$\displaystyle=$ $\displaystyle
1+\frac{1-r}{1+r}\sum_{n\in\mathbb{Z}}\exp(-k_{n}^{2}t/2)\int_{0}^{1}2ix\sin(k_{n}(x-y))dx=1,$

where we have used the fact that $\sin(k_{n}(x-y))$ is odd in $k_{n}$ and
$\exp(-k_{n}^{2}t/2)$ is even in the last equality. Since $f_{0}$ is a
probability measure, it follows that $\int_{0}^{1}f(x,t)dx=1$, i.e. part (1)
holds.

We next show that the integral kernel $K(r;x,y,t)$ is non-negative for $r\geq
0,t>0$ and $x,y\in[0,1]$. Suppose this were false for some
$(x_{0},y_{0})\in[0,1]^{2}$. The Poisson summation formula gives

$K(r;0,y,t)=\frac{2r}{1+r}K_{1}(y,t)>0,\quad
K(r;1,y,t)=\frac{2}{1+r}K_{1}(y,t)>0,$

and hence $(x_{0},y_{0})\in(0,1)^{2}$. Choose $f_{0}=u_{n}$ that integrates to
$1$ where $u_{n}(y)\geq 0$ and $u_{n}(y)=0$ unless $\|y-y_{0}\|\leq 1/n$. Then
for large $n$, $u_{n}$ satisfies the required boundary conditions (vanishes in
a neighborhood of the endpoints) and we must have that

$f_{n}(x_{0},t):=\int_{0}^{1}K(r;x_{0},y,t)u_{n}(y)dy\geq 0,$

by Proposition 1. But it clearly holds by continuity of the integral kernel
that $\lim_{n\rightarrow\infty}f_{n}(x_{0},t)=K(r;x_{0},y_{0},t)<0$, a
contradiction. This proves part (2) of the theorem. ∎

##  Appendix B Proof of Theorem 5

###### Proof.

We begin with the proof of 1. Recall that

$\textstyle\mathrm{Var}_{f_{X}}[f(x,t)]=\frac{\mathbb{E}_{f_{X}}[K(x,Y,t)^{2}]}{n}-\frac{\mathbb{E}_{f_{X}}[K(x,Y,t)]^{2}}{n},$

where $K$ is the kernel given by (40). The second of these terms is bounded by
a constant multiple of $1/n$ so we consider the first. Recall the
decomposition (40):

$\begin{split}\textstyle K(r;x,y,t)&=\textstyle
K_{1}(x-y,t)\big{[}1+(x-y)\frac{1-r}{1+r}\big{]}+K_{1}(x+y,t)(x+y-1)\frac{1-r}{1+r}\\\
&\quad\quad\quad\quad+\frac{t(1-r)}{1+r}\big{[}K_{1}^{\prime}(x+y,t)+K_{1}^{\prime}(x-y,t)\big{]},\end{split}$

where $K_{1}$ is the standard periodic heat kernel. For $x,y\in[0,1]$ we have
that

$\displaystyle K_{1}(x-y,t)$ $\displaystyle\textstyle\sim_{t\downarrow
0}\frac{1}{\sqrt{2\pi
t}}\big{[}\exp\big{(}-\frac{(x-y)^{2}}{2t}\big{)}+\exp\big{(}-\frac{(x-y-1)^{2}}{2t}\big{)}+\exp\big{(}-\frac{(x-y+1)^{2}}{2t}\big{)}\big{]}$
$\displaystyle K_{1}(x+y,t)$ $\displaystyle\textstyle\sim_{t\downarrow
0}\frac{1}{\sqrt{2\pi
t}}\big{[}\exp\big{(}-\frac{(x+y)^{2}}{2t}\big{)}+\exp\big{(}-\frac{(x+y-1)^{2}}{2t}\big{)}+\exp\big{(}-\frac{(x+y-2)^{2}}{2t}\big{)}\big{]},$

with the rest of the expansion exponentially small as $t\downarrow 0$ and the
asymptotics valid upon taking derivatives. Using this, it is straightforward
to show that we can write

$\textstyle K(r;x,y,t)=\frac{1}{\sqrt{2\pi t}}G(r;x,y,t),$

where $G$ is bounded. From the above asymptotic expansions, we can write

$\displaystyle G(r;x,y,t)$
$\displaystyle=\textstyle\exp\Big{(}-\frac{(x-y)^{2}}{2t}\Big{)}+h_{1}(x,y,t)\exp\Big{(}-\frac{(x-y-1)^{2}}{2t}\Big{)}$
$\displaystyle\quad\textstyle+h_{2}(x,y,t)\exp\Big{(}-\frac{(x-y+1)^{2}}{2t}\Big{)}+h_{3}(x,y,t)\exp\Big{(}-\frac{(x+y)^{2}}{2t}\Big{)}$
$\displaystyle\quad\quad\textstyle+h_{4}(x,y,t)\exp\Big{(}-\frac{(x+y-2)^{2}}{2t}\Big{)}+E(x,y,t),$

where the $h_{i}$ are bounded and the error term $E(x,y,t)$ is exponentially
small as $t\downarrow 0$ uniformly in $x,y$. Furthermore, we have

$\textstyle\lim_{t\downarrow
0}\int_{0}^{1}\int_{0}^{1}f_{X}(y)K(r;x,y,t)h_{1}(x,y,t)\exp\Big{(}-\frac{(x-y-1)^{2}}{2t}\Big{)}dydx=0$

by the dominated convergence theorem (by considering the inner integral as a
function of $x$). Similar results hold for the other $h_{i}$ multiplied by
their relative Gaussian functions. Similarly, we have

$\textstyle\lim_{t\downarrow
0}\frac{1}{\sqrt{t}}\int_{0}^{1}\int_{0}^{1}\exp\Big{(}-\frac{(x-y)^{2}}{2t}\Big{)}f_{X}(y)h_{1}(x,y,t)\exp\Big{(}-\frac{(x-y-1)^{2}}{2t}\Big{)}dydx=0$

and likewise for the other $h_{i}$ multiplied by their relative Gaussian
functions. The integral

$\textstyle\frac{1}{\sqrt{\pi
t}}\int_{0}^{1}f_{X}(y)\exp\Big{(}-\frac{(x-y)^{2}}{t}\Big{)}dy$

is bounded and converges pointwise for almost all $x\in[0,1]$ to $f_{X}(x)$.
It follows that

$\textstyle\int_{0}^{1}\frac{\mathbb{E}_{f_{X}}[K(x,Y,t)^{2}]}{n}dx=\frac{1}{2n\pi
t}\int_{0}^{1}\int_{0}^{1}f_{X}(y)\exp\Big{(}-\frac{(x-y)^{2}}{t}\Big{)}dydx+\underline{\text{o}}(\frac{1}{n\sqrt{t}}).$
(43)

The rate (15) now follows.

We now prove 2 and 3. Define the function $p_{0}(x)$ via
$f_{X}(x)=f_{X}(0)+x(1-r)f_{X}(1)+p_{0}(x)$, then the proof of Proposition 2
showed that

$\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)=\int_{0}^{1}K(r;x,y,t)[p_{0}(y)-p_{0}(x)]dy.$
(44)

Define the function

$w(x,y)=p_{0}(y)-p_{0}(x)+(x-y)\frac{1-r}{1+r}(p_{0}(y)+p_{0}(1-y)),$ (45)

then (44) and (40) imply that $\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)$ is equal
to

$\textstyle\int_{0}^{1}K_{1}(x-y,t)w(x,y)dy+\frac{t(1-r)}{1+r}\int_{0}^{1}K_{1}(x-y,t)[p_{0}^{\prime}(y)-p_{0}^{\prime}(1-y)]dy,$
(46)

where we have integrated by parts for the last term and used
$p_{0}(0)=p_{0}(1)=0$. Define the function

$\textstyle F(x,t)=\int_{0}^{1}K_{1}(x-y,t)w(x,y)dy,$ (47)

then taking the partial derivative with respect to time, integrating by parts
and using $w(x,0)=w(x,1)=0$ we have

$\textstyle\frac{\partial F}{\partial
t}(x,t)=K_{1}(x,t)[p_{0}^{\prime}(1)-p_{0}^{\prime}(0)]\frac{rx-
x-r}{1+r}+\frac{1}{2}\int_{0}^{1}K_{1}(x-y,t)\frac{\partial^{2}w}{\partial
y^{2}}(x,y)dy.$ (48)

First we assume that $p_{0}^{\prime}(0)\neq p_{0}^{\prime}(1)$. In this case
the above shows that

$\textstyle\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)=\frac{rx-
x-r}{1+r}[p_{0}^{\prime}(1)-p_{0}^{\prime}(0)]\int_{0}^{t}K_{1}(x,s)ds+\mathcal{O}(t),$
(49)

where the $\mathcal{O}(t)$ is uniform in $x$. Using the above asymptotics for
$K_{1}(x,t)$ and the dominated convergence theorem, it follows that

$\begin{split}&\textstyle\int_{0}^{1}\big{\\{}\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)\big{\\}}^{2}dx\\\
&\sim_{t\downarrow
0}\textstyle\frac{[p_{0}^{\prime}(1)-p_{0}^{\prime}(0)]^{2}}{2\pi(1+r)^{2}}\int_{0}^{1}(rx-
x-r)^{2}\Big{[}\int_{0}^{t}\frac{\exp\big{(}-\frac{x^{2}}{2s}\big{)}}{\sqrt{s}}+\frac{\exp\big{(}-\frac{(x-1)^{2}}{2s}\big{)}}{\sqrt{s}}ds\Big{]}^{2}dx.\end{split}$
(50)

Let
$\tau(x)=\exp(-x^{2})-\sqrt{\pi}\left|x\right|\mathrm{erfc}(\left|x\right|)$,
then we can perform the integral in the square brackets in terms of $\tau$ to
yield

$\displaystyle\textstyle\int_{0}^{1}\big{\\{}\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)\big{\\}}^{2}dx$
(51) $\displaystyle\textstyle\sim_{t\downarrow
0}t\Big{\\{}\frac{2[p_{0}^{\prime}(1)-p_{0}^{\prime}(0)]^{2}}{\pi(1+r)^{2}}\int_{0}^{1}(rx-
x-r)^{2}\big{[}\tau(\frac{x}{\sqrt{2t}})+\tau(\frac{x-1}{\sqrt{2t}})\big{]}^{2}dx\Big{\\}}$
(52) $\displaystyle\textstyle\sim_{t\downarrow
0}t^{3/2}\Big{\\{}\frac{2[p_{0}^{\prime}(1)-p_{0}^{\prime}(0)]^{2}}{\pi(1+r)^{2}}(r^{2}+1)\sqrt{2}\int_{0}^{\infty}\tau(y)^{2}dy\Big{\\}}.$
(53)

To finish the proof in this case, we have that

$\textstyle\int_{0}^{\infty}\tau(y)^{2}dy=\frac{1}{3}(\sqrt{2}-1)\sqrt{\pi}.$

Next suppose that $p_{0}^{\prime}(0)=p_{0}^{\prime}(1)$. In this case we have

$\displaystyle\textstyle\frac{\partial F}{\partial t}(x,t)$
$\displaystyle=\textstyle\frac{1}{2}\int_{0}^{1}K_{1}(x-y,t)\frac{\partial^{2}w}{\partial
y^{2}}(x,y)dy$ (54)
$\displaystyle=\textstyle\frac{1}{2}\frac{\partial^{2}w}{\partial
y^{2}}(x,x)+U(x,t),$ (55)

for some bounded function $U(x,t)$ which converges to $0$ as $t\downarrow 0$
for almost all $x\in[0,1]$. It follows that

$\textstyle F(x,t)=t\frac{1}{2}\frac{\partial^{2}w}{\partial
y^{2}}(x,x)+t\tilde{F}(x,t),$ (56)

for some bounded function $\tilde{F}(x,t)$ which converges to $0$ as
$t\downarrow 0$ for almost all $x\in[0,1]$. Similarly, we have

$\textstyle\frac{t(1-r)}{1+r}\int_{0}^{1}K_{1}(x-y,t)[p_{0}^{\prime}(y)-p_{0}^{\prime}(1-y)]dy=\frac{t(1-r)}{1+r}[p_{0}^{\prime}(x)-p_{0}^{\prime}(1-x)]+tV(x,t),$
(57)

for some bounded function $V(x,t)$ which converges to $0$ as $t\downarrow 0$
for almost all $x\in[0,1]$. It follows from (46) that

$\textstyle\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)=t\Big{\\{}\frac{1}{2}\frac{\partial^{2}w}{\partial
y^{2}}(x,x)+\frac{(1-r)}{1+r}[p_{0}^{\prime}(x)-p_{0}^{\prime}(1-x)]\Big{\\}}+tW(x,t),$
(58)

for some bounded function $W(x,t)$ which converges to $0$ as $t\downarrow 0$
for almost all $x\in[0,1]$. But we have

$\textstyle\frac{\partial^{2}w}{\partial
y^{2}}(x,x)=f_{X}^{{}^{\prime\prime}}(x)-2\frac{1-r}{1+r}[p_{0}^{\prime}(x)-p_{0}^{\prime}(1-x)].$

The dominated convergence theorem then implies that

$\textstyle\int_{0}^{1}\\{\mathbb{E}_{f_{X}}[f(x,t)]-f_{X}(x)\\}^{2}dx\sim_{t\downarrow
0}t^{2}\int_{0}^{1}\frac{1}{4}[f_{X}^{{}^{\prime\prime}}(x)]^{2}dx.$ (59)

∎

##  Appendix C Four Corners Matrix and Proof of Theorem 7

The ‘Four Corners Matrix’ (26), is a non-symmetric example of a ‘tridiagonal
Toeplitz matrix with four perturbed corners’ [70, 64]. Although we do not
pursue it further, one can also show (by extending the techniques of [64])
that all functions of (26) are the sum of (i) a Toeplitz part, which can be
thought of as the solution without boundary conditions; and (ii) a Hankel
part, which is precisely the correction due to the boundary conditions. Exact
and explicit formulas for the eigenvalues and eigenvectors are available and
we will use these to prove Theorem 7

There is a unique zero eigenvalue, corresponding to the stationary density as
$t\rightarrow\infty$. The stationary density is an affine function in the
continuous PDE setting. In the discrete setting the components of the
stationary eigenvector $\bm{v}$ are equally-spaced, in the sense that $\forall
i,j,\;\;v_{i}-v_{i+1}=v_{j}-v_{j+1}=\textrm{constant}$. All non-zero
eigenvalues of $\mathbf{A}$ are positive and we are in the setting of [70,
Theorem 3.2 (i)]. In the case that $r\neq 1$, we can group the spectral data
into two classes with eigenvalues:

$\lambda_{k}=2-2\cos\theta_{k},\qquad k=1,\ldots,m,$

where

$\theta_{k}=\begin{cases}k\frac{2\pi}{m}&\mbox{ if }1\leq
k\leq\left\lfloor\frac{m-1}{2}\right\rfloor\\\
(k-\left\lfloor\frac{m-1}{2}\right\rfloor-1)\frac{2\pi}{m+1}&\mbox{ if
}\left\lfloor\frac{m-1}{2}\right\rfloor+1\leq k\leq m.\end{cases}$

The zero eigenvalue, when $k=\left\lfloor\frac{m-1}{2}\right\rfloor+1$, has
already been discussed. Other eigenvalues correspond to eigenvectors with
components (listed via subscripts)

$\begin{cases}v^{k}_{j}=r\sin((j-1)\theta_{k})-\sin(j\theta_{k})&\mbox{ if
}1\leq k\leq\left\lfloor\frac{m-1}{2}\right\rfloor\\\
w^{k-\lfloor\frac{m-1}{2}\rfloor-1}_{j}=\sin(j\theta_{k})&\mbox{ if
}\left\lfloor\frac{m-1}{2}\right\rfloor+2\leq k\leq m.\end{cases}$ (60)

Some properties of the discrete model and its links to the continuous model
are:

  * •

All eigenvalues of the Four Corners Matrix are purely real. Also, the
eigenvalues of the operator in the continuous model are likewise purely real.
This is perhaps surprising since the matrix is not symmetric, and the operator
is not self-adjoint.

  * •

The Four Corners Matrix $\mathbf{A}$ is diagonalizable. In contrast, the
operator for the continuous PDE is not diagonalizable, and instead, the
analogy of the Jordan Normal Form from linear algebra applies to the operator.
Despite this, the following still hold:

    1. 1.

The eigenvalues of the discrete model matrix $\mathbf{A}$ converge to that of
the continuous model (including algebraic multiplicities). This holds, for
example, in the Attouch–Wets topology - the convergence is locally uniform.

    2. 2.

The eigenvectors converge to the generalized eigenfunctions of the continuous
operator. Letting $j=\lfloor(m+1)x\rfloor$ we have ($k\neq 0$)

$\displaystyle\lim_{m\rightarrow\infty}w^{k}_{j}=\sin(2\pi kx)$
$\displaystyle\lim_{m\rightarrow\infty}\frac{m}{4\pi^{2}k^{2}}\big{[}(r-1)w^{k}_{j}-v^{k}_{j}\big{]}=\phi_{k}(x).$

We prove Theorem 7 by invoking the celebrated Lax Equivalence Theorem, which
states that ‘stability and consistency implies convergence’ [52]. We will take
consistency for granted. Typically when proofs in the literature use the Lax
Equivalence Theorem, it is also taken for granted that the PDE is well-posed.
Fortunately, we have already established that the PDE is indeed well-posed in
Theorem 1. It remains only to show stability. Even though the matrix
$\mathbf{A}$ has non-negative eigenvalues, this does not immediately imply
stability of the backward Euler method since $\mathbf{A}$ is not normal, i.e.
$\mathbf{A}$ does not commute with $\mathbf{A}^{*}$. We establish stability
for our problem by showing that bounds for the continuous model in Proposition
1 have corresponding bounds in the discrete model as follows, where we use a
subscript $l^{p}$ to denote the operator $l^{p}$ norm. In particular, Lemma 3
shows the discrete model is stable in the maximum norm. Convergence then
follows from the Lax Theorem.

######  Lemma 3 (Stability and Well-posedness of Discrete Approximation).

Let $m\geq 2$, then the backward Euler method (27) preserves probability
vectors and satisfies the bound

$\left\|\left(\mathbf{I}+\mathbf{A}\right)^{-K}\right\|_{l^{\infty}}\leq\max\Big{\\{}\frac{2r}{1+r},\frac{2}{1+r}\Big{\\}},\quad\forall
K\in\mathbb{N}.$

As a result of this lemma, we also gain stability in any $p$–norm via
interpolation.

###### Proof of Lemma 3.

Since the sums of each column $\mathbf{A}$ are zero, it follows that

$\sum_{j=1}^{m}u_{j}^{k+1}=\sum_{j=1}^{m}u_{j}^{k}.$

Hence to prove the first part, it suffices to show that $\bm{u}^{k+1}$ is non-
negative if $\bm{u}^{k}$ is. Suppose this were false, and let
$j\in\\{1,..,m\\}$ be such that $u^{k+1}_{j}<0$ is the smallest component of
$\bm{u}^{k+1}$. We have that

$u_{j}^{k}=(1+\mathbf{A}_{jj})u_{j}^{k+1}+\sum_{l\neq
j}\mathbf{A}_{j,l}u_{l}^{k+1}.$

By choice of $j$ and the fact that $\mathbf{A}_{jj}$ is positive, the off-
diagonals of $\mathbf{A}$ are negative and the sum of the $j$th column of
$\mathbf{A}$ is zero, it follows that

$\mathbf{A}_{jj}u_{j}^{k+1}+\sum_{l\neq j}\mathbf{A}_{j,l}u_{l}^{k+1}\leq 0.$

But this then implies that $u_{j}^{k}\leq u_{j}^{k+1}$, the required
contradiction.

To prove the second part, let $\bm{u}\in\mathbb{R}_{\geq 0}^{m}$ be any
initial vector with $\|\bm{u}\|_{\infty}\leq 1$ and let $\mathbbm{1}$ denote
the vector with $1$ in all entries. The eigenvector in the kernel of
$\mathbf{A}$ is a linear multiple of $w^{0}$ defined by

$\textstyle w^{0}_{j}=1+\frac{1-r}{1+rm}(j-1).$

Define the vector $x$ via
$x(1-r)/(1+r)=\mathbbm{1}-2(1+rm)/[(m+1)(r+1)]w^{0}$. This has components

$\textstyle x_{j}=\frac{m+1-2j}{m+1}.$

Extend this vector to have $x_{0}=0$, then an application of the discrete
Fourier transform implies that we can write for $j\neq 0$

$\textstyle x_{j}=\frac{1}{m+1}\sum_{k=1}^{m}G_{m}(k)\exp\left(\frac{2\pi
ikj}{m+1}\right),$

where

$\textstyle G_{m}(k)=\frac{\exp(4\pi ik/(m+1))-1}{(\exp(2\pi
ik/(m+1))-1)^{2}}=\frac{-i}{2}\frac{\sin(2\pi
k/(m+1))}{\sin^{2}(k\pi/(m+1))}.$

Hence we have that

$\displaystyle x_{j}$
$\displaystyle=\textstyle\frac{1}{m+1}\sum_{k=1}^{m-\lfloor\frac{m-1}{2}\rfloor-1}\left[G_{m}(k)\exp\left(\frac{2\pi
ikj}{m+1}\right)-\overline{G_{m}(k)}\exp\left(-\frac{2\pi
ikj}{m+1}\right)\right]$
$\displaystyle=\textstyle\frac{1}{m+1}\sum_{k=1}^{m-\lfloor\frac{m-1}{2}\rfloor-1}\frac{\sin(2\pi
k/(m+1))}{\sin^{2}(k\pi/(m+1))}\sin\left(\frac{2\pi kj}{m+1}\right).$

This implies that we can write $1$ as a linear combination of eigenvectors:

$\textstyle
1=\frac{2(1+rm)w^{0}_{j}}{(m+1)(1+r)}+\frac{1-r}{(m+1)(1+r)}\sum_{k=1}^{m-\lfloor\frac{m-1}{2}\rfloor-1}\frac{\sin(2\pi
k/(m+1))}{\sin^{2}(k\pi/(m+1))}w^{k}_{j}.$

Define $\bm{Q}^{K}=\left(\mathbf{I}+\mathbf{A}\right)^{-K}\mathbbm{1}$, and
$\bm{q}^{K}=\left(\mathbf{I}+\mathbf{A}\right)^{-K}\bm{u}$. In particular,
using the eigenvalue decomposition we have

$\displaystyle Q_{j}^{K}$
$\displaystyle=\textstyle\frac{2(1+rm)w^{0}_{j}}{(m+1)(1+r)}$
$\displaystyle+\frac{1-r}{(m+1)(1+r)}\textstyle\sum_{k=1}^{m-\lfloor\frac{m-1}{2}\rfloor-1}\frac{\sin(2\pi
k/(m+1))}{\sin^{2}(k\pi/(m+1))}\sin\left(\frac{2\pi
kj}{m+1}\right)\left(3-2\cos\left(\frac{2\pi k}{m+1}\right)\right)^{-K}.$

Using similar arguments to the first part of the proof, it is easy to prove
the Discrete Maximum Principle:

$\sup_{K\in\mathbb{N}\cup\\{0\\}}\left\|\bm{Q}^{K}\right\|_{\infty}=\max\left\\{\sup_{K\in\mathbb{N}\cup\\{0\\}}Q_{1}^{K},\sup_{K\in\mathbb{N}\cup\\{0\\}}Q_{m}^{K},1\right\\}.$

Explicitly, we have that

$\displaystyle Q_{1}^{K}$ $\displaystyle=\textstyle\frac{2(1+rm)}{(m+1)(1+r)}$
$\displaystyle+\frac{1-r}{(m+1)(1+r)}\textstyle\sum_{k=1}^{m-\lfloor\frac{m-1}{2}\rfloor-1}\frac{\sin^{2}(2\pi
k/(m+1))}{\sin^{2}(k\pi/(m+1))}\left(3-2\cos\left(\frac{2\pi
k}{m+1}\right)\right)^{-K}.$

This is monotonic in $K$ with limit $2(1+rm)/[(m+1)(1+r)]$. Similarly, we have

$\displaystyle Q_{m}^{K}$ $\displaystyle=\textstyle\frac{2(m+r)}{(m+1)(1+r)}$
$\displaystyle-\frac{1-r}{(m+1)(1+r)}\textstyle\sum_{k=1}^{n-\lfloor\frac{m-1}{2}\rfloor-1}\frac{\sin^{2}(2\pi
k/(m+1))}{\sin^{2}(k\pi/(m+1))}\left(3-2\cos\left(\frac{2\pi
k}{m+1}\right)\right)^{-K},$

which is monotonic in $K$ with limit $2(m+r)/[(m+1)(1+r)]$. Now, we must have
that each entry of $\bm{Q}^{K}\pm\bm{q}^{K}$ is non-negative since this is
true for $K=0$. It follows that

$\left\|\bm{q}^{K}\right\|_{\infty}\leq\left\|\bm{Q}^{K}\right\|_{\infty}=\max\left\\{\frac{2(1+rm)}{(m+1)(1+r)},\frac{2(m+r)}{(m+1)(1+r)}\right\\}.$

Since the $l^{\infty}$ operator norm of a real matrix is independent of
whether the underlying field is $\mathbb{R}$ or $\mathbb{C}$, the lemma now
follows by taking suprema over $m$. ∎

## References

  * [1] N. Agarwal and N. R. Aluru.  A data-driven stochastic collocation approach for uncertainty quantification in MEMS.  International Journal for Numerical Methods in Engineering, 83(5):575–597, 2010.
  * [2] A. H. Al-Mohy and N. J. Higham.  A new scaling and squaring algorithm for the matrix exponential.  SIAM Journal on Matrix Analysis and Applications, 31(3):970–989, 2010.
  * [3] G. D. Birkhoff.  Boundary value and expansion problems of ordinary linear differential equations.  Transactions of the American Mathematical Society, 9(4):373–395, 1908.
  * [4] C. Bolley and M. Crouzeix.  Conservation de la positivité lors de la discrétisation des problèmes d’évolution paraboliques.  RAIRO. Analyse numérique, 12(3):237–245, 1978.
  * [5] Z. I. Botev, J. F. Grotowski, and D. P. Kroese.  Kernel density estimation via diffusion.  The Annals of Statistics, 38(5):2916–2957, 2010.
  * [6] S. X. Chen.  Beta kernel estimators for density functions.  Computational Statistics & Data Analysis, 31(2):131–145, 1999\.
  * [7] E. Chevallier, E. Kalunga, and J. Angulo.  Kernel density estimation on spaces of Gaussian distributions and symmetric positive definite matrices.  SIAM Journal on Imaging Sciences, 10(1):191–215, 2017.
  * [8] E. A. Coddington and N. Levinson.  Theory of ordinary differential equations.  Tata McGraw-Hill Education, 1955.
  * [9] M. J. Colbrook.  Extending the unified transform: curvilinear polygons and variable coefficient PDEs.  IMA Journal of Numerical Analysis, 40(2):976–1004, 2020.
  * [10] M. J. Colbrook and L. J. Ayton.  A spectral collocation method for acoustic scattering by multiple elastic plates.  Journal of Sound and Vibration, 461:114904, 2019.
  * [11] M. J. Colbrook, L. J. Ayton, and A. S. Fokas.  The unified transform for mixed boundary condition problems in unbounded domains.  Proceedings of the Royal Society A, 475(2222):20180605, 2019.
  * [12] M. J. Colbrook, N. Flyer, and B. Fornberg.  On the Fokas method for the solution of elliptic problems in both convex and non-convex polygonal domains.  Journal of Computational Physics, 374:996–1016, 2018.
  * [13] M. J. Colbrook, A. S. Fokas, and P. Hashemzadeh.  A hybrid analytical-numerical technique for elliptic PDEs.  SIAM Journal on Scientific Computing, 41(2):A1066–A1090, 2019.
  * [14] J. Dai and S. Sperlich.  Simple and effective boundary correction for kernel densities and regression with an application to the world income and Engel curve estimation.  Computational Statistics & Data Analysis, 54(11):2487–2497, 2010\.
  * [15] F. P. J. de Barros, M. J. Colbrook, and A. S. Fokas.  A hybrid analytical-numerical method for solving advection-dispersion problems on a half-line.  International Journal of Heat and Mass Transfer, 139:482–491, 2019\.
  * [16] B. Deconinck, Q. Guo, E. Shlizerman, and V. Vasan.  Fokas’s uniform transform method for linear systems.  arXiv preprint arXiv:1705.00358, 2017.
  * [17] B. Deconinck, B. Pelloni, and N. E. Sheils.  Non-steady-state heat conduction in composite walls.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 470(2165):20130605, 2014.
  * [18] B. Deconinck, T. Trogdon, and V. Vasan.  The method of Fokas for solving linear partial differential equations.  SIAM Review, 56(1):159–186, 2014.
  * [19] D. Devroye, J. Beirlant, R. Cao, R. Fraiman, P. Hall, M. C. Jones, G. Lugosi, E. Mammen, J. S. Marron, C. Sánchez-Sellero, J. de Una, F. Udina, and L. Devroye.  Universal smoothing factor selection in density estimation: theory and practice.  Test, 6(2):223–320, 1997.
  * [20] L. Dümbgen, A. Hüsler, and K. Rufibach.  Active set and EM algorithms for log-concave densities based on complete and censored data.  arXiv preprint arXiv:0707.4643, 2007.
  * [21] L. Dümbgen and K. Rufibach.  logcondens: Computations related to univariate log-concave density estimation.  Journal of Statistical Software, 39, 2010.
  * [22] L. Dümbgen and K. Rufibach.  logcondens: Computations related to univariate log-concave density estimation.  Journal of Statistical Software, 39(6):1–28, 2011.
  * [23] L. Dümbgen, K. Rufibach, and D. Schuhmacher.  Maximum-likelihood estimation of a log-concave density based on censored data.  Electronic Journal of statistics, 8(1):1405–1437, 2014.
  * [24] L. C. Evans.  Partial differential equations.  American Mathematical Society, 2010.
  * [25] A. S. Fokas.  A unified transform method for solving linear and certain nonlinear PDEs.  In Proc. R. Soc. A, volume 453, pages 1411–1443. The Royal Society, 1997.
  * [26] A. S. Fokas.  Integrable nonlinear evolution equations on the half-line.  Communications in mathematical physics, 230(1):1–39, 2002.
  * [27] A. S. Fokas.  A unified approach to boundary value problems.  SIAM, 2008.
  * [28] A. S. Fokas and B. Pelloni.  A transform method for linear evolution PDEs on a finite interval.  IMA journal of applied mathematics, 70(4):564–587, 2005.
  * [29] G. Geenens.  Probit transformation for kernel density estimation on the unit interval.  Journal of the American Statistical Association, 109(505):346–358, 2014.
  * [30] D. Gilbarg and N. S. Trudinger.  Elliptic partial differential equations of second order.  Springer, 2015.
  * [31] N. J. Higham.  The scaling and squaring method for the matrix exponential revisited.  SIAM Journal on Matrix Analysis and Applications, 26(4):1179–1193, 2005.
  * [32] Y. Hu and C. Scarrott.  evmix: An R package for extreme value mixture modeling, threshold estimation and boundary corrected kernel density estimation.  Journal of Statistical Software, 84(5):1–27, 2018.
  * [33] M. C. Jones and D. A. Henderson.  Kernel-type density estimation on the unit interval.  Biometrika, 94(4):977–984, 2007.
  * [34] M. C. Jones, J. S. Marron, and S. J. Sheather.  Progress in data-based bandwidth selection for kernel density estimation.  Department of Statistics [University of North Carolina at Chapel Hill], 1992.
  * [35] M. C. Jones, J. S. Marron, and S. J. Sheather.  A brief survey of bandwidth selection for density estimation.  Journal of the american statistical association, 91(433):401–407, 1996.
  * [36] R. Kafri, J. Levy, M. B. Ginzberg, S. Oh, G. Lahav, and M. W. Kirschner.  Dynamics extracted from fixed cells reveal feedback linking cell growth to cell cycle.  Nature, 494(7438):480–483, 2013.
  * [37] R. J. Karunamuni and T. Alberts.  On boundary correction in kernel density estimation.  Statistical Methodology, 2(3):191–212, 2005.
  * [38] K. Kuritz, D. Stöhr, D. S. Maichl, N. Pollak, M. Rehm, and F. Allgöwer.  Reconstructing temporal and spatial dynamics from single-cell pseudotime using prior knowledge of real scale cell densities.  Nature Scientific Reports, 10(1):3619, 2020.
  * [39] K. Kuritz, D. Stöhr, N. Pollak, and F. Allgöwer.  On the relationship between cell cycle analysis with ergodic principles and age-structured cell population models.  Journal of Theoretical Biology, 414:91–102, 2017.
  * [40] J. Locker.  Spectral theory of non-self-adjoint two-point differential operators.  American Mathematical Soc., 2000.
  * [41] M. Machover.  A generalized eigenfunction expansion of the Green’s function.  Proceedings of the American Mathematical Society, 16(3):348–352, 1965.
  * [42] P. Malec and M. Schienle.  Nonparametric kernel density estimation near the boundary.  Computational Statistics & Data Analysis, 72:57–76, 2014.
  * [43] D. Mantzavinos and A. S. Fokas.  The unified method for the heat equation: I. non-separable boundary conditions and non-local constraints in one dimension.  European Journal of Applied Mathematics, 24(6):857–886, 2013.
  * [44] J. S. Marron and D. Ruppert.  Transformations to reduce boundary bias in kernel density estimation.  Journal of the Royal Statistical Society. Series B (Methodological), pages 653–671, 1994.
  * [45] R. Mehmood, G. Zhang, R. Bie, H. Dawood, and H. Ahmad.  Clustering by fast search and find of density peaks via heat diffusion.  Neurocomputing, 208:210–217, 2016.
  * [46] R. Mennicken and M. Möller.  Non-self-adjoint boundary eigenvalue problems, volume 192.  Elsevier, 2003.
  * [47] P. D. Miller and D. A. Smith.  The diffusion equation with nonlocal data.  Journal of Mathematical Analysis and Applications, 466(2):1119–1143, 2018.
  * [48] C. Moler and C. Van Loan.  Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later.  SIAM review, 45(1):3–49, 2003.
  * [49] M. A. Naimark.  Linear differential operators, harrap, london, 1967.  Trans ER Dawson from Russian original, 1952.
  * [50] P. Olver, N. E. Sheils, and D. Smith.  Revivals and fractalisation in the linear free space schrödinger equation.  Quarterly of Applied Mathematics, 2019.
  * [51] B. Pelloni and D. A. Smith.  Nonlocal and multipoint boundary value problems for linear evolution equations.  Studies in Applied Mathematics, 141(1):46–88, 2018.
  * [52] R. D. Richtmyer and K. W. Morton.  Difference methods for initial-value problems.  Malabar, Fla.: Krieger Publishing Co.,— c1994, 2nd ed., 1994.
  * [53] W. Saelens, R. Cannoodt, H. Todorov, and Y. Saeys.  A comparison of single-cell trajectory inference methods.  Nature Biotechnology, 37:547–554, 2019.
  * [54] R. J. Samworth.  Recent progress in log-concave density estimation.  arXiv preprint arXiv:1709.03154, 2017.
  * [55] D. Santhosh and V. V. Srinivas.  Bivariate frequency analysis of floods using a diffusion based kernel density estimator.  Water Resources Research, 49(12):8328–8343, 2013.
  * [56] O. Scaillet.  Density estimation using inverse and reciprocal inverse gaussian kernels.  Nonparametric statistics, 16(1-2):217–226, 2004.
  * [57] S. J. Sheather and M. C. Jones.  A reliable data-based bandwidth selection method for kernel density estimation.  Journal of the Royal Statistical Society: Series B (Methodological), 53(3):683–690, 1991.
  * [58] N. E. Sheils.  Interface Problems using the Fokas Method.  PhD thesis, 2015.
  * [59] N. E. Sheils and B. Deconinck.  Heat conduction on the ring: Interface problems with periodic boundary conditions.  Applied Mathematics Letters, 37:107–111, 2014.
  * [60] N. E. Sheils and B. Deconinck.  Interface problems for dispersive equations.  Studies in Applied Mathematics, 134(3):253–275, 2015.
  * [61] N. E. Sheils and B. Deconinck.  The time-dependent schrödinger equation with piecewise constant potentials.  European Journal of Applied Mathematics, 31(1):57–83, 2020.
  * [62] B. W. Silverman.  Density estimation for statistics and data analysis.  Routledge, 2018.
  * [63] J. S. Simonoff.  Smoothing methods in Statistics.  Springer Science & Business Media, 2012.
  * [64] G. Strang and S. MacNamara.  Functions of difference matrices are Toeplitz plus Hankel.  SIAM Review, 56(3):525–546, 2014.
  * [65] J. Tamarkin.  Some general problems of the theory of ordinary linear differential equations and expansion of an arbitrary function in series of fundamental functions.  Mathematische Zeitschrift, 27(1):1–54, 1928.
  * [66] T. Trogdon and B. Deconinck.  The solution of linear constant-coefficient evolution PDEs with periodic boundary conditions.  Applicable Analysis, 91(3):529–544, 2012.
  * [67] A. B. Tsybakov.  Introduction to nonparametric estimation.  Springer Science & Business Media, 2008.
  * [68] M. P. Wand and M. C. Jones.  Kernel smoothing.  Chapman and Hall/CRC, 1994.
  * [69] X. Xu, Z. Yan, and S. Xu.  Estimating wind speed probability distribution by diffusion-based kernel density method.  Electric Power Systems Research, 121:28–37, 2015.
  * [70] W.-C. Yueh and S. S. Cheng.  Explicit eigenvalues and inverses of tridiagonal Toeplitz matrices with four perturbed corners.  ANZIAM Journal, 49(03):361–387, 2008.