arXiv:1709.02695v1 [] 8 Sep 2017

Applications of an algorithm for solving
Fredholm equations of the first kind∗
Minwoo Chae
Department of Mathematics, Applied Mathematics and Statistics
Case Western Reserve University
Ryan Martin
Department of Statistics, North Carolina State University
and
Stephen G. Walker
Department of Mathematics, University of Texas at Austin
September 11, 2017

Abstract
In this paper we use an iterative algorithm for solving Fredholm equations of the first kind.
The basic algorithm is known and is based on an EM algorithm when involved functions are
non–negative and integrable. With this algorithm we demonstrate two examples involving the
estimation of a mixing density and a first passage time density function involving Brownian
motion. We also develop the basic algorithm to include functions which are not necessarily
non–negative and again present illustrations under this scenario. A self contained proof of
convergence of all the algorithms employed is presented.

Keywords: Brownian motion first passage time; convergence; expectation–maximization; iterative
algorithm; mixture model.

1

Introduction

An important problem in statistics and applied mathematics is the solution to a so-called Fredholm
equation of the first kind, i.e., given a probability density function f (x) on X ⊂ R and a nonnegative kernel k(x, θ) which is a probability density function in x for each θ, find the probability
∗

Martin and Walker acknowledge NSF support with grants DMS-1611791 and DMS-1612891, respectively.

1

density function p(θ) such that
Z
f (x) =

k(x, θ) p(θ) dθ.

(1)

Θ

Such equations have a wide range of applications across a variety of fields, including signal processing, physics, and statistics. See Ramm (2005) for a general introduction, and Vangel (1992)
for statistical applications. We will also consider problems where the density and non-negativity
constraints are relaxed.
Given the importance of this problem, it should be no surprise that there is a substantial
literature on the theoretical and computational aspects of solving (1); see, for example, Corduneanu (1994); Groetsch (2007); Hansen (1999); Morozov (1984); Wing (1990). Existence of a
unique solution to (1) is a relevant question, though we will not speak directly on this issue here.
One important case for which there is a well-developed existence theory is when the operator
R
p 7→ k(·, θ) p(θ) dθ is compact, a consequence of square-integrability of k(x, θ) with respect to
dx × dθ. In this case, there exists a singular value system (σj , uj , vj ) for which (uj (x)) and (vj (θ))
are orthogonal bases in L2 (X) and L2 (Θ), respectively, and
Z
k(x, θ) vj (θ) dθ = σj uj (x).
Θ

Then f can be expressed as f (x) =

P

j

λj uj (x) for some (λj ), and it follows that

p(θ) =

X

(λj /σj ) vj (θ),

j

which exists and is square-integrable if

2
j (λj /σj )

P

< ∞. Of course, the above formula for p is

not of direct practical value since computing all the eigenvalues and eigenfunctions is not feasible.
Computationally efficient approximations are required.
The classic iterative algorithm for finding p is given by
Z

pm (θ) = pm−1 (θ) + k(x, θ) f (x) − fm−1 (x) dx,

(2)

X

where
Z
fm (x) =

k(x, θ) pm (θ) dθ.
Θ

Convergence properties of (2) are detailed in Landweber (1951). An issue here is that, under
the constraint that p and k(·, θ) are density functions, there is no guarantee that the sequence of
approximations (pm ) coming from (2) are density functions.
2

An alternative is to apply a discretization–based method. That is, specify grid points (θj , xi )
and approximate the original problem (1) via a discrete system
f (xj ) =

X

wi p(θi ) k(xj , θi ),

i

where the (wi ) are weighting coefficients for a quadrature formula. One then solves the linear
system of equations to get an approximation for p. See for example, Phillips (1962). Of course,
there is no guarantee here either that the solution p will be a density.
As an alternative to the additive updates in (2), in this paper we focus on a multiplicative
version, described in Section 2, that can guarantee the sequence of approximations are density
functions. Moreover, the theoretical convergence analysis of this multiplicative algorithm turns
out to be rather straightforward for the case where there exists a unique solution to equation (1).
For the general case, where a solution may not exist, the asymptotic behavior of the algorithm
can still be determined, and, in Section 3, we provide a characterization of this algorithm as
an alternating minimization algorithm, and employ the powerful tools developed in Csiszzár and
Tusnády (1984) to study its convergence. The remainder of the paper focuses on three applications
of this algorithm. The first, in Section 4, is estimating a smooth mixing density based on samples
X1 , . . . , Xn from the density f in (1). The second, in Section 5, is computing the density of the
first passage time for Brownian motion hitting an upper boundary, where the boundary need not
be concave. The third, in Section 6, is solving general Fredholm equations where the density
function constraints are relaxed. Finally, some concluding remarks are given in Section 7.

2

The algorithm

Following Vardi and Lee (1993), an alternative to the algorithms described above for iteratively
solving (1) is to fix an initial guess p0 and then repeat
Z
k(x, θ) f (x)
dx,
pm (θ) = pm−1 (θ)
fm−1 (x)
X

m ≥ 1.

(3)

This algorithm was also studied in Shyamalkumr (1996). First note that, if f0 is well-defined, then
it follows immediately from the multiplicative structure of the algorithm, and Fubini’s theorem,
that pm is also a density for m ≥ 1. This overcomes the difficulty faced with using (2). Second,

3

for a comparison with (2), we see that (3) operates in roughly the same way, but on the log–scale:
Z
k(x, θ) f (x)
dx.
(4)
log pm (θ) = log pm−1 (θ) + log
fm−1 (x)
X
This logarithmic version will be helpful when we discuss convergence in Section 3.
The sequence (3) also has a Bayesian interpretation. Indeed, if, at iteration m, we had seen
an observation X from the distribution with density f , then the Bayes update of “prior” pm−1 (·)
to “posterior” pm (· | X) would be given by
pm (θ | X) = pm−1 (θ)

k(X, θ)
.
fm−1 (X)

But without such an X, yet knowing X comes from f , the natural choice now is to use the “average
update”
Z
pm (θ | x) f (x) dx,

pm (θ) =
X

which is exactly (3).
For a quick proof–of–concept, consider a Pareto density f (x) = a(x + 1)−(a+1) , supported
on (0, ∞), with a > 0. It is easy to check that f is a gamma mixture of exponentials, i.e.,
R
f (x) = k(x, θ)p(θ) dθ, where k(x, θ) = θ e−θx is an exponential density and p(θ) = Γ(a)−1 θa−1 e−θ
is a gamma density. We can apply algorithm (3) to solve the above equation for p. Figure 1 shows
a plot of the estimate from (3), with a = 5, based on 200 iterations and p0 a half-Cauchy starting
density, along with the true gamma density p(x). Clearly, the approximation is quite accurate
over most of the support.
Before formally addressing the convergence properties of (3), it will help to provide some
intuition as to why it should work. The argument presented in Vardi and Lee (1993) proceeds by
considering i.i.d. samples X1 , . . . , Xn from the distribution with density f . Replacing the f (x) dx
in (3) with dFb(x), where Fb is the empirical distribution based on X1 , . . . , Xn , gives the algorithm
Z
pbm (θ) = pbm−1 (θ)

X

k(x, θ) b
dF (x).
fm−1 (x)

(5)

It turns out that this is precisely an EM algorithm to compute the nonparametric maximum
likelihood estimator of p, see e.g., Laird (1978). This connection to likelihood–based estimation
gives the algorithm (5) some justification for a fixed sample X1 , . . . , Xn . The dependence on a
particular sample can be removed by allowing n → ∞ and applying the law of large numbers to
4

0.20
0.15
0.10
0.00

0.05

p(θ)

0

5

10

15

20

θ

Figure 1: Plot of the mixing density p200 (black) from (3) and the true gamma density (gray).

get convergence of (5) to (3), hence motivation for the latter. However, despite the identification
of (5) as an EM algorithm, no formal convergence proof has been given that the iterates pbm in (5)
converge to the nonparametric maximum likelihood estimator, but see Chae et al. (2017).
The above argument giving intuitive support for algorithm (3) is not fully satisfactory because
it does not give any indication that the algorithm will converge to a solution of (1). For a more
satisfactory argument, note that, if the algorithm converges to a limit p∞ , then we must have
Z
k(x, θ)
f (x) dx = 1 ∀ θ ∈ supp(p∞ ),
X f∞ (x)
R
where f∞ (x) = k(x, θ) p∞ (θ) dθ. According to Lemma 2.1 in Patilea (2001) or Lemma 2.3 in
Kleijn and van der Vaart (2006), the above condition implies that
D(f, f∞ ) = inf D(f, fP ),
P

where D(f, g) =

(6)

R

log(f /g) f dx is the Kullback–Leibler divergence and the infimum is over all
R
densities of the form k(x, θ) dP (θ) for some probability measure P on Θ. So, if there exists a
solution to (1), then the limit p∞ would have to be one of them. Even if there is no solution to
(1), p∞ will be such that the corresponding f∞ is closest to f in the Kullback–Leibler sense. To
make this argument fully rigorous, we need to establish that algorithm (3) does indeed converge.

5

3

Convergence properties

In Section 3.4 of Vardi and Lee (1993), the authors consider (1), but with some restrictions. The
first is that p is a finite measure, i.e., point masses on a finite number of atoms, and the second is
that the solution p is piecewise constant. Our arguments here do not require such restrictions.
To assess the properties of (2), let p be a solution to (1). Multiply (4) by p(θ) throughout,
and then integrate over θ to get
Z
Z
Z
Z
k(x, θ) f (x)
dx dθ.
p(θ) log pm (θ) dθ = p(θ) log pm−1 (θ) dθ + p(θ) log
fm−1 (x)
By Jensen’s inequality, the last term is lower bounded by D(f, fm−1 ), the Kullback–Leibler divergence of fm−1 from f . Since D is non–negative, we can deduce that
D(p, pm ) ≤ D(p, pm−1 ) − D(f, fm−1 ).
This implies that D(p, pm ) is a non–negative and non–increasing sequence, hence has a limit, say,
c ≥ 0, which, in turn, implies fm → f strongly in the D or L1 sense. Therefore, in agreement with
Landweber (1951) for the additive algorithm (2), we have that pm is converging to a set S where
D(p, pe) = c for all pe ∈ S. Of course, if p is a unique solution to (1) then pm → p strongly. A similar
argument for convergence of (3) is presented in Shyamalkumr (1996), with some generalizations.
On the other hand, if (1) does not have a solution, then, as discussed above, we expect that
pm in algorithm (3) will converge to a density with limit p∞ such that the corresponding mixture
f∞ satisfies (6). This can be proved by considering (3) as an alternating minimization procedure;
see, for example, Csiszár (1975); Csiszzár and Tusnády (1984); Dykstra (1985).
In the following, define the joint densities on X × Θ,
f (x)qm−1 (x, θ)
.
qm (x, θ) = pm (θ)k(x, θ) and πm (x, θ) = R
qm−1 (x, θ0 )dθ0

(7)

Theorem 3.1. For an initial solution p0 > 0, let (pm ) be obtained via (3). Assume there exists a
sequence (p∗s )s≥1 of densities such that
D(f, fs∗ ) → inf D(f, fP )
P

and D(πs∗ , qm ) < ∞ for some m ≥ 0, where fs∗ (x) =
πs∗ (x, θ)

R

k(x, θ)p∗s (θ)dθ and

f (x)p∗s (θ)k(x, θ)
R
.
=
p∗s (θ0 )k(x, θ0 )dθ0

Then, D(f, fm ) decreases to inf P D(f, fP ).
6

Proof. Let P be the set of all bivariate densities π on X × Θ with x-marginal f , i.e., such that
R
π(x, θ) dθ = f (x). Similarly, let Q be the set of all bivariate densities q on X × Θ such that
q(x, θ) = k(x, θ) p(θ) for some density p. Note that qm and πm in (7) satisfy qm ∈ Q and πm ∈ P.
We first claim that (qm ) and (πm ) are obtained by the alternating minimization procedure of
Csiszzár and Tusnády (1984) with the objective function D(π, q), where π and q ranges over P
and Q, respectively. To see this, note that D(πm , qm−1 ) = D(f, fm−1 ) ≤ D(π, qm−1 ) for every
π ∈ P; the inequality holding since f and fm−1 are marginal densities of π and qm−1 , respectively.
It follows that πm = argminπ∈P D(π, qm−1 ). Also, note that
ZZ
πm (x, θ)
dθ dx
D(πm , q) =
πm (x, θ) log
q(x, θ)
Z
Z
= C − log p(θ) πm (x, θ) dx dθ
Z
Z
f (x)qm−1 (x, θ)
= C − log p(θ) R
dx dθ
qm−1 (x, θ0 )dθ0
Z
Z
k(x, θ)pm−1 (θ)
f (x) dx dθ
= C − log p(θ) R
k(x, θ0 )pm−1 (θ0 ) dθ0
Z
= C − log p(θ) dPm (θ),
where C does not depend on p. The last integral is maximized at pm , so qm = argminq∈Q D(πm , q).
Since P and Q are convex, we have
D(f, fm−1 ) = D(πm , qm−1 ) &

inf

p∈P0 , q∈Q

D(π, q),

by Theorem 3 in Csiszzár and Tusnády (1984), with P0 = {π ∈ P : D(π, qm ) < ∞ for some m}.
We have D(π, q) = D(f, fp ), and and the assumption D(πs∗ , qm ) < ∞ for some m implies
that inf p∈P0 inf q∈Q D(π, q) ≤ inf P D(f, fP ). Since D(f, fm ) ≥ inf P D(f, fP ), we conclude that
D(f, fm ) & inf P D(f, fP ).
Two remarks are in order. First, the integrability condition D(πs∗ , qm ) < ∞ for some m ≥ 0 is
not especially strong. For example, assume there exists P ∗ which minimizes D(f, fP ) and has a
density p∗ . Such a minimizer exists under a mild identifiability condition; see Lemma 3.1 of Kleijn
and van der Vaart (2006). Then D(π ∗ , qm ) < ∞ for some m implies the required integrability
condition, where
f (x)p∗ (θ)k(x, θ)
π ∗ (x, θ) = R ∗ 0
.
p (θ )k(x, θ0 )dθ0
7

Second, if Θ is compact and the minimizer P ∗ is unique, every subsequence of (pm ) has a further
subsequence weakly converging to P ∗ , so (pm ) converges to P ∗ .
It is also worth noting that the monotonicity property of Dm = D(f, fm ) from the theorem
can be used to define a stopping rule. For example, one could terminate algorithm (3) when Dm
itself and/or the difference Dm − Dm−1 falls below a certain user–specified tolerance.

4

Smooth mixing density estimation

Suppose we observe X1 , . . . , Xn as i.i.d. from a distribution with density f as in (1), and we wish
to estimate a smooth version of p. The nonparametric maximum likelihood estimator of p is known
to be discrete, so is not fully satisfactory for our purposes. One idea for a smooth estimate of p
is to maximize a penalized likelihood. The proposal in Liu et al. (2009) is to define a penalized
log–likelihood function
−1

`(η) = n

n
X

Z
log

k(Xi , θ)e

η(θ)

Z
dθ − log

η(θ)

e

Z
dθ − λ

[η 00 (θ)]2 dθ,

i=1

where λ > 0 is a smoothing parameter, and η determines the mixing density
p(θ) = R

eη(θ)
.
eη(θ) dθ

The right-most integral in the expression for `(η) measures the curvature of η, so maximizing `
will encourage solutions which are “less curved,” i.e., more smooth. An EM algorithm is available
to produce an estimate of η and, hence, of p. However, each iteration requires solving a non–trivial
functional differential equation.
Here we propose a more direct approach, namely, to use algorithm (3) with the true density f
replaced by, say, a kernel density estimate of f . For some h > 0, define the kernel estimate
n

1 X
f (x) =
φ
nh i=1
h



x − Xi
h



√
where φ(x) = exp(− 12 x2 )/ 2π is the standard normal density function. The bandwidth h can be
selected in a variety of ways. One option is the traditional approach (Sheather and Jones (1991))
of minimizing asymptotic mean integrated square error, which gives h of order n−1/5 . Another

8

0.15

0.20

f(x)

0.10

0.15

0.05

0.10

0.00

0.00

0.05

p(θ)

5

10

15

20

25

30

10

35

15

20

25

30

35

x

θ

(a) Mixing density

(b) Data and mixtures

Figure 2: Mixing and mixture density estimates based on algorithm (3) with a kernel density plugin for the galaxy data. In Panel (b), solid line is the kernel estimate and dashed is the mixture
corresponding to the mixing density in Panel (b).

idea is to get the nonparametric maximum likelihood estimator Pb of P , define the corresponding
mixture
Z
fb(x) =

k(x, θ) dPb(θ)

and then choose the bandwidth b
h = arg min d(f h , fb), where d is, say, the L1 distance. Finally, our
proposal is to estimate p by running algorithm (3) with f replaced by f ĥ . Since f ĥ is smooth, so
too will be our estimate of p.
As a quick real–data illustration, consider the well–known galaxy data set, available in the
MASS package in R (R Core Team (2015)). We estimate the mixture density f with the Gaussian
kernel method described above, using the default settings in density. Then we estimate the
mixing density via the procedure just described above. In this case, we used 25 iterations of (3)
and the results are displayed in Figure 2. Panel (a) shows p25 from (3), and Panel (b) shows the
data histogram, the kernel density estimator (solid), and the mixture corresponding to p25 . Both
densities in Panel (b) fit the data well, and the fact that the two are virtually indistinguishable
suggests that the density p in Panel (a) is indeed a solution to the mixture equation.
Next, for further investigation into the performance of the proposed mixing density estimation
procedure, we present some examples with simulated data. We start by considering a deconvolution
problem, where k(·, θ) is the normal density with mean θ and standard deviation σ = 0.05. Two
9

mixing densities on the unit interval [0,1] are illustrated;
p1 (θ) ∝ θ4 (1 − θ)4




θ − 0.3
θ − 0.7
p2 (θ) ∝ φ
+ 2φ
.
0.1
0.1
With a sample of size n = 300, we first obtained a kernel density estimator f h , using the Gaussian
kernel, with bandwidth as in Sheather and Jones (1991), and then ran (3). As suggested above,
monotonicity of Dm = D(f h , fm ) suggests a stopping rule, and we terminated the algorithm when
Dm − Dm−1 < 10−5 . The number of iterations used were m = 8 for p1 and m = 14 for p2 . The
estimates of the mixing and corresponding mixture densities are depicted in Figure 3. As can be
seen in the right figures, the fm and fb are indistinguishable; that is, Dm are effectively zero after
a few iterations.
Next, we consider scale mixtures of normal distributions where k(·, θ) is the centered normal
density with variance θ. For the mixing density, two well–known distributions on the positive real
line, the inverse-gamma and gamma densities, are considered, i.e.,
p1 (θ) ∝ θ−3 e−1/θ

and p2 (θ) ∝ e−5θ .

With a sample of size n = 300, the estimator of the mixing density is obtained as previously.
The number of iterations are m = 57 for p1 and m = 78 for p2 , and the estimates of p and f are
illustrated in Figure 4. Although the density estimates, f h and fm , are close to the true f , there
are deviations between p and pm , in particular for the inverse–gamma case. This is mainly due to
the ill–posedness of the problem.

5

First passage time for Brownian motion

The second example is computing the stopping time density for Brownian motion hitting a boundary function. For example, since Breiman (1966) there has been substantial interest in the first
passage time distribution of Brownian motion passing a square root boundary. To set the scene,
denote (B(t))t≥0 as standard Brownian motion started at 0, let h(t) : (0, ∞) → R be a continuous
function with h(0) > 0, the boundary function, and define
τ = inf{t > 0 : B(t) ≥ h(t)}.
10

3.0
2.0

f(x)

0.0

1.0

3.0
2.0
1.0
0.0

p(θ)

0.0

0.4

0.8

0.0

0.4

θ

3.0
2.0
1.0
0.0

1.0

f(x)

2.0

3.0

x

0.0

p(θ)

0.8

0.0

0.4

0.8

0.0

θ

0.4

0.8
x

Figure 3: Estimate of mixing density with location mixtures of normals. Left plots are the true
(black) and estimated (red) mixing density. Right plots are the true (black), kernel estimator
(red) and the one obtained by smooth NPMLE (green).

11

0.0

0.4

f(x)

0.8

2.0
1.0
0.0

p(θ)

0

1

2

3

4

−4 −2

θ

0

2

4

2

4

0.0

0.4

f(x)

1.0
0.0

p(θ)

0.8

2.0

x

0

1

2

3

4

−4 −2

θ

0
x

Figure 4: Estimate of mixing density with scale mixtures of normals. Left plots are the true
(black) and estimated (red) mixing density. Right plots are the true (black), kernel estimator
(red) and the one obtained by smooth NPMLE (green).

12

Then Peskir et al. (2002) provides an equation satisfied by the density p of τ , in particular,



Z t 
h(t) − h(s)
h(t)
√
Φ
dP (s),
= F (t) −
1−Φ √
t−s
t
0
where Φ and P are the probability measures/cumulative distribution functions corresponding to
φ and p, respectively. This leads to a Volterra equation and Peskir et al. (2002) demonstrates
√
solutions when h is constant and linear. When h is increasing and concave and h(t) ≤ h(0) + b t
then Peskir et al. (2002) proves that
Z

t

R(t, s) m(s) ds

P (t) = m(t) +
0

where

m(t) = 2 − 2Φ

h(t)
√
t


and R(t, s) =

∞
X

Km (t, s),

m=1

with
Z

t

K1 (t, r) Km−1 (r, s) dr,

Km (t, s) =
s

and



h(t) − h(s)
1
h(t) − h(s)
0
√
φ
.
K1 (t, s) = √
2h (s) −
t−s
t−s
t−s
√
Clearly h(t) = a + b t meets the requirements for this result.
On the other hand, we can solve using (3) after setting up a suitable Fredholm equation for
√
h(t) increasing and bounded by a + b t, but not necessarily concave. Define, for any x > 0, the
martingale

Xt = exp x B(t) − 12 x2 t .
With θ = Ta,b = inf{t > 0 : B(t) ≥ a + b h(t)} with a, b > 0, h(0) = 0 and h(t) ≤

√

t, we have that

Xt∧Ta,b is a martingale and, from the optional stopping theorem, that E(Xt∧Ta,b ) = E(X0 ) = 1.
Now it is easy to show that Xt∧Ta,b ≤ exp{xa + 21 b2 h2 (t)/t}, which is bounded and, hence, from
the dominated convergence theorem,
E(XTa,b ) = lim E(Xt∧Ta,b ) = E(X0 ) = 1.
t→∞

Effectively, we are also using the result that Ta,b < ∞ almost surely. Hence,
Z

exp xa + xbh(θ) − 12 x2 θ p(θ) dθ = 1
13

which after some algebra yields the Fredholm equation
Z ∞
k(x, θ) pe(θ) dθ = a e−xa = f (x),

(8)

0

where k(x, θ) is the normal density with mean b h(θ)/θ and variance 1/θ, constrained on R+ , and
√
√
√
1 2 2
a 2π e 2 b h (θ)/θ Ψ(−bh(θ)/ θ) p(θ) = θ pe(θ)
and Ψ(·) = 1 − Φ(·). Now (8) can be readily solved by algorithm (3).
The approach here is closely related to work appearing in Valov (2009). This considers first
passage times which are almost surely finite and of the type τ = inf{t > 0 : B(t) ≤ bα (t) =
b(t) + α t}, where bα (t) > −c for some c > 0. While Valov (2009) constructs a Fredholm equation
using Girsanov’s theorem, as we do, this is quite different and each equation depends on the type
of boundary, whereas ours is always based on the half–normal and exponential distributions.
Here we demonstrate the algorithm for solving the Fredholm equation for the hitting time
√
density with a square root boundary; specifically a + b t with a = 1 and b = 0.1. We start
with a grid of points θj = hj, where h = 0.05 and j = 1, . . . , 1000. The starting density is
p0 (θ) = 0.01 e−0.01θ and at each iteration we evaluate pm (θj ). Due to the ease of sampling from an
exponential density, we evaluate
Z

∞

0

k(x, θj ) f (x)
dx
fm (x)

using Monte Carlo methods, i.e.,
N

−1

N
X
k(Xi , θj )
i=1

fm (Xi )

where the (Xi ) are iid from the exponential distribution with mean 1/a, and we take N = 5000.
The fm (Xi ) is evaluated using the trapezoidal rule using the pem (θj ) values, so
M

fm (Xi ) ≈

hX
[k(Xi , θj−1 ) pem (θj−1 ) + k(Xi , θj ) pem (θj )] ,
2 j=1

where θ0 = 0. Consequently, we have numerically,
pem+1 (θj ) = N

−1

N
X
i=1

h

PM

j=1

2 k(Xi , θj ) pem (θj )
[k(Xi , θj−1 ) pem (θj−1 ) + k(Xi , θj ) pem (θj )]

.

The estimated p(θ), obtained by transforming pe(θ), is given in Figure 5, based on 200 iterations
of algorithm (3).
14

0.5
0.4
0.3
density
0.2
0.1
0.0
0

10

20

30

40

50

theta

Figure 5: Estimate of hitting time density for square root boundary

6

General Fredholm equation

Algorithm (3) can be applied to solve general Fredholm equations of the first kind, where f, p and
k(·, θ) are not necessarily probability density functions. Assume first that f, p and k(·, θ) in (1)
are non-negative functions, but not necessarily probability densities. Then equation (1) can be
rewritten as
Z
f (x) =

e
k(x, θ) q(θ) dθ,
Θ

R
R
where e
k(x, θ) = k(x, θ)/ k(x, θ)dx and q(θ) = p(θ) k(x, θ)dx. Thus, the update (3) can be
replaced with
pm−1 (θ)
pm (θ) = R
k(x, θ) dx
X

Z
X

k(x, θ) f (x)
dx,
fm−1 (x)

which was also considered in Vardi and Lee (1993). Here it is assumed that supx
R
supθ k(x, θ) dx are both finite.

(9)
R

k(x, θ)dθ and

Next, assume that the solution p and f may not be necessarily non-negative functions. In this
case, instead of the original equation (1), we solve an equivalent non-negative equation
Z
e
f (x) =
k(x, θ) pe(θ) dθ,

(10)

Θ

where pe(θ) = p(θ) + t, fe(x) = f (x) + t

R
Θ

k(x, θ) dθ, and t > 0 is a constant to be specified. Since
15

k is non-negative, so are both fe and pe for a sufficiently large t. As illustrated below, the value of t
rarely affects the convergence rate in practice. Therefore, t can be chosen as a very large number.
Finally, assume that k is not necessarily non-negative, and write k = k + −k − , where k + , k − ≥ 0.
This case can also be solved by transforming the original equation to a non-negative one. For
convenience, assume that Θ = [0, 1], then the original equation (1) can be written as
Z 1
Z 1
+
k (x, θ)p(θ)dθ −
k − (x, θ)p(θ)dθ.
f (x) =
0

Note that
Z

0

1

2

Z

−

k − (x, θ − 1)p(θ − 1)dθ = 0,

k (x, θ)p(θ)dθ −
0

(11)

1

and therefore, by adding the last two display equations, we have
Z 1
Z 2
+
f (x) =
k (x, θ)p(θ)dθ +
k − (x, θ − 1){−p(θ − 1)} dθ.
0

1

Let

 k + (x, θ)
if θ ∈ [0, 1]
e
k(x, θ) =
 k − (x, θ − 1) if θ ∈ (1, 2]
and

 p(θ)
if θ ∈ [0, 1]
pe(θ) =
 −p(θ − 1) if θ ∈ (1, 2],
then we have the Fredholm equation
Z
f (x) =

2

e
k(x, θ)e
p(θ)dθ

(12)

0

with a non-negative kernel, which can be solved as previously described.
Note that e
k and pe may have discontinuities; but this does not cause any problem once (12)
has the unique solution. If there is another solution to (12), say p, and (11) is not satisfied, the
restriction of p on [0,1] may not be a solution of the original equation. In this case we may apply
another decomposition of k. In many examples, however, the simple approach (12) works well.
For illustration, consider the equation (1), where k(·, θ) is the normal density with mean θ and
standard deviation σ = 0.05, and p has both positive and negative components. In particular, we
consider two examples
p1 (θ) = b2,5 (θ) − b4,1 (θ)
p2 (θ) = b10,1 (θ) − b1,10 (θ),
16

where ba,b (·) is the density of Beta(a, b) distribution. In both cases, a transformed equation (10)
is solved with t = 50 and the results are illustrated in Figure 6. It can be easily seen that the
true solution and pm are nearly the same. Different values of t (5 × 10k with k ≤ 4) have been
tried, and in any case, the algorithm has been stopped in 10 iterations yielding nearly the same

−3

−1

f(x)

0
−4

p(θ)

2

1 2

solution.

0.0

0.4

0.8

−0.2

0.2

θ

1.0

0 2 4
−4

0

f(x)

5

10

x

−10

p(θ)

0.6

0.0

0.4

0.8

−0.2

θ

0.2

0.6

1.0

x

Figure 6: Illustration with a positive kernel k(x, θ) = φσ (x − θ) and general solutions p1 (·) (top)
and p2 (·) (bottom).

17

1.0
−1.0

0.0

f(x)

1.0
0.0
−1.0

p(θ)

0.0

0.4

0.8

−1.5

−0.5

θ

1.5

0.5

1.5

−3

−1

1

f(x)

1

2

3

3

x

−1 0

p(θ)

0.5

0.0

0.4

0.8

−1.5

θ

−0.5
x

Figure 7: Illustration with a general kernel k(x, θ) = φσ (x − θ) − φσ (x + θ) for p1 (·) (top) and p2 (·)
(bottom).

18

Next, we consider a general kernel k(x, θ) = φσ (x − θ) − φσ (x + θ) with σ = 0.05 and
p1 (θ) = b2,3 (θ) − b3,2 (θ)
p2 (θ) = b2,7 (θ) + b3,4 (θ) − 1,
where φσ (x) = φ(x/σ)/σ. The equations are solved by setting k + (x, θ) = φσ (x − θ), k − (x, θ) =
φσ (x + θ) and t = 50. For both examples, the algorithm stopped in 5 iterations. Results are
illustrated in Figure 7.

7

Conclusion

In this paper, we focused on an algorithm for solving Fredholm integral equations of the first kind,
its properties, and some applications. For the mixing density estimation application described in
Section 4, we did not address the question of whether the estimate based on plugging in a kernel
density estimator for f in (3) would be consistent in the statistical estimation sense. The predictive
recursion method of Newton (2002) can also quickly produce a smooth estimator of the mixing
density, and it was shown in Tokdar et al. (2009) and Martin et al. (2009) that the estimator is
consistent, but non-standard arguments are needed because of its dependence on the data order.
We are optimistic that the estimator described in Section 4, through the simple formula (3) for
the updates, and the well-known behavior of the kernel density estimator, can have even stronger
convergence properties than those demonstrated for predictive recursion.
In Section 5 we were able to extend the class of boundary function for which the hitting time
density can be solved using a novel Fredholm equation. Future work would consider upper and
lower boundaries. Section 6 we were able to extend the basic algorithm, which was set up for
density functions, to non–negative functions.

SUPPLEMENTAL MATERIALS
Rcpp code: It contains code to perform the iterative algorithm proposed in this paper and
numerical experiments.

19

References
Breiman, L. (1966). First exit times from a square root boundary. In Proceedings of 5th Berkeley
Symposium on Mathematical Statistics & Probability, volume 2, pages 9–16.
Chae, M., Martin, R., and Walker, S. G. (2017). Fast nonparametric near–maximum likelihood
estimation of a mixing density. In Preparation.
Corduneanu, C. (1994). Integral Equations and Applications. Cambridge University Press, Cambridge.
Csiszár, I. (1975). I-divergence geometry of probability distributions and minimization problems.
Annals of Probability, 3:146–158.
Csiszzár, I. and Tusnády, G. (1984). Information geometry and alternating minimization procedures. Statistics & Decisions, Supplemental Issue No. 1, pages 205–237.
Dykstra, R. L. (1985). An iterative procedure for obtaining I-projections onto the intersection of
convex sets. Annals of Probability, 13:975–984.
Groetsch, C. (2007). Integral equations of the first kind, inverse problems and regularization.
73:1–32.
Hansen, P. C. (1999). Linear Intergal Equations. 2nd Edition, Volume 82 of Applied Mathematical
Sciences. Springer–Verlag, NY.
Kleijn, B. J. and van der Vaart, A. W. (2006). Misspecification in infinite-dimensional Bayesian
statistics. Annals of Statistics, 34:837–877.
Laird, N. (1978). Nonparametric maximum likelihood estimation of a mixing distribution. Journal
of the American Statistical Association, 73:805–811.
Landweber, L. (1951). An iteration formula for Fredholm integral equations of the first kind.
American Journal of Mathematics, 73:615–624.
Liu, L., Levine, M., and Zhu, Y. (2009). A functional EM algorithm for mixing density estimation
via nonparametric penalized likelihood maximization. Journal of Computational and Graphical
Statistics, 18:481–504.
20

Martin, R., Tokdar, S. T., et al. (2009). Asymptotic properties of predictive recursion: robustness
and rate of convergence. Electronic Journal of Statistics, 3:1455–1472.
Morozov, V. A. (1984). Methods of Solving Incorrectly Posed Problems. Springer–Verlag, NY.
Newton, M. A. (2002). On a nonparametric recursive estimator of the mixing distribution. Sankhyā
A, 64:306–322.
Patilea, V. (2001). Convex models, MLE and misspecification. Annals of Statistics, 29:94–123.
Peskir, G. et al. (2002). On integral equations arising in the first-passage problem for Brownian
motion. Journal Integral Equations and Applications, 14:397–423.
Phillips, D. L. (1962). A technique for the numerical solution of certain integral equations of the
first kind. Journal of the Association for Computing Machinery, 9:84–97.
R Core Team (2015). R: A language and environment for statistical computing. R Foundation for
Statistical Computing, Vienna, Austria.
Ramm, A. G. (2005). Inverse Problems. Springer–Verlag, NY.
Sheather, S. J. and Jones, M. C. (1991). A reliable data-based bandwidth selection method for
kernel density estimation. Journal of the Royal Statistical Society, Series B, 53:683–690.
Shyamalkumr, N. D. (1996). Cyclic I0 projections and its applications in statistics. Technical
Report #96–24, Department of Statistics, Purdue University.
Tokdar, S. T., Martin, R., and Ghosh, J. K. (2009). Consistency of a recursive estimate of mixing
distributions. Annals of Statistics, 37:2502–2522.
Valov, A. V. (2009). First passage times: Integral equations, randomization and analytical approximations. PhD Thesis, Department of Statistics, University of Toronto.
Vangel, M. (1992). Iterative algorithms for integral equations of the first kind with applications
to statistics. Technical report, PHD Thesis, Harvard University. Technical Report ONR-C-12.

21

Vardi, Y. and Lee, D. (1993). From image deblurring to optimal investments: Maximum likelihood
solutions for positive linear inverse problems. Journal of the Royal Statistical Society, Series B,
55:569–612.
Wing, G. M. (1990). A Primer on Integral Equations of the First Kind: The Problem of Deconvolution and Unfolding. SIAM, Philadelphia.

22

