# Self-Consistent Density Estimation in the Presence of Errors-in-Variables

Authors: Junhua Zhang; Yuping Hu; Sanying Feng
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper considers the estimation of the common probability density of independent and identically distributed variables observed with additive measurement errors. The self-consistent estimator of the density function is constructed when the error distribution is known, and a modification of the self-consistent estimation is proposed when the error distribution is unknown. The consistency properties of the proposed estimators and the upper bounds of the mean square error and mean integrated square error are investigated under some suitable conditions. Simulation studies are carried out to assess the performance of our proposed method and compare with the usual deconvolution kernel method. Two real datasets are analyzed for further illustration.


---

## Body

## 1. Introduction

The statistical methodology of density estimation has been widely studied in the literature (see [1–5]) and still is a hot issue in nonparametric statistics. For example, Park et al. [6] studied the kernel-based local likelihood density estimation. Jones and Henderson [7] proposed a Gaussian copula kernel estimator. Botev et al. [8] introduced an adaptive kernel density estimation method based on the smoothing properties of linear diffusion processes. The most commonly used nonparametric method is kernel density estimation. Usually, the choice of kernel function is not crucial, whereas the bandwidth parameter, which controls the degree of smoothing, must be chosen carefully. However, it is well known that the choice of bandwidth is difficult. Cross-validation techniques have been previously applied for this, but they are computationally expensive. Furthermore, it becomes difficult to estimate the density and choose the bandwidth when there exists the measurement error.In this paper, we consider the density estimation in the presence of additive measurement error. Suppose we haven i.i.d. available observations Y
1
,
…
,
Y
n which have the same distribution as that of Y, to estimate the unknown density f
(
x
) of a random variable X, where
(1)
Y
=
X
+
ɛ
,
with a measurement error ɛ, and ɛ is independent of X. The density function of ɛ is denoted as f
ɛ.When the distribution ofɛ is known, the statistical methodology for estimating the unknown density f
(
x
) has been extensively discussed in the literature. For instance, Carroll and Hall [2] and Fan [9, 10] discussed the optimal rates of convergence over a class of functions whose derivatives are Lipschitz continuous using kernel method. Zhang [11] studied the deconvolution kernel density estimation of the mixing densities and distributions and derived the optimal rates of convergence. See also Koo [12], Pensky and Vidakovic [13], Fan and Koo [14], and Comte et al. [15], among others, for earlier contributions. Recently, Butucea [4] and Butucea and Tsybakov [5] evaluated the minimax rate of convergence of the pointwise risk using the kernel method and computed upper bounds for the L
2 risk of the estimator. However, most of those papers deal only with theoretical aspects of the estimation, and very few focus attention on the yet important issue of choosing the bandwidth in practice. Delaigle and Gijbels [16] proposed a bandwidth selection procedure based on bootstrap techniques and proved that the mean integrated squared error (MISE) of the bootstrap estimator for the density function converges to the exact mean integrated squared error. However, the algorithm is complex and the quality of the estimator of the density depends strongly on the choice of the pilot bandwidth which must be chosen before the bootstrap procedure.Indeed, in most practical applications, the distribution ofɛ cannot be perfectly known. Meister [17] studied the effect of the misspecification of the error density on the MISE of the deconvolution estimator. He pointed out that the limit of MISE can be infinite in some cases. Sometimes, this problem can be solved by repeated observations of the same variable of interest, each time with an independent error; see, for example, Li and Vuong [18], Delaigle et al. [19], or Neumann [20] and references therein. However, there are also many application fields where it is not possible to do repeated measurements of the same variable. In that case, the information of the error distribution can be drawn from an additional experiment; this means that knowledge of f
ɛ can be replaced by observations ɛ
-
1
,
…
,
ɛ
-
M, which is a noise sample with distribution f
ɛ, and is independent of Y
1
,
…
,
Y
n. Neumann [21] and Kerkyacharian et al. [22] replaced the characteristic function of the error by its empirical version. Johannes [23] studied the density deconvolution with unknown (but observed) noise and showed that the proposed estimators are asymptotically optimal in a minimax sense over a Sobolev space. The resulting estimators depend on two bandwidth-type parameters, but the data-driven selection of these bandwidths was not done. Recently, Comte and Lacour [24] proposed an adaptive estimator of the density based on a data-driven model selection strategy and discussed the convergence rates of the estimator. In addition, they studied the link between M and n if one wants to preserve the rate that is found in the case where the distribution of ɛ is known.In the present paper, we study model (1), completed with a known or unknown distribution of error ɛ, and construct the self-consistent estimators of the density f
(
x
) by searching for the optimal kernel when the error distribution is ordinary smooth or super smooth. The method of constructing the estimator of density function by optimal kernel was first proposed by Watson and Leadbetter [25]. But the result can not be used directly, since the Fourier transform of the true density function is unknown. Bernacchia and Pigolotti [26] derived the exact expression of the density estimator based on optimal kernel and defined it as self-consistent estimate. They showed that the self-consistent method has better performance than existing methods for all examples that they studied. But the data in their paper did not include measurement errors. The self-consistent method shares some desired features: the choice of the bandwidth-type parameter is more convenient than that of kernel or spline method; thus the computing speed can be improved significantly; the resulting estimators are consistent; the proposed method is preferable for applications, especially for the large datasets. Thus, it is of great significance to extend the self-consistent estimate method to more general case.The paper is organized as follows. In Section2, we propose the self-consistent method for density estimation with the known distribution of ɛ, or with the noise sample ɛ
-
1
,
…
,
ɛ
-
M when f
ɛ is unknown, and give the asymptotic properties. In Section 3, some simulations are carried out to illustrate the efficacy of the proposed method. Two real data examples are used for illustration in Section 4. The proofs of the main results are included in Appendix.


## 2. Methodology and the Main Results

### 2.1. Error Distribution Is Known

For the sake of descriptive convenience, we first introduce some notations. For two real numbersa and b, we denote a
∧
b
=
min
⁡
(
a
,
b
) and a
+
=
max
⁡
(
a
,
0
). Let z be a complex number, let z
¯ denote its conjugate, and let |
z
| denote its modulus. Let g be the L
2-norm of g; that is, g
=
∫
R
‍
|
g
(
x
)
|
2
d
x. The Fourier transform of g is defined by
(2)
ϕ
g
t
=
∫
‍
exp
⁡
i
t
x
g
x
d
x
.
Similarly, we denote the characteristic function as ϕ
ɛ
(
t
) for the known density function f
ɛ of ɛ.The smoothness off
ɛ is described by the following assumption. There exist some positive constants d
0, d
1, β, and b and constants β
0
∈
R (β
0
>
0 if β
=
0), such that the characteristic function ϕ
ɛ
(
t
) of the error distribution satisfies
(3)
d
0
t
2
+
1
-
β
0
/
2
exp
⁡
-
b
t
β
≤
ϕ
ɛ
t
≤
d
1
t
2
+
1
-
β
0
/
2
exp
⁡
-
b
t
β
as
t
⟶
∞
.
When β
>
0 we call the distribution of ɛ as super smooth of order β. For example, standard normal and Cauchy distributions are super smooth with β
=
2 and 1, respectively. When β
=
0 we call the distribution of ɛ as ordinary smooth of order β
0. For example, gamma distribution G
a
(
α
,
λ
) and double exponential distribution are ordinary smooth with β
0
=
α and 2, respectively.Assume that the density functionf
(
x
) of X belongs to the following type of smoothness spaces:
(4)
F
r
,
a
,
δ
m
=
f
density:
∫
‍
ϕ
f
x
2
x
2
+
1
δ
exp
⁡
2
a
x
r
d
x
≤
m
,
where m
>
0, r
≥
0, a
>
0, δ
∈
R, and δ
>
1
/
2 if r
=
0. If r
>
0, the density function f
(
x
) is called super smooth and ordinary smooth otherwise.Letf
Y denote the common density of the Y
j’s. By (1) and the independence assumption between X and ɛ, we have
(5)
f
Y
=
f
*
f
ɛ
,
ϕ
Y
=
ϕ
f
ϕ
ɛ
,
where * denotes the convolution and ϕ
Y denotes the characteristic function of Y. Let
(6)
f
^
Y
y
=
1
n
∑
j
=
1
n
‍
K
y
-
Y
j
be a kernel estimator of f
Y with the kernel function K
(
·
). Then, the Fourier transform of f
^
Y is defined by ϕ
^
Y
(
t
)
=
ϕ
K
(
t
)
Δ
(
t
), where Δ
(
t
)
=
(
1
/
n
)
∑
j
=
1
n
exp
⁡
(
i
t
Y
j
) is the empirical characteristic function of Y. Therefore, under the classical assumption ϕ
ɛ
(
t
)
≠
0, for any t
∈
R, we can obtain an estimator of ϕ
f as follows:
(7)
ϕ
^
f
t
=
ϕ
K
t
Δ
t
ϕ
ɛ
t
.
By the inverse Fourier transform, we have
(8)
f
^
x
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
x
ϕ
K
t
Δ
t
ϕ
ɛ
t
d
t
.
Note that (8) does not depend on any bandwidth h, compared to the deconvolution kernel density estimation methods in the literature. Instead of looking for an arbitrary shape for the kernel K
(
·
) and choosing an optimal bandwidth (see [1]), we rather look for an optimal shape of the kernel K
(
·
) such that the estimator f
^
(
x
) minimizes the mean integrated square error
(9)
MISE
=
E
∫
-
∞
+
∞
‍
f
^
x
-
f
x
2
d
x
.
By means of Parseval’s theorem and (7), we have
(10)
MISE
=
1
2
π
E
∫
-
∞
+
∞
‍
ϕ
^
f
t
-
ϕ
f
t
2
d
t
=
1
2
π
E
∫
-
∞
+
∞
‍
ϕ
K
t
Δ
t
ϕ
ɛ
t
-
ϕ
f
t
2
d
t
=
1
2
π
E
∫
-
∞
+
∞
‍
ϕ
^
Y
t
-
ϕ
Y
t
2
ϕ
ɛ
t
2
d
t
.
Note that E
(
Δ
(
t
)
)
=
ϕ
Y
(
t
) and E
(
|
Δ
(
t
)
|
2
)
=
|
ϕ
Y
(
t
)
|
2
+
n
-
1
(
1
-
|
ϕ
Y
(
t
)
|
2
), and the MISE can be rewritten as
(11)
MISE
=
1
2
π
∫
-
∞
+
∞
‍
ϕ
ɛ
t
2
-
1
n
-
1
ϕ
K
t
2
1
-
ϕ
Y
t
2
+
ϕ
Y
t
2
×
1
-
ϕ
K
t
2
×
ϕ
ɛ
t
2
-
1
d
t
=
1
2
π
∫
-
∞
+
∞
‍
1
+
n
-
1
ϕ
Y
t
2
n
ϕ
ɛ
t
2
×
ϕ
K
(
t
)
-
n
ϕ
Y
t
2
1
+
n
-
1
ϕ
Y
t
2
2
d
t
+
1
2
π
∫
-
∞
+
∞
‍
ϕ
Y
t
2
1
-
ϕ
Y
t
2
ϕ
ɛ
t
2
1
+
n
-
1
ϕ
Y
t
2
d
t
.
Then we get the optimal kernel, which in Fourier space reads
(12)
ϕ
K
opt
t
=
n
n
-
1
+
ϕ
Y
t
-
2
.
Replacing ϕ
K
(
t
) in (7) by the optimal kernel (12), we have
(13)
ϕ
^
f
t
=
ϕ
K
opt
t
Δ
t
ϕ
ɛ
t
=
n
Δ
t
n
-
1
+
ϕ
Y
t
-
2
ϕ
ɛ
t
=
n
Δ
(
t
)
n
-
1
+
ϕ
f
t
ϕ
ɛ
t
-
2
ϕ
ɛ
t
.
Similar to the discussion of Bernacchia and Pigolotti [26], the self-consistent estimator of ϕ
f is equivalent to solving the equation
(14)
ϕ
^
f
sc
t
=
n
Δ
t
n
-
1
+
ϕ
^
f
sc
t
ϕ
ɛ
t
-
2
ϕ
ɛ
t
.
Simple calculation yields
(15)
ϕ
^
f
sc
(
t
)
=
n
Δ
t
2
n
-
1
ϕ
ɛ
t
1
+
1
-
4
n
-
1
n
2
Δ
t
2
I
A
t
,
where I
A
(
·
) is the indicator function (I
A
(
t
)
=
1 if t
∈
A; I
A
(
t
)
=
0 if t
∉
A) and A
⊆
B, B
=
{
t
:
|
Δ
(
t
)
|
2
≥
4
(
n
-
1
)
/
n
2
}, is a set of accepted frequencies (i.e., the frequencies giving a nonzero contribution to the estimate). In practical applications, when the sample Y
1
,
…
,
Y
n has been observed, we can choose a bounded interval
(16)
A
=
B
∩
[
-
t
~
,
t
~
]
,
where t
~ is a truncation parameter and t
~ can be chosen as follows.Step 1.
Give an initial valuet
~
0. In our simulation, we take the initial value t
~
0
=
π
/
(
Y
(
n
)
-
Y
(
1
)
), where Y
(
1
) and Y
(
n
) are the minimum and maximum sample quantiles of {
Y
1
,
…
,
Y
n
}.Step 2.
Sett
~
=
t
~
k, if the inequality condition Δ
t
~
2
≥
4
n
-
1
/
n
2 holds in one half of the interval -
t
~
k
,
t
~
k. Otherwise, set t
~
k
+
1
=
c
t
~
k, where c
>
1 is a positive constant, in our simulation we take c
=
1.2.Step 3.
Repeat Step2 until t
~
k satisfies the inequality condition.From the simulation studies in Section3, we can see that the choice of t
~ is convenient and fast, and the density estimator is not sensitive to the choice of t
~.By (15) and the inverse Fourier transform, the self-consistent estimate of f
(
x
) is defined by
(17)
f
^
sc
x
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
x
ϕ
^
f
sc
t
d
t
.
Let
(18)
K
*
u
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
u
n
2
n
-
1
ϕ
ɛ
t
×
1
+
1
-
4
n
-
1
n
2
Δ
t
2
I
A
t
d
t
.
Then, f
^
sc
(
x
) can be rewritten in the kernel form
(19)
f
^
sc
x
=
1
n
∑
j
=
1
n
‍
K
*
x
-
Y
j
.We state the asymptotic behavior of the self-consistent estimator in the following theorem.Theorem 1.
In addition to the regularity conditions (A1)–(A3) listed in the Appendix, assume that either of the following conditions holds.(a)
The error distribution is super smooth andt
~
=
o
(
1
/
2
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β.
(b)
The error distribution is ordinary smooth andt
~
=
o
n
1
/
2
(
β
0
+
1
),where γ
=
(
1
+
2
β
0
-
β
)
+
/
β.
Then the self-consistent density estimatorf
^
s
c
(
x
) defined by (17) is consistent; that is, f
^
s
c
(
x
)
→
P
f
(
x
), as n
→
∞.Remark 2.
As many other methods, for example, higher-order kernel estimators, spline estimators, wavelet estimators, orthogonal expansion estimators, and so forth, the resulting self-consistent estimators are not nonnegative, but those can be corrected without any error cost by translating the estimate downwards until the positive part is normalized to 1 and setting to 0 the negative part. For example, the modified estimator is defined byf
ˇ
sc
(
x
)
=
max
⁡
{
f
^
sc
(
x
)
-
d
,
0
}, where d is chosen in such a way that ∫
-
∞
+
∞
‍
f
ˇ
sc
(
x
)
d
x
=
1; see more details in Efromovich [27], Glad et al. [28], and so on.Theorem 3.
Assume that conditions (A1)–(A3) hold andf
∈
F
r
,
a
,
δ
(
m
).(a)
When the error distribution is super smooth andt
~
=
O
(
2
/
5
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β, one has
(20)
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
1
-
r
-
2
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
,
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
.
(b)
When the error distribution is ordinary smooth andt
~
=
O
n
3
/
5
(
2
β
0
+
1
), one has
(21)
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
1
-
r
-
2
δ
exp
⁡
-
2
a
t
~
r
+
n
-
2
/
5
,
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
2
/
5
.


### 2.2. Error Distribution Is Unknown

Suppose that, in addition to a sampleY
1
,
…
,
Y
n from f
Y, we observe a noise sample ɛ
-
1
,
…
,
ɛ
-
M from the unknown distribution f
ɛ. We assume that information about f
ɛ is obtained by the preliminary calibration measures before the work or is provided by an independent experiment.When the characteristic functionϕ
ɛ in (15) is unknown, we then cannot estimate the density function using (17) directly. Thus, we need to estimate the characteristic function ϕ
ɛ. Here, we can estimate ϕ
ɛ based on the preliminary noise sample as follows:
(22)
ϕ
^
ɛ
t
=
1
M
∑
j
=
1
M
‍
exp
⁡
i
t
ɛ
-
j
.
Similar to Neumann [21], we introduce the following truncated estimator:
(23)
1
ϕ
~
ɛ
t
=
1
{
|
ϕ
^
ɛ
(
t
)
|
≥
M
-
1
/
2
}
ϕ
^
ɛ
t
=
1
ϕ
^
ɛ
t
,
if
ϕ
^
ɛ
t
≥
M
-
1
/
2
,
0
,
otherwise
.
Then we can obtain the self-consistent estimator of the density function
(24)
f
~
sc
x
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
x
ϕ
~
f
sc
t
d
t
,
where ϕ
~
f
sc
(
t
)
=
(
n
Δ
t
/
2
n
-
1
ϕ
~
ɛ
t
)
1
+
1
-
(
4
n
-
1
/
n
2
|
Δ
(
t
)
|
2
)
I
A
(
t
).Theorem 4.
In addition to the regularity conditions (A1)–(A3) listed in the Appendix, assume that either of the following conditions holds.(a)
The error distribution is super smooth,β
≥
r, M
≥
n
log
⁡
(
n
)
-
(
1
+
2
(
β
0
∧
δ
)
-
β
)
+
/
β, and t
~
=
o
(
1
/
2
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β.
(b)
The error distribution is super smooth,β
<
r, M
→
∞ as n
→
∞, and t
~
=
o
(
1
/
2
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β.
(c)
The error distribution is ordinary smooth,β
0
>
δ, M
≥
n
2
(
β
0
-
δ
)
/
(
2
β
0
+
1
), and t
~
=
o
(
n
1
/
(
2
β
0
+
1
)
).
(d)
The error distribution is ordinary smooth,β
0
≤
δ, M
→
∞ as n
→
∞, and t
~
=
o
(
n
1
/
(
2
β
0
+
1
)
).Then the self-consistent density estimator f
~
s
c
(
x
) defined by (24) is consistent; that is, f
~
s
c
(
x
)
→
P
f
(
x
), as n
→
∞.Remark 5.
In the cases0
<
β
<
r or β
0
≤
δ and β
=
0, from the proof of theorem in Appendix, one can check that the estimator defined by (24) is consistent as long as the noise sample size M tends to infinity as n increases, while the rate of M tending to infinity is not restricted. Therefore, in practice, the noise sample size M does not need to be very large in these cases.Theorem 6.
Assume that conditions (A1)–(A3) hold andf
∈
F
r
,
a
,
δ
(
m
).(a)
When the error distribution is super smooth,β
≥
r, M
≥
n
log
⁡
(
n
)
-
(
1
+
2
(
β
0
∧
δ
)
-
β
)
+
/
β, and t
~
=
O
(
2
/
5
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β, then one has
(25)
E
f
~
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
.
(b)
When the error distribution is super smooth,β
<
r, M
≥
n
1
/
5
log
⁡
(
n
)
-
γ
/
5, and t
~
=
O
(
2
/
5
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β, then one has
(26)
E
f
~
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
.
(c)
The error distribution is ordinary smooth,M
≥
n
(
2
/
5
)
+
(
6
(
β
0
-
δ
)
+
/
5
2
β
0
+
1
), and t
~
=
O
n
3
/
5
(
2
β
0
+
1
); then one has
(27)
E
f
~
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
2
/
5
.


## 2.1. Error Distribution Is Known

For the sake of descriptive convenience, we first introduce some notations. For two real numbersa and b, we denote a
∧
b
=
min
⁡
(
a
,
b
) and a
+
=
max
⁡
(
a
,
0
). Let z be a complex number, let z
¯ denote its conjugate, and let |
z
| denote its modulus. Let g be the L
2-norm of g; that is, g
=
∫
R
‍
|
g
(
x
)
|
2
d
x. The Fourier transform of g is defined by
(2)
ϕ
g
t
=
∫
‍
exp
⁡
i
t
x
g
x
d
x
.
Similarly, we denote the characteristic function as ϕ
ɛ
(
t
) for the known density function f
ɛ of ɛ.The smoothness off
ɛ is described by the following assumption. There exist some positive constants d
0, d
1, β, and b and constants β
0
∈
R (β
0
>
0 if β
=
0), such that the characteristic function ϕ
ɛ
(
t
) of the error distribution satisfies
(3)
d
0
t
2
+
1
-
β
0
/
2
exp
⁡
-
b
t
β
≤
ϕ
ɛ
t
≤
d
1
t
2
+
1
-
β
0
/
2
exp
⁡
-
b
t
β
as
t
⟶
∞
.
When β
>
0 we call the distribution of ɛ as super smooth of order β. For example, standard normal and Cauchy distributions are super smooth with β
=
2 and 1, respectively. When β
=
0 we call the distribution of ɛ as ordinary smooth of order β
0. For example, gamma distribution G
a
(
α
,
λ
) and double exponential distribution are ordinary smooth with β
0
=
α and 2, respectively.Assume that the density functionf
(
x
) of X belongs to the following type of smoothness spaces:
(4)
F
r
,
a
,
δ
m
=
f
density:
∫
‍
ϕ
f
x
2
x
2
+
1
δ
exp
⁡
2
a
x
r
d
x
≤
m
,
where m
>
0, r
≥
0, a
>
0, δ
∈
R, and δ
>
1
/
2 if r
=
0. If r
>
0, the density function f
(
x
) is called super smooth and ordinary smooth otherwise.Letf
Y denote the common density of the Y
j’s. By (1) and the independence assumption between X and ɛ, we have
(5)
f
Y
=
f
*
f
ɛ
,
ϕ
Y
=
ϕ
f
ϕ
ɛ
,
where * denotes the convolution and ϕ
Y denotes the characteristic function of Y. Let
(6)
f
^
Y
y
=
1
n
∑
j
=
1
n
‍
K
y
-
Y
j
be a kernel estimator of f
Y with the kernel function K
(
·
). Then, the Fourier transform of f
^
Y is defined by ϕ
^
Y
(
t
)
=
ϕ
K
(
t
)
Δ
(
t
), where Δ
(
t
)
=
(
1
/
n
)
∑
j
=
1
n
exp
⁡
(
i
t
Y
j
) is the empirical characteristic function of Y. Therefore, under the classical assumption ϕ
ɛ
(
t
)
≠
0, for any t
∈
R, we can obtain an estimator of ϕ
f as follows:
(7)
ϕ
^
f
t
=
ϕ
K
t
Δ
t
ϕ
ɛ
t
.
By the inverse Fourier transform, we have
(8)
f
^
x
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
x
ϕ
K
t
Δ
t
ϕ
ɛ
t
d
t
.
Note that (8) does not depend on any bandwidth h, compared to the deconvolution kernel density estimation methods in the literature. Instead of looking for an arbitrary shape for the kernel K
(
·
) and choosing an optimal bandwidth (see [1]), we rather look for an optimal shape of the kernel K
(
·
) such that the estimator f
^
(
x
) minimizes the mean integrated square error
(9)
MISE
=
E
∫
-
∞
+
∞
‍
f
^
x
-
f
x
2
d
x
.
By means of Parseval’s theorem and (7), we have
(10)
MISE
=
1
2
π
E
∫
-
∞
+
∞
‍
ϕ
^
f
t
-
ϕ
f
t
2
d
t
=
1
2
π
E
∫
-
∞
+
∞
‍
ϕ
K
t
Δ
t
ϕ
ɛ
t
-
ϕ
f
t
2
d
t
=
1
2
π
E
∫
-
∞
+
∞
‍
ϕ
^
Y
t
-
ϕ
Y
t
2
ϕ
ɛ
t
2
d
t
.
Note that E
(
Δ
(
t
)
)
=
ϕ
Y
(
t
) and E
(
|
Δ
(
t
)
|
2
)
=
|
ϕ
Y
(
t
)
|
2
+
n
-
1
(
1
-
|
ϕ
Y
(
t
)
|
2
), and the MISE can be rewritten as
(11)
MISE
=
1
2
π
∫
-
∞
+
∞
‍
ϕ
ɛ
t
2
-
1
n
-
1
ϕ
K
t
2
1
-
ϕ
Y
t
2
+
ϕ
Y
t
2
×
1
-
ϕ
K
t
2
×
ϕ
ɛ
t
2
-
1
d
t
=
1
2
π
∫
-
∞
+
∞
‍
1
+
n
-
1
ϕ
Y
t
2
n
ϕ
ɛ
t
2
×
ϕ
K
(
t
)
-
n
ϕ
Y
t
2
1
+
n
-
1
ϕ
Y
t
2
2
d
t
+
1
2
π
∫
-
∞
+
∞
‍
ϕ
Y
t
2
1
-
ϕ
Y
t
2
ϕ
ɛ
t
2
1
+
n
-
1
ϕ
Y
t
2
d
t
.
Then we get the optimal kernel, which in Fourier space reads
(12)
ϕ
K
opt
t
=
n
n
-
1
+
ϕ
Y
t
-
2
.
Replacing ϕ
K
(
t
) in (7) by the optimal kernel (12), we have
(13)
ϕ
^
f
t
=
ϕ
K
opt
t
Δ
t
ϕ
ɛ
t
=
n
Δ
t
n
-
1
+
ϕ
Y
t
-
2
ϕ
ɛ
t
=
n
Δ
(
t
)
n
-
1
+
ϕ
f
t
ϕ
ɛ
t
-
2
ϕ
ɛ
t
.
Similar to the discussion of Bernacchia and Pigolotti [26], the self-consistent estimator of ϕ
f is equivalent to solving the equation
(14)
ϕ
^
f
sc
t
=
n
Δ
t
n
-
1
+
ϕ
^
f
sc
t
ϕ
ɛ
t
-
2
ϕ
ɛ
t
.
Simple calculation yields
(15)
ϕ
^
f
sc
(
t
)
=
n
Δ
t
2
n
-
1
ϕ
ɛ
t
1
+
1
-
4
n
-
1
n
2
Δ
t
2
I
A
t
,
where I
A
(
·
) is the indicator function (I
A
(
t
)
=
1 if t
∈
A; I
A
(
t
)
=
0 if t
∉
A) and A
⊆
B, B
=
{
t
:
|
Δ
(
t
)
|
2
≥
4
(
n
-
1
)
/
n
2
}, is a set of accepted frequencies (i.e., the frequencies giving a nonzero contribution to the estimate). In practical applications, when the sample Y
1
,
…
,
Y
n has been observed, we can choose a bounded interval
(16)
A
=
B
∩
[
-
t
~
,
t
~
]
,
where t
~ is a truncation parameter and t
~ can be chosen as follows.Step 1.
Give an initial valuet
~
0. In our simulation, we take the initial value t
~
0
=
π
/
(
Y
(
n
)
-
Y
(
1
)
), where Y
(
1
) and Y
(
n
) are the minimum and maximum sample quantiles of {
Y
1
,
…
,
Y
n
}.Step 2.
Sett
~
=
t
~
k, if the inequality condition Δ
t
~
2
≥
4
n
-
1
/
n
2 holds in one half of the interval -
t
~
k
,
t
~
k. Otherwise, set t
~
k
+
1
=
c
t
~
k, where c
>
1 is a positive constant, in our simulation we take c
=
1.2.Step 3.
Repeat Step2 until t
~
k satisfies the inequality condition.From the simulation studies in Section3, we can see that the choice of t
~ is convenient and fast, and the density estimator is not sensitive to the choice of t
~.By (15) and the inverse Fourier transform, the self-consistent estimate of f
(
x
) is defined by
(17)
f
^
sc
x
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
x
ϕ
^
f
sc
t
d
t
.
Let
(18)
K
*
u
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
u
n
2
n
-
1
ϕ
ɛ
t
×
1
+
1
-
4
n
-
1
n
2
Δ
t
2
I
A
t
d
t
.
Then, f
^
sc
(
x
) can be rewritten in the kernel form
(19)
f
^
sc
x
=
1
n
∑
j
=
1
n
‍
K
*
x
-
Y
j
.We state the asymptotic behavior of the self-consistent estimator in the following theorem.Theorem 1.
In addition to the regularity conditions (A1)–(A3) listed in the Appendix, assume that either of the following conditions holds.(a)
The error distribution is super smooth andt
~
=
o
(
1
/
2
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β.
(b)
The error distribution is ordinary smooth andt
~
=
o
n
1
/
2
(
β
0
+
1
),where γ
=
(
1
+
2
β
0
-
β
)
+
/
β.
Then the self-consistent density estimatorf
^
s
c
(
x
) defined by (17) is consistent; that is, f
^
s
c
(
x
)
→
P
f
(
x
), as n
→
∞.Remark 2.
As many other methods, for example, higher-order kernel estimators, spline estimators, wavelet estimators, orthogonal expansion estimators, and so forth, the resulting self-consistent estimators are not nonnegative, but those can be corrected without any error cost by translating the estimate downwards until the positive part is normalized to 1 and setting to 0 the negative part. For example, the modified estimator is defined byf
ˇ
sc
(
x
)
=
max
⁡
{
f
^
sc
(
x
)
-
d
,
0
}, where d is chosen in such a way that ∫
-
∞
+
∞
‍
f
ˇ
sc
(
x
)
d
x
=
1; see more details in Efromovich [27], Glad et al. [28], and so on.Theorem 3.
Assume that conditions (A1)–(A3) hold andf
∈
F
r
,
a
,
δ
(
m
).(a)
When the error distribution is super smooth andt
~
=
O
(
2
/
5
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β, one has
(20)
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
1
-
r
-
2
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
,
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
.
(b)
When the error distribution is ordinary smooth andt
~
=
O
n
3
/
5
(
2
β
0
+
1
), one has
(21)
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
1
-
r
-
2
δ
exp
⁡
-
2
a
t
~
r
+
n
-
2
/
5
,
E
f
^
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
2
/
5
.


## 2.2. Error Distribution Is Unknown

Suppose that, in addition to a sampleY
1
,
…
,
Y
n from f
Y, we observe a noise sample ɛ
-
1
,
…
,
ɛ
-
M from the unknown distribution f
ɛ. We assume that information about f
ɛ is obtained by the preliminary calibration measures before the work or is provided by an independent experiment.When the characteristic functionϕ
ɛ in (15) is unknown, we then cannot estimate the density function using (17) directly. Thus, we need to estimate the characteristic function ϕ
ɛ. Here, we can estimate ϕ
ɛ based on the preliminary noise sample as follows:
(22)
ϕ
^
ɛ
t
=
1
M
∑
j
=
1
M
‍
exp
⁡
i
t
ɛ
-
j
.
Similar to Neumann [21], we introduce the following truncated estimator:
(23)
1
ϕ
~
ɛ
t
=
1
{
|
ϕ
^
ɛ
(
t
)
|
≥
M
-
1
/
2
}
ϕ
^
ɛ
t
=
1
ϕ
^
ɛ
t
,
if
ϕ
^
ɛ
t
≥
M
-
1
/
2
,
0
,
otherwise
.
Then we can obtain the self-consistent estimator of the density function
(24)
f
~
sc
x
=
1
2
π
∫
-
∞
+
∞
‍
exp
⁡
-
i
t
x
ϕ
~
f
sc
t
d
t
,
where ϕ
~
f
sc
(
t
)
=
(
n
Δ
t
/
2
n
-
1
ϕ
~
ɛ
t
)
1
+
1
-
(
4
n
-
1
/
n
2
|
Δ
(
t
)
|
2
)
I
A
(
t
).Theorem 4.
In addition to the regularity conditions (A1)–(A3) listed in the Appendix, assume that either of the following conditions holds.(a)
The error distribution is super smooth,β
≥
r, M
≥
n
log
⁡
(
n
)
-
(
1
+
2
(
β
0
∧
δ
)
-
β
)
+
/
β, and t
~
=
o
(
1
/
2
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β.
(b)
The error distribution is super smooth,β
<
r, M
→
∞ as n
→
∞, and t
~
=
o
(
1
/
2
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β.
(c)
The error distribution is ordinary smooth,β
0
>
δ, M
≥
n
2
(
β
0
-
δ
)
/
(
2
β
0
+
1
), and t
~
=
o
(
n
1
/
(
2
β
0
+
1
)
).
(d)
The error distribution is ordinary smooth,β
0
≤
δ, M
→
∞ as n
→
∞, and t
~
=
o
(
n
1
/
(
2
β
0
+
1
)
).Then the self-consistent density estimator f
~
s
c
(
x
) defined by (24) is consistent; that is, f
~
s
c
(
x
)
→
P
f
(
x
), as n
→
∞.Remark 5.
In the cases0
<
β
<
r or β
0
≤
δ and β
=
0, from the proof of theorem in Appendix, one can check that the estimator defined by (24) is consistent as long as the noise sample size M tends to infinity as n increases, while the rate of M tending to infinity is not restricted. Therefore, in practice, the noise sample size M does not need to be very large in these cases.Theorem 6.
Assume that conditions (A1)–(A3) hold andf
∈
F
r
,
a
,
δ
(
m
).(a)
When the error distribution is super smooth,β
≥
r, M
≥
n
log
⁡
(
n
)
-
(
1
+
2
(
β
0
∧
δ
)
-
β
)
+
/
β, and t
~
=
O
(
2
/
5
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β, then one has
(25)
E
f
~
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
.
(b)
When the error distribution is super smooth,β
<
r, M
≥
n
1
/
5
log
⁡
(
n
)
-
γ
/
5, and t
~
=
O
(
2
/
5
b
)
log
⁡
n
(
log
⁡
(
n
)
)
-
γ
1
/
β, then one has
(26)
E
f
~
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
1
/
5
log
⁡
n
γ
/
5
.
(c)
The error distribution is ordinary smooth,M
≥
n
(
2
/
5
)
+
(
6
(
β
0
-
δ
)
+
/
5
2
β
0
+
1
), and t
~
=
O
n
3
/
5
(
2
β
0
+
1
); then one has
(27)
E
f
~
s
c
x
-
f
x
2
≤
O
t
~
2
+
1
-
δ
exp
⁡
-
2
a
t
~
r
+
n
-
2
/
5
.


## 3. Simulation Studies

In this section, we report on the simulation studies to illustrate the finite sample performances of the proposed self-consistent (SC) estimation method and compare it with the kernel (KN) method. For comparison, we compute the estimators for different signal densities and different types of noise. The following densities are considered:(a)
gamma distributionf
(
x
)
=
10.125
x
4
e
-
3
x with shape parameter α
=
5 and scale parameter λ
=
3;
(b)
standard normal distributionf
(
x
)
=
(
1
/
2
π
)
exp
⁡
(
-
x
2
/
2
).Two kinds of error distributions are considered to study their effects on the MISE of the estimators.(c)
Double exponential error: the density ofɛ is given by
(28)
f
ɛ
x
=
2
σ
0
-
1
exp
⁡
-
2
x
σ
0
for
σ
0
2
=
1
16
.
(d)
Normal error: the density ofɛ is given by
(29)
f
ɛ
x
=
1
2
π
σ
0
exp
⁡
-
x
2
2
σ
0
2
for
σ
0
2
=
1
16
.For kernel method, we choose the Gaussian kernel for double exponential error, and for normal error we suppose that the kernel has a Fourier transformϕ
K
(
t
)
=
(
1
-
t
2
)
+
2. For the unknown error distributions, we suppose that the noise sample ɛ
-
1
,
…
,
ɛ
-
M comes from the above two error distributions. Throughout the simulations, we set n
=
100,200,400 and M
=
50,100. The MISE of the estimators are computed with 1000 random samples, and the results are reported in Tables 1 and 2. When f
ɛ is known or unknown, the gamma density estimator through the double exponential noise and the normal density estimator through the normal noise are shown in Figures 1 and 2, respectively.Table 1
MISE for the self-consistent estimator and the kernel estimator when the error distribution is double exponential.
n
Case
Gamma
Normal
SC
KN
SC
KN
100
f
ε known
0.0021
0.0036
0.0105
0.0102
M
=
50
0.0022
0.0041
0.0105
0.0106
M
=
100
0.0021
0.0039
0.0102
0.0103
200
f
ε known
0.0011
0.0019
0.0056
0.0065
M
=
50
0.0011
0.0025
0.0063
0.0068
M
=
100
0.0013
0.0021
0.0057
0.0067
400
f
ε known
0.0006
0.0011
0.0031
0.0037
M
=
50
0.0007
0.0012
0.0036
0.0035
M
=
100
0.0006
0.0012
0.0032
0.0033Table 2
MISE for the self-consistent estimator and the kernel estimator when the error distribution is normal.
n
Case
Gamma
Normal
SC
KN
SC
KN
100
f
ε known
0.0028
0.0039
0.0102
0.0109
M
=
50
0.0030
0.0042
0.0106
0.0111
M
=
100
0.0029
0.0038
0.0106
0.0110
200
f
ε known
0.0012
0.0023
0.0061
0.0073
M
=
50
0.0015
0.0023
0.0069
0.0075
M
=
100
0.0012
0.0024
0.0059
0.0071
400
f
ε known
0.0007
0.0012
0.0035
0.0039
M
=
50
0.0007
0.0013
0.0038
0.0039
M
=
100
0.0007
0.0011
0.0036
0.0038Figure 1
Estimates of gamma density through double exponential noise forn
=
400. The thick solid line is the true density, the thin solid line and the dashed line are the kernel estimator and the SC estimator when f
ɛ is known, respectively, and the dot-dashed line and the dotted line are the SC estimators when the noise sample sizes M
=
50 and M
=
10, respectively.Figure 2
Estimates of normal density through normal noise forn
=
400. The thick solid line is the true density, the thin solid line and the dashed line are the kernel estimator and the SC estimator when f
ɛ is known, respectively, and the dot-dashed line and the dotted line are the SC estimators when the noise sample sizes M
=
50 and M
=
10, respectively.Tables1 and 2 indicate the following simulation results.(1)
When the error distribution is unknown, we can see that the estimator of the characteristic function of the error does not spoil the procedure so much. It even happens that the estimation with unknown error distribution works better, which is probably due to the truncation (23).
(2)
For the same observation sample, the self-consistent method performs slightly better than the kernel method. It is worth mentioning that the self-consistent method can save more time than the kernel method, because the truncation parametert
~ can be chosen conveniently.
(3)
The MISE obtained by self-consistent and kernel methods decrease asn increases whether the error distribution is known or unknown.Figures1 and 2 clearly show that the self-consistent estimators perform well, no matter whether the error distribution is known or unknown. Moreover, for M
=
10, the results show that our method is still very satisfactory for small noise sample size.Next, we evaluate the sensitivity of the self-consistent density estimation procedure on the choice of truncation parametert
~. In this simulation, we only consider the case when the error distribution is double exponential, and the value of t
~ is fixed at t
~
=
t
~
*, t
~
*
/
1.25, and 1.25
t
~
*, respectively, where t
~
* is obtained by the iterative algorithm given in Section 2.1. The MISE of the self-consistent estimator is reported in Table 3. From Table 3 we can see that, for given n, the performance of the proposed method does not depend sensitively on the choice of the value of t
~.Table 3
MISE for the self-consistent estimator with different values oft
~ when the error distribution is double exponential.
n
Case
Gamma
Normal
t
~
*
t
~
*
/
1.25
1.25
t
~
*
t
~
*
t
~
*
/
1.25
1.25
t
~
*
100
f
ε known
0.0020
0.0021
0.0020
0.0105
0.0105
0.0107
M
=
50
0.0022
0.0024
0.0020
0.0108
0.0106
0.0106
M
=
100
0.0021
0.0022
0.0022
0.0104
0.0103
0.0108
200
f
ε known
0.0011
0.0014
0.0012
0.0057
0.0054
0.0057
M
=
50
0.0013
0.0015
0.0013
0.0060
0.0059
0.0058
M
=
100
0.0012
0.0014
0.0012
0.0059
0.0058
0.0062
400
f
ε known
0.0006
0.0007
0.0007
0.0032
0.0034
0.0032
M
=
50
0.0007
0.0008
0.0007
0.0035
0.0034
0.0033
M
=
100
0.0006
0.0006
0.0007
0.0032
0.0035
0.0032


## 4. Application to Real Data

We now use two real datasets to illustrate the proposed method. The first dataset is the density of direction chosen by an ant to an evenly illuminated black target. The experiment with 100 ants was described by Fisher [29], who concluded that the ants tend to run toward the target with some moderate variation. For this particular example, Y, a point where an ant intersects a circle, can be treated as indirect observation of the direction X chosen by an ant. The estimation method of (17) is used to analyze the dataset. In Figure 3, the solid line denotes the estimate based on the assumption that the data can be accurately measured (i.e., Y
=
X), and the dashed and dot-dashed lines denote the case of N
(
0
,
σ
2
) measurement errors with σ
=
0.3 and σ
=
0.5, respectively. We compare our results with Efromovich [30] who already used the orthogonal series approach to estimate the density.Figure 3
Estimates of density of movement of ants. The solid line is the estimate based on direct data, and the dashed and dot-dashed lines are the estimates under assumption on normal error with zero mean and standard deviationsσ
=
0.3 and σ
=
0.5, respectively.From Figure3, we find that when the data are accurately measured (σ
=
0), the corresponding density estimate is rather smooth and has a pronounced mode (direction of the target) with a relatively large deviation. When the standard deviation σ increases, it makes the density more picky interestingly and makes two additional local modes more significant. Therefore, the assumption about measurement error reveals the presence of three distinct groups (strata) of ants. A majority clearly tends to move toward the target with a rather small variation, and a minority tends to move in two other directions about the target also with small variation. These findings basically agree with those that were discovered by Efromovich [30]. But the two directions are not symmetric here, and the variation of the right side is slightly less than that of the left side.The second example involves estimating density of the magnitudes of Alaskan earthquakes for the period from 1969 to 1978. The data comes from the National Oceanic and Atmospheric Administration’s Hypocenter Data File [31]. In this dataset, Y denotes the logarithm of the seismogram amplitude of longitudinal body waves for 62 Alaskan earthquakes. There is a measurement error associated with the observations, which includes errors made in determining the amplitude of ground motion arising from such things as the orientation of a limited number of observation stations to the fault plane of the earthquake. In Figure 4, the solid line denotes the estimate based on assumption that the body waves are accurately measured (i.e., Y
=
X), and the dashed and dot-dashed lines denote the case of Laplace measurement errors with position parameter α
=
0 and scale parameters λ
=
0.3 and λ
=
0.5, respectively.Figure 4
Estimates of density of body waves. The solid line is the estimate based on direct data, and the dashed and dot-dashed lines are the estimates under assumption on Laplace error with position parameterα
=
0 and scale parameters λ
=
0.3 and λ
=
0.5, respectively.Figure4 clearly shows the nature of the indirect data. Under the assumption that the data are accurately measured (λ
=
0), the corresponding density estimate is smooth and has two pronounced modes. For the case of λ
≠
0, the distinctions between the density estimates become more noticeable. It makes the data become a unimodal distribution with λ increasing. It is clear that the presence of the measurement error allows one to gain some insight into the data.



---

*Source: 958702-2014-12-10.xml*